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

by Trent Guidry10. December 2011 09:20
In this post I create a Gauss Newton 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 GaussNewton : INonLinearRegressor{private Matrix _jacobian;private Matrix _residuals;private Matrix _regressionParameters0;private Derivatives _derivatives;private Parameter[] _regressionParameters;private Parameter[] _observedParameters;private Func<double>[] _regressionFunctions;private Matrix _weightMatrix;double[] _weights;private double[,] _data;public GaussNewton(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);_regressionParameters0 = new Matrix(numberOfParameters, 1);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 GaussNewton(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 k = 0; k < numberOfPoints; k++){for (int l = 0; l < _observedParameters.Length; l++){_observedParameters[l].Value = _data[k, l];}double functionValue = _regressionFunctions[i]();double residual = _data[k, _observedParameters.Length + i] - functionValue;_residuals[k + i * numberOfPoints, 0] = residual;currentResidual += residual * residual;for (int l = 0; l < numberOfParameters; l++){_jacobian[k + i * numberOfPoints, l] = _derivatives.ComputePartialDerivative(_regressionFunctions[i], _regressionParameters[l], 1, functionValue);}}}for (int i = 0; i < numberOfParameters; i++){_regressionParameters0[i, 0] = _regressionParameters[i];}Matrix matJacobianTrans = _jacobian.Transpose();Matrix jacobianTransposeWeighted = matJacobianTrans * _weightMatrix;Matrix matJTranResid = jacobianTransposeWeighted * _residuals;Matrix matJTranJ = jacobianTransposeWeighted * _jacobian;Matrix matDelta = matJTranJ.SolveFor(matJTranResid); ;Matrix matNewRegPars = _regressionParameters0 + matDelta;for (int i = 0; i < numberOfParameters; i++){_regressionParameters[i].Value = matNewRegPars[i, 0];}}}
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 GaussNewton(regressionFunctions, weights, regressionParameters, observedParameters, data, 5);for (int i = 0; i < 50; i++){regressor.Iterate();}