Solving Systems of Non-Linear Equations

Newton's method is a generalization of the Newton-Raphson method for solving equations in one variable to the multi-dimensional case. The Jacobian of the system takes the place of the derivative.

Newton's method is implemented by the NewtonRaphsonSystemSolver class

Powell's dogleg method, also called Powell's hybrid method, attempts to minimize the sum of the squares of the function values. It does this using a combination of Newton's method and the steepest descent method. This is a so-called trust region method. This means that every step moves the current point to within a finite region. This makes the method more stable than Newton's method.

On the other hand, the fact that the method is, in essence, a specialized minimizer means that the algorithm can get stuck in a local minimum that does not correspond to a solution of the system of equations.

The DoglegSystemSolver class implements the dog leg algorithm. This class uses derivatives when available, and uses an optimized variant using numerical derivatives when they are not. It has one additional property, TrustRegionRadius, which allows you to set the initial trust region radius. The default is 1.

The equations that make up the system of equations can be defined in two ways. Viewed as a set of equations, the system is a collection of multivariate functions that are each set equal to a real number (usually 0). The functions must be supplied as an array of FuncT, TResult delegates that map a VectorT to a real number. These can be supplied to the constructor, or set later using the SetTargetFunctions method. The example below defines the system of equations that is the basis for the Rosenbrock function:

1 - x = 0

10(y - x²) = 0

Func<Vector<double>, double>[] rosenbrock =
{
    x => 1 - x[0],
    x => 10 * (x[1] - x[0] * x[0])
};

Dim rosenbrock() As Func(Of Vector(Of Double), Double) =
{
    Function(x) 1 - x(0),
    Function(x) 10 * (x(1) - x(0) * x(0))
}

let rosenbrock : Func<Vector<float>, float>[] = 
  [|
    Func<_,_>(fun x -> 1.0 - x.[0]);
    Func<_,_>(fun x -> 10.0 * (x.[1] - x.[0] * x.[0]))
  |]

