[R] solve.QP with box and equality constraints

[Selwyn McCracken]

Berwin,

Many thanks for your reply - your solution is exactly what I was
looking for and the additional advice is also much appreciated.

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