Constrained quadratic optimization with the quadProg library

2.7k views Asked by At

I have a vector A of length N. Also I have N*N matrix C. I want to maximize following equation :

minimize (- (w_transpose * A) + p * w_transpose * C * w)

Where w is a vector of length N, with constraints that each w is non-negative and sum of all w is 1.

I have seen a package called quadProg. There I need to specify :

Dmat = C, dvec = A, and bvec = w

but not sure how to apply above mentioned constraints there.

I suppose I could provide Amat to be an identity matrix, which will keep all w non-negative. But not sure how to keep w normalized (sum equal to zero). Actually I could normalize them later as well, but still wondering if I can do it here itself.

1

There are 1 answers

0
josliber On BEST ANSWER

You can do this with the solve.QP function from quadprog. From ?solve.QP, we read that solve.QP solves systems of the form min_b {-d'b + 0.5 b'Db | A'b >= b0}. You are solving a problem of the form min_w {-A'w + pw'Cw | w >= 0, 1'w = 1}. Thus, mapping between the forms:

  • d = A (this is called dvec in the arguments to solve.QP)
  • D = 2pC (this is called Dmat in the arguments to solve.QP)
  • For the first set of constraints, you have I'w >= 0. The final constraint can be reformulated as 1'w >= 1 and -1'w >= -1. Therefore your A matrix of constraints (Amat in the arguments to solve.QP) is the identity matrix with a 1 vector and a -1 vector appended on the right, and the right-hand side b0 (bvec in the arguments to solve.QP) is the 0 vector with a 1 and -1 appended.

You can put it all together in R pretty easily:

library(quadprog)
solve.my.QP <- function(A, p, C) {
  solve.QP(Dmat=2*p*C,
           dvec=A,
           Amat=cbind(diag(1, length(A)), rep(1, length(A)), rep(-1, length(A))),
           bvec=c(rep(0, length(A)), 1, -1))$solution
}

You can test it on some simple 2-dimensional problems:

# Even penalty
solve.my.QP(c(0, 0), 1, diag(1, 2))
# [1] 0.5 0.5

# Encourage inclusion of first variable
solve.my.QP(c(0.1, 0), 1, diag(1, 2))
# [1] 0.525 0.475