[R] Is it solve.QP or is it me?
Talbot Katz
topkatz at msn.com
Fri Sep 21 18:38:05 CEST 2007
Hi.
Here are three successive examples of simple quadratic programming problems
with the same structure. Each problem has 2*N variables, and should have a
solution of the form (1/N,0,1/N,0,...,1/N,0). In these cases, N=4,5,6. As
you will see, the N=4 and 6 cases give the expected solution, but the N=5
case breaks down.
>cm8
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0
[5,] 0 0 0 0 1 0 0 0
[6,] 0 0 0 0 0 1 0 0
[7,] 0 0 0 0 0 0 1 0
[8,] 0 0 0 0 0 0 0 1
>dv8
[1] 0 0 0 0 0 0 0 0
>Am8
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 1 -1 0 0 0 0 0 0 0
[2,] -1 1 0 1 0 0 0 0 0 0
[3,] 1 1 0 0 -1 0 0 0 0 0
[4,] -1 1 0 0 0 1 0 0 0 0
[5,] 1 1 0 0 0 0 -1 0 0 0
[6,] -1 1 0 0 0 0 0 1 0 0
[7,] 1 1 0 0 0 0 0 0 -1 0
[8,] -1 1 0 0 0 0 0 0 0 1
>bv8
[1] 1 1 -1 0 -1 0 -1 0 -1 0
>meq
[1] 2
>liSM8<-solve.QP(cm8,dv8,Am8,bv8,meq)
>liSM8$solution
[1] 2.500000e-01 2.858899e-53 2.500000e-01 0.000000e+00 2.500000e-01
0.000000e+00 2.500000e-01 1.387779e-17
>cma
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 0 0 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0 0 0
[5,] 0 0 0 0 1 0 0 0 0 0
[6,] 0 0 0 0 0 1 0 0 0 0
[7,] 0 0 0 0 0 0 1 0 0 0
[8,] 0 0 0 0 0 0 0 1 0 0
[9,] 0 0 0 0 0 0 0 0 1 0
[10,] 0 0 0 0 0 0 0 0 0 1
>dva
[1] 0 0 0 0 0 0 0 0 0 0
>Ama
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] 1 1 -1 0 0 0 0 0 0 0 0 0
[2,] -1 1 0 1 0 0 0 0 0 0 0 0
[3,] 1 1 0 0 -1 0 0 0 0 0 0 0
[4,] -1 1 0 0 0 1 0 0 0 0 0 0
[5,] 1 1 0 0 0 0 -1 0 0 0 0 0
[6,] -1 1 0 0 0 0 0 1 0 0 0 0
[7,] 1 1 0 0 0 0 0 0 -1 0 0 0
[8,] -1 1 0 0 0 0 0 0 0 1 0 0
[9,] 1 1 0 0 0 0 0 0 0 0 -1 0
[10,] -1 1 0 0 0 0 0 0 0 0 0 1
>bva
[1] 1 1 -1 0 -1 0 -1 0 -1 0 -1 0
>meq
[1] 2
>liSMa<-solve.QP(cma,dva,Ama,bva,meq)
Error in solve.QP(cma, dva, Ama, bva, meq) :
constraints are inconsistent, no solution!
>cmc
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] 1 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 1 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 1 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 1 0 0 0 0 0
[8,] 0 0 0 0 0 0 0 1 0 0 0 0
[9,] 0 0 0 0 0 0 0 0 1 0 0 0
[10,] 0 0 0 0 0 0 0 0 0 1 0 0
[11,] 0 0 0 0 0 0 0 0 0 0 1 0
[12,] 0 0 0 0 0 0 0 0 0 0 0 1
>dvc
[1] 0 0 0 0 0 0 0 0 0 0 0 0
>Amc
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[,14]
[1,] 1 1 -1 0 0 0 0 0 0 0 0 0 0
0
[2,] -1 1 0 1 0 0 0 0 0 0 0 0 0
0
[3,] 1 1 0 0 -1 0 0 0 0 0 0 0 0
0
[4,] -1 1 0 0 0 1 0 0 0 0 0 0 0
0
[5,] 1 1 0 0 0 0 -1 0 0 0 0 0 0
0
[6,] -1 1 0 0 0 0 0 1 0 0 0 0 0
0
[7,] 1 1 0 0 0 0 0 0 -1 0 0 0 0
0
[8,] -1 1 0 0 0 0 0 0 0 1 0 0 0
0
[9,] 1 1 0 0 0 0 0 0 0 0 -1 0 0
0
[10,] -1 1 0 0 0 0 0 0 0 0 0 1 0
0
[11,] 1 1 0 0 0 0 0 0 0 0 0 0 -1
0
[12,] -1 1 0 0 0 0 0 0 0 0 0 0 0
1
>bvc
[1] 1 1 -1 0 -1 0 -1 0 -1 0 -1 0 -1 0
>meq
[1] 2
>liSMc<-solve.QP(cmc,dvc,Amc,bvc,meq)
>liSMc$solution
[1] 1.666667e-01 3.703906e-18 1.666667e-01 3.703906e-18 1.666667e-01
0.000000e+00 1.666667e-01 3.703906e-18 1.666667e-01 3.703906e-18
1.666667e-01
[12] 2.220988e-17
The problem, of course, is rounding error. A small jitter in the
constraints takes care of it:
>(bvan<-c(bva[1:11],-0.0000000001))
[1] 1e+00 1e+00 -1e+00 0e+00 -1e+00 0e+00 -1e+00 0e+00 -1e+00 0e+00
-1e+00 -1e-10
>liSMan<-solve.QP(cma,dva,Ama,bvan,meq)
>liSMan$solution
[1] 2.000000e-01 0.000000e+00 2.000000e-01 8.002963e-53 2.000000e-01
0.000000e+00 2.000000e-01 0.000000e+00 2.000000e-01 -1.664250e-17
>
When I was testing, I ran into the bad case right away, and it knocked me
for a loop, because I thought I'd made a mistake in my code. My problem is
that I'm trying to write a program that will run significantly larger
problems (hundreds of variables) with a variety of possible constraints
allowed, and it might be difficult to determine whether a problem is
legitimately infeasible or just being kicked out of the feasible region due
to rounding errors. I was wondering whether anyone has any tricks to share
for mitigating these kind of problems and still generating feasible
solutions. Thanks!
-- TMK --
212-460-5430 home
917-656-5351 cell
More information about the R-help
mailing list