Nonlinearconstraint in AugmentedLagrangian

[scott_g...@trimble.com]

Firstly, thank you for providing this fantastic library.

I've been experimenting with AugmentedLagrangian. I started with a underdetermined linear system, A * x = b that could be solved using A.Solve(b). In the underdetermined case I believe .Solve will minimise the 2-norm of x subject to A * x = b.

I wanted to extended this to minimise a function subject to a system of linear constraints using the method of AugmentedLagrangian. For the function I chose the two norm of x, so that I could compare it to .Solve. The AugmentedLagrangian code works and satisfies the constraints to 1e-8. When I try to tighten the tolerance on the constraint, constraintsTolerance below, it doesn't seem to have any affect.

Is this the correct way to tighten the tolerance on the constraint? Will there be support for AugmentedLagrangian to take LinearConstraints?

// TEST
int nVariablesTest = 4; // number of variables
int nConstraintsTest = 2; // number of constraints
double constraintsTolerance = 1e-100;
double[,] ATest = new double[,] { { 1, 2, 3, 4 }, { 0, 4, 3, 1 } }; // arbitary A matrix. A*X = b
double[,] bTest = new double[,] { { 0 }, { 2 } }; // arbitary A matrix. A*X = b

double[,] XSolve = ATest.Solve(bTest); // uses the pseudoinverse to minimise norm(X) subject to A*X = b

// recreate Solve function using AugmentedLagrangian
var fTest = new NonlinearObjectiveFunction(nVariablesTest, ds => ds.InnerProduct(ds), ds => ds.Multiply(2.0)); // minimise norm(ds)

var nonlinearConstraintsTest = new List<NonlinearConstraint>(nConstraintsTest); // linear constraints A*X = b
for (int i = 0; i < nConstraintsTest; i++)
{
int j = i; // http://blogs.msdn.com/b/ericlippert/archive/2009/11/12/closing-over-the-loop-variable-considered-harmful.aspx
nonlinearConstraintsTest.Add(new NonlinearConstraint(fTest, ds => ATest.GetRow(j).InnerProduct(ds) - (double)bTest.GetValue(j, 0), ConstraintType.EqualTo, 0.0, ds => ATest.GetRow(j), constraintsTolerance));
}

var innerSolverTest = new BroydenFletcherGoldfarbShanno(nVariablesTest);
var solverTest = new Accord.Math.Optimization.AugmentedLagrangian(innerSolverTest, fTest, nonlinearConstraintsTest);
bool didMinimise = solverTest.Minimize();

var errorConstraintRelative = XSolve.Subtract(solverTest.Solution, 1).ElementwiseDivide(XSolve); // relative error between .Solve and .Minimize
var errorConstraintAbsolute = XSolve.Subtract(solverTest.Solution, 1); // absolute error between .Solve and .Minimize

double[] errorConstraintsTest = new double[nConstraintsTest];
for (int i = 0; i < nConstraintsTest; i++)
{
errorConstraintsTest[i] = nonlinearConstraintsTest[i].Function(solverTest.Solution);
}

[César]

Hi there!

I have been trying to improve the accuracy of the specified constraints, but I still didn't have success with it. I am starting to believe this might be a limitation of the Augmented-Lagragian method itself and the finite precision of floating-point numbers. I am wondering if you wouldn't get better results by applying some transformation to your problem, such as transforming it to the log-domain. If your function is indeed only a inner product, this might be feasible and it might indeed help. I am taking a look on it!

I am also registering your code in an unit test so I can continue checking if it can be improved. Please let me know if you also find something!

Best regards,

Cesar

[scott_g...@trimble.com]

Hey Cesar,

Thanks for looking into this. I'm currently unable to debug Accord.net because I'm using MSVS2012 and need to change the toolset to v120 somehow. I'll try and get that working soon.

The stopping condition would be the first place that I would look. I think bool feasible should be false in the case above because the extremely tight tolerance hasn't been met, line 597 of AugmentedLagrangian.cs. So the algorithm must be stopping because ICM == 0 or noProgressCounter > maxCount, line 668 and line 673.

My function isn't an inner product. I was just using that as a test so I could compare the AugmentedLagrangian solution to the .Solve solution.