bugGNU Octave - Bugs: bug #47314, optim package: quadprog() equality...

 
 

bug #47314: optim package: quadprog() equality constraints have opposite sign

Submitter:  Ray Zimmerman <rdzman>
Submitted:  Tue 01 Mar 2016 03:35:39 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  None
Originator Name:  Ray Zimmerman Open/Closed:  * Closed
Release:  * 4.0.1-rc1 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 27 Apr 2016 06:30:15 PM UTC, comment #22: 

@Olaf:
There you are.

Philip Nienhuis <philipnienhuis>
Group Member
Wed 27 Apr 2016 01:53:48 PM UTC, comment #21: 

Could someone close this bug?

Olaf Till <i7tiol>
Group Member
Tue 29 Mar 2016 01:32:33 PM UTC, comment #20: 

I ran a few tests and all seems fine to me.

Thanks.

Ray Zimmerman <rdzman>
Mon 28 Mar 2016 09:28:25 PM UTC, comment #19: 

Sorry, cut-and-paste error, the cloning command line should be:

hg clone http://hg.code.sf.net/p/octave/optim fixed-optim

Olaf Till <i7tiol>
Group Member
Mon 28 Mar 2016 09:19:28 PM UTC, comment #18: 

Done, with respective code comments and docs change. Could you do:

hg clone http://i7tiol@hg.code.sf.net/p/octave/optim fixed-optim

cd fixed-optim/inst/

