optimatika / ojAlgo

oj! Algorithms
http://ojalgo.org
MIT License
459 stars 207 forks source link

Dual variables for convex optimization solution not returned #280

Closed Programmer-Magnus closed 3 years ago

Programmer-Magnus commented 3 years ago

I need the Lagrangian multipliers for a Quadratic Programming solution but they are not returned by result.getMultipliers() in this example using version 48.2.0 of ojAlgo:

package com.mycompany.testojalgo;

import java.util.Optional;
import org.ojalgo.matrix.store.MatrixStore;
import org.ojalgo.matrix.store.Primitive64Store;
import org.ojalgo.optimisation.Optimisation;
import org.ojalgo.optimisation.convex.ConvexSolver;
import org.ojalgo.structure.Access1D;

public class ojAlgoQP {

   public static void main(String[] args) {
      testOjAlgoQuadraticProgramming2();
   }

   public static void testOjAlgoQuadraticProgramming2() {
//  QP Example 16.2 p453 in 'Numerical Optimization', 2ed, (2006), Jorge Nocedal and Stephen J. Wright.
//  minimize function F(x1,x2,x3) = 3*x1*x1 + 2*x1*x2 + x1*x3 + 2.5*x2*x2 + 2*x2*x3 + 2*x3*x3 - 8*x1 - 3*x2 - 3*x3
//  x = [x1, x2, x3]'
//  F(x) = 1/2*x'*H*x + x'*g
//  constraints x1 + x3 = 3, x2 + x3 = 0
//  A*x = b

//objectiveGradient
      Primitive64Store g = Primitive64Store.FACTORY.rows(new double[][]{
         {-8}, {-3}, {-3}
      });
//objectiveHessian
      Primitive64Store H = Primitive64Store.FACTORY.rows(new double[][]{
         {6, 2, 1},
         {2, 5, 2},
         {1, 2, 4}
      });
// constraint equations
      Primitive64Store A = Primitive64Store.FACTORY.rows(new double[][]{
         {1, 0, 1},
         {0, 1, 1}
      });
// required constraint values
      Primitive64Store b = Primitive64Store.FACTORY.rows(new double[][]{
         {3}, {0}
      });
      ConvexSolver.Builder builder = ConvexSolver.getBuilder();
      builder.equalities(A, b);
      builder.objective(H, g.negate());
      ConvexSolver solver = builder.build();
      Optimisation.Result result = solver.solve();

// How do I get the multipliers?  multipliers = Optional.empty
      Optional<Access1D<?>> multipliers = result.getMultipliers();
// value1 = -3.5
      double value1 = result.getValue();

// Verify result:
// x= [2.0, -0.9999999999999996, 0.9999999999999997]';
// value = -3.5
// residual =[-4.440892098500626E-16, 1.1102230246251565E-16]'
      Primitive64Store x = Primitive64Store.FACTORY.column(result.toRawCopy1D());
      MatrixStore<Double> value = x.transpose().multiply(H.multiply(0.5)).multiply(x).add(x.transpose().multiply(g));
      MatrixStore<Double> residual = A.multiply(x).subtract(b);

   }
}
Programmer-Magnus commented 3 years ago

This issue refere to https://stackoverflow.com/questions/63951773/how-to-get-multipliers-after-solving-a-quadratic-program-in-ojalgo

apete commented 3 years ago

Try with this 1ca62fa83ee9a41a10c82b4ea55cd22b485008cd from the https://github.com/optimatika/ojAlgo/tree/ConvexLagrangeMultipliers branch.

Can you verify that it works? Perhaps contribute a couple of unit tests.

It's best if you first fork ojAlgo to your own account. Then you clone that to a local copy on your development machine, and import that project/module to your IDE.

Programmer-Magnus commented 3 years ago

I might have made my first git and github pull request to contribute a test for multipliers. My experiments say multipliers have the wrong sign. I look forward to your comments.

apete commented 3 years ago

