Thursday 23 July 2015

Week 8 and 9: The ordering of lambda in quadprog

I was trying to dig through __qp__.cc to figure out why the ordering of lambda in [1] does not match that of quadprog in Matlab.
The example shows how the values were different:

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 = 1*ones(4,1);
H = C ' * C;

f = -C ' * d;

[x, obj_qp, INFO, lambda] = qp ([],H,f,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


Reordering lambda based on the length of constraints resulted in
lambda.eqlin = 
   0.01654
lambda.lower =
   0.06743
   0.00000
   0.24993
   0.00000

lambda.upper =
   0.00000
   0.00000
   0.00000
   0.00000

lambda.ineqlin =
   0.00000
   0.49819
   0.00000


Matlab however,  gave lambda.eqlin= -0.0165 for this example and lambda.lower =
    0.0674
    0.2499
         0
         0

There were two issues with the ordering:
  1.  the sign for Lagrange multipliers corresponding to linear equality constraint is always different from Matlab's
  2. The multipliers corresponding to the bound constraints as underlined above are swapped.
I tried several different examples to understand what is going on. For all the examples, the sign for lambda.eqlin was consistently different. Although, I still can't pinpoint why but for now I am just multiplying lambda.eqlin by -1.

For the swapping issue, I tried the same example with just the bound constraints:

[x, obj_qp, INFO, lambda] = qp ([],H,f,[],[],lb,ub)
lambda= 
   0.07377
   0.00000
   0.29791
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000

but only considering lower bound constraints gave:
[x, obj_qp, INFO, lambda] = qp ([],H,f,[],[],lb,[])
lambda=
   0.07377
   0.29791
   0.00000
   0.00000


which is how it is supposed to be. Tracing back, I found out that the ordering of lambda vector in qp.m [2] is not [equality constraints; lower bounds; upper bounds; other inequality constraints] like I previously assumed. From lines 287 and 288 in qp.m [2], the bound constraints are added to the inequality constraint matrix alternatively. So the issue wasn't swapping but understanding how the constraints are passed to __qp__.

In my code in [1], I had to make significant changes in the original qp.m code such as:
- Inequality constraint matrix has the order: [Linear inequality constraints; lower bounds; upper bounds]
- Check for too close bounds forming an equality constraint- This brings indexing issues as now the Lagrange multiplier value corresponding to bounds is in place of multipliers corresponding to linear equality constraints.
Also, Matlab only accepts too close bounds when using a medium scale algorithm and since the lower bound is approximately equal to the upper bound and it is considered as a single equality constraint, the single Lagrange multiplier is placed in the corresponding lambda.upper field while the corresponding lambda.lower value is zero.

Continuing the above example:

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


[x, obj_qp, INFO, lambda] = qp ([],H,f,[],[],lb,ub)
lambda =

  -0.04828
   0.04154
   0.00000
   0.28615
   0.00000
   0.00000
   0.01546


Here, lb(4) = ub(4) and hence the constraint is treated as an equality constraint so the value for corresponding Lagrange multipliers is present on top (underlined)
I added checks for such cases and now my code in [1] gives same results as Matlab:

[x,obj,flag,op,lambda]=quadprog(H,f,[],[],[],[],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


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

Plans for the next weeks:
- Get feedback from my mentor on the changes in quadprog.
- Begin intial work on optimoptions.

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

No comments:

Post a Comment