Re: RFC: quadprog/lsqlin with __qp__ (Re: GSoC 2015, optim)

[Asma Afzal]

Hi Olaf,

> The updated code still accepts (H, f, A, [], ...) and (H, f, [], A, ...).

I have changed it.

> > The values are different from Matlab's, but I was more concerned about
> > arranging multipliers from qp in the right order. Running this example
> > with qp I get:
> >
> > [x, obj_qp, INFO, lambda]  = qp ([], C'*C,-C'*d,Aeq,beq,lb,ub,[],A,b)
> >
> > lambda =
> >
> >    0.01654
> >    0.06743
> >    0.00000
> >    0.24993
> >    0.00000
> >    0.00000
> >    0.00000
> >    0.00000
> >    0.00000
> >    0.00000
> >    0.49819
> >    0.00000
> >
> >  I am not sure why the values are swapped (third and fourth values)
>
> I'd rather say the 2nd and 3rd values are swapped.
>
> > and the sign is different for lambda.eqlin.
> >  lambda.eqlin in matlab is -0.0165 for this example and lambda.lower =
> >     0.0674
> >     0.2499
> >          0
> >          0
> >
> > Can you explain why this is so?
>
> I don't think I've currently the time to dig into __qp__.cc to explain
> the swapping. But it's an indication that nobody ever used the lambda
> returned by __qp__.cc.

I spent some time to dig into __qp__.cc, but it turns out that the
bound constraints are places alternatively in qp.m (lines 287-288)
[1]. I figured this out when I removed all constraints in this example
and just used lower bounds:

C = [0.9501    0.7620    0.6153    0.4057
    0.2311    0.4564    0.7919    0.9354
    0.6068    0.0185    0.9218    0.9169
    0.4859    0.8214    0.7382    0.4102
    0.8912    0.4447    0.1762    0.8936];
d = [0.0578
    0.3528
    0.8131
    0.0098
    0.1388];
A =[0.2027    0.2721    0.7467    0.4659
    0.1987    0.1988    0.4450    0.4186
    0.6037    0.0152    0.9318    0.8462];
b =[0.5251
    0.2026
    0.6721];
Aeq = [3 5 7 9];
beq = 4;
lb = -0.1*ones(4,1);
ub = ones(4,1);

[x,obj,flag,op,lambda]=quadprog(C'*C,-C'*d,[],[],[],[],lb,ub)

lambda =

   0.07377
   0.00000
   0.29791
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000

[x, obj_qp, INFO, lambda] = qp ([],C'*C,-C'*d,[],[],lb,[])

lambda=

   0.07377
   0.29791
   0.00000
   0.00000


So I changed this in the updated code in [2].  For the sign of
lambda.eqlin, it is always the case for all the examples I tried so
I'm just using -1* lambda.eqlin for now.

Few more things that affect the ordering of lambda:

- Too close bound constraints:
qp treats it as an equality constraint and hence the resulting
Lagrange multiplier is placed next to the multipliers corresponding to
the equality constraints in lambda generated by qp.
In Matlab, the resulting multiplier is placed in the lambda.upper
field while a zero is placed in the respective position in
lambda.lower.
Continuing the same example:

lb(4) = 0.3;
ub = 0.3*ones(4,1);

[x, obj_qp, INFO, lambda] = qp ([],C'*C,-C'*d,[],[],lb,ub)

lambda =

  -0.04828
   0.04154
   0.00000
   0.28615
   0.00000
   0.00000
   0.01546

This will be reordered (as in Matlab):

[x,obj,flag,op,lambda]=quadprog(C'*C,-C'*d,[],[],[],[],lb,ub)

lambda =
  scalar structure containing the fields:
    lower =
       0.04154
       0.28615
       0.00000
       0.00000
    upper =

       0.000000
       0.000000
       0.015456
       0.048279

- -Inf bound constraints:
qp.m strips off the -Inf constraints before passing them to __qp__. I
am doing the same in quadprog. I have added further checks to make
sure the multipliers are placed in the right positions in their
respective fields.

Kind Regards,
Asma.

[1] http://fossies.org/dox/octave-4.0.0/qp_8m_source.html
[2] https://github.com/AsmaAfzal/octave_workspace/blob/master/quadprog.m