Using the Levenberg-Marquardt to solve a set of non-linear equations

[profhsgill]

I would like to use the Levenberg-Marquardt to solve a set of non linear equations.

As an example:

x − y = −1
y = x*x + 1
Solutions are (x=1, y=2) and (x=0, y=1)

The above can be written as two functions to be minimized

f(1) = x - y + 1 = 0
f(2) = x*x - y + 1 = 0

I have read the source codes:
LevenbergMarquardtMinimizer.cs
ObjectiveFunction.cs

And I have also searched the tests files, however, I am not sure how to formulate the code to use the minimizer.

While the above equations are analytical and their Jacobian can be easily computed, my final goal is to simply use the objective function with a set of guess values to find the solution of the equations.

Any help on this will be appreciated.

[diluculo]

[EDIT] I've revised the comment to express the Levenberg-Marquardt method and its underlying assumptions more clearly.

Rosenbrock_LM_Der() or Rosenbrock_LM_Dif() can be of assistance. As you know, the Levenberg-Marquardt method seeks to determine parameters that minimize the squared difference between the model and the measured values.

For your optimization, assume that the measured values, represented by the Y vectors, are zeros, and consider the following model values:

// model: 
//       f(x; a, b) = f1^2 + f2^2 = (a - b + 1)^2 + (a^2 - b + 1)^2
// derivatives:
//       df/da = 4 a^3 + a (6 - 4 b) - 2 b + 2
//       df/db = -2 (a^2 + a - 2 b + 2)

Please note that the results may vary depending on the initial guess.

[profhsgill]

Thank you.

[diluculo]

Sample code:

// model: 
//       f(x; a, b) = (a - b + 1)^2 + (a^2 - b + 1)^2
// derivatives:
//       df/da = 4 a^3 + a(6 - 4 b) - 2 b + 2
//       df/db = -2 (a^2 + a - 2 b + 2)

Vector<double> Model(Vector<double> p, Vector<double> x)
{
    var y = CreateVector.Dense<double>(x.Count);
    for (var i = 0; i < x.Count; i++)
    {
        y[i] = (p[0] - p[1] + 1) * (p[0] - p[1] + 1) + (p[0] * p[0] - p[1] + 1) * (p[0] * p[0] - p[1] + 1);
        //y[i] = (p[0] - p[1] + 1) + (p[0] * p[0] - p[1] + 1) ;
    }
    return y;
}

Matrix<double> Prime(Vector<double> p, Vector<double> x)
{
    var prime = Matrix<double>.Build.Dense(x.Count, p.Count);
    for (var i = 0; i < x.Count; i++)
    {
        prime[i, 0] = 4 * p[0] * p[0] * p[0] + p[0] * (6 - 4 * p[1]) - 2 * p[1] + 2;
        prime[i, 1] = -2 * (p[0] * p[0] + p[0] - 2 * p[1] + 2);
    }
    return prime;
}

var X = Vector<double>.Build.Dense(2); // any positive size is available
var Y = Vector<double>.Build.Dense(2);
var obj = ObjectiveFunction.NonlinearModel(Model, Prime, X, Y);

var start1 = Vector<double>.Build.DenseOfArray(new[] { 0.0, 0.0 });            
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, start1);
// Minimizing Points: [ 0.0, 1.0 ]

var start2 = Vector<double>.Build.DenseOfArray(new[] { 1.0, 1.0 });
result = solver.FindMinimum(obj, start2);
// Minimizing points: [ 1.0, 2.0 ]

[profhsgill]

This is awesome, thank you very much. This helps a lot.

[profhsgill]

As a follow up, I was able to:

1. Not square the equations
2. Keep them separate
3. Use Numerical Jacobian instead of analytically computing it

Here is the full code

using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.Optimization;
using MathNet.Numerics.Differentiation;
using MathNet.Numerics.LinearAlgebra.Double;

// This is an example from https://www.mathworks.com/help/optim/ug/fsolve.html
//
// Equations:
//
// x − y = −1
// y = x*x + 1
// Solutions are (x=1, y = 2) and (x=0, y = 1)
// 

//===================================================================
// USING LevenbergMarquardtMinimizer in MathNet.Numerics.Optimization
//===================================================================
// Does NOT square the equations/functions to be solved
//===================================================================
// Uses Numerical Jacobian 
//===================================================================

// Define x, y data set for model
// x data is immaterial. It can be a vector of any random numbers
var X = Vector<double>.Build.Dense(2); // By default fills with 0 
//var X = Vector<double>.Build.Dense(2, 11.123); // For example fills with 11.123 => ALSO WORKS

