The Affine Covariant Newton Method

The affine covariant Newton method presented here is an implementation of the NLEQ_ERR algorithm [Deu06], and is available as a method option to (RefLink[FindRoot, paclet:ref/FindRoot]). To illustrate how the algorithm works and its advantages, some simple examples of the affine covariant Newton solver will be shown. Its main application area is for solving discretized nonlinear partial differential equations. These systems of equations can become large, and it is important to pay attention to evaluation costs both in the computational speed and the amount of memory needed to solve these equations.

We start by specifying an example system of equations for which the roots are sought. Consider the nonlinear system of equations:

({{(e^(x-2))-y, =, 0}, {(y^2)-x, =, 0}})

Given the variables and starting values, (RefLink[FindRoot, paclet:ref/FindRoot]) will find the roots for this system of equations.

Use (RefLink[FindRoot, paclet:ref/FindRoot]) to find the zero crossings:

FindRoot[{Exp[x-2]-y,y^2-x},{{x,1},{y,1}}]

To make things a bit easier and more uniform throughout this tutorial, set up a vector-valued function f that takes a vector argument and evaluates all equations.

Set up a function:

f[X_]:=Block[{x,y},{x,y}=X;{Exp[x-2]-y,y^2-x}]

Group variable names and starting vector:

vars={x,y};
start={1,1};

Call (RefLink[FindRoot, paclet:ref/FindRoot]) and replace the variables with values:

vars/.FindRoot[f[vars],Evaluate[Transpose[{vars,start}]]]

For comparisons in the following text, it is useful to have a reference solution that is known to be very precise.

Use higher precision to compute a reference solution:

referenceSolution=vars/.FindRoot[f[vars],Evaluate[Transpose[{vars,start}]],WorkingPrecision->32]

Here is a helper function that makes calling (RefLink[FindRoot, paclet:ref/FindRoot]) easier for the type of root finding this tutorial is presenting.

A wrapper function for (RefLink[FindRoot, paclet:ref/FindRoot]):

callFindRoot[f_,vars_,start_,opts___]:=vars/.FindRoot[f[vars],Evaluate[Transpose[{vars,start}]],opts]

Find the roots:

callFindRoot[f,vars, start]

Newton methods, such as the default method of (RefLink[FindRoot, paclet:ref/FindRoot]), will need to evaluate the function f and compute and evaluate the Jacobian of the function f. The evaluation of the Jacobian can be computationally expensive. The fact that the Jacobian needs to be evaluated is not something that can be avoided in the type of Newton methods discussed here; however, there exist algorithms that minimize the amount of the Jacobian evaluations. To investigate the evaluations needed by (RefLink[FindRoot, paclet:ref/FindRoot]), monitors such as (RefLink[EvaluationMonitor, paclet:ref/EvaluationMonitor]) and (RefLink[StepMonitor, paclet:ref/StepMonitor]) can be made use of.

It is instructive to monitor the number of function evaluations, the number of steps needed to get to a solution and the number of Jacobian evaluations during the root-finding process.

Inspect various evaluation counts in (RefLink[FindRoot, paclet:ref/FindRoot]):

Block[{s = 0, e = 0, j = 0}, 
sol=callFindRoot[f,vars,start,EvaluationMonitor:> e++, StepMonitor:> s++, Jacobian->{Automatic, EvaluationMonitor:> j++}];
{"Steps"->s,"Function Evaluations"-> e,"Jacobian Evaluations"-> j}]//TableForm

Seven steps and 10 function evaluations were needed to find a solution. Sometimes a trial step is made, which increases the number of function evaluations, but that trial step is then rejected. A rejected step does not increase the step count but increases the function evaluation count. We also see that seven Jacobian evaluations were needed to complete the solution process.

The default solution found by (RefLink[FindRoot, paclet:ref/FindRoot]) is very close to the reference solution:

referenceSolution-sol

(RefLink[FindRoot, paclet:ref/FindRoot]) has several method options, among them "AffineCovariantNewton". The affine covariant Newton method is explained in detail in [Deu06].

Inspect various evaluation counts in the "AffineCovariantNewton" method of (RefLink[FindRoot, paclet:ref/FindRoot]):

