Implementing a Levenberg Marquardt algorithm for the nonlinear least squares regression of multiple weighted simultaneous equations in C#

by Trent Guidry10. December 2011 12:04
In this post I create a Levenberg Marquardt class in C# that can be used for the nonlinear least squares regression of multiple simultaneous weighted equations.
public interface INonLinearRegressor{void Iterate();}public class LevenbergMarquardt : INonLinearRegressor{private Matrix _jacobian;private Matrix _residuals;private Matrix _regressionParameters0;private Derivatives _derivatives;private Parameter[] _regressionParameters;private Parameter[] _observedParameters;private Func<double>[] _regressionFunctions;private double[] _weights;private Matrix _weightMatrix;private double[,] _data;private double _l0 = 100.0;private double _v = 10.0;public LevenbergMarquardt(Func<double>[] regressionFunctions, double[] weights, Parameter[] regressionParameters, Parameter[] observedParameters, double[,] data, int numberOfDerivativePoints){Debug.Assert(data.GetLength(1) == observedParameters.Length + regressionFunctions.Length);_data = data;_observedParameters = observedParameters;_regressionParameters = regressionParameters;_regressionFunctions = regressionFunctions;_weights = weights;int numberOfParameters = _regressionParameters.Length;int numberOfPoints = data.GetLength(0);_derivatives = new Derivatives(numberOfDerivativePoints);int numberOfFunctions = _regressionFunctions.Length;_jacobian = new Matrix(numberOfFunctions * numberOfPoints, numberOfParameters);_residuals = new Matrix(numberOfFunctions * numberOfPoints, 1);_regressionParameters0 = new Matrix(numberOfParameters, 1);_weightMatrix = new Matrix(numberOfPoints * numberOfFunctions, numberOfPoints * numberOfFunctions);for (int i = 0; i < numberOfFunctions * numberOfPoints; i++){for (int j = 0; j < numberOfFunctions * numberOfPoints; j++){if (i == j){int weightNumber = i / numberOfPoints;_weightMatrix[i, i] = _weights[weightNumber];}else{_weightMatrix[i, j] = 0.0;}}}}public LevenbergMarquardt(Func<double>[] regressionFunctions, double[] weights, Parameter[] regressionParameters, Parameter[] observedParameters, double[,] data) :this(regressionFunctions, weights, regressionParameters, observedParameters, data, 3){}public void Iterate(){int numberOfPoints = _data.GetLength(0);int numberOfParameters = _regressionParameters.Length;int numberOfFunctions = _regressionFunctions.Length;double currentResidual = 0.0;for (int i = 0; i < numberOfFunctions; i++){for (int j = 0; j < numberOfPoints; j++){for (int k = 0; k < _observedParameters.Length; k++){_observedParameters[k].Value = _data[j, k];}double functionValue = _regressionFunctions[i]();double residual = _data[j, _observedParameters.Length + i] - functionValue;_residuals[j + i * numberOfPoints, 0] = residual;currentResidual += residual * residual;for (int k = 0; k < numberOfParameters; k++){_jacobian[j + i * numberOfPoints, k] = _derivatives.ComputePartialDerivative(_regressionFunctions[i], _regressionParameters[k], 1, functionValue);}}}for (int i = 0; i < numberOfParameters; i++){_regressionParameters0[i, 0] = _regressionParameters[i];}Matrix jacobianTranspose = _jacobian.Transpose();Matrix jacobianTransposeWeighted = jacobianTranspose * _weightMatrix;Matrix jacobianTransposeResiduals = jacobianTransposeWeighted * _residuals;Matrix jacobianTransposeJacobian = jacobianTransposeWeighted * _jacobian;Matrix jacobianTransposeJacobianDiagnol = new Matrix(jacobianTransposeJacobian.RowCount, jacobianTransposeJacobian.RowCount);for (int i = 0; i < jacobianTransposeJacobian.RowCount; i++){jacobianTransposeJacobianDiagnol[i, i] = jacobianTransposeJacobian[i, i];}Matrix delta = null;Matrix newRegressionParameters = null;double newResidual = currentResidual + 1.0;_l0 /= _v;while (newResidual > currentResidual){newResidual = 0.0;_l0 *= _v;Matrix matLHS = jacobianTransposeJacobian + _l0 * jacobianTransposeJacobianDiagnol;delta = matLHS.SolveFor(jacobianTransposeResiduals); ;newRegressionParameters = _regressionParameters0 + delta;for (int i = 0; i < numberOfParameters; i++){_regressionParameters[i].Value = newRegressionParameters[i, 0];}for (int i = 0; i < numberOfFunctions; i++){for (int j = 0; j < numberOfPoints; j++){for (int k = 0; k < _observedParameters.Length; k++){_observedParameters[k].Value = _data[j, k];}double functionValue = _regressionFunctions[i]();double residual = _data[j, _observedParameters.Length + i] - functionValue;newResidual += residual * residual;}}}_l0 /= _v;}
Example
Func<double, double>[] actualFunctions = new Func<double, double>[] {(t) => -4 * Math.Exp(t) * Math.Sin(2.0 * t),(t) => 4 * Math.Exp(t) * Math.Cos(2.0 * t)};Parameter A = new Parameter(-2.0);Parameter B = new Parameter(3.0);Parameter C = new Parameter(4.4);Parameter D = new Parameter(3.0);Parameter time = new Parameter(0);Func<double>[] regressionFunctions = new Func<double>[] {() => A * Math.Exp(time) * Math.Sin(B * time),() => C * Math.Exp(time) * Math.Cos(D * time)};double[] weights = new double[] { 1, 1 };int numberOfFunctions = actualFunctions.Length;int numberOfPoints = 10;double[,] data = new double[numberOfPoints, numberOfFunctions + 1];for (int i = 0; i < numberOfPoints; i++){double t = ((double)i) / 10.0;data[i, 0] = t;for (int j = 0; j < numberOfFunctions; j++){data[i, j + 1] = actualFunctions[j](t);}}Parameter[] regressionParameters = new Parameter[] { A, B, C, D };Parameter[] observedParameters = new Parameter[] { time };INonLinearRegressor regressor = new LevenbergMarquardt(regressionFunctions, weights, regressionParameters, observedParameters, data, 5);for (int i = 0; i < 50; i++){regressor.Iterate();}