R help - Jun 2006 - Problems using quadprog for solving quadratic programming problem

Hi,

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

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


b_1+b_3=1
b_2+b_4=1

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

In my opinion the solution should be b_1=b_2=0.2 und b_3=b_4=0.8.

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

I would be very happy, if anyone could help and tell me, where's my mistake.
Thank you very much. Fabian

My R-code also containing the sampel of quadprog starts here:

#sample from quadprog package
library(quadprog)

Dmat       <- matrix(0,3,3)
diag(Dmat) <- 1
dvec       <- c(0,5,0)
Amat       <- matrix(c(-4,-3,0,2,1,0,0,-2,1),3,3)
bvec       <- c(-8,2,0)
erg<-solve.QP(Dmat,dvec,Amat,bvec=bvec)

print(erg)


erg<-erg$solution
-dvec%*%erg+.5*erg^T%*%Dmat%*%erg


# my "non working" sample

n<-2
k<-2
Dmatj<-diag((n-1),n)
Dmatj[lower.tri(Dmatj)]<--2
Dmat<-matrix(0,nrow=n*k,ncol=n*k)
for(j in 1:k){
  Dmat[((j-1)*n+1):(j*n),((j-1)*n+1):(j*n)]<-Dmatj
}
print(Dmat)
Amat<-matrix(ncol=n,nrow=k*n,0)
ind<-seq(1,((k-1)*n+1),n)
for(i in 1:n){
  Amat[(ind+(i-1)),i]<-1
}
Amat<-cbind(Amat,diag(n*k),-diag(n*k))
print(Amat)
bvec<-c(rep(1,n),c(0.1,0.2,0.8,0.6,-0.2,-0.4,-0.9,-0.8))
print(bvec)
dvec=rep(0,(n*k))
  
erg<-solve.QP(Dmat,dvec,Amat,bvec,meq=n)

print(erg)


erg<-erg$solution
-dvec%*%erg+.5*erg^T%*%Dmat%*%erg

G'day Fabian,
>>>>> "FB" == Fabian Barth <Fabian.Barth at
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 at maths.uwa.edu.au
Australia                        http://www.maths.uwa.edu.au/~berwin