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

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();
			}