Block[{s = 0, e = 0, j = 0},
sol=callFindRoot[f,vars,start,Method->{"AffineCovariantNewton"},
EvaluationMonitor:> e++, StepMonitor:> s++, Jacobian->{Automatic, EvaluationMonitor:> j++}];
{"Steps"->s,"Function Evaluations"-> e,"Jacobian Evaluations"-> j}]//TableForm

The affine covariant Newton method needed more steps and function evaluations but fewer Jacobian evaluations. In a case such as solving nonlinear partial differential equations, where evaluating the Jacobian is computationally expensive, not having to evaluate the Jacobian as often can make a big difference in being able to find a solution to the system of equations efficiently.

Compare to the reference solution:

referenceSolution-sol

Note that the precision of the affine covariant solution is less than what the default (RefLink[FindRoot, paclet:ref/FindRoot]) method found, but it is still within the default settings for (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]). Before we look at that, however, we inspect the "AffineCovariantNewton" method option "BroydenUpdates" that can be set to (RefLink[True, paclet:ref/True]) or (RefLink[False, paclet:ref/False]). Broyden updates are a relatively inexpensive way of estimating updated factors of the Jacobian without having to refactor the Jacobian. The affine covariant Newton solver estimates when to use Broyden updates, and computes them efficiently. If during the solution process the nonlinearity increases and Broyden updates can no longer be made useful, the algorithm switches back to compute the Jacobians. The Broyden update algorithm implemented is the QNERR algorithm [Deu06].

Find the root of the system without Broyden updates:

Block[{s = 0, e = 0, j = 0}, 
sol=callFindRoot[f,vars,start,Method->{"AffineCovariantNewton","BroydenUpdates"->False},
EvaluationMonitor:> e++, StepMonitor:> s++, Jacobian->{Automatic, EvaluationMonitor:> j++}];
{"Steps"->s,"Function Evaluations"-> e,"Jacobian Evaluations"-> j}]//TableForm

With Broyden updates switched off, finding a solution takes fewer steps and function evaluations but more Jacobian evaluations than when Broyden updates are done.

Compare to the reference solution:

referenceSolution-sol

Coming back to the previous issue of increasing the precision of a sought solution, the precision increase can be done by increasing the precision goal.

Specify a (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]) to reduce the difference to the Newton solution:

Block[{s = 0, e = 0, j = 0}, 
sol=callFindRoot[f,vars,start,Method->{"AffineCovariantNewton"},
PrecisionGoal->14,
EvaluationMonitor:> e++, StepMonitor:> s++, Jacobian->{Automatic, EvaluationMonitor:> j++}];
{"Steps"->s,"Function Evaluations"-> e,"Jacobian Evaluations"-> j}]//TableForm

Compare to the reference solution:

referenceSolution-sol

The affine covariant Newton method is an error norm–controlled algorithm, and as such, an (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) is not monitored. In other words, the algorithm only monitors the error and not the residual. The error is associated with the (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]), while the residual is associated with the accuracy goal.

Method Options

The affine covariant Newton method takes a number of options in addition to the (RefLink[FindRoot, paclet:ref/FindRoot]) options. These method options will be discussed. First, a list of available options is presented. Then these options are illustrated with examples, and typical scenarios where the options may be beneficial are presented.

The following method options can be given for the (RefLink[Method, paclet:ref/Method])->"AffineCovariantNewton" in (RefLink[FindRoot, paclet:ref/FindRoot]).

Method options for (RefLink[Method, paclet:ref/Method])->"AffineCovariantNewton" in (RefLink[FindRoot, paclet:ref/FindRoot]).

"BroydenUpdates"

Broyden updates are a way to reuse a previously computed Jacobian. Reusing a previously computed Jacobian is much more efficient than recomputing a new Jacobian. The downside is that the Broyden update may not be as precise as a newly computed Jacobian. The affine covariant Newton method algorithm has a mechanism implemented that switches to Broyden updates when it is safe to do so, and may also switch back to compute the Jacobians when Broyden updates are no longer usable. Nevertheless, Broyden updates can be disabled.

