CODEWITHSUNDEEP
roptimizationconstraintsquadprogquadratic-programming
vonjd
vonjd

Reputation: 4380

Solve Constrained Quadratic Programming with R

I really love R but from time to time it really gives me a headache...

I have the following simple quadratic minimization problem which can be formulated and solved within no time in Excel (click on picture to enlarge):

enter image description here

and

![enter image description here

The problem itself is pretty straightforward: I want to minimize (w1^2+w2^2)/2 by finding the best combination of w1,w2 and b under the constraints that for all Y*(w1*X1+w2*X2+b) >= 1

I know that there is the quadprog package for solving these kinds of problems but I find it so unintuitive that I am not able to specify the problem correctly :-( I hate to say it but Excel seems to be better suited for specifying optimization problems like these :-(((

My question
How to formulate the above problem correctly so that it can be solved with R (no matter which package) and the program arrives at the correct values for w1, w2 and b (as can be seen in the picture above). Please don't just post links but please give actual code that works. It would be great if you could comment your code so it becomes clear why you do the things you do. Thank you!

The necessary data is here:

data <- matrix(c(2.947814,6.626878, 1,
                 2.530388,7.785050, 1,
                 3.566991,5.651046, 1,
                 3.156983,5.467077, 1,
                 2.582346,4.457777,-1,
                 2.155826,6.222343,-1,
                 3.273418,3.520687,-1),ncol=3,byrow=T)
colnames(data) <- c("X1","X2","y")

Addendum
Some people took offense at my request to provide code (and not simply links). I apologized for that and gave my reasons that I did not find any good approaches in the answers so far on SO. The deeper reason for that is that the problem is unusual in the sense that b is only in the constraint and not in the objective function. So I still think that this question is a good fit for SO.

Upvotes: 1

Views: 1515

Answers (2)

vonjd
vonjd

Reputation: 4380

For future reference: this problem can also be solved quite intuitively using the CVXR package, which provides an interface for convex optimization problems in R:

# Install and load the CVXR package
if (!requireNamespace("CVXR", quietly = TRUE)) {
  install.packages("CVXR")
}
library(CVXR)
# Define the data
data <- matrix(c(2.947814,6.626878, 1,
                 2.530388,7.785050, 1,
                 3.566991,5.651046, 1,
                 3.156983,5.467077, 1,
                 2.582346,4.457777,-1,
                 2.155826,6.222343,-1,
                 3.273418,3.520687,-1), ncol = 3, byrow = T)
colnames(data) <- c("X1", "X2", "y")

# Define the variables and parameters
w <- Variable(2)
b <- Variable()

# Define the constraints
constraints <- list()
for (i in 1:nrow(data)) {
  constraints[[i]] <- data[i, "y"] * (sum(data[i, 1:2] * w) + b) >= 1
}

# Define the objective function
objective <- Minimize(sum(w^2)/2)

# Solve the problem
problem <- Problem(objective, constraints)
result <- solve(problem)

# Print the results
cat("w1:", result$getValue(w)[1], "\nw2:", result$getValue(w)[2], "\nb:", result$getValue(b))
## w1: 2.90391 
## w2: 1.201258 
## b: -14.73496

Upvotes: 0

cryo111
cryo111

Reputation: 4474

Actually, the problem is a little tricky because b is only present in the inequality constraint matrix but not in the objective function. Therefore the matrix in the quadratic programming problem is only positive semidefinite but not positive definite.

My approach is therefore to set the matrix entry corresponding to b to a very small value - in my case 1e-9. Someone else more familiar with such optimization problems might know how to solve the problem properly...

Calculate solve.QP input

c1=data[,"X1"]*data[,"y"]
c2=data[,"X2"]*data[,"y"]

#I use 1e-9 for the b entry
Dmat=matrix(`[<-`(numeric(9),c(1,5,9),c(1,1,1e-9)),3,3)
dvec=rep(0,3)
Amat=cbind(c1,c2,data[,"y"])
bvec=rep(1,nrow(Amat))

Solve with solve.QP

library(quadprog)
sol=solve.QP(Dmat=Dmat,dvec=dvec,Amat=t(Amat),bvec=bvec)$solution
sol
#[1]   2.903910   1.201258 -14.734964

Same as excel.

Upvotes: 5

Related Questions