// We want y = 0 since this is where the 'roots' are
var Y = Vector<double>.Build.Dense(2); // By default fills with 0, and we want this to be 0

// Define the equations to solve
double[] Equations(double[] x)
{
    var y = new double[2];
    y[0] = x[0] - x[1] + 1; // Write equations/functions as f(x) = 0
    y[1] = x[0] * x[0] - x[1] + 1;
    return y;
}


// Create Function of type that can be called by MathNet
Vector<double> FuncToSolve(Vector<double> p)
{
    var p1 = p.ToArray();
    var pComputed = Equations(p1);
    return CreateVector.Dense<double>(pComputed);
}

// Define the Model
Vector<double> Model(Vector<double> p, Vector<double> x)
{
    var y = FuncToSolve(p);
    return y;
}

// Define the Numerical form of Jacobian
// Reference: https://github.com/mathnet/mathnet-numerics/blob/master/src/Numerics.Tests/DifferentiationTests/NumericalJacobianTests.cs
// I would like to not have to rewrite the equations
// I have not figured out how to do that yet
// I defined a separate 'Equations' function to use in creating the 'Model' via FuncToSolve
// Perhaps 'Equations' could be used below as well rather than having to rewrite them

Matrix<double> Jacobian(Vector<double> p, Vector<double> x)
{
    // Here a delegate is being used on lambda function to create the type required by 'Evaluate'
    Func<double[], double>[] fJacobian =
        {
            (p) => (p[0] - p[1] + 1),
            (p) => (p[0] * p[0] - p[1] + 1)
        };
    var jacobian = new NumericalJacobian();
    var jeval1 = jacobian.Evaluate(fJacobian, p.ToArray());
    return DenseMatrix.OfArray(jeval1);
}

// Use LevenbergMarquardt 
var obj = ObjectiveFunction.NonlinearModel(Model, Jacobian, X, Y);
var guessValue1 = Vector<double>.Build.DenseOfArray(new[] { 0.0, 0.0 }); 
var solver = new LevenbergMarquardtMinimizer();
var results = solver.FindMinimum(obj, guessValue1);
Console.WriteLine("Solution with LevenbergMarquardtMinimizer is \n" + results.MinimizingPoint);
// Gives solution
// x = -5.01042E-09 , y = 1 
// Expected solution:
// x = 0, y = 1

var guessValue2 = Vector<double>.Build.DenseOfArray(new[] { 1.0, 1.0 }); 
var solver2 = new LevenbergMarquardtMinimizer();
var results2 = solver2.FindMinimum(obj, guessValue2);
Console.WriteLine("Solution with LevenbergMarquardtMinimizer is \n" + results2.MinimizingPoint);
// Gives solution
// x = 1 , y = 2
// Expected solution:
// x = 1, y = 2

[diluculo]

@profhsgill, great work! Your model is more accurate and straightforward.

1. Analytical Jacobian can be implemented as follows:

// model: 
//    f(x; a, b) = [ f(x0; a, b)   f(x1; a, b) ]
//               = [   a - b + 1   a^2 - b + 1 ]
// derivatives:
//    ∂f/∂a = [  1   2*a ]
//    ∂f/∂b = [ -1    -1 ]
//    i.e., Jacobian = [ ∂f(x0)/∂a  ∂f(x0)/∂b ] = [   1  -1 ]
//                     [ ∂f(x1)/∂a  ∂f(x1)/∂b ]   [ 2*a  -1 ]

2. There's no need for an additional numerical Jacobian function since the NonlinearModel handles the numerical calculation.
   Here, the calculation of f'(x) is performed with an accuracy order of 2 using the expression {f(x + h) - f(x - h)} / 2h + O(h^2), where h represents the step size.

var obj = ObjectiveFunction.NonlinearModel(Model, X, Y, accuracyOrder: 2);

[profhsgill]

That's good to know that the Jacobian is builtin. It saves many lines of code if we simply use :

var obj = ObjectiveFunction.NonlinearModel(Model, X, Y, accuracyOrder: 2);

@diluculo Thank you for your help.

I checked and the same code also works if we simply replace

var solver = new LevenbergMarquardtMinimizer();

with MathNet.Numerics.Optimization.TrustRegion

TrustRegionNewtonCGMinimizer

var solver = new TrustRegionNewtonCGMinimizer();

or TrustRegionDogLegMinimizer

var solver = new TrustRegionDogLegMinimizer();