A system of equations can also be looked at as a single vector equation, where a multivariate function returns a vector that is set equal to a vector (usually all 0's). This function is supplied as one FuncT1, T2, TResult delegate. The first argument is a VectorT that contains the x-values. The second argument is a vector that is to hold the result. The delegate can be passed to the constructor. It can also be set and retrieved through the TargetFunction property. This is the default representation of the system of equations. The example below illustrates how to define the same system of equations using a single delegate that evaluates all equations simultaneously.

Vector<double> Rosenbrock(Vector<double> x, Vector<double> f) {
    if (f == null)
        f = Vector.Create<double>(2);
    f[0] = 1 - x[0];
    f[1] = 10 * (x[1] - x[0] * x[0]);
    return f;
}

Shared Function Rosenbrock(x As Vector(Of Double), f As Vector(Of Double)) As Vector(Of Double)
    If (f Is Nothing) Then f = Vector.Create(Of Double)(2)
    f(0) = 1 - x(0)
    f(1) = 10 * (x(1) - x(0) * x(0))
    Return f
End Function

let Rosenbrock (x : Vector<double>) (f : Vector<double>) =
    f.[0] <- 1.0 - x.[0]
    f.[1] <- 10.0 * (x.[1] - x.[0] * x.[0])
    f

The RightHandSide property is a VectorT that contains the right-hand sides of the equations. If no value is specified, all right-hand sides are assumed to be zero.

Some algorithms require the Jacobian (matrix of partial derivatives) of the system. Each row of the Jacobian corresponds to the gradient of the corresponding target function. There are two ways to supply this function. The first is through the JacobianFunction property. This property is a FuncT1, T2, TResult. This is a method that takes a VectorT argument and returns the matrix result in a second MatrixT. The return value is also a reference to this matrix.

Matrix<double> RosenbrockJacobian(Vector<double> x, Matrix<double> J)
{
    if (J == null)
        J = Matrix.Create<double>(2, 2);
    J[0, 0] = -1;
    J[0, 1] = 0;
    J[1, 0] = -20 * x[0];
    J[1, 1] = 10;
    return J;
}

Function RosenbrockJacobian(x As Vector(Of Double), J As Matrix(Of Double)) As Matrix(Of Double)
    If (J Is Nothing) Then J = Matrix.Create(Of Double)(2, 2)
    J(0, 0) = -1
    J(0, 1) = 0
    J(1, 0) = -20 * x(0)
    J(1, 1) = 10
    Return J
End Function

let RosenbrockJacobian (x : Vector<double>) (J : Matrix<double>) =
    J.[0, 0] <- -1.0
    J.[0, 1] <- 0.0
    J.[1, 0] <- -20.0 * x.[0]
    J.[1, 1] <- 10.0
    J

Alternatively, the individual gradients can be specified using the SetGradientFunctions property. This method takes an array of either FuncT, TResult. or FuncT1, T2, TResult delegates. If the delegate has a second argument, it represents the vector that is to hold the result.

Func<Vector<double>, Vector<double>>[] rosenbrockGradients =
{
    x => Vector.Create(-1.0, 0.0),
    x => Vector.Create(-20.0 * x[0], 10)
};

Dim rosenbrockGradients() As Func(Of Vector(Of Double), Vector(Of Double)) =
{
    Function(x) Vector.Create(-1.0, 0.0),
    Function(x) Vector.Create(-20.0 * x(0), 10)
}

let rosenbrockGradients : Func<Vector<float>, Vector<float>>[] = [|
    Func<_,_>(fun x -> Vector.Create(-1.0, 0.0) * 1.0);
    Func<_,_>(fun x -> Vector.Create(-20.0 * x.[0], 10.0) * 1.0) |]

The following code puts this all together. It creates a DogLegSystemSolver and defines the system of equations using the functions defined in earlier examples:

DoglegSystemSolver dogleg = new DoglegSystemSolver();
dogleg.InitialGuess = Vector.Create(0.5, 0.5);
dogleg.TargetFunction = Rosenbrock;
dogleg.SetGradientFunctions(rosenbrockGradients);

Dim dogleg As DoglegSystemSolver = New DoglegSystemSolver()
dogleg.InitialGuess = Vector.Create(0.5, 0.5)
dogleg.TargetFunction = AddressOf Rosenbrock
dogleg.SetGradientFunctions(rosenbrockGradients)

let dogleg = DoglegSystemSolver()
dogleg.InitialGuess <- Vector.Create(0.5, 0.5)
dogleg.TargetFunction <- Func<_,_,_>(Rosenbrock)
dogleg.SetGradientFunctions(rosenbrockGradients)

If the functions are of a simple enough form, then automatic differentiation can be used. This completely eliminates the need to define derivatives, gradients, or Jacobians by hand. The code is also very concise:

dogleg.SetSymbolicTargetFunctions(
    x => 1 - x[0],
    x => 10 * (x[1] - x[0] * x[0]));

dogleg.SetSymbolicTargetFunctions(
    Function(x) 1 - x(0),
    Function(x) 10 * (x(1) - x(0) * x(0)))

dogleg.SetSymbolicTargetFunctions(
    (fun (x : Vector<float>) -> 1.0 - x.[0]),
    (fun (x : Vector<float>) -> 10.0 * (x.[1] - x.[0] * x.[0])))

Algorithms for solving systems of nonlinear equations typically have two criteria that may each signal that a solution has been found.

The ValueTest property returns a VectorConvergenceTestT object that represents the convergence test based on the function values. The test is successful when some norm of the vector of function values is less than the tolerance. By default, the function value with the largest absolute value is used. The test returns Divergent if one of the function values is infinite, and BadFunction when the value is NaN.

The SolutionTest property returns a VectorConvergenceTestT object that represents the convergence test based on the estimated solution. The test is successful if the change in the approximate solution is less than the tolerance. The test returns a value of Divergent if one of the elements of the solution is infinite, and BadFunction when one of the values is NaN.

The VectorConvergenceTestT has many options to specify the details of the convergence test. By default, the error is calculated as the maximum of the errors in each of the elements of the change of the solution from the previous iteration compared to the corresponding elements of the solution.

The following example specifies that the algorithm is to end when the sum of the absolute values of the function values is less than 0.00001. The solution test is not used:

dogleg.ValueTest.Tolerance = 0.00001;
dogleg.ValueTest.Norm = VectorConvergenceNorm.SumOfAbsoluteValues;
dogleg.SolutionTest.Enabled = false;

dogleg.ValueTest.Tolerance = 0.00001
dogleg.ValueTest.Norm = VectorConvergenceNorm.SumOfAbsoluteValues
dogleg.SolutionTest.Enabled = False

dogleg.ValueTest.Tolerance <- 0.00001
dogleg.ValueTest.Norm <- VectorConvergenceNorm.SumOfAbsoluteValues
dogleg.SolutionTest.Enabled <- false

You must specify an initial guess for the solution by setting the InitialGuess property to a suitable VectorT value.

The Solve method does the actual work of solving the system of equations. When called without parameters, it returns the approximation to the zero of the target function. When called with a single VectorT argument, it attempts to find a point where the target functions equal the corresponding components of the vector.

The Solve method always returns the best estimate for the solution of the system. Successive calls to the Result property will also return this value, until the next call to Solve. An optional parameter lets you specify the right-hand side of the equation TargetFunction(x) = rightHandSide that is to be solved. You can also set this value through the RightHandSide property.

If the ThrowExceptionOnFailure property is set to true, an exception is thrown if the algorithm has failed to converge to a solution within the desired accuracy. If false, the Solve method returns the best approximation to the zero, regardless of whether it is within the requested tolerance.

The Status property indicates how the algorithm terminated. Its possible values and their meaning are listed below.

The normal status for successful completion is Converged. Of particular interest in this instance is the value ConvergedToFalseSolution. This status is returned when the algorithm converges, but the solution that was found is not a solution of the system of equations. This can happen when an algorithm that solves the system by minimizing the sum of the squares of the function values converges to a local minimum that does not coincide with the global minimum of the sum of squares.

Several properties provide more information about the algorithm. The IterationsNeeded property returns the number of iterations the algorithm needed to complete. The EvaluationsNeeded property returns the number of times the target functions were evaluated. The convergence tests, ValueTest and SolutionTest, each have an Error property that gives an estimate for the error in the solution.

In the next example, we print a summary of the results of the algorithm:

Console.WriteLine("  Initial guess: {0:F2}", dogleg.InitialGuess);
Console.WriteLine("  Status: {0}", dogleg.Status);
Console.WriteLine("  Solution: {0}", dogleg.Result);
Console.WriteLine("  Estimated error: {0}", dogleg.ValueTest.Error);
Console.WriteLine("  # iterations: {0}", dogleg.IterationsNeeded);
Console.WriteLine("  # evaluations: {0}", dogleg.EvaluationsNeeded);

Console.WriteLine("  Initial guess: 0:F2", dogleg.InitialGuess)
Console.WriteLine("  Status: {0}", dogleg.Status)
Console.WriteLine("  Solution: {0}", dogleg.Result)
Console.WriteLine("  Estimated error: {0}", dogleg.ValueTest.Error)
Console.WriteLine("  # iterations: {0}", dogleg.IterationsNeeded)
Console.WriteLine("  # evaluations: {0}", dogleg.EvaluationsNeeded)

Console.WriteLine("  Initial guess: 0:F2", dogleg.InitialGuess)
Console.WriteLine("  Status: 0", dogleg.Status)
Console.WriteLine("  Solution: 0", dogleg.Result)
Console.WriteLine("  Estimated error: 0", dogleg.ValueTest.Error)
Console.WriteLine("  # iterations: 0", dogleg.IterationsNeeded)
Console.WriteLine("  # evaluations: 0", dogleg.EvaluationsNeeded)

J.J. Moré, M.Y. Cosnard, Numerical Solution of Nonlinear Equations", ACM Transactions on Mathematical Software , Vol 5, No 1, 1979.

M. J. D. Powell, A Hybrid Method for Nonlinear Equations. Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz, editor. Gordon and Breach, 1970.