Find the root of the system without Broyden updates:

Block[{s = 0, e = 0, j = 0}, 
callFindRoot[f,vars,start,Method->{"AffineCovariantNewton","BroydenUpdates"->False},
EvaluationMonitor:> e++, StepMonitor:> s++, Jacobian->{Automatic, EvaluationMonitor:> j++}];
{"Steps"->s,"Function Evaluations"-> e,"Jacobian Evaluations"-> j}]//TableForm

"MinimalDampingFactor"

During step k of the solution process, this equation is solved: F'((x^k))((δx)^k)=-F((x^k)). The Jacobian F' is inverted to obtain a new increment ((δx)^k). This increment is then added to the previous solution to obtain an updated solution (x^(k+1))=(x^k)+((λδx)^k). λ is a damping factor between (Subscript[λ, min])<λ≤1 that restricts the step size taken. For extreme nonlinearities, the default minimal damping factor (Subscript[λ, min]) can be reduced. The default value is (10^(-4)). Often, this should not be necessary.

"InitialDampingFactor"

When it is known that the nonlinearity is not too extreme, a larger initial damping factor can be chosen. The default initial damping factor is 1/100. A chosen value should be larger than the "MinimalDampingFactor" but smaller or equal to 1: (Subscript[λ, min])<(Subscript[λ, init])≤1.

Specify a larger "InitialDampingFactor":

Block[{s = 0, e = 0, j = 0}, 
callFindRoot[f,vars,start,Method->{"AffineCovariantNewton","InitialDampingFactor"->1},
EvaluationMonitor:> e++, StepMonitor:> s++, Jacobian->{Automatic, EvaluationMonitor:> j++}];
{"Steps"->s,"Function Evaluations"-> e,"Jacobian Evaluations"-> j}]//TableForm

A larger initial damping factor can save Jacobian evaluations.

"InverseJacobian"

The Newton algorithm actually needs the inverse of the Jacobian to do its job. This inverse is computed by creating a (RefLink[LinearSolveFunction, paclet:ref/LinearSolveFunction]). In some cases, the inverse Jacobian may be known or easy to compute. In such a scenario, the inverse Jacobian can be passed to the "AffineCovariantNewton" method via a method option.

Construct a function that returns the evaluated inverse Jacobian at the evaluation point ep:

makeInverseJacobianFunction[f_,vars_]:=With[{inv=Inverse[D[f[vars],{vars}]],v=vars},
Function[{X},
Block[v,v=X;With[{invEval=inv},
Function[ep,invEval.ep]]]
]
]

Construct the Jacobian inverse for the problem function f:

Jinverse=makeInverseJacobianFunction[f,vars]

Find the roots with a specified inverse of the Jacobian:

Block[{s = 0, e = 0, j = 0}, 
sol=callFindRoot[f,vars,start,Method->{"AffineCovariantNewton","InverseJacobian"->Jinverse},
EvaluationMonitor:> e++, StepMonitor:> s++, Jacobian->{Automatic, EvaluationMonitor:> j++}];
{"Steps"->s,"Function Evaluations"-> e,"Jacobian Evaluations"-> j}]//TableForm

sol

Specifying the inverse of the Jacobian does not have any affect on the number of Jacobian evaluations. The fact that the count of Jacobian evaluations are shown to be 0 is an unfortunate shortcoming of the way the inverse Jacobian interacts with (RefLink[FindRoot, paclet:ref/FindRoot]).

"LinearSolver"

It is possible change the linear solver or the options used by (RefLink[LinearSolve, paclet:ref/LinearSolve]) to construct the inverse of the Jacobian. This is done with the "LinearSolver" method option.

Specify the iterative "Krylov" method for (RefLink[LinearSolve, paclet:ref/LinearSolve]) to use:

callFindRoot[f,vars,start,Method->{"AffineCovariantNewton","LinearSolver"->{Automatic,Method->"Krylov"}}]

When an iterative method like "Krylov" is specified as a method to (RefLink[LinearSolve, paclet:ref/LinearSolve]), then the factorization of the Jacobian is not stored when possible, but computed anew every time it is needed.