From: Berwin A Turlach <berwin_at_maths.uwa.edu.au>
Date: Tue 06 Jun 2006 - 21:05:29 EST

G'day Fabian,

>>>>> "FB" == Fabian Barth <Fabian.Barth@web.de> writes:

FB> I'm using the package quadprog to solve the following     FB> quadratic programming problem.

```    FB> I want to minimize the function
FB> (b_1-b_2)^2+(b_3-b_4)^2
FB> by the following constraints b_i, i=1,...,4:

```

FB> b_1+b_3=1
FB> b_2+b_4=1

```    FB> 0.1<=b_1<=0.2
FB> 0.2<=b_2<=0.4
FB> 0.8<=b_3<=0.9
FB> 0.6<=b_4<=0.8

```

FB> In my opinion the solution should be b_1=b_2=0.2 und b_3=b_4=0.8. This could well be the correct solution. For these values for the b_i the function evaluates to zero and you definitely won't get anything smaller than that. :)

```    FB> Unfortunately R doesn't find this solution and what's
FB> surprising to me, evaluation the solution of solve.QP with my
FB> function doesn't lead the minimal "value" calculated by
FB> solve.QP

```

FB> I would be very happy, if anyone could help and tell me,     FB> where's my mistake.
Mmh, the first code that you are sending seems to involve only three unknowns. Did you try to eliminate one by hand? If so, which one. Also, in that part of the code you do not tell solve.QP that there are equality constraints. So I am a bit confused what problem the first part is supposed to solve. :-)

The second code, labelled `my "non working" sample' seems to implement the above problem directly. The problem with that code is that the matrix Dmat that you construct is not symmetric, there is the (apparently undocumented and unchecked) assumption that Dmat is symmetric. (Note if the matrix D in the quadratic form is not summetric, then you can always replace it by (D+D^T)/2 without changing the objective function. Thus, without loss of generality, the matrix D in the objective function is always assumed to be symmetric). So adding code like

Dmat <- (Dmat+t(Dmat))/2
in front of your `print(Dmat)' statement would be necessary.

Two minor observations. First, note that the help page for solve.QP states that the function to be minimised is

-d^T b + 1/2 b^T D b
where D is passed in Dmat and d is passed in dvec. In your case it might not matter, but strictly speaking, you should multiply your Dmat by 2 to take care of the 1/2 factor in the definition of the objective function. Secondly, the FORTRAN code that is used by solve.QP only uses the upper triangular part of the matrix Dmat that is passed to solve.QP. Thus, essentially you were minimizing

min 1/2 (b_1^2 + b_2^2 + b_3^2 + b_4^2) under all the conditions stated. And, for that problem solve.QP gave the correct answer.

Once you correct the problem, and a faster way of coding your problem is probabably the following:

> Dmat1 <- matrix(c(2, -2, -2, 2), 2, 2)
> Dmat2 <- matrix(0, 2, 2)
> Dmat <- rbind(cbind(Dmat1, Dmat2), cbind(Dmat2, Dmat1))
> dvec <- rep(0,4)
> Amat <- cbind(c(1,0,1,0), c(0,1,0,1), diag(4), -diag(4))
> bvec <- c(1, 1, 0.1, 0.2, 0.8, 0.6, -0.2, -0.4, -0.9, -0.8)

you will see that there is still a problem:

> solve.QP(Dmat,dvec,Amat,bvec=bvec, meq=2)
Error in solve.QP(Dmat, dvec, Amat, bvec = bvec, meq = 2) :

matrix D in quadratic function is not positive definite!

The Goldfarb-Idnani algorithm starts off by calculating the unconstrained solution. Thus, it requires that the matrix D in the objective function is positive definite. In your problem the matrix is indefinite.

Hope this helps.

Cheers,

Berwin

• Full address ============================ Berwin A Turlach Tel.: +61 (8) 6488 3338 (secr) School of Mathematics and Statistics +61 (8) 6488 3383 (self) The University of Western Australia FAX : +61 (8) 6488 1028 35 Stirling Highway Crawley WA 6009 e-mail: berwin@maths.uwa.edu.au Australia http://www.maths.uwa.edu.au/~berwin

R-help@stat.math.ethz.ch mailing list