LevenbergMarquardt Class

public class LevenbergMarquardt : BaseLeastSquaresMethod, 
	ILeastSquaresMethod, IConvergenceLearning

Public Class LevenbergMarquardt
	Inherits BaseLeastSquaresMethod
	Implements ILeastSquaresMethod, IConvergenceLearning

// Suppose we would like to map the continuous values in the
// second column to the integer values in the first column.
double[,] data =
{
    { -40,    -21142.1111111111 },
    { -30,    -21330.1111111111 },
    { -20,    -12036.1111111111 },
    { -10,      7255.3888888889 },
    {   0,     32474.8888888889 },
    {  10,     32474.8888888889 },
    {  20,      9060.8888888889 },
    {  30,    -11628.1111111111 },
    {  40,    -15129.6111111111 },
};

// Extract inputs and outputs
double[][] inputs = data.GetColumn(0).ToJagged();
double[] outputs = data.GetColumn(1);

// Create a Nonlinear regression using 
var nls = new NonlinearLeastSquares()
{
    NumberOfParameters = 3,

    // Initialize to some random values
    StartValues = new[] { 4.2, 0.3, 1 },

    // Let's assume a quadratic model function: ax² + bx + c
    Function = (w, x) => w[0] * x[0] * x[0] + w[1] * x[0] + w[2],

    // Derivative in respect to the weights:
    Gradient = (w, x, r) =>
    {
        r[0] = w[0]* w[0]; // w.r.t a: a²  // https://www.wolframalpha.com/input/?i=diff+ax²+%2B+bx+%2B+c+w.r.t.+a
        r[1] = w[0];       // w.r.t b: b   // https://www.wolframalpha.com/input/?i=diff+ax²+%2B+bx+%2B+c+w.r.t.+b
        r[2] = 1;          // w.r.t c: 1   // https://www.wolframalpha.com/input/?i=diff+ax²+%2B+bx+%2B+c+w.r.t.+c
    },

    Algorithm = new LevenbergMarquardt()
    {
        MaxIterations = 100,
        Tolerance = 0
    }
};


var regression = nls.Learn(inputs, outputs);

// Use the function to compute the input values
double[] predict = regression.Transform(inputs);

' Suppose we would Like to map the continuous values in the
' second column to the integer values in the first column.
Dim data(,) As Double =
{
    {-40, -21142.1111111111},
    {-30, -21330.1111111111},
    {-20, -12036.1111111111},
    {-10, 7255.3888888889},
    {0, 32474.8888888889},
    {10, 32474.8888888889},
    {20, 9060.8888888889},
    {30, -11628.1111111111},
    {40, -15129.6111111111}
}

' Extract inputs And outputs
Dim inputs As Double()() = data.GetColumn(0).ToJagged()
Dim outputs As Double() = data.GetColumn(1)

' Create a Nonlinear regression using 
Dim nls As NonlinearLeastSquares = New NonlinearLeastSquares
With nls
    .NumberOfParameters = 3

    ' Initialize to some random values
    .StartValues = {4.2, 0.3, 1}

    ' Let's assume a quadratic model function: ax² + bx + c
    .Function = Function(w, x) w(0) * x(0) * x(0) + w(1) * x(0) + w(2)

    ' Derivative in respect to the weights
    .Gradient = Sub(w, x, r)
                    r(0) = 2 * w(0) ' w.r.t a:     2a  
                    r(1) = w(1)     ' w.r.t b: b
                    r(2) = w(2)     ' w.r.t c: 0
                End Sub
End With

Dim algorithm As LevenbergMarquardt = New LevenbergMarquardt
With algorithm
    .MaxIterations = 100
    .Tolerance = 0
End With

nls.Algorithm = algorithm

Dim regression As NonlinearRegression = nls.Learn(inputs, outputs)

' Use the function to compute the input values
Dim predict As Double() = regression.Transform(inputs)

// Example from https://en.wikipedia.org/wiki/Gauss%E2%80%93Newton_algorithm

// In this example, the Gauss–Newton algorithm will be used to fit a model to 
// some data by minimizing the sum of squares of errors between the data and 
// model's predictions.

// In a biology experiment studying the relation between substrate concentration [S]
// and reaction rate in an enzyme-mediated reaction, the data in the following table
// were obtained:

double[][] inputs = Jagged.ColumnVector(new [] { 0.03, 0.1947, 0.425, 0.626, 1.253, 2.500, 3.740 });
double[] outputs = new[] { 0.05, 0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317 };

// It is desired to find a curve (model function) of the form
// 
//   rate = \frac{V_{max}[S]}{K_M+[S]}
// 
// that fits best the data in the least squares sense, with the parameters V_max
// and K_M to be determined. Let's start by writing model equation below:

LeastSquaresFunction function = (double[] parameters, double[] input) =>
{
    return (parameters[0] * input[0]) / (parameters[1] + input[0]);
};

// Now, we can either write the gradient function of the model by hand or let
// the model compute it automatically using Newton's finite differences method:

LeastSquaresGradientFunction gradient = (double[] parameters, double[] input, double[] result) =>
{
    result[0] = -((-input[0]) / (parameters[1] + input[0]));
    result[1] = -((parameters[0] * input[0]) / Math.Pow(parameters[1] + input[0], 2));
};

// Create a new Levenberg-Marquardt algorithm
var gn = new LevenbergMarquardt(parameters: 2)
{
    Function = function,
    Gradient = gradient,
    Solution = new[] { 0.9, 0.2 } // starting from b1 = 0.9 and b2 = 0.2
};

// Find the minimum value:
gn.Minimize(inputs, outputs);

// The solution will be at:
double b1 = gn.Solution[0]; // will be 0.362
double b2 = gn.Solution[1]; // will be 0.556