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

In this post I create a Levenberg 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 Levenberg : 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 Levenberg(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 Levenberg(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 identity = Matrix.Identity(jacobianTransposeJacobian.RowCount);

			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 * identity;
				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 Levenberg(regressionFunctions, weights, regressionParameters, observedParameters, data, 5);
			for (int i = 0; i < 50; i++)
			{
				regressor.Iterate();
			}