and test it? (Yes, I've tested also, but ...)

Anyway, please comment on whether you also think it is as it should be now.

Olaf Till <i7tiol>
Group Member
Mon 28 Mar 2016 08:10:12 PM UTC, comment #17: 

You are right ... the sign of lambda for bounds should be kept non-negative even for too close bounds, as in the other bounds, which are treated as inequality constraints and therefore have non-negative lambdas. So if negative, it should be negated again and therefore assigned to lower bounds instead of upper bounds (the latter would otherwise naturally resemble equality constraints).

This complicates the code changes a bit, I'll tell you when I'm ready with it.

Olaf Till <i7tiol>
Group Member
Mon 28 Mar 2016 06:38:25 PM UTC, comment #16: 

I'm not sure I understand what you mean by "Matlab chooses to assign lambda of too close upper and lower bounds to the lower bound". It would be incorrect to always assign it to lower or always to upper. That depends on the problem.

I assume when the bounds are "too close" they are being replaced by an equality constraint. In that case it is the the sign of the lambda from that constraint that tells you which bound is binding, and therefore which bound to assign the lambda too.

Here is a trivial example, run in Matlab, that illustrates.


>> [x,obj,flag,op,lambda]=quadprog(1,0,[],[],[],[],1,1+eps);

Solution found during presolve.

Some combination of the bounds, linear constraints, and linear terms
in the objective function immediately lead to the solution.

>> lambda

lambda =

    ineqlin: [0x1 double]
      eqlin: [0x1 double]
      lower: 1
      upper: 0

>> [x,obj,flag,op,lambda]=quadprog(1,0,[],[],[],[],-1-eps,-1);

Solution found during presolve.

Some combination of the bounds, linear constraints, and linear terms
in the objective function immediately lead to the solution.

>> lambda

lambda =

    ineqlin: [0x1 double]
      eqlin: [0x1 double]
      lower: 0
      upper: 1.0000


Ray Zimmerman <rdzman>
Mon 28 Mar 2016 06:11:52 PM UTC, comment #15: 

Thanks for the tests and the theory discussion.

So Matlab chooses to assign lambda of too close upper and lower bounds to the lower bound (this is was Asma Afzal told me, but I wanted to be sure). This means that the sign of this lambda has to be opposite to the sign returned (by Matlab and now by us) for lambda in equality constraints. Matlab does this correctly.

As it turns out, this is what we currently do for lambda of too close upper and lower bounds, so I only inserted a further comment and made no further code changes.

This treatment of lambda is now documented in the user docs.

A new release will probably be requested anyway during this week.

(What I still have to do is to notify Asma that I told her something wrong ...)

Thanks again, and if nobody has further comments this bug-report can be closed.

Olaf Till <i7tiol>
Group Member
Mon 28 Mar 2016 06:00:35 PM UTC, comment #14: 

Adding Rays comment here (replying per E-Mail does not reach the bug tracker ...):

On Mon, Mar 28, 2016 at 11:27:43AM -0400, Ray Zimmerman wrote:

> I just posted the results from Matlab.
>
> Another “justification” for the convention chosen by Matlab (and numerous other solvers, such as Octave’s own LP solver, GLPK) is that a positive shadow price on an inequality corresponds to the sensitivity of the objective to relaxing the constraint (i.e. to an increase in b for the standard A*x <= b formulation). Similarly, for an inequality constraint, it is the sensitivity of the objective function to an increase in the righthand side of the equality constraint A * x  = b.
>
> I hadn’t thought about the fact that tight variable bounds might be translated into an equality constraint. But it is true, that in such cases, you’ll need to be careful to apply the absolute value of the lambda to the correct bound.
>
> Thanks,
>
>    Ray

Olaf Till <i7tiol>
Group Member
Sun 27 Mar 2016 02:18:01 PM UTC, comment #13: 

On an auxiliary help page, Matlab indeed implicitly defines lambda for equality constraints, similar to '.5 gradient_of_func + lambda gradient_constraints = 0', i.e. opposite to our current lambda. This has to do with their treating of _in_equlality constraints with '<=', while we use '>='.

Since it's defined, we indeed must negate our lambda for equality constraints. I've applied your changeset, and a further with an explaining comment. The story is not ready, because I need to make our treatment of too close lower and upper bounds consistent with this, and have to wait for some Matlab test requested at the maintainers mailing list.

Olaf Till <i7tiol>
Group Member
Sat 26 Mar 2016 08:09:12 PM UTC, comment #12: 

Forgot to paste 'A' in last example, it was A = [-3, 5];. Also, the second call to quadprog should be with '-A', and in the output 'lambda2' (which has opposite sign).

Thinking further about what is in the last post, this argument for Octaves variant of lambda seems quite strong, reformulated:

- quadprog specifies, in the doc, Aeq x = b for equality constraints (not the -- equivalent -- -Aeq x = -b).

- the usual definition of lambda in equality constraints is that it satisfies the equation ('H' and 'f' as used by the doc): '.5Hx + f = transpose_of_given_constraint_matrix * lambda' (It seems not usual to negate the given constraint matrix, and formally it would be an unnecessary additional step.)

- so someone using lambda should expect (if he doesn't check it) that the matrix is not (implicitely) negated for the returned gamma -- should expect what Octaves quadprog does, not Matlabs.

So I'm in favour of keeping our quadprog's current behaviour, but I still hope for some input.

Olaf Till <i7tiol>
Group Member
Sat 26 Mar 2016 03:17:46 PM UTC, comment #11: 

My memory was wrong, quadprog does not seem to return always positve lambda. But take a look at the example anyway:


T = [3, 1; 1, 4];
H = T * T';
F = [1, -2];
[x, ~, ~, ~, lambda1] = quadprog (H, F, [], [], A, 0)
x =

   0.0081566
   0.0048940

lambda1 =

  scalar structure containing the fields:

    lower = [](0x0)
    upper = [](0x0)
    eqlin = -0.37194
    ineqlin = [](0x0)

[x, ~, ~, ~, lambda1] = quadprog (H, F, [], [], A, 0)
x =

   0.0081566
   0.0048940

lambda1 =

  scalar structure containing the fields:

    lower = [](0x0)
    upper = [](0x0)
    eqlin = -0.37194
    ineqlin = [](0x0)

F' + H * x
ans =

   1.1158
  -1.8597

A' * lambda1.eqlin
ans =

   1.1158
  -1.8597


The last two commands compare the function derivative with constraints derivative times lambda. That they are equal even in sign seems more intuitive (though not specified) for using in an application than Matlabs behaviour (though also correct). (You probably checked the same for Matlab, since you use lambda?)

There is also the aspect that there are other frontends in Octave for that algorythm. I'm still not sure what to do.


Olaf Till <i7tiol>
Group Member
Sat 26 Mar 2016 01:33:50 PM UTC, comment #10: 

I'd been fooled by a certain solver implementation where all lambdas were kept positive, even those of equality constraints.

I see now that you're right, for equality constraints lambdas of both signs are equivalent in some sense.

But since you (Ray) say the sign matters for your application:

As I see it now, what you call convention may not apply to matrix A as a whole -- it would be legitimate to arbitrarily choose the sign of individual rows of A (and corresponding elements of b), and so the sign of individual elements of lambda. If the sign is significant, don't you also have to know whether the corresponding row of A had been negated by the algorithm or not? I haven't deeply studied this particular algorithm, but it seems to always return positive lambdas -- so I'm afraid it implicitly changes the sign of some rows of A. Does that not destroy the usefulness of lambda for you?

Olaf Till <i7tiol>
Group Member
Sat 26 Mar 2016 06:29:46 AM UTC, comment #9: 

Olaf, I can confirm Ray's description of the mathematics.

The mnemonic I use to remember that equality constraints can be either positive or negative whereas for inequalities they can be chosen to be always positive, is that inequalities have an intrinsic "direction" whereas equalities do not.

Lachlan Andrew <lachlan>
Fri 25 Mar 2016 02:18:45 PM UTC, comment #8: 

I don't know why you say that lambda's should be inherently non-negative. That is certainly true for inequality constraints, but if you have a set of equality constraints you can have some multipliers negative and others positive ... and yes, the sign is significant.

For example, in my application, in power systems, these lambdas refer to the nodal prices of energy in a power grid. These prices can be positive at some nodes, while negative at others ... and it matters that I be able to know which ones are which. So, to answer your question ... yes, I do need these lambdas for something and their signs do matter.

For an equality constraint, such as A*x = b, it can be implemented as either ...

 A * x =  b, or
-A * x = -b

Both are correct and yield lambdas of opposite signs. One simply needs to pick a convention. Currently, quadprog in Octave and quadprog in Matlab have chosen opposite conventions. So, currently I have to check whether I'm calling quadprog from Matlab or quadprog from Octave and either swap the sign on the constraints or swap the sign on the lambda to get consistent answers.

In my view, this is purely a question of Matlab compatibility and not a question of correctness. Both are correct. There is no value to having them behave differently. There is value to having them behave identically.

My 2 cents.

Ray Zimmerman <rdzman>
Fri 25 Mar 2016 01:57:49 PM UTC, comment #7: 

I quote myself from

http://lists.gnu.org/archive/html/octave-maintainers/2015-07/msg00250.html

Lambda should be inherently non-negative. Don't change it's sign to get the same value as Matlab. I don't know why they return a negative value.

(end quote)

Please explain if you have a different opinion. ('Matlab makes it different' is no explanation.)

I was never lucky with returning lambda at all. Do you indeed need it for something, or have you just seen that Matlab returns a different value?

Olaf Till <i7tiol>
Group Member
Fri 25 Mar 2016 08:44:35 AM UTC, comment #6: 

Here is a Mercurial changeset for the change.

I'm adding Olaf to the Cc list, as he maintains this package.

(file #36760)

Lachlan Andrew <lachlan>
Fri 25 Mar 2016 04:03:06 AM UTC, comment #5: 

Ray, thanks for the bug report.

The fix looks very simple and is presumably correct, but could you please test it on a few other QPs just to make sure that it is simply the sign convention, rather than something more obscure and sinister?

Thanks!

Lachlan Andrew <lachlan>
Wed 02 Mar 2016 09:08:37 PM UTC, comment #4: 

Yes, in my opinion this is the final fix. Any other approach would require modifying _qp_ which may have compatibility implications for other wrappers. So, yes, I think adding that single negative sign is the correct fix. The only effect should be to make it compatible with Matlab's quadprog(), which I presume is the primary reason for the existence of the quadprog() wrapper in the first place.

Ray Zimmerman <rdzman>
Wed 02 Mar 2016 08:50:32 PM UTC, comment #3: 

So, Is this the final fix or Something can be done to overcome this?

Ankush Khandelwal <ank2206>
Wed 02 Mar 2016 01:21:02 PM UTC, comment #2: 

The sign of the shadow prices (lambdas) on the equality constraints depends on the convention used to implement the constraints. Apparently, the Octave implementation of quadprog() uses the opposite convention from the Matlab version for the multipliers on the equality constraints. This results in an incompatibility with Matlab.

As mentioned in my original bug report, the fix is trivial. It is simply to add this negative sign to line 382 in quadprog.m:


      lambda.eqlin = -qp_lambda(lm_idx:lm_idx + n_eq - count_not_ineq - 1);
                     ^
                     |


Ray Zimmerman <rdzman>
Wed 02 Mar 2016 07:34:18 AM UTC, comment #1: 

Can you elaborate about what needs to be fix in this issue? I am willing to work on this.

Ankush Khandelwal <ank2206>
Tue 01 Mar 2016 03:35:39 PM UTC, original submission:  

quadprog() from optim-1.5.0 implements the equality constraints with the opposite sign from Matlab. That is, the multipliers on the equality constraints have the opposite sign.

I have attached an example that gives the following respective output:

for Octave:

ans =  2.1933

for Matlab:

ans =

   -2.1933

The fix is to simply add a negative sign in line 382:


      lambda.eqlin = -qp_lambda(lm_idx:lm_idx + n_eq - count_not_ineq - 1);


Ray Zimmerman <rdzman>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #36521:  quadprogtest.m added by rdzman (456B - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by i7tiol (Posted a comment)
  • -email is unavailable- added by lachlan
  • -email is unavailable- added by lachlan (Posted a comment)
  • -email is unavailable- added by ank2206 (Posted a comment)
  • -email is unavailable- added by rdzman (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-02-26 mtmiller Carbon-CopyRemoved 80942 -
    2016-04-27 philipnienhuis StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2016-03-25 lachlan Attached File- Added bug_47314_quadprog_eq_sign.cset, #36760
        Carbon-Copy- Added i7tiol
    2016-03-03 mtmiller StatusNone Patch Submitted
        Operating SystemMac OS Any
        Summaryquadprog() equality constraints have opposite sign optim package: quadprog() equality constraints have opposite sign
    2016-03-01 rdzman Attached File- Added quadprogtest.m, #36521

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code