I implement a Sequential quadratic programming (SQP) optimizer and use ojAlgo for the quadratic programming (QP) subproblem.
My question is: How do I get hold of the "Lagrange multipliers" for the QP solution?
In the attached example code that solve an QP result.getMultipliers() only return an empty Optional.
package com.mycompany.testojalgo;
import java.math.BigDecimal;
import java.util.ArrayList;
import java.util.List;
import java.util.Optional;
import org.ojalgo.matrix.Primitive64Matrix;
import org.ojalgo.optimisation.Expression;
import org.ojalgo.optimisation.ExpressionsBasedModel;
import org.ojalgo.optimisation.Optimisation;
import org.ojalgo.optimisation.Variable;
import org.ojalgo.structure.Access1D;
import org.ojalgo.type.StandardType;
import org.ojalgo.type.context.NumberContext;
public class ojAlgoQP {
public static void main(String[] args) {
testOjAlgoQuadraticProgramming();
}
public static void testOjAlgoQuadraticProgramming() {
// 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
Primitive64Matrix g = Primitive64Matrix.FACTORY.rows(new double[][]{
{-8}, {-3}, {-3}
});
//objectiveHessian
Primitive64Matrix H = Primitive64Matrix.FACTORY.rows(new double[][]{
{6, 2, 1},
{2, 5, 2},
{1, 2, 4}
});
Variable x1 = new Variable("x1");
Variable x2 = new Variable("x2");
Variable x3 = new Variable("x3");
// constraint equations
Primitive64Matrix A = Primitive64Matrix.FACTORY.rows(new double[][]{
{1, 0, 1},
{0, 1, 1}
});
// required constraint values
Primitive64Matrix b = Primitive64Matrix.FACTORY.rows(new double[][]{
{3}, {0}
});
List<Variable> variables = new ArrayList<>();
variables.add(x1);
variables.add(x2);
variables.add(x3);
ExpressionsBasedModel model = new ExpressionsBasedModel(variables);
Expression energy = model.addExpression("Energy");
energy.setLinearFactors(variables, g);
//divide by two to express function using hessian.
energy.setQuadraticFactors(variables, H.divide(2));
energy.weight(BigDecimal.ONE);
//create constraint equations
for (int i = 0; i < A.countRows(); i++) {
Expression expression = model.addExpression("Constraint#"+i);
for (int j = 0; j < A.countColumns(); j++) {
expression.set(variables.get(j), A.get(i, j));
}
expression.level(b.get(i));
}
Optimisation.Result result = model.minimise();
NumberContext accuracy = StandardType.PERCENT.withPrecision(1);
boolean ok = model.validate(result, accuracy);
Optimisation.State v = result.getState();
// How do I get the multipliers
Optional<Access1D<?>> multipliers = result.getMultipliers();
double value1 = result.getValue();
// Get result and check value and constraint
Primitive64Matrix x = Primitive64Matrix.FACTORY.rows(new double[][]{
{x1.getValue().doubleValue()}, {x2.getValue().doubleValue()}, {x3.getValue().doubleValue()}
});
//divide by two to express function using hessian, again.
Primitive64Matrix value = x.transpose().multiply(H.divide(2)).multiply(x).add(x.transpose().multiply(g));
Primitive64Matrix residual= A.multiply(x).subtract(b);
}
}
Update 1: Here is my reworked example using org.ojalgo.optimisation.convex.ConvexSolver.getBuilder();
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 gStore = Primitive64Store.FACTORY.rows(new double[][]{
{-8}, {-3}, {-3}
});
//objectiveHessian
Primitive64Store HStore = Primitive64Store.FACTORY.rows(new double[][]{
{6, 2, 1},
{2, 5, 2},
{1, 2, 4}
});
// constraint equations
Primitive64Store AStore = Primitive64Store.FACTORY.rows(new double[][]{
{1, 0, 1},
{0, 1, 1}
});
// required constraint values
Primitive64Store bStore = Primitive64Store.FACTORY.rows(new double[][]{
{3}, {0}
});
ConvexSolver.Builder builder = ConvexSolver.getBuilder();
builder.equalities(AStore, bStore);
builder.objective(HStore, gStore.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(HStore.multiply(0.5)).multiply(x).add(x.transpose().multiply(gStore));
MatrixStore<Double> residual = AStore.multiply(x).subtract(bStore);
}
}
I believe that is an Optional because it was (sometimes) too messy to map the Lagrange multipliers from the solver to the constraints of the model.
If you're implementing an SQP solver may I suggest that you don't implement it in terms of ExpressionsBasedModel, but delegate to the convex solvers directly. Build something that implements org.ojalgo.optimisation.Optimisation.Solver and delegate to the various classes in the org.ojalgo.optimisation.convex package. Then you code more directly with the matrices, vectors and multipliers.
To make that solver usable by ExpressionsBasedModel you also implement an org.ojalgo.optimisation.Optimisation.Integration and register that by calling ExpressionsBasedModel.addPreferredSolver(myIntegeration) or ExpressionsBasedModel.addFallbackSolver(myIntegeration).
Implementing a solver and making it usable from the modelling tool are two separate things.