[Rd] Solve.QP

[Berwin A Turlach]

G'day Serge,

On Thu, 6 Dec 2007 13:34:58 +0100
"de Gosson de Varennes Serge (4100)"
<serge.de.gosson.de.varennes at forsakringskassan.se> wrote:

> I have a major problem (major for me that is) with solve.QP and I'm
> new at this. You see, to solve my quadratic program I need to have
> the lagrange multipliers after each iteration. Solve.QP gives me the
> solution, the unconstrained solution aswell as the optimal value.
> Does anybody have an idea for how I could extract the multipliers?

You could calculate them.  For the quadratic program 

min  1/2 x'Dx - d'x   such that A'x >= b

the KKT conditions are:

( D -A) (x) = (d)
( A' 0) (l) = (b)

plus the complementary conditions.  solve.QP tells you the solution x*
and which constraints are active.  If A* is the submatrix of A that
contains only those columns of A that correspond to active constraints
at the solution, then the first line in the above set of equations
imply that the corresponding Lagrange multiplier l* fulfill the
equations:

D x* - A* l* = d  -->  l* = (A*' A*)^{-1} A*'(D x* - d)

all other Lagrange multiplier would be zero.  

Thus, using the example from the help page and expanding it, the
following code should calculate the Lagrange multipliers:

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)
res <- solve.QP(Dmat,dvec,Amat,bvec=bvec)
xst <- res$solution
tmp <- Dmat %*% xst -dvec
Ast <- Amat[,res$iact]
ll <- solve(crossprod(Ast,Ast), crossprod(Ast, tmp))
## a small check
cbind(tmp, Ast %*% ll)
lagr <- rep(0, ncol(Amat))
lagrange[res$iact] <- ll
lagrange
[1] 0.0000000 0.2380952 2.0952381


Alternatively, somewhere down in the FORTRAN code the Lagrange
multipliers are actually calculated.  About 4 years ago somebody asked
me about this and he could locate where they are calculated.  He
modified the FORTRAN code and the R code such that the Lagrange
multipliers would be returned too.  Curiously, he sent the modified
code to me and not to the package maintainer, but I had no time at that
moment to check his modification, so I never passed anything on to the
package maintainer either.  But if you are interested in modifying the
FORTRAN and R code and recompile the package yourself, I can see if I
can find that code.

I still think that this package should be worked over by someone with a
better understanding of the kind of fudges that do not come back to
bite and of finite precision arithmetic than the original author's
appreciation of such issues when the code was written. ;-))  Given
Bill's recent comments on r-help, I wonder whether this package is one
of those on his list of downright dangerous packages.  LOL.

Hope this helps.

Cheers,

	Berwin

=========================== Full address =============================
Berwin A Turlach                            Tel.: +65 6516 4416 (secr)
Dept of Statistics and Applied Probability        +65 6516 6650 (self)
Faculty of Science                          FAX : +65 6872 3919       
National University of Singapore     
6 Science Drive 2, Blk S16, Level 7          e-mail: statba at nus.edu.sg
Singapore 117546                    http://www.stat.nus.edu.sg/~statba