[R] solve.QP with box and equality constraints

[Berwin A Turlach]

G'day Selwyn,

On Mon, 16 Feb 2009 10:11:20 +0000
Selwyn McCracken <selwyn.mccracken at gmail.com> wrote:

> I am trying to follow an example that estimates a 2x2 markov
> transition matrix across several periods from aggregate data using
> restricted least squares.
> 
> I seem to be making headway using solve.QP(quadprog) as the
> unrestricted solution matches the example I am following, and I can
> specify simple equality and inequality constraints. However, I cannot
> correctly specify a constraint matrix (Amat) with box constraints per
> cell and equality constraints that span multiple cells. Namely the
> solution matrix I am aiming for needs to respect the following
> conditions:
>   - each row must sum to 1
>   - each cell must >= 0
>   - each cell must <= 1

Note that this set of constraints contains some redundancy.  If each
row has to sum to one and the summands are all negative then,
necessarily, each summand is at most one.  So the last set of
constraints is not necessary.

[...]
> (b2 = c(1,1,rep(0:1,times=4)))

This should be 
   (b2 <- c(1,1, rep(0:(-1), times=4)))

Since x <= 1 iff -x >= -1

> solve.QP(Dmat = XX,dvec=Xy,Amat = a2,bvec=b2,meq=2) #complains that
> Amat and dvec are incompatible

The constraints of the QP are A^T b >= b0 and b0 is passed to b2 while
A is passed to Amat.  Your matrix a2 contains A^T, so you have to pass
t(a2) to solve.QP:

R> solve.QP(Dmat = XX,dvec=Xy,Amat = t(a2),bvec=b2,meq=2)
Error in solve.QP(Dmat = XX, dvec = Xy, Amat = t(a2), bvec = b2, meq = 2) : 
  constraints are inconsistent, no solution!

While this might not be very helpful, the problem now is that the
entries in your matrix XX (and Xy) are quite large and this can lead to
numerical problems in the FORTRAN code that quadprog relies on.  But
this is easily fixed.

R> XX <- XX*1e-9
R> Xy <- Xy*1e-9
R> solve.QP(Dmat = XX,dvec=Xy,Amat = t(a2),bvec=b2,meq=2)
$solution
[1] 0.9367934 0.0000000 0.0632066 1.0000000

$value
[1] -63.38694

$unconstrainted.solution
[1] 1.004181694 0.002401266 0.012673485 1.012848987

$iterations
[1] 4 0

$iact
[1] 10  1  2

But, as pointed out earlier, your constraints contain some 
redundancies, so it would be shorter to code them as:

R> a2 = rbind(
   c(1,0,1,0), #specify the two cells in the 'row' I want to sum to one(???)
   c(0,1,0,1), # likewise
   diag(rep(1,4))
   )
R> b2 <- c(1,1, rep(0, times=4))
R> solve.QP(Dmat = XX,dvec=Xy,Amat = t(a2),bvec=b2,meq=2)
$solution
[1] 0.9367934 0.0000000 0.0632066 1.0000000

$value
[1] -63.38694

$unconstrainted.solution
[1] 1.004181694 0.002401266 0.012673485 1.012848987

$iterations
[1] 4 0

$iact
[1] 1 2 4


HTH.

Cheers,

	Berwin

=========================== Full address =============================
Berwin A Turlach                            Tel.: +65 6515 4416 (secr)
Dept of Statistics and Applied Probability        +65 6515 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