I see that you forked ojAlgo and created a branch test1, but I do not see any commits nor a pull request.

If you're new to git (GitHub) it takes a while to get used to. Perhaps try the GitHub Desktop GUI. It hides a bit of the complexity.

Are the multipliers the wrong sign, are you sure? With or without inequality constraints different solvers are used. Are the multipliers the wrong sign in both cases?

Programmer-Magnus commented 3 years ago

I made a second attempt. There should now be some actual code in the pull request.

About the sign of the multipliers, I am not 100% sure. But I added a calculation in the end of the test that convinced me they have a minus to much. Please check the expressions and tell me if my assumptions are wrong.

I assume this is what is solved:

| [Q] -[A]^T| |x| = |-c|
| [A]   [0] | |L|   |b|

And check result of [L] using: [Q]*[x] - [A]^T*[L] = -c

also check [A]*[x] = b

apete commented 3 years ago

Looking at the code this is what is solved:

| [Q]  [A]^T| |x| = |c|
| [A]  [0] |  |L|   |b|

In the ConstrainedSolver class look at the getIterationKKT() and getIterationRHS() methods.

Remember the solver is designed to have c negated (so that negating it doesn't have to be part of the algorithm).

Are the Lagrange multipliers given in that text book example?

Programmer-Magnus commented 3 years ago

Yes multipliers are given in the book at p453. The matrix equation I wrote (16.4) is on p451. In the test I wrote Q and C are the Hessian and Gradient of the function at x=[0.0, 0.0, 0.0]'

apete commented 3 years ago

If we leave the ojAlgo solvers out for a while, and just take the numbers from that Nocedal and Wright example and form the KKT equation system

| [Q]  [A]^T| |x| = |-c|
| [A]  [0] |  |L|   |b|

Here's a small pdf using the same example that outlines how to solve using KKT: http://www.lendek.net/teaching/opt_en/qp.pdf

I do that and get the same solution as ojAlgo does.

Programmer-Magnus commented 3 years ago

I believe its the minus sign in front of [A]^T in the KKT equation system. It affect the sign of |L|, and thereby the interpretation of a negative/positive value and nothing else as I understand it. I am fine with any choice. And believe its your call. I update the test to make them pass with the present implementation. In the text book by Nocedal and Wright there is a minus in-front of the [A]^T, i do not find their reason for it.

apete commented 3 years ago

I've been reading up on the subject this morning (it's been a while). I believe for equality constraints you have a choice on how to define the augmented Lagrangian, is that true? With inequalities it's well defined...

If there's a problem I want to fix it, but if it's a matter of choice/definition then it stays the way it is.

Programmer-Magnus commented 3 years ago

I think it would be confusing if multiplier sign had different meaning for equality vs inequality constraints. I would like to avoid that for ojAlgo users. Lagrangian function L(x, lambda) = F(x) - lambdaG(x) could as well be written L(x, lambda) = F(x) + (-lambda)G(x). Convincing me that the sign is not important, but I would expect the same for inequality constraints. What is the reasoning there?

apete commented 3 years ago

I constructed another test based on an example found here: https://people.duke.edu/~hpgavin/cee201/LagrangeMultipliers.pdf

I pushed to your PR. Please verify that I didn't make any mistake.

The ojAlgo solvers mimic the results of that example using unconstrained, equality constrained and inequality constrained problem formulations (which results in different solvers being used). The results are consistent.

Haven't done this (now) but something like this:

  1. Split each equality constraint in to 2 inequalities (resulting in 2 Lagrange multipliers)
  2. Using those inequalities formulate the Lagrangian as a relaxation of the original problem (that means there are restrictions on the signs)
  3. ...
Programmer-Magnus commented 3 years ago

I checked both test codes in LagrangeTest.java and the only mistake I found was a '-' that should be a '+' in a comment. I made a commit for that, do i also need to make push request for that, or can you already see it? Have a nice week end!