# [R] solve.QP with box and equality constraints

Berwin A Turlach berwin at maths.uwa.edu.au
Mon Feb 16 16:13:09 CET 2009

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.
>
> 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