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

Selwyn McCracken selwyn.mccracken at gmail.com
Mon Feb 16 16:42:35 CET 2009

Berwin,

Have a good day :-)

Selwyn

On Mon, Feb 16, 2009 at 3:13 PM, Berwin A Turlach
<berwin at maths.uwa.edu.au> wrote:
> 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
>