2.1 Local convergence of Newton’s method

The reader may know Theorem 2.2 well. The simple statement is that if the standard assumptions hold and the initial iterateFootnote ¹ is sufficiently near $\mathbf{x}^{\ast }$ (hence the term local), then the Newton iteration $\{\mathbf{x}_{n}\}$ exists (i.e.  $\mathbf{F}^{\prime }(\mathbf{x}_{n})$ is non-singular for all $n\geq 0$ ) and converges quadratically to $\mathbf{x}^{\ast }$ . The Newton iterates are, of course,

$$\begin{eqnarray}\mathbf{x}_{n+1}=\mathbf{x}_{n}-\mathbf{F}(\mathbf{x}_{n})^{-1}\mathbf{F}(\mathbf{x}_{n}),\end{eqnarray}$$

for $n\geq 0$ . This is exactly (1.2) with $\mathbf{x}_{c}$ replaced by $\mathbf{x}_{n}$ and $\mathbf{x}_{+}$ by $\mathbf{x}_{n+1}$ . The advantage of the $\mathbf{x}_{c}$ and $\mathbf{x}_{+}$ notation is that the transition from $\mathbf{x}_{c}$ to $\mathbf{x}_{+}$ is central and the iteration counter is generally irrelevant to the convergence analysis. Quadratic convergence means that

(2.2) $$\begin{eqnarray}\Vert \mathbf{e}_{+}\Vert =O(\Vert \mathbf{e}_{c}\Vert ^{2}),\end{eqnarray}$$

where the error in $\mathbf{x}$ is $\mathbf{e}=\mathbf{x}-\mathbf{x}^{\ast }$ . Quadratic convergence says that the number of significant figures in the result roughly doubles with each iteration.

We will begin with a precise statement of this result which will, among other things, exhibit the constant in the $O$ -term. We will also begin to develop a taxonomy of convergence types. We say that $\mathbf{x}_{n}\rightarrow \mathbf{x}^{\ast }$ q-quadratically if (2.2) holds and r-quadratically if there is a real sequence $\{\unicode[STIX]{x1D709}_{n}\}$ which converges q-quadratically to $0$ such that $\Vert \mathbf{e}_{n}\Vert \leq \unicode[STIX]{x1D709}_{n}$ . We say the convergence is q-linear if there is $\unicode[STIX]{x1D6FC}\in [0,1)$ such that $\Vert \mathbf{e}_{n+1}\Vert \leq \unicode[STIX]{x1D6FC}\Vert \mathbf{e}_{n}\Vert$ for $n$ sufficiently large; $\unicode[STIX]{x1D6FC}$ is called the q-factor. The convergence is q-superlinear if

$$\begin{eqnarray}\lim _{n\rightarrow \infty }\displaystyle \frac{\Vert \mathbf{e}_{n+1}\Vert }{\Vert \mathbf{e}_{n}\Vert }=0.\end{eqnarray}$$

Finally, we will quantify ‘sufficiently near $\mathbf{x}^{\ast }$ ’. At a minimum, all the local convergence results in this section require

(2.3) $$\begin{eqnarray}|\mathbf{x}-\mathbf{x}^{\ast }\Vert \leq \min \biggl(\displaystyle \frac{\Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })^{-1}\Vert ^{-1}}{2\unicode[STIX]{x1D6FE}},\unicode[STIX]{x1D70C}^{\ast }\biggr),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}$ is the Lipschitz constant for $\mathbf{F}^{\prime }$ and $\unicode[STIX]{x1D70C}^{\ast }$ is the radius of the ball about $\mathbf{x}^{\ast }$ in which the Lipschitz continuity assumption on $\mathbf{F}^{\prime }$ holds.

2.2 Classical Newton’s method

We will prove Theorem 2.2 in detail. Not only is the proof illuminating in its own right, but some of the details lead to algorithmic insights.

Theorem 2.2. Let Assumption 2.1 hold and assume that $\mathbf{x}_{c}$ satisfies (2.3). Then $\mathbf{F}(\mathbf{x}_{c})$ is non-singular,

(2.4) $$\begin{eqnarray}\Vert \mathbf{F}^{\prime }(\mathbf{x}_{c})\Vert ^{-1}\leq 2\Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })\Vert ^{-1},\end{eqnarray}$$

and

(2.5) $$\begin{eqnarray}\Vert \mathbf{e}_{+}\Vert \leq \Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })^{-1}\Vert \unicode[STIX]{x1D6FE}\,\Vert \mathbf{e}_{c}\Vert ^{2}\leq \Vert \mathbf{e}_{c}\Vert /2.\end{eqnarray}$$

Proof. The standard assumptions and (2.3) imply that

$$\begin{eqnarray}\Vert \mathbf{F}^{\prime }(\mathbf{x}_{c})-\mathbf{F}^{\prime }(\mathbf{x}^{\ast })\Vert \leq \unicode[STIX]{x1D6FE}\Vert \mathbf{e}_{c}\Vert \leq \displaystyle \frac{\Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })^{-1}\Vert }{2},\end{eqnarray}$$

and hence

$$\begin{eqnarray}\Vert \mathbf{I}-\mathbf{F}^{\prime }(\mathbf{x}^{\ast })^{-1}\mathbf{F}^{\prime }(\mathbf{x}_{c})\Vert \leq \Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })^{-1}\Vert \Vert \mathbf{F}^{\prime }(\mathbf{x}_{c})-\mathbf{F}^{\prime }(\mathbf{x}^{\ast })\Vert \leq 1/2,\end{eqnarray}$$

so $\mathbf{F}^{\prime }(\mathbf{x}^{\ast })^{-1}$ is an approximate inverse of $\mathbf{F}^{\prime }(\mathbf{x}_{c})$ , and

$$\begin{eqnarray}\Vert \mathbf{F}^{\prime }(\mathbf{x}_{c})^{-1}\Vert \leq \displaystyle \frac{\Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })^{-1}\Vert }{2},\end{eqnarray}$$

as asserted.

The remainder of the proof follows from the fundamental theorem of calculus, which implies that

(2.6) $$\begin{eqnarray}\displaystyle \mathbf{F}(\mathbf{x}_{c}) & = & \displaystyle \int _{0}^{1}\mathbf{F}^{\prime }(\mathbf{x}^{\ast }+t\mathbf{e}_{c})\mathbf{e}_{c}\,\text{d}t\nonumber\\ \displaystyle & = & \displaystyle \mathbf{F}^{\prime }(\mathbf{x}_{c})\mathbf{e}_{c}+\int _{0}^{1}(\mathbf{F}^{\prime }(\mathbf{x}^{\ast }+t\mathbf{e}_{c})-\mathbf{F}^{\prime }(\mathbf{x}_{c}))\mathbf{e}_{c}\,\text{d}t.\end{eqnarray}$$

Hence,

$$\begin{eqnarray}\displaystyle \mathbf{e}_{+} & = & \displaystyle \mathbf{e}_{c}-\mathbf{F}^{\prime }(\mathbf{x}_{c})^{-1}\biggl(\mathbf{F}^{\prime }(\mathbf{x}_{c})\mathbf{e}_{c}+\int _{0}^{1}(\mathbf{F}^{\prime }(\mathbf{x}^{\ast }+t\mathbf{e}_{c})-\mathbf{F}^{\prime }(\mathbf{x}_{c}))\mathbf{e}_{c}\,\text{d}t\biggr)\nonumber\\ \displaystyle & = & \displaystyle -\mathbf{F}^{\prime }(\mathbf{x}_{c})^{-1}\biggl(\int _{0}^{1}(\mathbf{F}^{\prime }(\mathbf{x}^{\ast }+t\mathbf{e}_{c})-\mathbf{F}^{\prime }(\mathbf{x}_{c}))\mathbf{e}_{c}\,\text{d}t\biggr).\nonumber\end{eqnarray}$$

Note that

$$\begin{eqnarray}\biggl\|\int _{0}^{1}(\mathbf{F}^{\prime }(\mathbf{x}^{\ast }+t\mathbf{e}_{c})-\mathbf{F}^{\prime }(\mathbf{x}_{c}))\mathbf{e}_{c}\,\text{d}t\biggr\|\leq \int _{0}^{1}\unicode[STIX]{x1D6FE}\Vert \mathbf{e}_{c}\Vert ^{2}(1-t)\,\text{d}t=\unicode[STIX]{x1D6FE}\Vert \mathbf{e}_{c}\Vert ^{2}/2.\end{eqnarray}$$

Hence,

(2.7) $$\begin{eqnarray}\displaystyle \Vert \mathbf{e}_{+}\Vert & \leq & \displaystyle \displaystyle \frac{\Vert \mathbf{F}^{\prime }(\mathbf{x}_{c})^{-1}\Vert \unicode[STIX]{x1D6FE}}{2}\Vert \mathbf{e}_{c}\Vert ^{2}\leq \displaystyle \frac{2\Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })^{-1}\Vert \unicode[STIX]{x1D6FE}}{2}\Vert \mathbf{e}_{c}\Vert ^{2}\nonumber\\ \displaystyle & \leq & \displaystyle \Vert \mathbf{e}_{c}\Vert /2.\end{eqnarray}$$

This completes the proof and shows that the constant in the $O$ -term for quadratic convergence is no larger than $\Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })^{-1}\Vert \unicode[STIX]{x1D6FE}$ .☐

Our bounds for the distance of the initial iterate from the root and for the convergence rate depend only on the norm of the inverse of the Jacobian at the root and the Lipschitz constant of the Jacobian near the root. This observation leads directly to the Kantorovich theorem (Kantorovich and Akilov Reference Kantorovich and Akilov1982: see Section 3) and the implicit function theorem (Keller Reference Keller1987) in Section 4, which are the basis for the pseudo-arclength continuation method we discuss in Section 4.

2.3 Termination criteria

The only obvious ways to terminate the Newton iteration are to examine the norm of the residual $\Vert \mathbf{F}(\mathbf{x}_{n})\Vert$ or the norm of the step $\Vert \mathbf{x}_{n+1}-\mathbf{x}_{n}\Vert$ . Either way is fine if one computes the Newton step exactly.

The norm of the step is a very good surrogate for the norm of the error at the previous iteration. To see this, note that quadratic convergence implies that

$$\begin{eqnarray}\mathbf{x}_{n+1}-\mathbf{x}_{n}=\mathbf{e}_{n+1}-\mathbf{e}_{n}=-\mathbf{e}_{n}+O(\Vert \mathbf{e}_{n}\Vert ^{2}).\end{eqnarray}$$

So, as the iteration converges, the norm of the step is asymptotically equal to the norm of the previous error. Suppose, for example, one wishes to terminate the iteration when $\Vert \mathbf{e}_{n}\Vert \leq \unicode[STIX]{x1D70F}$ . One could very safely stop the iteration when $\Vert \mathbf{x}_{n+1}-\mathbf{x}_{n}\Vert \leq \unicode[STIX]{x1D70F}$ and return $\mathbf{x}_{n+1}$ as the solution. Alternatively, one could terminate the iteration when $\Vert \mathbf{x}_{n}-\mathbf{x}_{n-1}\Vert \leq \unicode[STIX]{x1D6FC}\sqrt{\unicode[STIX]{x1D70F}}$ , where $\unicode[STIX]{x1D6FC}$ is a small constant. Then quadratic convergence would imply that

$$\begin{eqnarray}\Vert \mathbf{e}_{n}\Vert =O(\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D70F}),\end{eqnarray}$$

which would suffice if $\unicode[STIX]{x1D6FC}$ were small enough to balance the constant in the $O$ -term.

If the iteration converges q-superlinearly, then

$$\begin{eqnarray}\Vert \mathbf{s}\Vert =\Vert \mathbf{x}_{+}-\mathbf{x}_{c}\Vert =\Vert \mathbf{e}_{c}\Vert +o(\Vert \mathbf{e}_{c}\Vert )\end{eqnarray}$$

and the step is still an excellent surrogate for the error in the previous iteration, but one has less information than in the quadratically convergent case and cannot use $\Vert \mathbf{s}\Vert$ to estimate $\Vert \mathbf{e}_{n+1}\Vert$ . One can apply similar logic if one has an accurate upper bound for the q-factor in a q-linearly convergent iteration. If one knows that

$$\begin{eqnarray}\Vert \mathbf{e}_{+}\Vert \leq \unicode[STIX]{x1D6FC}\Vert \mathbf{e}_{c}\Vert ,\end{eqnarray}$$

then

$$\begin{eqnarray}(1-\unicode[STIX]{x1D6FC})\Vert \mathbf{e}_{c}\Vert \leq \Vert \mathbf{s}\Vert .\end{eqnarray}$$

From this we can recover an estimate for the error in terms of the step,

$$\begin{eqnarray}\Vert \mathbf{e}_{+}\Vert \leq \unicode[STIX]{x1D6FC}\Vert \mathbf{e}_{c}\Vert \leq \displaystyle \frac{\unicode[STIX]{x1D6FC}}{1-\unicode[STIX]{x1D6FC}}\Vert \mathbf{s}\Vert .\end{eqnarray}$$

See Petzold (Reference Petzold and Stepleman1983) or Ascher and Petzold (Reference Ascher and Petzold1998) for examples of how this can be used in an initial value problem integration, and Tocci, Kelley and Miller (Reference Tocci, Kelley and Miller1997) for an example of the limitations of this idea.

The relation of the residual norm to the norm of the error is very similar to that for the linear case. In the linear case the equation is $\mathbf{A}\mathbf{x}=\mathbf{b}$ , the residual is $\mathbf{r}=\mathbf{b}-\mathbf{A}\mathbf{x}$ , and the error is $\mathbf{e}=\mathbf{x}-\mathbf{A}^{-1}b$ . The standard result is

$$\begin{eqnarray}\unicode[STIX]{x1D705}(\mathbf{A})^{-1}\displaystyle \frac{\Vert \mathbf{r}\Vert }{\Vert \mathbf{b}\Vert }\leq \displaystyle \frac{\Vert \mathbf{e}\Vert }{\Vert \mathbf{A}^{-1}b\Vert }\leq \unicode[STIX]{x1D705}(\mathbf{A})\displaystyle \frac{\Vert \mathbf{r}\Vert }{\Vert \mathbf{b}\Vert }\end{eqnarray}$$

(Kelley Reference Kelley1995, Demmel Reference Demmel1997, Golub and Van Loan Reference Golub and Van Loan1996). This familiar estimate compares the residual at $\mathbf{x}$ to the residual at the zero vector. In the nonlinear case it is not generally useful to use the zero vector as a reference, so we will use the initial iterate $\mathbf{x}_{0}$ instead.

Suppose the standard assumptions hold and $\mathbf{x}_{0}$ and $\mathbf{x}$ both satisfy (2.3). Then, using the fundamental theorem of calculus, as we did in (2.6),

(2.8) $$\begin{eqnarray}\displaystyle \mathbf{F}(\mathbf{x}) & = & \displaystyle \int _{0}^{1}\mathbf{F}^{\prime }(\mathbf{x}^{\ast }+t\mathbf{e})\mathbf{e}\,\text{d}t\nonumber\\ \displaystyle & = & \displaystyle \mathbf{F}^{\prime }(\mathbf{x}^{\ast })\mathbf{e}+\int _{0}^{1}(\mathbf{F}^{\prime }(\mathbf{x}^{\ast }+t\mathbf{e})-\mathbf{F}^{\prime }(\mathbf{x}^{\ast }))\mathbf{e}\,\text{d}t.\end{eqnarray}$$

Therefore, using (2.3),

$$\begin{eqnarray}\displaystyle \Vert \mathbf{F}(\mathbf{x})\Vert & \leq & \displaystyle \Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })\mathbf{e}\Vert +\unicode[STIX]{x1D6FE}\Vert \mathbf{e}\Vert ^{2}/2,\nonumber\\ \displaystyle & \leq & \displaystyle \Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })\mathbf{e}\Vert +\Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })^{-1}\Vert ^{-1}\Vert \mathbf{e}\Vert /4\leq \displaystyle \frac{5}{4}\Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })\mathbf{e}\Vert .\nonumber\end{eqnarray}$$

Similarly

$$\begin{eqnarray}\Vert \mathbf{F}(\mathbf{x})\Vert \geq \displaystyle \frac{3}{4}\Vert \mathbf{F}(\mathbf{x}^{\ast })\mathbf{e}\Vert ,\end{eqnarray}$$

and the same inequalities hold for $\mathbf{x}_{0}$ . Hence

(2.9) $$\begin{eqnarray}\displaystyle \frac{3}{5}\unicode[STIX]{x1D705}(\mathbf{F}^{\prime }(\mathbf{x}^{\ast }))^{-1}\displaystyle \frac{\Vert \mathbf{e}\Vert }{\Vert \mathbf{e}_{0}\Vert }\leq \displaystyle \frac{\Vert \mathbf{F}(\mathbf{x})\Vert }{\Vert \mathbf{F}(\mathbf{x}_{0})\Vert }\leq \displaystyle \frac{5}{3}\unicode[STIX]{x1D705}(\mathbf{F}^{\prime }(\mathbf{x}^{\ast }))\displaystyle \frac{\Vert \mathbf{e}\Vert }{\Vert \mathbf{e}_{0}\Vert }.\end{eqnarray}$$

There is nothing magic about the numbers $3/5$ and $5/3$ . Both are artifacts of the fraction $1/2$ in (2.3). As $\mathbf{x}$ and $\mathbf{x}_{0}$ approach $\mathbf{x}^{\ast }$ , both coefficients will approach $1$ . So, the inequality (2.9) is satisfyingly consistent with the linear case (where $\unicode[STIX]{x1D6FE}=0$ and $\unicode[STIX]{x1D70C}^{\ast }=\infty$ , so any $\mathbf{x}$ satisfies (2.3)).

Most implementations of Newton’s method do not attempt to compute the step with high accuracy, as we will see in the following sections. Instead one accepts low accuracy in the Jacobian, the linear solve for the step, or even the residual itself. In these cases it is usually unwise to terminate on small steps, and one must terminate on small residuals and accept the effects of ill-conditioning. In the descriptions of algorithms, we will terminate when

(2.10) $$\begin{eqnarray}\Vert \mathbf{F}(\mathbf{x})\Vert \leq \unicode[STIX]{x1D70F}_{a}+\unicode[STIX]{x1D70F}_{r}\Vert \mathbf{F}(\mathbf{x}_{0})\Vert .\end{eqnarray}$$

2.4 Implementation: LU factorization of $\mathbf{F}^{\prime }$

The outline of a Newton iteration is simple. One evaluates the residual, computes the step, and continues until a termination criterion is satisfied. A broad outline of the Newton iteration is shown in algorithm newton.

Algorithm newton leaves out all the important details. We will generally terminate the Newton iteration when the residual norm is small using (2.10). However, that alone is not enough. One must limit the number of iterations to avoid an infinite loop when, for example, the equation has no solution. One must also decide how to solve the linear equation for the step. If the Jacobian $\mathbf{F}^{\prime }$ is small, dense and unstructured, the natural implementation of algorithm newton is to use Gaussian elimination and compute an LU factorization of $\mathbf{F}^{\prime }$ . The resulting algorithm, newton_LU, is now quite specific and the reader should be able to implement it easily.

Algorithm newton_LU works well and is widely used. However, there is more to consider. One important matter is how one computes $\mathbf{F}^{\prime }$ . The best way, if one can do it, is to compute the Jacobian analytically. Analytic Jacobians are usually less expensive computationally and avoid any possible problems with differencing. Computing analytic Jacobians is not possible in a general-purpose code, however, and a typical approach is to compute $\mathbf{F}^{\prime }$ with a forward difference. One way to do this is to approximate the $j$ th column of $\mathbf{F}^{\prime }(\mathbf{x})$ with the difference

(2.11) $$\begin{eqnarray}\displaystyle \frac{\mathbf{F}(\mathbf{x}+h\mathbf{u}_{j})-\mathbf{F}(\mathbf{x})}{h},\end{eqnarray}$$

where $\mathbf{u}_{j}$ is the unit vector in the $j$ th coordinate direction. A finite-difference Jacobian $\unicode[STIX]{x1D6FB}_{h}\mathbf{F}(\mathbf{x})$ , therefore, has an $O(h)$ error. The computational cost is $N$ additional function evaluations, one for each direction. If the cost of an evaluation of $\mathbf{F}$ is $O(N^{2})$ , as it would be for a linear equation, then the cost of computing the finite-difference Jacobian is $O(N^{3})$ , which is the same order as LU factorization. Hence, the construction and factorization of the Jacobian can dominate the cost of the solve. We will show later that analytic Jacobians and finite-difference Jacobians usually produce very similar Newton iterations. The difference is in the expense of computing the finite-difference Jacobian. If the Jacobian is sparse, there are ways to do the differencing much more efficiently. Curtis, Powell and Reid (Reference Curtis, Powell and Reid1974) describe one such method.

The chord method modifies the Newton iteration by moving the Jacobian evaluation and factorization out of the loop. The coefficient matrix in the linear equation for the step is $\mathbf{F}(\mathbf{x}_{0})$ for every nonlinear iteration. For example, algorithm newton_LU becomes an implementation of the chord method by moving one line, shown as algorithm chord_LU.

The chord iteration is

$$\begin{eqnarray}\mathbf{x}_{+}=\mathbf{x}_{c}-\mathbf{F}^{\prime }(\mathbf{x}_{0})^{-1}\mathbf{F}(\mathbf{x}_{c}),\end{eqnarray}$$

so one has simply replaced the Jacobian at the current point with an approximation. The error in that approximation is

$$\begin{eqnarray}\Vert \mathbf{F}^{\prime }(\mathbf{x}_{c})-\mathbf{F}^{\prime }(\mathbf{x}_{0})\Vert \leq \unicode[STIX]{x1D6FE}\Vert \mathbf{x}_{c}-\mathbf{x}_{0}\Vert =\unicode[STIX]{x1D6FE}\Vert \mathbf{e}_{c}-\mathbf{e}_{0}\Vert .\end{eqnarray}$$

One can prove local convergence with the standard assumptions if the initial iterate is sufficiently near $\mathbf{x}^{\ast }$ . However, (2.3) may not be good enough. The next section looks at a longer list of approximations one can apply to Newton’s method and their effects on the iteration.

2.5 Errors in $\mathbf{F}$ and $\mathbf{F}^{\prime }$

Suppose one approximates Newton’s method by

$$\begin{eqnarray}\mathbf{x}_{+}=\mathbf{x}_{c}+\mathbf{s},\end{eqnarray}$$

where

(2.12) $$\begin{eqnarray}\Vert \mathbf{J}_{c}\mathbf{s}+(\mathbf{F}(x_{c})+\unicode[STIX]{x1D716}(\mathbf{x}_{c}))\Vert \leq \unicode[STIX]{x1D702}_{c}\Vert \mathbf{F}(x_{c})+\unicode[STIX]{x1D716}(\mathbf{x}_{c})\Vert\end{eqnarray}$$

and

(2.13) $$\begin{eqnarray}\Vert \mathbf{J}_{c}-\mathbf{F}^{\prime }(\mathbf{x}_{c})\Vert \leq \unicode[STIX]{x1D6E5}_{c}.\end{eqnarray}$$

We allow for errors in every possible way in this approximation. The evaluation of $\mathbf{F}$ has an error $\unicode[STIX]{x1D716}$ . We have an approximate Jacobian $\mathbf{J}$ for the linear equation for the Newton step. Finally, we do not even solve that incorrect equation for the step exactly, rather we take a step $\mathbf{s}$ which satisfies the inexact Newton condition

(2.14) $$\begin{eqnarray}\Vert \mathbf{J}_{c}\mathbf{s}+\mathbf{F}(\mathbf{x}_{c})+\unicode[STIX]{x1D716}(\mathbf{x}_{c})\Vert \leq \unicode[STIX]{x1D702}_{c}\Vert \mathbf{F}(\mathbf{x}_{c})+\unicode[STIX]{x1D716}(\mathbf{x}_{c})\Vert\end{eqnarray}$$

(Dembo, Eisenstat and Steihaug Reference Dembo, Eisenstat and Steihaug1982). One way to interpret the inexact Newton condition is as the termination criterion for an iterative linear solver (small relative residuals).

One should expect the iteration to converge if the standard assumptions hold and the errors are sufficiently small. Theorem 2.3 (Kelley Reference Kelley1995) quantifies that.

Theorem 2.3. Let Assumption 2.1, (2.3), and (2.14) hold. Then

(2.15) $$\begin{eqnarray}\Vert \mathbf{e}_{+}\Vert =O(\Vert \mathbf{e}_{c}\Vert ^{2}+(\Vert \unicode[STIX]{x1D702}_{c}\Vert +\unicode[STIX]{x1D6E5}_{c})\Vert \mathbf{e}_{c}\Vert +\Vert \unicode[STIX]{x1D716}(\mathbf{x}_{c})\Vert ).\end{eqnarray}$$

Theorem 2.3 is very satisfying and explains most of the algorithms used in practice. We will now apply (2.15) to several examples.

For the chord method, $\unicode[STIX]{x1D716}=0$ , $\unicode[STIX]{x1D702}=0$ and $\mathbf{J}=\mathbf{F}^{\prime }(\mathbf{x}_{0})$ . Hence

$$\begin{eqnarray}\unicode[STIX]{x1D6E5}_{c}=\Vert \mathbf{F}^{\prime }(\mathbf{x}_{0})-\mathbf{F}^{\prime }(\mathbf{x}_{c})\Vert \leq \unicode[STIX]{x1D6FE}\Vert \mathbf{x}_{0}-\mathbf{x}_{c}\Vert \leq \unicode[STIX]{x1D6FE}(\Vert \mathbf{e}_{0}\Vert +\Vert \mathbf{e}_{c}\Vert ).\end{eqnarray}$$

If (2.3) holds, then $\Vert \mathbf{e}_{1}\Vert \leq \Vert \mathbf{e}_{0}\Vert /2$ , because the first chord iteration is a Newton iteration. However, one needs a better initial iterate to compensate for the error in the Jacobian. If the initial iterate is sufficiently good, then

$$\begin{eqnarray}\Vert \mathbf{e}_{n+1}\Vert =O(\Vert \mathbf{e}_{n}\Vert ^{2}+\Vert \mathbf{e}_{0}\Vert \Vert \mathbf{e}_{n}\Vert )=O(\Vert \mathbf{e}_{0}\Vert \Vert \mathbf{e}_{n}\Vert )<\Vert \mathbf{e}_{n}\Vert .\end{eqnarray}$$

Hence, the convergence of the chord method is not q-quadratic, but rather q-linear, with a q-factor proportional to $\Vert \mathbf{e}_{0}\Vert$ .

Theorem 2.3 is also the tool one needs to understand the effects of approximating the Jacobian with finite differences. For this we will assume that $\unicode[STIX]{x1D716}$ is independent of $\mathbf{x}$ , $\unicode[STIX]{x1D702}=0$ , and the error in the Jacobian is first-order in the difference increment. The statement

$$\begin{eqnarray}\Vert \unicode[STIX]{x1D6FB}_{h}\mathbf{F}(\mathbf{x})-\mathbf{F}^{\prime }(\mathbf{x})\Vert =O(h)\end{eqnarray}$$

hides the prefactor of $\unicode[STIX]{x1D6FE}/2$ in the $O$ -term. If $\unicode[STIX]{x1D6FE}$ , the Lipschitz constant of $\mathbf{F}^{\prime }$ , is not too large, then a finite-difference Jacobian can be used safely. Most of the time this approximation is fine, but there are exceptions (Kerkhoven and Jerome Reference Kerkhoven and Jerome1990, Coughran and Jerome Reference Coughran, Jerome and Bank1990).

Theorem 2.3 tells us that in this case

$$\begin{eqnarray}\Vert \mathbf{e}_{+}\Vert =O(\Vert \mathbf{e}_{c}\Vert ^{2}+h\Vert \mathbf{e}_{c}\Vert +\Vert \unicode[STIX]{x1D716}\Vert ).\end{eqnarray}$$

The estimate implies that we cannot hope to reduce $\Vert \mathbf{e}\Vert$ to any less than $\Vert \unicode[STIX]{x1D716}\Vert$ and that the error terms balance when $h=O(\sqrt{\Vert \unicode[STIX]{x1D716}\Vert })$ (which is a standard lesson in numerical analysis about finite-difference derivatives: Kelley Reference Kelley1995). The more subtle message in the estimate is that if $h=O(\sqrt{\Vert \unicode[STIX]{x1D716}\Vert })$ , the iteration is indistinguishable from the Newton iteration with an exact derivative until $\Vert \mathbf{e}\Vert \approx \sqrt{\Vert \unicode[STIX]{x1D716}\Vert }$ .

As a final example, we will consider the secant method for scalar equations (i.e. equations for one variable) $f(x)=0$ . Here the model derivative is

$$\begin{eqnarray}j_{c}=\displaystyle \frac{f(x_{c})-f(x_{-})}{x_{c}-x_{-}},\end{eqnarray}$$

where $x_{-}$ is the iterate before $x_{c}$ . One must, of course, decide what $x_{-1}$ should be. One good choice is $1.01\times x_{0}$ , which we will use in the examples. Similarly to the analysis of the chord method, we have

$$\begin{eqnarray}|j_{c}-f^{\prime }(x_{c})|=O(|e_{c}|+|e_{-}|).\end{eqnarray}$$

Theorem 2.3 with $\unicode[STIX]{x1D716}=0$ and $\unicode[STIX]{x1D702}=0$ says that if $|e_{c}|\leq |e_{-}|$ are sufficiently small,

(2.16) $$\begin{eqnarray}|e_{+}|=O(|e_{c}||e_{-}|).\end{eqnarray}$$

So, if the initial iterations are sufficiently good, the secant iteration converges and $|e_{n+1}|=O(|e_{n}||e_{n-1}|)$ , which implies that the convergence is q-superlinear, that is,

(2.17) $$\begin{eqnarray}\lim _{n\rightarrow \infty }\displaystyle \frac{|e_{n+1}|}{|e_{n}|}=0.\end{eqnarray}$$

The secant method is limited to scalar equations. In fact, the secant method converges q-superlinearly with q-order $\unicode[STIX]{x1D6FC}=(1+\sqrt{5})/2$ , that is,

$$\begin{eqnarray}|e_{+}|=O(|e_{c}|^{\unicode[STIX]{x1D6FC}}).\end{eqnarray}$$

The quasi-Newton methods, which we discuss in Section 2.8, extend the secant method to systems of equations. Scalar equations are no longer an active topic of research, but at one time there was considerable activity (Traub Reference Traub1964). We will have very little to say about them in this article aside from a few examples.

We can illustrate these results with a simple example. The residual for the scalar equation

$$\begin{eqnarray}f(x)=x-\text{e}^{-x}\cos (x)=0\end{eqnarray}$$

can be evaluated to about $15$ figures of accuracy. The Newton iteration is

$$\begin{eqnarray}x_{+}=x_{c}-f(x_{c})/(1+\text{e}^{-x_{c}}(\sin (x_{c})+\cos (x_{c}))).\end{eqnarray}$$

Using the initial iterate $x_{0}=1$ , we compare the iteration histories of Newton’s method, Newton’s method with a forward difference derivative, and the chord method. We will use both a table and a plot, and invite the reader to pick her or his own favourite way to present such data.

To illustrate the effects of the error in $f$ , we tabulate and plot $|f|$ for the three methods. We set the termination criteria to $\unicode[STIX]{x1D70F}_{a}=\unicode[STIX]{x1D70F}_{r}=10^{-20}$ . Since $\unicode[STIX]{x1D716}\approx 10^{-15}$ in this example, the theory does not imply that it is possible to drive the residual to a value as small as $10^{-20}$ , and the computation confirms that. We used a difference increment of $h=10^{-7}$ .

Table 2.1.

Iteration histories: $|f(x_{n})|$ .

Those of you accustomed to looking at columns of figures may have noticed that, as the theory predicts, there is very little difference between the finite-difference Newton method and the version with analytic derivatives until the iteration stagnates at roughly the level of machine precision. One can also notice that the residuals for the chord method decay more slowly, by a factor of 4–5 with each iteration. The secant method converges faster than chord, but not as fast as Newton’s method. One weakness of the secant method is exposed by the $10$ th iteration, where there is a floating-point exception. The problem is that the $f(x_{n})=f(x_{n+1})$ and $x_{n}=x_{n+1}$ at this point, so one gets $j_{c}=0/0$ , which is reported in IEEE arithmetic (Overton Reference Overton2001) as NaN (Not a Number). In this example we have run the iteration far beyond any sensible termination point. The NaN and the stagnation in the Newton iteration are signs of that.

It is more illuminating, at least in the author’s opinion, to visualize iteration histories, and Table 2.1 is one of the very few tables we will use for that. In Figure 2.1 we visualize the data from the table in a semi-log plot. It is very clear that there is little difference between the two realizations of Newton’s method in terms of the number of nonlinear iterations needed to converge to the limiting level of precision. One can also see the signature of superlinear convergence in downward concavity of the residual history for Newton’s method and the secant method. The NaN is missing from the plot by convention. The q-linear convergence of the chord method appears as a linear residual history when plotted in this way.

Figure 2.1.

Visual iteration histories.

As a general rule, if residuals are accurate to machine unit roundoff, the iterations will stagnate at that level. There are exceptions. For example, if the floating-point implementation of a function has a root in the floating-point number system, then one can continue the iteration for much longer. The equation $\arctan (x)=0$ is an example of this phenomenon.

Our final example in this section is the Shamanskii method (Shamanskii Reference Shamanskii1967). This method is a hybrid between Newton’s method and the chord method. The idea is to update the Jacobian every $m\geq 1$ iterations. Clearly $m=1$ is Newton’s method and $m=\infty$ is the chord method. If Gaussian elimination is the linear solver, the Shamanskii iteration for finite $m$ is shamanskii_LU.

In algorithm shamanskii_LU the iteration counter is incremented outside of the inner $j$ loop. Keeping this in mind, Theorem 2.3 states that

(2.18) $$\begin{eqnarray}\Vert \mathbf{e}_{n+1}\Vert =O(\Vert \mathbf{e}_{n}\Vert ^{m+1}).\end{eqnarray}$$

If $m>1$ the convergence rate is faster than the q-quadratic rate for Newton’s method. We say the convergence is q-superlinear with q-order $m+1$ . The Shamanskii method is less appealing than it appears because if the Jacobian is sufficiently near $\mathbf{F}^{\prime }(\mathbf{x}^{\ast })$ the modest reduction in the iterations is not worth the cost of computing and factoring the Jacobian. In many cases $m=\infty$ (the chord method) is a better choice (Brent Reference Brent1973).

This section has looked at the effects of $\unicode[STIX]{x1D6E5}$ and $\unicode[STIX]{x1D716}$ on the convergence of Newton’s method. Typically proofs set $\unicode[STIX]{x1D716}=0$ and proceed as if evaluations $\mathbf{F}$ were exact. The users of the results generally know that $\unicode[STIX]{x1D716}\neq 0$ and when the error in the function evaluation is important. We will follow that approach in the next section, where the inexactness in the linear solver, measured by $\unicode[STIX]{x1D702}$ , is the important part of the algorithm.

2.6 Inexact Newton methods and JFNK

When one solves a linear equation $\mathbf{A}\mathbf{x}=\mathbf{b}$ with an iterative method, one usually terminates the iteration when the relative residual $\Vert \mathbf{A}\mathbf{x}-\mathbf{b}\Vert /\Vert \mathbf{b}\Vert$ is sufficiently small. An inexact Newton method (Dembo et al. Reference Dembo, Eisenstat and Steihaug1982) approximates Newton’s method by using a step $\mathbf{s}$ that satisfies the inexact Newton condition

(2.19) $$\begin{eqnarray}\Vert \mathbf{F}^{\prime }(\mathbf{x}_{c})\mathbf{s}+\mathbf{F}(\mathbf{x}_{c})\Vert \leq \unicode[STIX]{x1D702}_{c}\Vert \mathbf{F}(\mathbf{x}_{c})\Vert .\end{eqnarray}$$

Here $\unicode[STIX]{x1D702}$ is called the forcing term.

While the theory does not say how one realizes (2.19), in practice it is usually the outcome of an iterative method for solving $\mathbf{F}^{\prime }(\mathbf{x}_{c})\mathbf{x}=-\mathbf{F}(\mathbf{x}_{c})$ . The term Newton-iterative method is used in that case. The iteration for $\mathbf{x}$ is called the outer or nonlinear iteration. The iterative method for the linear equation is called the inner or linear iteration. A specific linear solver or class of solvers is often specified. For example Newton–Krylov and Newton–GMRES are common choices. Preconditioners can also be part of the name. Newton–Krylov–Schwarz (Cai, Gropp, Keyes and Tidriri Reference Cai, Gropp, Keyes, Tidriri and Rannacher1994) methods use a Krylov linear solver and a Schwarz domain decomposition method as a preconditioner.

A straightforward application of Theorem 2.3 with $\unicode[STIX]{x1D6E5}=0$ and $\unicode[STIX]{x1D716}=0$ leads to

(2.20) $$\begin{eqnarray}\Vert \mathbf{e}_{+}\Vert =O(\Vert \mathbf{e}_{c}\Vert ^{2}+\unicode[STIX]{x1D702}_{c}\Vert \mathbf{e}_{c}\Vert ).\end{eqnarray}$$

The estimate (2.20) leads to a convergence theorem.

There are, as one might suspect from Theorem 2.4, many approaches to managing $\unicode[STIX]{x1D702}$ as the iteration progresses. Similar to the distinction between Newton’s method and the chord method, it is rarely a good idea to make $\unicode[STIX]{x1D702}$ very small, especially at the beginning of the iteration when only limited accuracy is needed to get the same reduction in error that one would get with Newton’s method itself. While one could make $\unicode[STIX]{x1D702}$ small once the residuals are small, it is not clear that the additional cost in the linear solve makes the reduction in nonlinear iterations worthwhile. There are useful discussions of this issue in Kelley (Reference Kelley1995) and Eisenstat and Walker (Reference Eisenstat and Walker1996). The author of this article has had success with $\unicode[STIX]{x1D702}_{n}\equiv 1/10$ .

Theorem 2.4 does not specify any particular norm. If one uses the weighed norm

$$\begin{eqnarray}\Vert \mathbf{x}\Vert _{\ast }=\Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })\mathbf{x}\Vert ,\end{eqnarray}$$

then the theory no longer needs a small $\bar{\unicode[STIX]{x1D702}}$ . Any $\bar{\unicode[STIX]{x1D702}}<1$ will do.

Theorem 2.5. Let Assumption 2.1 hold. Then if $\mathbf{x}_{0}$ is sufficiently near $\mathbf{x}^{\ast }$ and $0\leq \unicode[STIX]{x1D702}_{n}\leq \bar{\unicode[STIX]{x1D702}}<\tilde{\unicode[STIX]{x1D702}}<1$ ,

$$\begin{eqnarray}\Vert \mathbf{e}_{n+1}\Vert _{\ast }\leq \tilde{\unicode[STIX]{x1D702}}\Vert \mathbf{e}_{n}\Vert _{\ast }\end{eqnarray}$$

and the other conclusions of Theorem 2.4 hold.

Assuming that $\mathbf{e}_{n}\neq 0$ for all $n$ , Theorem 2.5 implies that

$$\begin{eqnarray}\limsup \displaystyle \frac{\Vert \mathbf{e}_{n+1}\Vert _{\ast }}{\Vert \mathbf{e}_{n}\Vert _{\ast }}\leq \bar{\unicode[STIX]{x1D702}}.\end{eqnarray}$$

If the linear solver is a Krylov method which only needs Jacobian-vector products, it is not necessary or desirable to compute and store a Jacobian matrix. For example, one can approximate the Jacobian-vector product with a forward difference. Methods of this type are called Jacobian-free Newton–Krylov (JFNK) methods. JFNK methods are the most common choice for those large-scale nonlinear equations which come from differential and integral equations. Knoll and Keyes (Reference Knoll and Keyes2004) provide an excellent account of JNFK methods and applications.

The linear solver in a JFNK method often requires preconditioning to work well enough to be useful. In most implementations preconditioning is done at the level of the linear solver. Preconditioning can also be encoded in the nonlinear map itself and it is interesting to examine that. We will let $\mathbf{M}$ be the preconditioner for the linear equation for the Newton step, and assume that $\mathbf{M}$ does not depend on the $\mathbf{x}$ , the nonlinear iteration. If we precondition from the left, the equation for the Newton step is transformed into

$$\begin{eqnarray}\mathbf{M}\mathbf{F}^{\prime }(\mathbf{x}_{c})\mathbf{s}=-\mathbf{M}\mathbf{F}(\mathbf{x}_{c}).\end{eqnarray}$$

This is exactly the Newton step for the equation

$$\begin{eqnarray}\mathbf{M}\mathbf{F}(\mathbf{x})=0.\end{eqnarray}$$

So, in the case of left preconditioning, one can place the preconditioning in the definition of the nonlinear map, replacing $\mathbf{F}$ by $\mathbf{M}\mathbf{F}$ , or apply it to the linear equation for the step. In either case, the steps and the iterations will be the same as in the unpreconditioned case. Most of the production codes put the preconditioning in the linear solve and measure the unpreconditioned residual $\mathbf{F}$ when terminating the linear iteration. When one does this, however, the termination criterion for the linear iteration would be

$$\begin{eqnarray}\Vert \mathbf{M}\mathbf{F}^{\prime }(\mathbf{x}_{c})\mathbf{s}+\mathbf{M}\mathbf{F}(\mathbf{x}_{c})\Vert \leq \unicode[STIX]{x1D702}\Vert \mathbf{M}\mathbf{F}(\mathbf{x}_{c})\Vert ,\end{eqnarray}$$

which is the inexact Newton condition for $\mathbf{M}\mathbf{F}$ . This does not change the theory if $\unicode[STIX]{x1D702}$ is sufficiently small. To summarize, replacing $\mathbf{F}$ with $\mathbf{M}\mathbf{F}$ does not change the iteration or the steps, but does change the norm of the Jacobian at the solution and the Lipschitz constant for the Jacobian.

Similarly, if one preconditions the linear equation from the right, the equation is

$$\begin{eqnarray}\mathbf{F}^{\prime }(\mathbf{x}_{c})\mathbf{M}\mathbf{z}=-\mathbf{F}(\mathbf{x}_{c}),\quad \mathbf{s}=\mathbf{M}\mathbf{z},\end{eqnarray}$$

and the corresponding nonlinear system is

$$\begin{eqnarray}\mathbf{F}(\mathbf{M}\mathbf{y})=0,\quad \mathbf{x}=\mathbf{M}\mathbf{y}.\end{eqnarray}$$

In this case the residuals are unchanged and the inexact Newton condition has its original meaning if one puts the preconditioning in the linear solver.

2.7 Global convergence

The Newton iteration for $\arctan (x)=0$ with $x_{0}=1$ exhibits classic quadratic convergence. With a poor initial iterate, $x_{0}=10$ for example, the first five iterations are

$$\begin{eqnarray}10,\quad -138,\quad 2.9\times 10^{4},\quad -1.5\times 10^{9},\quad 9.9\times 10^{17}.\end{eqnarray}$$

This divergence is consistent with the theory because the initial iterate is so poor. The Armijo line search (Armijo Reference Armijo1966) is a wonderful and simple solution to this problem.

The idea is to observe that even though the Newton step moves farther from the root, the direction is correct. To clarify this we make a distinction between the Newton direction

$$\begin{eqnarray}\mathbf{d}=-\mathbf{F}^{\prime }(\mathbf{x}_{c})^{-1}\mathbf{F}(\mathbf{x}_{c})\end{eqnarray}$$

and the Newton step

$$\begin{eqnarray}\mathbf{s}=\mathbf{x}_{+}-\mathbf{x}_{c}.\end{eqnarray}$$

The simplest strategy is to find the least $\unicode[STIX]{x1D706}=2^{-m}$ for $m=0,1,\ldots \,$ so that

(2.21) $$\begin{eqnarray}\Vert \mathbf{F}(\mathbf{x}_{c}+\unicode[STIX]{x1D706}\mathbf{d})\Vert <\Vert \mathbf{F}(\mathbf{x}_{c})\Vert\end{eqnarray}$$

and use $\mathbf{s}=\unicode[STIX]{x1D706}\mathbf{d}$ . The simple decrease condition (2.21) is a bit too weak for a convergence analysis, but is close enough to save the $\arctan$ iteration, as one can see in Figure 2.2.

Figure 2.2.

The Armijo line search for $\arctan (x)=0$ , $x_{0}=10$ .

In the figure the circles are iterations for which the step length required reduction, and the number of stepsize reductions are indicated next to the circles. After the fourth iteration the iteration took full steps, and the local theory holds. Note that the residual reduction continued past the resolution of the floating-point system. The reason for this, as we mentioned in Section 2.5, is that $x^{\ast }=0$ is also a solution of the equation as implemented in MATLAB with IEEE floating-point arithmetic.

Convergence theory requires more than simple decrease as in (2.21). The sufficient decrease condition that one needs for theory is

(2.22) $$\begin{eqnarray}\Vert \mathbf{F}(\mathbf{x}_{c}+2^{-m}\mathbf{d})\Vert <(1-\unicode[STIX]{x1D6FC}2^{-m})\Vert \mathbf{F}(\mathbf{x}_{c})\Vert .\end{eqnarray}$$

In most of the literature and codes, $\unicode[STIX]{x1D6FC}=10^{-4}$ . In most cases sufficient decrease conditions such as (2.22) will lead to the same stepsize decisions as simple decrease, and the purpose is to enable theory. It is common to adaptively change the stepsize reduction factor. One way to do this (Kelley Reference Kelley1995, Dennis and Schnabel Reference Dennis and Schnabel1996) is to begin with a stepsize of $\unicode[STIX]{x1D706}=1$ , and with each failure of the sufficient decrease condition

(2.23) $$\begin{eqnarray}\Vert \mathbf{F}(\mathbf{x}_{c}+\unicode[STIX]{x1D706}\mathbf{d})\Vert <(1-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D706})\Vert \mathbf{F}(\mathbf{x}_{c})\Vert ,\end{eqnarray}$$

reduce $\unicode[STIX]{x1D706}$ by a factor $\unicode[STIX]{x1D70E}\in [0.1,0.5]$ . The standard way to do this is to use this history of failed steps to generate a polynomial approximation of $\unicode[STIX]{x1D719}(\unicode[STIX]{x1D706})=\Vert \mathbf{F}(\mathbf{x}_{c}+\unicode[STIX]{x1D706}\mathbf{d})\Vert$ and minimize that polynomial. This adaptivity is very useful in practice.

Algorithm newton_armijo is an inexact formulation of the Newton–Armijo method. It includes the essential ideas and makes the theory easy to state. It is critical that one uses $\mathbf{F}^{\prime }(\mathbf{x})$ and not an approximation. The Armijo rule is not, for example, theoretically supported for the chord method. Adding a line search to a Newton code is easy and requires only a few new lines.

Algorithm newton_armijo does not say how or if the forcing term $\unicode[STIX]{x1D702}$ changes with the iteration. Most of the codes do vary $\unicode[STIX]{x1D702}$ and use the ideas in Eisenstat and Walker (Reference Eisenstat and Walker1996). One example from Eisenstat and Walker (Reference Eisenstat and Walker1996) which is common is

(2.24) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{n}=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D702}_{\max }\quad & n=0,\\ \min (\unicode[STIX]{x1D702}_{\max },\unicode[STIX]{x1D702}_{n}^{A})\quad & n>0.\end{array}\right.\end{eqnarray}$$

Here

$$\begin{eqnarray}\unicode[STIX]{x1D702}_{n}^{A}=\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}}\Vert F(x_{n})\Vert ^{2}/\Vert F(x_{n-1})\Vert ^{2}\end{eqnarray}$$

and $\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}}$ is a parameter. The parameter $\unicode[STIX]{x1D702}_{\max }$ is an upper limit on the sequence $\{\unicode[STIX]{x1D702}_{n}\}$ . Eisenstat and Walker (Reference Eisenstat and Walker1996) use the choices $\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}}=0.9$ and $\unicode[STIX]{x1D702}_{\max }=0.9999$ . The author of this article, however, likes the choice $\unicode[STIX]{x1D702}_{n}\equiv 0.1$ .

Theorem 2.6 is very satisfying. For sufficiently smooth problems, the Newton–Armijo iteration has only three possible outcomes. One is convergence to a solution which satisfies the standard assumptions. In that case the stepsize $\unicode[STIX]{x1D706}$ will be one in the terminal phase of the iteration. The other two outcomes are failures, which are easy to detect numerically. One failure mode is that the iteration becomes unbounded. An example of such a problem is the scalar equation $e^{x}=0$ , where the Newton–Armijo iteration diverges to $-\infty$ . The second failure mode is that the Jacobian drifts to singularity. An example is the scalar equation $x^{2}+1=0$ .

The Newton–Armijo method does not solve all problems. Even in the successful case, there is no guarantee that the iteration converges to a useful solution. Nonlinear equations can have multiple solutions, and there are often constraints such as dynamic stability or correct signs for physical quantities to which the Newton iteration is oblivious. We will consider a few ways to address non-uniqueness in Sections 4 and 7.

Another approach to globalization is that of trust region methods (Powell Reference Powell and Rabinowitz1970). These methods are widely used in optimization (Conn, Gould and Toint Reference Conn, Gould and Toint2000), but less so for nonlinear equations. The idea is to model $\Vert \mathbf{F}(\mathbf{x}_{c})\Vert _{2}^{2}$ with a quadratic and minimize that quadratic in a bounded set, the trust region, centred at $x_{c}$ . We will not discuss these methods in detail here. Absil, Baker and Gallivan (Reference Absil, Baker and Gallivan2007) and Higham (Reference Higham1999) report on some interesting applications which are connected to the continuation methods in Section 7.

2.8 Broyden’s method

Quasi-Newton methods construct a model Jacobian from the history of the iteration. One maintains both an approximation $\mathbf{x}_{n}$ of the solution and an approximation $\mathbf{B}_{n}$ of the Jacobian. There are many of these methods (Kelley Reference Kelley1995, Dennis and Schnabel Reference Dennis and Schnabel1996, Dennis and Walker Reference Dennis and Walker1981) and they are widely used in optimization. They have largely been replaced by JFNK methods for nonlinear equations. We will briefly discuss Broyden’s method, the simplest of them and the one that is used in the NOX code from Trilinos (Heroux et al. Reference Heroux, Bartlett, Howle, Hoekstra, Hu, Kolda, Lehoucq, Long, Pawlowski, Phipps, Salinger, Thornquist, Tuminaro, Willenbring, Williams and Stanley2005).

The Broyden update is

(2.25) $$\begin{eqnarray}\mathbf{B}_{+}=\mathbf{B}_{c}+{\displaystyle \frac{(\mathbf{y}-\mathbf{B}_{c}\mathbf{s})\mathbf{s}^{T}}{\mathbf{s}^{T}\mathbf{s}}}.\end{eqnarray}$$

Here $\mathbf{y}=\mathbf{F}(\mathbf{x}_{+})-\mathbf{F}(\mathbf{x}_{c})$ and $\mathbf{s}=\mathbf{x}_{+}-\mathbf{x}_{c}$ . One can think of this as a generalization of the secant method for scalar equations. For the secant equation, the model derivative at the new point $x_{+}$ is

$$\begin{eqnarray}b_{+}=\displaystyle \frac{f(x_{+})-f(x_{c})}{x_{+}-x_{c}}=\displaystyle \frac{y}{s},\end{eqnarray}$$

giving $b_{+}s=y$ . For systems of equations, the secant equation

$$\begin{eqnarray}\mathbf{B}_{+}\mathbf{s}=\mathbf{y}\end{eqnarray}$$

is a system of $N$ equations in $N^{2}$ unknowns. This enables construction of updates that satisfy structural constraints such as sparsity or positivity. Dennis and Schnabel (Reference Dennis and Schnabel1996) and Kelley (Reference Kelley1995) discuss several kinds of secant updates for nonlinear equations. In optimization, for example, one often wants the model Hessian to be symmetric and positive definite. There are many quasi-Newton update which do that (Dennis and Schnabel Reference Dennis and Schnabel1996, Kelley Reference Kelley1999, Nocedal and Wright Reference Nocedal and Wright1999). One can also design updates to capture functional analytic properties of infinite-dimensional problems and their discretizations (Kelley and Sachs Reference Kelley and Sachs1993, Kelley and Sachs Reference Kelley and Sachs1995, Kelley and Sachs Reference Kelley and Sachs1987, Kelley and Sachs Reference Kelley and Sachs1989, Hart and Soul Reference Hart and Soul1973). The most general accounts of theory can be found in Dennis and Walker (Reference Dennis and Walker1981) and Dennis and Schnabel (Reference Dennis and Schnabel1979).

One can implement the Broyden update by storing two vectors for each iteration and using the Sherman–Morrison formula (Kelley Reference Kelley1995) to update the product of $\mathbf{B}^{-1}$ and a vector. The storage burden can be reduced to one vector per nonlinear iteration by using the dependence of $\mathbf{y}$ on $\mathbf{s}$ (Kelley Reference Kelley1995, Deuflhard, Freund and Walter Reference Deuflhard, Freund and Walter1990). JFNK methods place the storage burden on the linear iteration, and that seems to be best.

The formula for the update (2.25) allows for a line search in which

$$\begin{eqnarray}\mathbf{s}=-\unicode[STIX]{x1D706}\mathbf{B}_{c}^{-1}\mathbf{F}(\mathbf{x}_{c}).\end{eqnarray}$$

Theorem 2.6 does not apply to a Broyden–Armijo algorithm patterned after algorithm newton_armijo, but in practice such an algorithm often works fine.

The convergence theory is only local and requires the standard assumptions and accurate initial approximations to the solution and the Jacobian at the solution.

Theorem 2.7. If the standard assumptions hold and $\mathbf{x}_{0}$ and $\mathbf{B}_{0}$ are sufficiently near $\mathbf{x}^{\ast }$ and $\mathbf{F}^{\prime }(\mathbf{x}^{\ast })$ , then $\mathbf{x}_{n}\rightarrow \mathbf{x}^{\ast }$ q-superlinearly:

$$\begin{eqnarray}\lim _{n\rightarrow \infty }\displaystyle \frac{\Vert \mathbf{e}_{n+1}\Vert }{\Vert \mathbf{e}_{n}\Vert }=0.\end{eqnarray}$$

Broyden’s method and JFNK methods have remarkable similarities when applied to discretizations of infinite-dimensional problems. The preconditioning issues are very closely related; see, for example, Kelley and Sachs (Reference Kelley and Sachs1985), Kelley and Xue (Reference Kelley and Xue1996), Burmeister (Reference Burmeister1975) and Nevanlinna (Reference Nevanlinna1993).

2.9 Fréchet and Gâteaux derivatives

We have not fully explained how the results in this section map to the infinite-dimensional setting. The reason for this is that we have defined the derivative in the context of its matrix representation as a Jacobian matrix. We are now at a point where we must make a coordinate-free definition to consider our first infinite-dimensional example in the next section.

If $D\subset \mathbb{R}^{M}$ , we will let $D^{\circ }$ denote the interior of $D$ .

Definition 2.8. A function $\mathbf{F}:D\subset \mathbb{R}^{N}\rightarrow \mathbb{R}^{M}$ is Fréchet differentiable at $x\in D^{\circ }$ if there is a linear map $\mathbf{F}^{\prime }(\mathbf{x})$ from $\mathbb{R}^{M}$ to $\mathbb{R}^{N}$ such that

(2.26) $$\begin{eqnarray}\lim _{\mathbf{h}\rightarrow 0}\displaystyle \frac{\Vert \mathbf{F}(\mathbf{x}+\mathbf{h})-\mathbf{F}(\mathbf{x})-\mathbf{F}^{\prime }(\mathbf{x})\mathbf{h}\Vert }{\Vert \mathbf{h}\Vert }=0.\end{eqnarray}$$

$\mathbf{F}^{\prime }$ is called the Fréchet derivative (or simply the derivative) of $\mathbf{F}$ at $\mathbf{x}$ .

Fréchet differentiability simply means that the difference quotients converge to the derivative uniformly in the direction $\mathbf{h}/\Vert \mathbf{h}\Vert$ as $\Vert \mathbf{h}\Vert \rightarrow 0$ . The Jacobian matrix is simply the matrix representation of the Fréchet derivative in the basis of coordinate directions. So the $j$ th column of the Jacobian is $\mathbf{F}^{\prime }(\mathbf{x})\mathbf{u}_{j}$ as expressed in (2.11). The results on convergence of Newton’s method in this section do not change in the infinite-dimensional case. For quasi-Newton methods, such as Broyden’s method, superlinear convergence depends on compactness properties (Kelley and Sachs Reference Kelley and Sachs1985) of $\mathbf{I}-{\mathcal{F}}^{\prime }$ , which are trivially satisfied in finite-dimensional problems.

One useful way to compute the Fréchet derivative is to apply (2.26) to compute $\mathbf{F}^{\prime }(\mathbf{x})\mathbf{u}$ for an arbitrary $\mathbf{u}$ . Often one can easily extract $\mathbf{F}^{\prime }$ directly by looking at the results, since

(2.27) $$\begin{eqnarray}\mathbf{F}^{\prime }(\mathbf{x})\mathbf{u}=\displaystyle \frac{\text{d}}{\,\text{d}t}\mathbf{F}(\mathbf{x}+t\mathbf{u})\bigg|_{t=0}.\end{eqnarray}$$

As an example, suppose that

$$\begin{eqnarray}\mathbf{F}(\mathbf{x})=\mathbf{A}\mathbf{x}+f(\mathbf{x}),\end{eqnarray}$$

where $\mathbf{A}$ is a linear operator and $f(\mathbf{x})$ is the substitution operator

$$\begin{eqnarray}f(\mathbf{x})\equiv \left(\begin{array}{@{}c@{}}f(x_{1})\\ f(x_{2})\\ \vdots \\ f(x_{N})\end{array}\right)\!.\end{eqnarray}$$

Then (2.27) implies that

$$\begin{eqnarray}\mathbf{F}^{\prime }(\mathbf{x})=\mathbf{A}+\text{diag}(f^{\prime }(\mathbf{x})).\end{eqnarray}$$

If it exists, the limit

(2.28) $$\begin{eqnarray}\,\text{d}\mathbf{F}(\mathbf{x}:\mathbf{u})=\lim _{t\downarrow 0}\displaystyle \frac{\mathbf{F}(\mathbf{x}+t\mathbf{u})-\mathbf{F}(\mathbf{x})}{t}\end{eqnarray}$$

is called the directional derivative of $\mathbf{F}$ at $\mathbf{x}$ in the direction $\mathbf{u}$ . The scalar example $f(x)=|x|$ at $x=0$ has directional derivatives in all directions, but is not Fréchet differentiable. This example also shows that $\text{d}\mathbf{F}(\mathbf{x}:\mathbf{u})$ need not be linear in $\mathbf{u}$ . If the limit

(2.29) $$\begin{eqnarray}\,\text{d}\mathbf{F}(\mathbf{x}:\mathbf{u})=\lim _{t\rightarrow 0}\displaystyle \frac{\mathbf{F}(\mathbf{x}+t\mathbf{u})-\mathbf{F}(\mathbf{x})}{t}\end{eqnarray}$$

exists, then $\text{d}\mathbf{F}(\mathbf{x}:\mathbf{u})$ is called the Gâteaux derivative of $\mathbf{F}$ at $\mathbf{x}$ in the direction $\mathbf{u}$ . In this case

$$\begin{eqnarray}\,\text{d}\mathbf{F}(\mathbf{x}:\mathbf{u})=\displaystyle \frac{\text{d}}{\,\text{d}t}\mathbf{F}(\mathbf{x}+t\mathbf{u})\bigg|_{t=0},\end{eqnarray}$$

but still need not be linear in $\mathbf{u}$ . With the exception of Sections 6 and 7.1, all functions in this article will be Fréchet differentiable. If $\mathbf{F}$ is Fréchet differentiable at $\mathbf{x}$ , then $\text{d}\mathbf{F}(\mathbf{x}:\mathbf{u})=\mathbf{F}^{\prime }(\mathbf{x})\mathbf{u}$ for all $\mathbf{u}$ .

2.10 Example: Chandrasekhar $H$ -equation

This is the first example of an operator in a function space. We will denote such operators by script upper-case letters to distinguish them from their discretizations. The discretizations will be functions on $\mathbb{R}^{N}$ and will have the boldface fonts we have been using.

The Chandrasekhar $H$ -equation (Chandrasekhar Reference Chandrasekhar1960, Busbridge Reference Busbridge1960) is

(2.30) $$\begin{eqnarray}{\mathcal{F}}(H)(\unicode[STIX]{x1D707})=H(\unicode[STIX]{x1D707})-\biggl(1-{\displaystyle \frac{\unicode[STIX]{x1D714}}{2}}\displaystyle \int _{0}^{1}{\displaystyle \frac{\unicode[STIX]{x1D707}H(\unicode[STIX]{x1D708})\,\text{d}\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D707}+\unicode[STIX]{x1D708}}}\biggr)^{-1}=0.\end{eqnarray}$$

The equation arises in radiative transfer theory and is a very tractable example of a nonlinear equation in a function space. We will regard ${\mathcal{F}}$ as a map on $C[0,1]$ , the Banach space of continuous function on the interval $[0,1]$ with the $\Vert \cdot \Vert _{\infty }$ -norm. The unknown is a function $H\in C[0,1]$ ; $\unicode[STIX]{x1D714}\in [0,1]$ is a parameter.

We will discretize the equation with the composite midpoint rule

$$\begin{eqnarray}\displaystyle \int _{0}^{1}f(\unicode[STIX]{x1D707})\,\text{d}\unicode[STIX]{x1D707}\approx \displaystyle \frac{1}{N}\mathop{\sum }_{j=1}^{N}f(\unicode[STIX]{x1D707}_{j}),\end{eqnarray}$$

where $\unicode[STIX]{x1D707}_{i}=(i-1/2)/N$ for $1\leq i\leq N$ .

This leads to an equation in $\mathbb{R}^{N}$ :

(2.31) $$\begin{eqnarray}\mathbf{F}^{N}(\mathbf{h})_{i}=h_{i}-\biggl(1-{\displaystyle \frac{\unicode[STIX]{x1D714}}{2N}}\mathop{\sum }_{j=1}^{N}{\displaystyle \frac{\unicode[STIX]{x1D707}_{i}h_{j}}{\unicode[STIX]{x1D707}_{i}+\unicode[STIX]{x1D707}_{j}}}\biggr)^{-1}.\end{eqnarray}$$

We will discuss some important properties of this equation in the context of the infinite-dimensional problem. Our conclusions will be equally valid for the discrete problem because the midpoint rule integrates constant functions exactly.

2.10.1 The continuous problem

Our first task is to differentiate ${\mathcal{F}}$ . Applying (2.26), we see that for all $v\in C[0,1]$ ,

$$\begin{eqnarray}({\mathcal{F}}^{\prime }(H)w)=w-\displaystyle \frac{\unicode[STIX]{x1D714}Lw}{(1-\unicode[STIX]{x1D714}LH)^{2}},\end{eqnarray}$$

where $L$ is the integral operator defined by

$$\begin{eqnarray}Lw(\unicode[STIX]{x1D707})=\displaystyle \frac{1}{2}\int _{0}^{1}{\displaystyle \frac{\unicode[STIX]{x1D707}w(\unicode[STIX]{x1D708})\,\text{d}\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D707}+\unicode[STIX]{x1D708}}}.\end{eqnarray}$$

Hence,

$$\begin{eqnarray}{\mathcal{F}}^{\prime }(H)=I-\displaystyle \frac{\unicode[STIX]{x1D714}L}{(1-\unicode[STIX]{x1D714}LH)^{2}}.\end{eqnarray}$$

We will now check the standard assumptions. To see that there is a solution, we note that the sequence

$$\begin{eqnarray}H_{0}=1,\quad H_{n+1}=1+\unicode[STIX]{x1D714}H_{n}LH_{n}\end{eqnarray}$$

is increasing. We will show convergence by showing that the $L^{1}$ -norms are bounded:

(2.32) $$\begin{eqnarray}\Vert H_{n}\Vert _{1}=\int _{0}^{1}H_{n}(\unicode[STIX]{x1D708})\,\text{d}\unicode[STIX]{x1D708}\leq \unicode[STIX]{x1D6FD}\equiv \displaystyle \frac{1+\sqrt{1-\unicode[STIX]{x1D714}}}{2/\unicode[STIX]{x1D714}}.\end{eqnarray}$$

Clearly (2.32) holds for $n=1$ . To proceed by induction, suppose $\Vert H_{n}\Vert _{1}\leq \unicode[STIX]{x1D6FD}$ . We seek to show that

(2.33) $$\begin{eqnarray}\int _{0}^{1}H_{n+1}(\unicode[STIX]{x1D707})\,\text{d}\unicode[STIX]{x1D707}=1+\displaystyle \frac{\unicode[STIX]{x1D714}}{2}\int _{0}^{1}\displaystyle \frac{\unicode[STIX]{x1D707}H_{n}(\unicode[STIX]{x1D708})H_{n}(\unicode[STIX]{x1D707})}{\unicode[STIX]{x1D707}+\unicode[STIX]{x1D708}}\,\text{d}\unicode[STIX]{x1D708}\end{eqnarray}$$

for all $0\leq \unicode[STIX]{x1D707}\leq 1$ . The trick is to integrate both sides of (2.33) with respect to $\unicode[STIX]{x1D707}$ and note that

$$\begin{eqnarray}\displaystyle \int _{0}^{1}\int _{0}^{1}\displaystyle \frac{\unicode[STIX]{x1D707}H_{n}(\unicode[STIX]{x1D708})H_{n}(\unicode[STIX]{x1D707})}{\unicode[STIX]{x1D707}+\unicode[STIX]{x1D708}}\,\text{d}\unicode[STIX]{x1D708}\,\text{d}\unicode[STIX]{x1D707} & = & \displaystyle \int _{0}^{1}\int _{0}^{1}\displaystyle \frac{\unicode[STIX]{x1D708}H_{n}(\unicode[STIX]{x1D708})H_{n}(\unicode[STIX]{x1D707})}{\unicode[STIX]{x1D707}+\unicode[STIX]{x1D708}}\,\text{d}\unicode[STIX]{x1D708}\,\text{d}\unicode[STIX]{x1D707}\nonumber\\ \displaystyle & = & \displaystyle \Vert H_{n}\Vert _{1}^{2}/2\leq \unicode[STIX]{x1D6FD}^{2}/2.\nonumber\end{eqnarray}$$

We are done since

$$\begin{eqnarray}\unicode[STIX]{x1D714}\unicode[STIX]{x1D6FD}^{2}/4+1\leq \unicode[STIX]{x1D6FD}\end{eqnarray}$$

by the quadratic formula. So our sequence $H_{n}$ converges in $L^{1}$ to a function $H$ which satisfies $H=1+\unicode[STIX]{x1D714}HLH$ . This implies that $H$ satisfies (2.30), since if $1-\unicode[STIX]{x1D714}LH$ vanishes at any $\unicode[STIX]{x1D707}\in [0,1]$ , then $H(\unicode[STIX]{x1D707})=1+H(\unicode[STIX]{x1D707})$ implying that $H(\unicode[STIX]{x1D707})=0$ , which violates the equation since $H(\unicode[STIX]{x1D707})\geq 1$ . Also, since $H_{n}$ is an increasing sequence, $H$ is a positive function of $\unicode[STIX]{x1D707}$ for $\unicode[STIX]{x1D714}\in [0,1]$ .

So $H\in L^{1}[0,1]$ satisfies (2.30). This implies that $H$ is continuous since $L$ is a bounded operator from $L^{1}[0,1]$ to $C[0,1]$ .

Lipschitz continuity of ${\mathcal{F}}^{\prime }$ is easy to check. The most interesting point is non-singularity of ${\mathcal{F}}^{\prime }(H)$ . If $H$ is the solution of (2.30), then it is clear by the argument above that

(2.34) $$\begin{eqnarray}\int _{0}^{1}H(\unicode[STIX]{x1D707})\,\text{d}\unicode[STIX]{x1D707}=\displaystyle \frac{1+\sqrt{1-\unicode[STIX]{x1D714}}}{2/\unicode[STIX]{x1D714}}.\end{eqnarray}$$

This clearly shows that there is no real solution for $\unicode[STIX]{x1D714}>1$ .

Since ${\mathcal{F}}^{\prime }$ is the sum of a compact integral operator and the identity, singularity of the operator implies that there is a non-trivial null space. Suppose ${\mathcal{F}}^{\prime }(H)w=0$ . Then

$$\begin{eqnarray}w=\displaystyle \frac{\unicode[STIX]{x1D714}Lw}{(1-\unicode[STIX]{x1D714}LH)^{2}}.\end{eqnarray}$$

We may take $w\geq 0$ by the Perron–Frobenius theorem (Karlin Reference Karlin1959) and may therefore assume that $\int _{0}^{1}w(\unicode[STIX]{x1D707})\,\text{d}\unicode[STIX]{x1D707}>0$ . Since $H$ is a solution, we have $(1-\unicode[STIX]{x1D714}LH)^{-1}=H$ , and so

$$\begin{eqnarray}w(1-\unicode[STIX]{x1D714}wLH)=\unicode[STIX]{x1D714}HLw.\end{eqnarray}$$

So

$$\begin{eqnarray}\int _{0}^{1}w(\unicode[STIX]{x1D707})\,\text{d}\unicode[STIX]{x1D707}=\unicode[STIX]{x1D714}/2\int _{0}^{1}w(\unicode[STIX]{x1D707})\,\text{d}\unicode[STIX]{x1D707}\int _{0}^{1}H(\unicode[STIX]{x1D708})\,\text{d}\unicode[STIX]{x1D708}\end{eqnarray}$$

implies that $\unicode[STIX]{x1D714}\Vert H\Vert _{1}=2$ , which by (2.34) implies that $\unicode[STIX]{x1D714}=1$ .

So ${\mathcal{F}}^{\prime }(H)$ is non-singular unless $\unicode[STIX]{x1D714}=1$ . When $\unicode[STIX]{x1D714}=1$ , ${\mathcal{F}}^{\prime }(H)$ is indeed singular and, again by the Perron–Frobenius theorem, the null space has dimension one and is spanned by a non-negative function. In fact, that function is $w(\unicode[STIX]{x1D707})=\unicode[STIX]{x1D707}H(\unicode[STIX]{x1D707})$ . We see that the singularity structure of the $H$ -equation is quite simple and will use the $H$ -equation as an example again in Section 4.

2.10.2 The discrete problem

Now we return to the discrete version. The purpose of this section is to illustrate how exploiting problem structure can give very different cost estimates from the simple accounting of $O(N^{3})$ work for a linear solve and $N$ function evaluations for a Jacobian.

We begin with the cost of an evaluation. Using $\unicode[STIX]{x1D707}_{i}=(i-1/2)/N$ for $1\leq i\leq N$ and (2.31), we obtain

(2.35) $$\begin{eqnarray}\mathbf{F}^{N}(\mathbf{h})_{i}=h_{i}-\biggl(1-{\displaystyle \frac{\unicode[STIX]{x1D714}}{2N}}\mathop{\sum }_{j=1}^{N}{\displaystyle \frac{ih_{j}}{i+j-1}}\biggr)^{-1}.\end{eqnarray}$$

The approximate integral operator $\mathbf{L}$ , where

$$\begin{eqnarray}\mathbf{L}_{ij}=\displaystyle \frac{\unicode[STIX]{x1D714}i}{2N(i+j-1)},\end{eqnarray}$$

is the product of a diagonal and a Hankel matrix (Golub and Van Loan Reference Golub and Van Loan1996). In fact $\mathbf{L}=\mathbf{D}_{1}\mathbf{K}$ , where

$$\begin{eqnarray}\mathbf{D}_{1}=\text{diag}(i/2N)\end{eqnarray}$$

and $\mathbf{K}$ is the Hankel matrix

$$\begin{eqnarray}\mathbf{K}_{ij}=1/(i+j-1).\end{eqnarray}$$

So evaluation of $\mathbf{L}\mathbf{h}$ can be done at a cost of $O(N\log N)$ work if one computes the product of the Hankel matrix and a vector with a fast Fourier transform (Golub and Van Loan Reference Golub and Van Loan1996). Since the remaining cost of the evaluation of $\mathbf{F}^{N}$ reduces to simple binary operations, the cost of a function evaluation is $O(N\log N)$ work.

We may now proceed exactly as we did in the continuous case. The Jacobian matrix is

(2.36) $$\begin{eqnarray}(\mathbf{F}^{N})^{\prime }(\mathbf{h})_{ij}=\unicode[STIX]{x1D6FF}_{ij}-\displaystyle \frac{\unicode[STIX]{x1D714}\mathbf{L}_{ij}}{1-\unicode[STIX]{x1D714}(\mathbf{L}\mathbf{u})_{i}}.\end{eqnarray}$$

Hence the Jacobian can be expressed as the identity plus the product of a diagonal and a Hankel matrix. If we set

$$\begin{eqnarray}\mathbf{D}_{2}=\text{diag}(1-\unicode[STIX]{x1D714}(\mathbf{L}\mathbf{u})_{i})^{-1},\end{eqnarray}$$

then

$$\begin{eqnarray}(\mathbf{F}^{N})^{\prime }(\mathbf{h})=\mathbf{I}=\unicode[STIX]{x1D714}\mathbf{D}_{1}\mathbf{D}_{2}\mathbf{K},\end{eqnarray}$$

and the Jacobian matrix can be constructed and stored with $O(N^{2})$ work and the Jacobian-vector product computed with $O(N\log N)$ work.

One conclusion that may be unexpected is that, even if one does the matrix–vector product with $\mathbf{K}$ naively for $O(N^{2})$ work, the solution with Newton–GMRES is far faster than one using a direct method to factor  $\mathbf{F}^{\prime }$ . The integral equation structure can be exploited to show that the number of Krylov iterations per nonlinear iteration can be bounded independently of $N$ . Hence, as $N$ becomes large the advantage of Newton–Krylov over Newton–LU grows rapidly. The reader should try this and see. Chandrasekhar (Reference Chandrasekhar1960) tabulates the solution to several figures.

This problem has a rich structure and we will report the costs of a computation in some detail. Table 2.2 summarizes our discussion on the costs of function evaluations, Jacobian evaluations and Jacobian-vector products as functions of $N$ .

Table 2.2.

Iteration component costs.

We will now solve the $H$ -equation for a few values of $\unicode[STIX]{x1D714}$ . The problem becomes more difficult as $\unicode[STIX]{x1D714}$ increases. We will use $\unicode[STIX]{x1D714}=0.5$ (easy), $\unicode[STIX]{x1D714}=0.99$ (less easy) and $\unicode[STIX]{x1D714}=1$ (tricky). The initial iterate in all cases will be

$$\begin{eqnarray}\mathbf{h}_{0}=(1,1,\ldots ,1)^{T}.\end{eqnarray}$$

We will terminate the iteration when

$$\begin{eqnarray}\Vert \mathbf{F}^{N}(\mathbf{h}_{n})\Vert \leq 10^{-10}\Vert \mathbf{F}^{N}(\mathbf{h}_{0})\Vert ,\end{eqnarray}$$

so $\unicode[STIX]{x1D70F}_{r}=10^{-10}$ and $\unicode[STIX]{x1D70F}_{a}=0$ . We will use a mesh with $N=1000$ points to begin with and then demonstrate that the performance of all the methods does not vary much as the mesh is refined. We will do this numerically in this chapter and discuss the theory in Section 3.3.

We will first look at Newton’s method. The dominant cost for a Newton iteration is the matrix factorization. Therefore Newton’s method, even though the number of nonlinear iterations is small, is the most expensive approach for this problem when $\unicode[STIX]{x1D714}<1$ and the number of mesh points $N$ is large. In Figure 2.3 we plot the relative residual ( $\Vert \mathbf{F}^{N}(\mathbf{h}_{n})\Vert /\Vert \mathbf{F}^{N}(\mathbf{h}_{0})\Vert$ ) histories for the three values of $\unicode[STIX]{x1D714}$ .

Figure 2.3.

Newton’s method for the $H$ -equation example.

Figure 2.4.

The chord method for the $H$ -equation example.

Figure 2.3 reinforces the ideas from this section. The concave plots of residual histories for $\unicode[STIX]{x1D714}<1$ indicate superlinear convergence. The linear plot for $\unicode[STIX]{x1D714}=1$ , however, shows that the standard assumptions are violated in this case, as we demonstrated above. One can analyse this (Decker and Kelley Reference Decker and Kelley1980) and show that the simple scalar equation $x^{2}=0$ fully explains the q-/linear convergence rate. For any $x_{0}\neq 0$ , the Newton iteration for $x^{2}=0$ is simply

$$\begin{eqnarray}x_{n}=x_{n-1}/2.\end{eqnarray}$$

Hence the convergence is q-linear with q-factor $1/2$ . The residuals

$$\begin{eqnarray}x_{n}^{2}=x_{n-1}^{2}/4\end{eqnarray}$$

converge q-linearly to $0$ with q-factor $1/4$ .

In Figure 2.4 we plot the residual histories for the chord method. For $\unicode[STIX]{x1D714}=0.5$ and $\unicode[STIX]{x1D714}=0.99$ one can see the q-linear convergence predicted by theory. However, for $\unicode[STIX]{x1D714}=1$ something different is happening. The discussion for the continuous case applies without modification to the discrete problem and $\mathbf{F}^{\prime }(\mathbf{h}^{\ast })$ is singular for $\unicode[STIX]{x1D714}=1$ . There is theory for that case (Decker and Kelley Reference Decker and Kelley1983) and the results can be explained by considering the scalar equation $x^{2}=0$ . With the initial iterate $x_{0}=1$ the chord iterations are

$$\begin{eqnarray}x_{n}=x_{n-1}(1-x_{n-1}/2)=O(1/n),\end{eqnarray}$$

which is very slow convergence.

If one refines the mesh for this problem, one will have a more accurate approximation of the solution of the integral equation. However, mesh refinement will have a very small effect on the iteration statistics. Table 2.3 illustrates this by showing some relative residual norms for the Newton iteration with $\unicode[STIX]{x1D714}=0.5$ and a few values of $N$ . As you can see, the relative residual norms are identical to three figures until the terminal iteration. This effect is called mesh independence. We will return to this topic in Section 3.3. We were limited to $N=8000$ because we could not store or factor larger dense matrices in our computing environment.

Table 2.3.

Newton residuals and timings.

We also report computer times in seconds using the MATLAB tic and toc commands. We do this only to show that the run times do not, as you might expect from the $O(N^{3})$ work in the matrix factorizations, increase by a factor of eight as we double the size of the problem. The reason for this is that there is a great deal of $O(N^{2})$ work in the Jacobian evaluations, which is significant in this case, and the LU factorization in MATLAB is very fast. The effects of the $O(N^{2})$ work are also seen in the chord iterations. One would expect the chord iterations to be much faster if the $O(N^{3})$ work in the matrix factorization dominated the computation. In Table 2.4 we show the first seven scaled residual norms for the chord iteration and the timings for the entire sequence of 63 iterations. As is the case with Newton’s method, the iteration statistics agree. Figures 2.3 and 2.4 would not change with higher values of $N$ .

Table 2.4.

Chord residuals and timings.

The performance of the JFNK iteration is far better because there is no $O(N^{2})$ work or storage at all – only the $O(N\log N)$ cost of the function evaluations. We will illustrate this by repeating the computations above with Newton–GMRES, a finite-difference Jacobian-vector product, and a constant forcing term of $\unicode[STIX]{x1D702}=0.1$ . The results were not affected by switching to an analytic Jacobian-vector product. In Figure 2.5 we plot the iteration histories. The curves are not as smooth as those in Figure 2.3, which reflects the limited accuracy of the linear solves. However, the overall cost was, as you will see in Table 2.5, much less.

Figure 2.5.

Newton–GMRES for the $H$ -equation example.

The iteration is so fast that we must solve much larger problems to see the $O(N\log N)$ cost reflected in the timings. Table 2.5 shows the mesh independence of the nonlinear iteration. The performance of the linear iteration is also mesh-independent. For each value of $N$ the entire iteration needed $19$ calls to $\mathbf{F}^{N}$ . These $19$ calls include both the evaluation of the residual and the additional calls for the finite-difference Jacobian-vector product. One could also use an analytic Jacobian-vector product and both the cost and the results in Table 2.5 would be the same.

Table 2.5.

Newton–GMRES residuals and timings.

3.1 The theorem

The Kantorovich theorem (Kantorovich and Akilov Reference Kantorovich and Akilov1982) is a convergence result for Newton’s method that replaces the standard assumption that there is a solution with the statement that $\mathbf{F}(\mathbf{x}_{0})$ is sufficiently small (consistency). That coupled with a uniform bound on $\mathbf{F}^{\prime }(\mathbf{x})^{-1}$ in a neighbourhood of $\mathbf{x}_{0}$ (stability) implies that there is a unique solution $\mathbf{x}^{\ast }$ in a (smaller) neighbourhood of $\mathbf{x}_{0}$ and that $\mathbf{x}^{\ast }$ is near $\mathbf{x}_{0}$ (convergence!).

We use the formulation of the theorem from Ortega and Rheinboldt (Reference Ortega and Rheinboldt1970). There are several variations of the theorem. In addition to Ortega and Rheinboldt (Reference Ortega and Rheinboldt1970) and Kantorovich and Akilov (Reference Kantorovich and Akilov1982), one can find examples of alternative formulations in Kelley (Reference Kelley1995), Dennis (Reference Dennis1969), and Dennis (Reference Dennis and Rall1971). In Section 6.2 we will present a version from Qi and Sun (Reference Qi and Sun1993) for non-smooth functions.

The assumptions are similar to the standard assumptions.

Note that the assumptions on $\mathbf{x}_{0}$ cannot be stated in terms of $\Vert \mathbf{x}_{0}-\mathbf{x}^{\ast }\Vert$ because the existence of $\mathbf{x}^{\ast }$ is part of the conclusion of the theorem. Rather, one uses $\Vert \mathbf{F}^{\prime }(\mathbf{x}_{0})^{-1}\mathbf{F}(\mathbf{x}_{0})\Vert$ as a surrogate. However, the theorem will imply that

$$\begin{eqnarray}\unicode[STIX]{x1D702}\geq \Vert \mathbf{F}^{\prime }(\mathbf{x}_{0})^{-1}\mathbf{F}(\mathbf{x}_{0})\Vert =\Vert \mathbf{e}_{0}\Vert +O(\Vert \mathbf{e}_{0}\Vert ^{2}).\end{eqnarray}$$

Since $\unicode[STIX]{x1D6FD}=\Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })^{-1}\Vert +O(\Vert \mathbf{e}_{0}\Vert )$ , the assumption $\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D702}<1/2$ implies that

$$\begin{eqnarray}\Vert \mathbf{e}_{0}\Vert \leq \displaystyle \frac{\Vert \mathbf{F}^{\prime }(\mathbf{x}^{\ast })^{-1}\Vert ^{-1}}{2\unicode[STIX]{x1D6FE}}+O(\Vert \mathbf{e}_{0}\Vert ^{2}),\end{eqnarray}$$

which is close to (2.3) if $\Vert \mathbf{e}_{0}\Vert$ is small enough to permit one to neglect the second-order term. One remarkable fact about the Kantorovich theorem is that there is no need for $\Vert \mathbf{e}_{0}\Vert$ to be small to obtain convergence. In that sense the theorem is a global convergence result. This is why one cannot expect q-quadratic convergence for the entire iteration. The r-quadratic convergence one obtains from the theorem is remarkable.

Theorem 3.2. Let Assumption 3.1 hold. Then there is a unique root $x^{\ast }$ of $F$ in ${\mathcal{B}}(\mathbf{x}_{0},r_{+})$ , the Newton iteration with $\mathbf{x}_{0}$ as the initial iterate converges to $\mathbf{x}^{\ast }$ , and $x_{n}\in {\mathcal{B}}(\mathbf{x}_{0},r_{-})$ for all $n$ . The errors satisfy the estimate

(3.1) $$\begin{eqnarray}\Vert \mathbf{e}_{n}\Vert \leq \displaystyle \frac{(2\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FE})^{2^{n}}}{2^{n}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}},\end{eqnarray}$$

and hence the convergence is r-quadratic.

We refer to Ortega and Rheinboldt (Reference Ortega and Rheinboldt1970) for the proof and to Kelley (Reference Kelley1995) for a simpler proof of an analogous theorem for the chord method. The proof in Kelley (Reference Kelley1995) is based on directly showing that the chord map

$$\begin{eqnarray}\mathbf{G}_{C}(\mathbf{x})=\mathbf{x}-\mathbf{F}^{\prime }(\mathbf{x}_{0})^{-1}\mathbf{F}(\mathbf{x})\end{eqnarray}$$

is a contraction on ${\mathcal{B}}(\mathbf{x}_{0},r_{-})$ . The proof of Theorem 3.2 uses a recursion relation for a bound on the steps to show that the Newton map

$$\begin{eqnarray}\mathbf{G}_{N}(\mathbf{x})=\mathbf{x}-\mathbf{F}^{\prime }(\mathbf{x})^{-1}\mathbf{F}(\mathbf{x})\end{eqnarray}$$

is a contraction and to obtain the convergence rate estimate (3.1).

3.2 A boundary value problem

In this section we will illustrate in some detail how the Kantorovich theorem can be applied to show convergence of a finite-difference approximation to a nonlinear two-point boundary value problem. We will begin with an existence/uniqueness result for the continuous problem and then directly estimate $\unicode[STIX]{x1D702}$ , $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$ to apply the Kantorovich theorem.

The boundary problem is

(3.2) $$\begin{eqnarray}-u^{\prime \prime }(x)-\cos (u)=0,\quad u(0)=u(1)=0.\end{eqnarray}$$

We will use the contraction mapping theorem in the space $C[0,1]$ to prove existence and uniqueness of a solution of (3.2).

Let ${\mathcal{G}}$ be the inverse of the negative Laplacian in one space dimension with homogeneous Dirichlet boundary conditions, that is,

(3.3) $$\begin{eqnarray}{\mathcal{G}}(u)(x)=\int _{0}^{1}g(x,y)u(y)\,\text{d}y,\end{eqnarray}$$

where the Green’s function $g$ is

$$\begin{eqnarray}g(x,z)=\left\{\begin{array}{@{}ll@{}}x(1-z)\quad & 0<x<z,\\ z(1-x)\quad & z<x<1.\end{array}\right.\end{eqnarray}$$

If $f\in C[0,1]$ and $u={\mathcal{G}}(f)$ , then

$$\begin{eqnarray}-u^{\prime \prime }=f,\quad u(0)=u(1)=0.\end{eqnarray}$$

Hence (3.2) has a solution if and only if the equivalent integral equation

(3.4) $$\begin{eqnarray}{\mathcal{F}}(u)(x)=u(x)-\int _{0}^{1}g(x,y)\cos (u(y))\,\text{d}y=0\end{eqnarray}$$

has a solution. We express the integral equation as a fixed-point problem

$$\begin{eqnarray}u(x)={\mathcal{K}}(u)(x)\equiv {\mathcal{G}}(\cos (u))(x),\end{eqnarray}$$

and will show that ${\mathcal{K}}$ is a contraction on $C[0,1]$ .

Clearly

$$\begin{eqnarray}|\text{cos}(u)-\cos (v)|\leq |u-v|,\quad \text{for all }u,v.\end{eqnarray}$$

So

$$\begin{eqnarray}|{\mathcal{K}}(u)(x)-{\mathcal{K}}(v)(x)|\leq \Vert u-v\Vert _{\infty }\int _{0}^{1}g(x,y)\,\text{d}y.\end{eqnarray}$$

Note that

$$\begin{eqnarray}w(x)=\int _{0}^{1}g(x,y)\,\text{d}y=\displaystyle \frac{x(1-x)}{2}\end{eqnarray}$$

and hence $\Vert w\Vert _{\infty }=1/8$ . Thus

$$\begin{eqnarray}\Vert {\mathcal{K}}(u)-{\mathcal{K}}(v)\Vert _{\infty }\leq \Vert u-v\Vert _{\infty }/8.\end{eqnarray}$$

So ${\mathcal{K}}$ is a contraction on $C[0,1]$ . We now know that the boundary value problem has a unique solution $u^{\ast }\in C[0,1]$ . We have also derived a uniform bound for ${\mathcal{F}}^{\prime }(u)^{-1}$ . For any $w\in C[0,1]$ we have (using (2.27))

$$\begin{eqnarray}{\mathcal{F}}^{\prime }(u)w(x)=w(x)+\int _{0}^{1}g(x,y)\sin (u(y))w(y)\,\text{d}y\end{eqnarray}$$

and so

$$\begin{eqnarray}\Vert {\mathcal{F}}^{\prime }(u)w\Vert _{\infty }\geq (7/8)\Vert w\Vert _{\infty }.\end{eqnarray}$$

Hence

$$\begin{eqnarray}\Vert {\mathcal{F}}^{\prime }(u)^{-1}\Vert _{\infty }\leq 8/7\end{eqnarray}$$

for all $u\in C[0,1]$ . We can also conclude from the above that ${\mathcal{F}}^{\prime }$ is Lipschitz continuous with Lipschitz constant $\unicode[STIX]{x1D6FE}\leq 1/8$ .

We discretize the problem with central finite differences with a uniform mesh width $h$ . The discrete problem is

$$\begin{eqnarray}\mathbf{D}^{h}\mathbf{u}^{h}+\cos (\mathbf{u}^{h})=0.\end{eqnarray}$$

We precondition with $\mathbf{G}^{h}\equiv (\mathbf{D}^{h})^{-1}$ and obtain the analogue of (3.4):

(3.5) $$\begin{eqnarray}\mathbf{F}^{h}(\mathbf{u}^{h})=\mathbf{u}^{h}+\mathbf{G}^{h}\cos (\mathbf{u}^{h})\equiv \mathbf{u}^{h}+\mathbf{K}^{h}(\mathbf{u}^{h}).\end{eqnarray}$$

The smoothness of $u^{\ast }$ implies that there is $\unicode[STIX]{x1D702}_{0}>0$ such that

(3.6) $$\begin{eqnarray}\Vert \mathbf{F}^{h}({\mathcal{E}}^{h}u^{\ast })\Vert _{\infty }\leq \unicode[STIX]{x1D702}_{0}h^{2}\end{eqnarray}$$

and that

(3.7) $$\begin{eqnarray}\Vert (\mathbf{F}^{h})^{-1}(u)\Vert _{\infty }\leq 2\end{eqnarray}$$

uniformly in $u$ and $h$ . Clearly the Lipschitz constant $\unicode[STIX]{x1D6FE}^{h}$ of of $(\mathbf{F}^{h})^{\prime }$ is $\leq 1/4$ for $h$ sufficiently small.

We are can now apply the Kantorovich theorem to show that a solution $\mathbf{u}^{h}$ of (3.5) exists and that

(3.8) $$\begin{eqnarray}\Vert \mathbf{u}^{h}-{\mathcal{E}}^{h}u^{\ast }\Vert _{\infty }=O(h^{2}).\end{eqnarray}$$

We will apply the theorem with $\unicode[STIX]{x1D6FD}=2$ from (3.7), $\unicode[STIX]{x1D702}=2\unicode[STIX]{x1D702}_{0}h^{2}$ from (3.6) and (3.7), and $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}^{h}=1/4$ . Then

$$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D702}_{0}h^{2}<1/2\end{eqnarray}$$

for $h$ sufficiently small. So there is a solution $\mathbf{u}^{h}$ in ${\mathcal{B}}({\mathcal{E}}^{h}u^{\ast },r_{-})$ where

$$\begin{eqnarray}r_{-}=\displaystyle \frac{1-\sqrt{1-2\unicode[STIX]{x1D6FC}}}{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}}=\displaystyle \frac{1-\sqrt{1-2\unicode[STIX]{x1D702}_{0}h^{2}}}{1/2}=O(h^{2}).\end{eqnarray}$$

Summarizing, stability (3.7) and consistency (3.6) imply convergence.

3.3 Mesh independence results

The example in Section 3.2 made use of the fact that the parameters in the Kantorovich theorem could be estimated independently of the mesh $h$ . Hence the convergence estimate in (3.1) is also independent of $h$ . This is a simple example of a mesh independence theorem. In this section we will explore two variations of this idea. Classical mesh independence results are about the convergence of the iteration statistics to those of an underlying infinite-dimensional problem as a grid is refined. In this case the errors are deterministic and we look at some of the ideas in this section. After that we consider tracking theorems. These results describe convergence of the iteration statistics when the errors are stochastic and depend, for example, on the sample size.

3.3.1 Deterministic errors

The concept of mesh independence has its origin in Allgower, Böhmer, Potra and Rheinboldt (Reference Allgower, Böhmer, Potra and Rheinboldt1986). While the results may seem obvious, mesh independence is a very useful concept in both nonlinear equations (Kelley and Sachs Reference Kelley and Sachs1991, Ferng and Kelley Reference Ferng and Kelley2000) and optimization (Sachs Reference Sachs1990, Hintermüller and Ulbrich Reference Hintermüller and Ulbrich2003).

The idea is that one is approximating an infinite-dimensional equation

(3.9) $$\begin{eqnarray}{\mathcal{F}}(u)=0,\end{eqnarray}$$

defined on a Banach space $X$ , with a sequence of finite-dimensional problems

(3.10) $$\begin{eqnarray}\mathbf{F}^{h}(\mathbf{u})=0\end{eqnarray}$$

on $\mathbb{R}^{N}$ . The boundary value problem from Section 3.2 and the Chandrasekhar $H$ -equation from Section 2.10 are examples of this situation. We will be interested in a solution $u^{\ast }$ of (3.9) at which the standard assumptions hold.

In both of these examples the performance of Newton’s method is mesh-independent in the sense that the relative residuals not only converge to zero, but also converge for each iteration as the mesh is refined. In Section 3.2 we introduced one way to quantify this. Let ${\mathcal{E}}^{h}:X\rightarrow \mathbb{R}^{N}$ be the projection from the space $X$ to $\mathbb{R}^{N}$ that encodes the discretization. For a finite-difference approximation, for example, ${\mathcal{E}}^{h}$ could be evaluation at the grid points.

Convergence of the approximation means that the solution $\mathbf{u}^{h}$ of (3.10) satisfies

(3.11) $$\begin{eqnarray}\lim _{h\rightarrow 0}\Vert \mathbf{u}^{h}-{\mathcal{E}}^{h}u^{\ast }\Vert =0.\end{eqnarray}$$

In many cases, including the examples in Sections 3.2 and 2.10, one can show that the Newton iterations also converge in the sense that

(3.12) $$\begin{eqnarray}\lim _{h\rightarrow 0}\Vert \mathbf{u}_{n}^{h}-{\mathcal{E}}^{h}u_{n}\Vert =0,\end{eqnarray}$$

where $\mathbf{u}_{n}^{h}$ is the $n$ th Newton iteration for (3.10) and $u_{n}$ is the $n$ th Newton iteration for (3.9).

We will express the fact that the iteration statistics converge in two ways, both of which follow from (3.11) and (3.12). These results also hold for the JFNK methods we discussed in Section 2.6. We will denote the infinite-dimensional iteration

$$\begin{eqnarray}u_{n+1}=u_{n}-{\mathcal{F}}^{\prime }(u_{n})^{-1}{\mathcal{F}}(u_{n}).\end{eqnarray}$$

• ∙ Let $\unicode[STIX]{x1D716}>0$ and let $k^{h}$ be the least $k$ such that $\Vert \mathbf{F}^{h}(\mathbf{u}_{k}^{h})\Vert /\Vert \mathbf{F}^{h}(\mathbf{u}_{0}^{h})\Vert <\unicode[STIX]{x1D716}$ . Let $k^{0}$ be the least $k$ such that

  $$\begin{eqnarray}\Vert {\mathcal{F}}(u_{k})\Vert /\Vert {\mathcal{F}}(u_{0})\Vert <\unicode[STIX]{x1D716}.\end{eqnarray}$$
  Then, for all $h$ sufficiently small, $|k^{h}-k^{0}|\leq 1$ .

• ∙ Let $K>0$ and let $\unicode[STIX]{x1D716}>0$ be given. Then, for $h$ sufficiently small,

  $$\begin{eqnarray}\biggl|\displaystyle \frac{\Vert \mathbf{F}^{h}(\mathbf{u}_{k}^{h})\Vert }{\Vert \mathbf{F}^{h}(\mathbf{u}_{0}^{h})\Vert }-\displaystyle \frac{\Vert {\mathcal{F}}(u_{k}^{h})\Vert }{\Vert {\mathcal{F}}(u_{0})\Vert }\biggr|<\unicode[STIX]{x1D716}\end{eqnarray}$$
  for all $0\leq k\leq K$ .

3.3.2 Tracking theorems

The results in Willert, Chen and Kelley (Reference Willert, Chen and Kelley2015) consider the case where the function $\mathbf{F}$ , the Jacobian $\mathbf{F}^{\prime }$ and Jacobian-vector products are not evaluated directly, but are approximated using internal Monte Carlo computations. The work was motivated by problems in neutron transport (Willert, Kelley, Knoll and Park Reference Willert, Kelley, Knoll and Park2013, Knoll, Park and Smith Reference Knoll, Park and Smith2011), where a Monte Carlo simulation was embedded in the residual. The results are similar to mesh independence theorems in that one seeks to show that some finite subset of the iterations converges, but does not seek to drive the approximation error (or its variance) to zero. However, the assumptions and convergence theorems are more technical. We used the term ‘tracking’ in Willert et al. (Reference Willert, Chen and Kelley2015) rather than convergence.

We approximate functions, Jacobians, and Jacobian-vector projects with a Monte Carlo simulation with a sample size $N_{\text{MC}}$ . The notation from Willert et al. (Reference Willert, Chen and Kelley2015) is as follows.

• ∙ $N_{\text{MC}}$ is the sample size for the function and $N_{\text{MC}}^{J}$ is the sample size for the Jacobian or Jacobian-vector product.

• ∙ $\tilde{\mathbf{F}}(\mathbf{x},N_{\text{MC}})$ is an outcome of the simulation for the residual $\mathbf{F}(\mathbf{x})$ .

• ∙ $\mathbf{J}(\mathbf{x},N_{\text{MC}}^{J})$ is an outcome of the simulation for the Jacobian $\mathbf{F}^{\prime }(\mathbf{x})$ .

• ∙ $\mathbf{J}_{p}(\mathbf{x},\mathbf{v},N_{\text{MC}}^{J})$ is an outcome of the simulation for the Jacobian-vector product $\mathbf{F}^{\prime }(\mathbf{x})\mathbf{v}$ .

We assume that the evaluations of $\tilde{\mathbf{F}}$ , $\mathbf{J}$ and $\mathbf{J}_{p}$ are independent.

Suppose the standard assumptions hold and that Newton’s method converges for all $\mathbf{u}\in {\mathcal{B}}(\mathbf{x}^{\ast },\unicode[STIX]{x1D70C})$ . We make a consistency assumption for the function and Jacobian evaluations.

Assumption 3.3. There are functions $c_{F}$ , $c_{J}$ and $c_{Jv}$ and an open set ${\mathcal{B}}^{\prime }$ which contains ${\mathcal{B}}(\mathbf{x}^{\ast },\unicode[STIX]{x1D70C})$ such that, for all $\mathbf{x}\in {\mathcal{B}}^{\prime }$ , unit vectors $\mathbf{v}\in \mathbb{R}^{N}$ and $\unicode[STIX]{x1D6FF}>0$ ,

(3.13) $$\begin{eqnarray}\displaystyle \operatorname{Prob}\biggl(\Vert \mathbf{F}(\mathbf{x})-\tilde{\mathbf{F}}(\mathbf{x},N_{\text{MC}})\Vert >\displaystyle \frac{c_{F}(\unicode[STIX]{x1D6FF})}{\sqrt{N_{\text{MC}}}}\biggr) & {<} & \displaystyle \unicode[STIX]{x1D6FF},\end{eqnarray}$$

(3.14) $$\begin{eqnarray}\displaystyle \operatorname{Prob}\biggl(\Vert \mathbf{F}^{\prime }(\mathbf{x})-\mathbf{J}(\mathbf{x},N_{\text{MC}}^{J})\Vert >\displaystyle \frac{c_{J}(\unicode[STIX]{x1D6FF})}{\sqrt{N_{\text{MC}}^{J}}}\biggr) & {<} & \displaystyle \unicode[STIX]{x1D6FF}.\end{eqnarray}$$

In this article we will only summarize the results from Willert et al. (Reference Willert, Chen and Kelley2015) on methods which use full Jacobians. The assumptions on the quality of the residual evaluation and the Jacobian are not sufficient to guarantee quadratic convergence. The algorithm increases the number of samples as the iteration progresses, which sometimes reflects practice. One would never increase the number of samples rapidly enough to capture superlinear convergence, and the theory reflects that. The algorithm from Willert et al. (Reference Willert, Chen and Kelley2015) is given in algorithm newton_MC.

Theorem 3.4. Let (3.13) and (3.14) from Assumption 3.3 and the assumptions of Theorem 2.2 hold. Let $r_{\text{Newton}}\in (0,1)$ be given and assume that $\Vert \mathbf{e}_{0}\Vert \leq \unicode[STIX]{x1D70C}$ and $\hat{\unicode[STIX]{x1D702}}$ are small enough for the inexact Newton iteration to converge q-linearly with a q-factor $r_{\text{Newton}}$ . Let a positive integer $K$ , $r\in (r_{\text{Newton}},1)$ and $\unicode[STIX]{x1D714}\in (0,1)$ be given. Then there are $\hat{\unicode[STIX]{x1D702}}$ , $N_{\text{MC}}$ , $N_{\text{MC}}^{J}$ and $N_{\text{inc}}$ , such that, with probability $(1-\unicode[STIX]{x1D714})$ for all $1\leq n\leq K$ , the iteration produced by algorithm newton_MC satisfies

(3.15) $$\begin{eqnarray}\Vert \mathbf{e}_{n}\Vert \leq r^{n}\Vert \mathbf{e}_{0}\Vert ,\end{eqnarray}$$

If one approximates matrix–vector products rather than the Jacobian itself, then assumption (3.14) must be replaced by one on the matrix–vector product computation. The subtle, and important, problem with this is that there is no underlying matrix. This can (and does: see Willert et al. Reference Willert, Kelley, Knoll and Park2013) produce error accumulation in the Krylov method. Simoncini and Szyld (Reference Simoncini and Szyld2003a , Reference Simoncini and Szyld2003b , Reference Simoncini and Szyld2007) explain this in detail. Willert et al. (Reference Willert, Chen and Kelley2015) use that analysis to derive very technical tracking results, which we mercifully omit from this paper, for matrix-free methods.

4.1 The implicit function theorem

The implicit function theorem says that algorithm simple_continuation will successfully follow the solution path as long as $\mathbf{F}^{\prime }$ is safely non-singular and also provides an estimate of $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D706}}$ . We give a proof that is a nice application of the Kantorovich theorem.

We begin with some notation. If $\mathbf{F}(\mathbf{x},\unicode[STIX]{x1D706})$ is a function of $\mathbf{x}\in \mathbb{R}^{N}$ and $\unicode[STIX]{x1D706}\in R$ , then $\mathbf{F}_{\mathbf{x}}$ will denote the Jacobian in the $\mathbf{x}$ variable. Hence $\mathbf{F}_{\mathbf{x}}$ is an $N\times N$ matrix. Similarly $\mathbf{F}_{\unicode[STIX]{x1D706}}\in \mathbb{R}^{N}$ is the partial of $\mathbf{F}$ with respect to the scalar variable $\unicode[STIX]{x1D706}$ .

Proof. Similarly to the proof of Theorem 2.2, Lipschitz continuity of $\mathbf{F}_{\mathbf{x}}$ implies that $\mathbf{F}_{\mathbf{x}}(\mathbf{x},\unicode[STIX]{x1D706})$ is non-singular in the set

$$\begin{eqnarray}{\mathcal{D}}=\biggl\{(\mathbf{x},\unicode[STIX]{x1D706})\mid \Vert \mathbf{x}-\mathbf{x}_{0}\Vert +|\unicode[STIX]{x1D706}-\unicode[STIX]{x1D706}_{0}|\leq \displaystyle \frac{1}{2\unicode[STIX]{x1D6FE}\Vert \mathbf{F}_{\mathbf{x}}^{-1}(\mathbf{x}_{0},\unicode[STIX]{x1D706}_{0})\Vert }\biggr\},\end{eqnarray}$$

and that

$$\begin{eqnarray}\Vert \mathbf{F}_{\mathbf{x}}^{-1}(\mathbf{x},\unicode[STIX]{x1D706})\Vert \leq \unicode[STIX]{x1D6FD}\equiv 2\unicode[STIX]{x1D6FD}_{0}\end{eqnarray}$$

for all $(\mathbf{x},\unicode[STIX]{x1D706})\in {\mathcal{D}}$ .

We have $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$ in hand for an application of the Kantorovich theorem. We use Lipschitz continuity of $\mathbf{F}$ to obtain

$$\begin{eqnarray}\Vert \mathbf{F}(\mathbf{x}_{0},\unicode[STIX]{x1D706})\Vert =\Vert \mathbf{F}(\mathbf{x}_{0},\unicode[STIX]{x1D706})-\mathbf{F}(\mathbf{x}_{0},\unicode[STIX]{x1D706}_{0})\Vert \leq \unicode[STIX]{x1D6FE}|\unicode[STIX]{x1D706}-\unicode[STIX]{x1D706}_{0}|\equiv \unicode[STIX]{x1D702}.\end{eqnarray}$$

So if $\unicode[STIX]{x1D6FE}^{2}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6E5}<1/2$ and $|\unicode[STIX]{x1D706}-\unicode[STIX]{x1D706}_{0}|\leq \unicode[STIX]{x1D6E5}$ , we can apply the Kantorovich theorem with

$$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D702}=\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}^{2}\unicode[STIX]{x1D6E5}\leq 1/2,\end{eqnarray}$$

to complete the proof of existence with

$$\begin{eqnarray}r=\displaystyle \frac{1\pm \sqrt{1+2\unicode[STIX]{x1D6FC}}}{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}}.\end{eqnarray}$$

To show differentiability, we formally differentiate (4.1) with respect to $\unicode[STIX]{x1D706}$ and note that $\mathbf{x}^{\prime }(\unicode[STIX]{x1D706})$ is the solution $\mathbf{z}$ of

If $\mathbf{F}$ has higher-order derivatives, one can show that $\mathbf{x}$ has as many derivatives in $\unicode[STIX]{x1D706}$ as $\mathbf{F}$ does in $(\mathbf{x},\unicode[STIX]{x1D706})$ (Rabinowitz Reference Rabinowitz1971).

If $\mathbf{F}_{\mathbf{x}}(\mathbf{x}^{\ast },\unicode[STIX]{x1D706}^{\ast })$ is singular, then the implicit function cannot be used to assert that there are solutions near $\unicode[STIX]{x1D706}^{\ast }$ . In the case of the $H$ -equation, $\mathbf{F}_{\mathbf{h}}(\mathbf{h}^{\ast },1)$ is singular and (2.34) implies that there is no real solution for $\unicode[STIX]{x1D714}>1$ .

So, does the solution arc stop dead at the point $(H(1),1)$ ? The answer is that for many common singularities, including the one for the $H$ -equation, the solution arc does not terminate abruptly but either loops back or becomes unbounded (Crandall and Rabinowitz Reference Crandall and Rabinowitz1971, Rabinowitz Reference Rabinowitz1971, Keller Reference Keller1987). We will discuss some specific examples in the rest of this section.

4.2 Simple folds and pseudo-arclength continuation

We will begin this section with an example that is both simple and general. Consider the scalar equation

$$\begin{eqnarray}f(x,\unicode[STIX]{x1D706})=x^{2}-\unicode[STIX]{x1D706}.\end{eqnarray}$$

The function $f$ is Lipschitz continuous if $x$ is restricted to a bounded set. Since $f_{x}=2x$ the Lipschitz constant of the derivative is $2$ . There is a unique solution $x(0)=0$ when $\unicode[STIX]{x1D706}=0$ , no real solution when $\unicode[STIX]{x1D706}<0$ , and two solutions when $\unicode[STIX]{x1D706}>0$ . The singularity of $f_{x}$ at $\unicode[STIX]{x1D706}=0$ has, as we will see, exactly the same structure as the singularity of the $H$ -equation when $\unicode[STIX]{x1D714}=1$ .

Suppose we begin with $(x,\unicode[STIX]{x1D706})=(1,1)$ and try to use simple continuation to reach and pass $\unicode[STIX]{x1D706}=0$ . The continuation will fail because there is no solution for $\unicode[STIX]{x1D706}<0$ . However, the path of solutions continues beyond $\unicode[STIX]{x1D706}=0$ , the only difference being that the sign of $x$ changes. How can we modify the simple continuation algorithm to follow the path without getting stuck at the singularity?

One way would be to interchange the roles of $\unicode[STIX]{x1D706}$ and $x$ . This is the idea of the PITCON (Rheinboldt Reference Rheinboldt1986) code. If we do that, then $f_{\unicode[STIX]{x1D706}}=-1$ is never singular and we can continue with ease. The problem with this approach is that one must identify the variable to exchange with $\unicode[STIX]{x1D706}$ . Rheinboldt (Reference Rheinboldt1986) identifies many situations where this is readily done. There are many advanced continuation codes (Salinger et al. Reference Salinger, Bou-Rabee, Pawlowski, Wilkes, Burroughs, Lehoucq and Romero2002, Doedel and Kernévez Reference Doedel and Kernévez1986, Govaerts Reference Govaerts2000) and some very good books (Govaerts Reference Govaerts2000, Kuznetsov Reference Kuznetsov1998, Doedel Reference Doedel1997, Keller Reference Keller1987, Rheinboldt Reference Rheinboldt1986) on the topic.

The approach we describe in this section, pseudo-arclength continuation, attempts to parametrize the solution arc and then add an approximation to arclength as a new parameter (Keller Reference Keller1987).

To see how this would work, suppose that $x(s)$ and $\unicode[STIX]{x1D706}(s)$ are functions of an arclength parameter $s$ . Setting ${\dot{x}}=\text{d}x/\text{d}s$ and $\dot{\unicode[STIX]{x1D706}}=\text{d}\unicode[STIX]{x1D706}/\text{d}s$ , we use the arclength normalization

(4.2) $$\begin{eqnarray}{\dot{x}}^{2}+\dot{\unicode[STIX]{x1D706}}^{2}=1.\end{eqnarray}$$

We define an expanded equation in $\mathbf{z}=(x,\unicode[STIX]{x1D706})^{T}$ with $s$ as the parameter:

(4.3) $$\begin{eqnarray}\mathbf{G}(\mathbf{z},s)\equiv \left(\begin{array}{@{}c@{}}f(x,\unicode[STIX]{x1D706})\\ {\dot{x}}^{2}+\dot{\unicode[STIX]{x1D706}}^{2}-1\end{array}\right)=\left(\begin{array}{@{}c@{}}0\\ 0\end{array}\right)\!.\end{eqnarray}$$

If $\mathbf{g}_{\mathbf{z}}$ is non-singular and one can effectively approximate ${\dot{x}}$ and $\dot{\unicode[STIX]{x1D706}}$ , then one can apply simple continuation to the expanded system.

We will investigate the non-singularity of $\mathbf{G}_{\mathbf{z}}$ first. Clearly

$$\begin{eqnarray}\mathbf{G}_{\mathbf{z}}=\left(\begin{array}{@{}cc@{}}f_{x} & f_{\unicode[STIX]{x1D706}}\\ {\dot{x}} & \dot{\unicode[STIX]{x1D706}}\end{array}\right)=\left(\begin{array}{@{}cc@{}}2x & -1\\ {\dot{x}} & \dot{\unicode[STIX]{x1D706}}\end{array}\right)\!,\end{eqnarray}$$

so

$$\begin{eqnarray}\det (\mathbf{G}_{\mathbf{z}})=2x\dot{\unicode[STIX]{x1D706}}+{\dot{x}}.\end{eqnarray}$$

Use $f(x,\unicode[STIX]{x1D706})=0$ to get

$$\begin{eqnarray}0=\,\text{d}f(x,\unicode[STIX]{x1D706})/\,\text{d}s=f_{x}{\dot{x}}+f_{\unicode[STIX]{x1D706}}\dot{\unicode[STIX]{x1D706}}=2x{\dot{x}}-\dot{\unicode[STIX]{x1D706}}\end{eqnarray}$$

and so $2x{\dot{x}}=\dot{\unicode[STIX]{x1D706}}$

Now multiply $\det (\mathbf{G}_{\mathbf{z}})=2x\dot{\unicode[STIX]{x1D706}}+{\dot{x}}$ by ${\dot{x}}$ and

$$\begin{eqnarray}{\dot{x}}\det (\mathbf{G}_{\mathbf{z}})=2x{\dot{x}}\dot{\unicode[STIX]{x1D706}}+{\dot{x}}^{2}.\end{eqnarray}$$

Use $2x{\dot{x}}=\dot{\unicode[STIX]{x1D706}}$ and

$$\begin{eqnarray}{\dot{x}}\det (\mathbf{G}_{\mathbf{z}})=\dot{\unicode[STIX]{x1D706}}^{2}+{\dot{x}}^{2}=1.\end{eqnarray}$$

So $\det (\mathbf{G}_{\mathbf{z}})\neq 0$ .

One way to approximate ${\dot{x}}$ and $\dot{\unicode[STIX]{x1D706}}$ is to use simple continuation for two values of $\unicode[STIX]{x1D706}$ and then switch to continuation in $s$ . If one has two solutions $(x(\unicode[STIX]{x1D706}_{0}),\unicode[STIX]{x1D706}_{0})$ and $(x(\unicode[STIX]{x1D706}_{-1}),\unicode[STIX]{x1D706}_{-1})$ in hand, then one could approximate $ds$ by using

$$\begin{eqnarray}{\dot{x}}^{2}+\dot{\unicode[STIX]{x1D706}}^{2}=1\end{eqnarray}$$

to conclude that if

$$\begin{eqnarray}(x_{0}-x_{-1})^{2}+(\unicode[STIX]{x1D706}_{0}-\unicode[STIX]{x1D706}_{-1})^{2}\equiv ds^{2}\end{eqnarray}$$

then

$$\begin{eqnarray}{\dot{x}}_{0}\approx \displaystyle \frac{x_{0}-x_{-1}}{ds}\quad \text{and}\quad \dot{\unicode[STIX]{x1D706}}_{0}\approx \displaystyle \frac{\unicode[STIX]{x1D706}_{0}-\unicode[STIX]{x1D706}_{-1}}{ds}.\end{eqnarray}$$

Having ${\dot{x}}_{0}$ and $\dot{\unicode[STIX]{x1D706}}_{0}$ , we can use any increment $ds$ in $s$ we like, and replace ${\dot{x}}^{2}+\dot{\unicode[STIX]{x1D706}}^{2}-1=0$ with

$$\begin{eqnarray}{\dot{x}}_{0}(x-x_{0})+\dot{\unicode[STIX]{x1D706}}_{0}(\unicode[STIX]{x1D706}-\unicode[STIX]{x1D706}_{0})-ds=0.\end{eqnarray}$$

So, the equation $\mathbf{G}(\mathbf{z},s)=0$ changes as we increment $s$ because the approximation of ${\dot{z}}$ depends on the current point in the path. However, for sufficiently small $ds$ this approximation works very well, and the implicit function theorem applies at every stage of the continuation. The analysis is a bit trickier in several variables, but this simple problem captures all the essential ideas.

In several variables we must require that either $\mathbf{F}_{\mathbf{x}}$ be non-singular or that the singularity be a simple fold.

For our simple scalar example, the first two conditions in the definition are trivially true. Since $f_{x}(0,0)=0$ and $f_{\unicode[STIX]{x1D706}}(0,0)=-1$ , the third condition holds as well.

Pseudo-arclength continuation in the case of several variables proceeds in the same way as our scalar example. We seek to advance in arclength by $ds$ from a point $(\mathbf{x}_{0},\unicode[STIX]{x1D706}_{0})$ on the solution path. We do this by solving the expanded system for $\mathbf{z}=(\mathbf{x}^{T},\unicode[STIX]{x1D706})^{T}$ :

(4.4) $$\begin{eqnarray}\mathbf{G}(\mathbf{z},s)\equiv \left(\begin{array}{@{}c@{}}\mathbf{F}(\mathbf{x},\unicode[STIX]{x1D706})\\ \mathbf{N}(\mathbf{z},s)\end{array}\right)=\left(\begin{array}{@{}c@{}}0\\ 0\end{array}\right)\!.\end{eqnarray}$$

In (4.4),

$$\begin{eqnarray}\mathbf{N}(\mathbf{z},s)\approx \Vert \dot{\mathbf{x}}\Vert ^{2}+\dot{\unicode[STIX]{x1D706}}^{2}-1\end{eqnarray}$$

is a normalization term. One common choice is the secant normalization

(4.5) $$\begin{eqnarray}\displaystyle \mathbf{N}(\mathbf{z},s) & = & \displaystyle \unicode[STIX]{x1D703}\biggl(\displaystyle \frac{\mathbf{x}(s_{0})-\mathbf{x}(s_{-1})}{s_{0}-s_{-1}}\biggr)^{T}(\mathbf{x}(s)-\mathbf{x}(s_{0}))\nonumber\\ \displaystyle & & \displaystyle \quad +(1-\unicode[STIX]{x1D703})\biggl(\displaystyle \frac{\unicode[STIX]{x1D706}(s_{0})-\unicode[STIX]{x1D706}(s_{-1})}{s_{0}-s_{-1}}\biggr)(\unicode[STIX]{x1D706}(s)-\unicode[STIX]{x1D706}(s_{0}))-(s-s_{0}),\end{eqnarray}$$

which we can use when two points on the path are available. To start the continuation, a typical choice is the norm-based normalization

(4.6) $$\begin{eqnarray}\mathbf{N}(\mathbf{z},s)=\unicode[STIX]{x1D703}\Vert \mathbf{x}-\mathbf{x}(s_{0})\Vert ^{2}+(1-\unicode[STIX]{x1D703})|\unicode[STIX]{x1D706}-\unicode[STIX]{x1D706}(s_{0})|^{2}-(s-s_{0})^{2}.\end{eqnarray}$$

The parameter $\unicode[STIX]{x1D703}$ is a scaling parameter that balances the size of the vector norm with the scalar parameter $\unicode[STIX]{x1D706}$ .

Summarizing, algorithm arclength_continuation is a formal description of pseudo-arclength continuation.

We describe the continuation method as a theorem.

Theorem 4.3 asserts that if the parameter is changed from $\unicode[STIX]{x1D706}$ to $s$ , then the singularity has been eliminated and the path of solutions is homeomorphic to a line segment (Crandall and Rabinowitz Reference Crandall and Rabinowitz1971). In that event, all one has to do in order to follow the path of solutions is to apply algorithm arclength_continuation. In addition, one can prove a mesh independence result (Ferng and Kelley Reference Ferng and Kelley2000) for discretizations of continuous problems.

As you might imagine, there are a few details to resolve. The parameter $\unicode[STIX]{x1D703}$ in (4.5) and (4.6) must be used to maintain mesh independence if one is discretizing a continuous problem. The reason for this is that if one uses a solver which is based on the discrete $\ell ^{2}$ -norm, that norm does not converge to the integral $L^{2}$ -norm as the mesh is refined. In the case of the midpoint rule, for example,

$$\begin{eqnarray}\int _{0}^{1}f(\unicode[STIX]{x1D707})\,\text{d}\unicode[STIX]{x1D707}\approx \displaystyle \frac{1}{N}\mathop{\sum }_{i=1}^{N}f(\unicode[STIX]{x1D707}_{i}),\end{eqnarray}$$

and one would use $\unicode[STIX]{x1D703}=1/N$ , where $N$ is the number of spatial mesh points. Without this the discrete approximation would not be consistent with the continuous problem.

Another important detail is the continuation itself. The theorem limits the range of $s$ for a good reason. In some cases, the $H$ -equation being one of them, $\unicode[STIX]{x1D706}(s)\rightarrow \unicode[STIX]{x1D706}_{\infty }=0$ as $s\rightarrow \infty$ while $\Vert \mathbf{x}(s)\Vert \rightarrow \infty$ . In this case one needs to reduce $ds$ as the continuation progresses in order to ensure that the solutions stay on the path. Even with a line search, there is no guarantee that the solutions will stay on the path (as opposed to jumping to a different solution branch) if the initial iterate is poor. Hence we limit the range of $s$ in the theorem so that one choice of $ds$ will suffice.

Those readers who are familiar with initial value problems will not be surprised to hear that the initial iterate for the solve step in algorithm arclength_continuation is called the predictor. The algorithm uses the trivial predictor, that is, $\mathbf{z}(s)$ as the initial iterate for $\mathbf{z}(s+ds)$ . Since one has an estimate for $\dot{\mathbf{z}}(s)$ ,

$$\begin{eqnarray}\dot{\mathbf{z}}(s)\approx ((\mathbf{x}(s)-\mathbf{x}(s-ds))^{T},\unicode[STIX]{x1D706}(s)-\unicode[STIX]{x1D706}(s-ds)^{T}/ds,\end{eqnarray}$$

one could use linear extrapolation to form the linear predictor

(4.7) $$\begin{eqnarray}\mathbf{z}_{0}(s+ds)=\mathbf{z}(s)+\dot{\mathbf{z}}(s)ds,\end{eqnarray}$$

which generally performs better. Finally, similarly to numerical integration of initial value problems, it is a bad sign if the nonlinear solver requires too many iterations or a line search. Those events are signals that the predictor is poor and that one should reduce $ds$ .

In our numerical example we used $ds=0.002$ and let the continuation run until the line search failed because the limit on stepsize reductions had been exceeded. One would expect that to happen at some point because the solutions are moving farther apart and the predictor’s performance is therefore becoming worse as the continuation progresses.

We will illustrate the output for the $H$ -equation. In the example we used Newton–GMRES for the nonlinear solver and a forward-difference Jacobian-vector product. We used the forcing term from (2.24). We would expect from (2.34) that there are two solutions for $0<\unicode[STIX]{x1D714}<1$ , one for each choice of sign for the square root. The continuation computation confirms that.

Typically one plots the progress of a continuation with $\unicode[STIX]{x1D706}$ on the horizontal axis and a functional of $\mathbf{x}$ on the vertical axis. For the $H$ -equation and its discretization the function is strictly positive (in fact $\geq 1$ ) and increasing. One could use either $H(1)$ (or $h_{N}$ in the discrete case) or the $L^{1}$ -norm

$$\begin{eqnarray}\int _{0}^{1}H(\unicode[STIX]{x1D707})\,\text{d}\unicode[STIX]{x1D707}\quad \text{or}\quad \displaystyle \frac{1}{N}\mathop{\sum }_{i=1}^{N}h_{i}.\end{eqnarray}$$

We use $h_{N}=\Vert \mathbf{h}\Vert _{\infty }$ to illustrate the steep rise in the size of the solution. See Figure 4.1.

We used the secant normalization for all but the first point on the path and exploit the fact that if $\unicode[STIX]{x1D706}(0)=0$ then $H\equiv 1$ . This is also true for the discrete case. Hence there is no need to solve the equation at the first point on the path.

Figure 4.1.

Solution path for $H$ -equation.

4.3 The Bratu problem

In this section we illustrate how multiple solutions of a nonlinear equation can have very different properties when considered as steady-state solutions of a time-dependent problem. Consider the system of ordinary differential equations

(4.8) $$\begin{eqnarray}\dot{\mathbf{x}}=-\mathbf{F}(\mathbf{x}),\end{eqnarray}$$

where $\dot{\mathbf{x}}=\,\text{d}\mathbf{x}/\text{d}t$ . We chose the sign of $\mathbf{F}$ to be consistent with the standard practice in pseudo-transient continuation (see Section 7).

A solution $\mathbf{x}^{\ast }$ is a steady-state solution if it is independent of time. In that case (4.8) implies that

$$\begin{eqnarray}0=\dot{\mathbf{x}}^{\ast }=-\mathbf{F}(\mathbf{x}^{\ast }).\end{eqnarray}$$

Here $\mathbf{x}^{\ast }$ is dynamically stable if the solution $\mathbf{x}$ of the initial value problem

$$\begin{eqnarray}\dot{\mathbf{x}}=-\mathbf{F}(\mathbf{x}),\quad \mathbf{x}(0)=\mathbf{x}_{0}\end{eqnarray}$$

converges to $\mathbf{x}^{\ast }$ as $t\rightarrow \infty$ for all $\mathbf{x}_{0}$ sufficiently close to $\mathbf{x}^{\ast }$ . It is often easier to check for linear stability, which is sufficient (but not necessary) for stability. We say $\mathbf{x}^{\ast }$ is linearly stable if the eigenvalues of $\mathbf{F}^{\prime }(\mathbf{x}^{\ast })$ all have positive real part. Not all steady-state solutions are stable, and algorithms such as newton_armijo are not aware of dynamics and can (and sometimes do) converge to unstable steady states. We will give an example of this in Section 7.

The Bratu problem (Bratu Reference Bratu1914) is an example of this. The steady-state problem is the two-point boundary value problem

(4.9) $$\begin{eqnarray}-u_{xx}=\unicode[STIX]{x1D706}\,\text{e}^{u},\quad u(0)=u(1)=0.\end{eqnarray}$$

We would express the boundary value problem as $-{\mathcal{F}}(u)=0$ in the function space. The Fréchet derivative of ${\mathcal{F}}(u)$ is clearly symmetric and positive definite if $\unicode[STIX]{x1D706}$ satisfies

$$\begin{eqnarray}\unicode[STIX]{x1D70B}^{2}-\unicode[STIX]{x1D706}\,\text{e}^{u}>0.\end{eqnarray}$$

If $\unicode[STIX]{x1D706}=0$ , the unique solution of (4.9) is $u\equiv 0$ . The Fréchet derivative ${\mathcal{F}}_{u}$ at $(u,\unicode[STIX]{x1D706})=(0,0)$ is $-\text{d}^{2}/\text{d}x^{2}$ with homogeneous Dirichlet boundary conditions and is positive definite. Hence the implicit function implies that there is a solution of (4.9) for sufficiently small $\unicode[STIX]{x1D706}>0$ and that the solution will be linearly stable as long as the smallest eigenvalue of ${\mathcal{F}}_{x}$ is positive. As the continuation progresses, the eigenvalue will change sign at the singularity of ${\mathcal{F}}_{x}$ and the solution will lose linear stability (if it exists at all). The structure of the Bratu problem (and the $H$ -equation as well) implies that the eigenvalue will change sign when one passes the singularity on the continuation path (Keller Reference Keller1987).

One can solve (4.9) explicitly (Ascher, Mattheij and Russell Reference Ascher, Mattheij and Russell1995). The solution is

$$\begin{eqnarray}u(x)=-2\ln \biggl(\displaystyle \frac{\cosh [(x-1/2)\unicode[STIX]{x1D703}/2]}{\cosh (\unicode[STIX]{x1D703}/4)}\biggr),\end{eqnarray}$$

where $\unicode[STIX]{x1D703}$ is the solution of the scalar equation

(4.10) $$\begin{eqnarray}\unicode[STIX]{x1D703}=\sqrt{2\unicode[STIX]{x1D706}}\cosh (\unicode[STIX]{x1D703}/4).\end{eqnarray}$$

Similarly to (2.34), (4.10) has two solutions for $0<\unicode[STIX]{x1D706}<\unicode[STIX]{x1D706}^{\ast }$ and no real solution for $\unicode[STIX]{x1D706}>\unicode[STIX]{x1D706}^{\ast }\approx 3.52$ . The continuation plot for the Bratu problem is similar to that for the $H$ -equation and we invite the reader to do that computation.

5.1 Algorithmic description

Algorithm anderson is a formal description of the method and is the one we will use for analysis. Implementation is a different matter, and there are many examples of efficient implementations (Walker and Ni Reference Walker and Ni2011, Toth and Kelley Reference Toth and Kelley2015, Collier et al. Reference Collier, Hindmarsh, Serban and Woodward2015, Hindmarsh et al. Reference Hindmarsh, Brown, Grant, Lee, Serban, Shumaker and Woodward2005, Toth and Pawlowski Reference Toth and Pawlowski2015, Toth Reference Toth2016).

The parameter $\unicode[STIX]{x1D6FD}$ is called the mixing parameter. This is the same as the damping parameter in Picard iteration. In many cases one must damp Picard iteration to secure convergence. The damped Picard iteration is

$$\begin{eqnarray}\mathbf{x}_{n+1}=(1-\unicode[STIX]{x1D6FD})\mathbf{x}_{n}+\unicode[STIX]{x1D6FD}\mathbf{G}(\mathbf{x}_{n})\equiv \mathbf{G}_{\unicode[STIX]{x1D6FD}}(\mathbf{x}_{n}).\end{eqnarray}$$

Anderson acceleration with mixing parameter $\unicode[STIX]{x1D6FD}$ is the same as applying the algorithm with $\unicode[STIX]{x1D6FD}=1$ to the map $\mathbf{G}_{\unicode[STIX]{x1D6FD}}$ . Hence there is no loss of generality in setting $\unicode[STIX]{x1D6FD}\equiv 1$ for analysis.

Anderson maintains a limited history of the iteration of size $m+1$ ; $m$ is called the depth. The iteration uses the most recent $m+1$ residuals $\mathbf{F}(\mathbf{x}_{j})$ for $k-m_{k}\leq j\leq k$ where $m_{k}\leq \min (k,m)$ . The key step in the iteration is solving the optimization problem

(5.5) $$\begin{eqnarray}\min \biggl\|\mathop{\sum }_{j=0}^{m_{k}}\unicode[STIX]{x1D6FC}_{j}^{k}\mathbf{F}(\mathbf{x}_{k-m_{k}+j})\biggr\|\quad \text{subject to }\mathop{\sum }_{j=0}^{m_{k}}\unicode[STIX]{x1D6FC}_{j}^{k}=1,\end{eqnarray}$$

for the coefficients $\{\unicode[STIX]{x1D6FC}_{j}^{k}\}$ .

Any vector norm can be used in the optimization problem with no change in the theory. The optimization problem is easier to solve if one uses the $\ell ^{2}$ -norm, and that is standard practice. In this case the optimization problem for the coefficients can be expressed as a linear least-squares problem and solved very inexpensively. One way to do this is to solve the linear least-squares problem

(5.6) $$\begin{eqnarray}\min \biggl\|\mathbf{F}(\mathbf{x}_{k})-\mathop{\sum }_{j=0}^{m_{k}-1}\unicode[STIX]{x1D6FC}_{j}^{k}(\mathbf{F}(\mathbf{x}_{k-m_{k}+j})-\mathbf{F}(\mathbf{x}_{k}))\biggr\|_{2}^{2},\end{eqnarray}$$

for $\{\unicode[STIX]{x1D6FC}_{j}^{k}\}_{j=0}^{m_{k}-1}$ . Then one recovers $\unicode[STIX]{x1D6FC}_{m_{k}}^{k}$ by

$$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{m_{k}}^{k}=1-\mathop{\sum }_{j=0}^{m_{k}-1}\unicode[STIX]{x1D6FC}_{j}^{k}.\end{eqnarray}$$

Toth and Kelley (Reference Toth and Kelley2015) point out that other norms could be used. In particular, the optimization problem for the coefficients in either the $\ell ^{1}$ - or $\ell ^{\infty }$ -norms can be formulated as a linear programming problem and solved directly with many codes (CVX Research 2012). However, the least-squares approach using (5.6) is more efficient, and we will use the $\ell ^{2}$ -norm only in this article.

The choice of $m_{k}$ is, in the original form, simply $\min (m,k)$ . One can adapt $m_{k}$ as the iteration progresses to, for example, enforce well-conditioning of the linear least-squares problem (5.6) (Walker and Ni Reference Walker and Ni2011, Toth Reference Toth2016, An et al. Reference An, Jia and Walker2017).

One can show (Fang and Saad Reference Fang and Saad2009, Saad, Chelikowsky and Shontz Reference Saad, Chelikowsky and Shontz2010, Rohwedder and Schneider Reference Rohwedder and Schneider2011, Walker and Ni Reference Walker and Ni2011, Potra and Engler Reference Potra and Engler2013) that Anderson acceleration is related to multisecant quasi-Newton methods or, in the case of linear problems, GMRES. None of these results lead to a convergence proof, even in the linear case, unless the available storage is large enough to allow GMRES to take a number of iterations equal to the dimension of the problem.

One result from Walker and Ni (Reference Walker and Ni2011) illustrates the power of unlimited storage. While not the case seen in practice, this result does illustrate why Anderson acceleration may perform better than Picard iteration in some cases. Unlike the remainder of the results we present in this section, contractivity is not necessary.

This result says that, under most circumstances, Anderson acceleration with $m\geq N$ performs exactly as well as GMRES and that the analysis of preconditioning is the same as for GMRES. GMRES is well known (Nevanlinna Reference Nevanlinna1993, Campbell, Ipsen, Kelley and Meyer Reference Campbell, Ipsen, Kelley and Meyer1996a , Campbell et al. Reference Campbell, Ipsen, Kelley, Meyer and Xue1996b ) to converge rapidly for discretizations of second-kind Fredholm integral equations, for example.

5.2 Convergence theory

Toth and Kelley (Reference Toth and Kelley2015) has the first convergence analysis of the method as used in practice. The central idea in Toth and Kelley (Reference Toth and Kelley2015) was to show that Anderson acceleration does no harm rather than to prove a general convergence result. The results are consistent with the observations from computational chemistry (Foresman and Frisch Reference Foresman and Frisch1996). For example, a good initial iterate is needed for convergence in many cases.

5.2.1 Linear problems

We will begin with an analysis of the linear case. The convergence theory, at least as it stands today, only shows that Anderson acceleration does not degrade the convergence of Picard iteration. However, in practice Anderson acceleration is often (but not always) much better. There is at present no satisfactory characterization of problems for which Anderson is better.

Theorem 5.2. Let $\mathbf{M}$ be an $N\times N$ matrix with $\Vert \mathbf{M}\Vert =c<1$ . Let $m\geq 0$ . Then Anderson( $m$ ) acceleration, when applied to $\mathbf{G}(\mathbf{x})=\mathbf{M}\mathbf{x}+\mathbf{b}$ , converges to the solution $\mathbf{x}^{\ast }=(\mathbf{I}-\mathbf{M})^{-1}\mathbf{b}$ . Moreover, the residuals $\mathbf{F}(\mathbf{x})=\mathbf{b}-(\mathbf{I}-\mathbf{M})\mathbf{x}$ converge to zero with a q-factor no larger than $c$ .

Proof. In this proof the optimization problem is used in an important way. Given $\mathbf{x}_{k}$ , we note that since $\sum \unicode[STIX]{x1D6FC}_{j}^{k}=1$ , the new residual is

$$\begin{eqnarray}\displaystyle \mathbf{F}(\mathbf{x}_{k+1}) & = & \displaystyle b-(I-\mathbf{M})\mathbf{x}_{k+1}\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{j=0}^{m_{k}}\unicode[STIX]{x1D6FC}_{j}^{k}[b-(I-\mathbf{M})(b+\mathbf{M}\mathbf{x}_{k-m_{k}+j})]\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{j=0}^{m_{k}}\unicode[STIX]{x1D6FC}_{j}^{k}\mathbf{M}[b-(I-\mathbf{M})\mathbf{x}_{k-m_{k}+j}]\nonumber\\ \displaystyle & = & \displaystyle \mathbf{M}\mathop{\sum }_{j=0}^{m_{k}}\unicode[STIX]{x1D6FC}_{j}^{k}\mathbf{F}(\mathbf{x}_{k-m_{k}+j}).\nonumber\end{eqnarray}$$

We take norms and use $\Vert \mathbf{M}\Vert =c$ to obtain

$$\begin{eqnarray}\Vert \mathbf{F}(\mathbf{x}_{k+1})\Vert \leq c\biggl\|\mathop{\sum }_{j=0}^{m_{k}}\unicode[STIX]{x1D6FC}_{j}\mathbf{F}(\mathbf{x}_{k-m_{k}+j})\biggr\|.\end{eqnarray}$$

Optimality implies that

$$\begin{eqnarray}\biggl\|\mathop{\sum }_{j=0}^{m_{k}}\unicode[STIX]{x1D6FC}_{j}\mathbf{F}(\mathbf{x}_{k-m_{k}+j})\biggr\|\leq \Vert \mathbf{F}(\mathbf{x}_{k})\Vert .\end{eqnarray}$$

Hence

(5.7) $$\begin{eqnarray}\Vert \mathbf{F}(\mathbf{x}_{k+1})\Vert \leq c\Vert \mathbf{F}(\mathbf{x}_{k})\Vert ,\end{eqnarray}$$

as asserted. ☐

One might think that the analysis could proceed like that for Newton’s method in that the result for the linear problem (convergence in one iteration for Newton) would imply a result for the nonlinear problem after a Taylor expansion if the initial iterate were accurate enough to neglect the high-order terms. In fact, that analogy is correct. We will illustrate the point with two theorems from Toth and Kelley (Reference Toth and Kelley2015).

In the special case $m=1$ , one can solve the optimization problem analytically. One can use this to show that (5.7) holds for Anderson(1) in the nonlinear case if the initial iterate is sufficiently near the solution. The assumption of continuous differentiability is weaker than the one in Toth and Kelley (Reference Toth and Kelley2015) and is also used in Chen and Kelley (Reference Chen and Kelley2017).

5.2.2 The special case $m=1$

If $m=1$ we need only assume that $\mathbf{G}$ is a continuously differentiable contraction to obtain q-linear convergence of the residuals.

Assumption 5.3. $\mathbf{G}$ has a fixed point $\mathbf{x}^{\ast }$ .

• ∙ $\mathbf{G}$ is continuously differentiable in the ball ${\mathcal{B}}(\mathbf{x}^{\ast },\hat{\unicode[STIX]{x1D70C}})=\{\mathbf{x}\mid \Vert \mathbf{e}\Vert \leq \hat{\unicode[STIX]{x1D70C}}\}$ for some $\hat{\unicode[STIX]{x1D70C}}>0$ .

• ∙ There is $c\in (0,1)$ such that for all $\mathbf{x},\mathbf{y}\in {\mathcal{B}}(\mathbf{x}^{\ast },\hat{\unicode[STIX]{x1D70C}})$ , $\Vert \mathbf{G}(\mathbf{x})-\mathbf{G}(\mathbf{y})\Vert \leq c\Vert \mathbf{x}-\mathbf{y}\Vert$ .

Theorem 5.4 is a generalization of a result from Toth and Kelley (Reference Toth and Kelley2015) with stronger convergence and slightly weaker assumptions.

Theorem 5.4. Assume that Assumption 5.3 holds. Then if $\Vert e_{0}\Vert$ is sufficiently small, the Anderson(1) residuals with $\ell ^{2}$ -optimization converge q-linearly. Moreover,

(5.8) $$\begin{eqnarray}\limsup _{k\rightarrow \infty }\displaystyle \frac{\Vert \mathbf{F}(\mathbf{x}_{k+1})\Vert }{\Vert \mathbf{F}(\mathbf{x}_{k})\Vert }\leq c.\end{eqnarray}$$

Proof. We will express the iteration as

(5.9) $$\begin{eqnarray}\mathbf{x}_{k+1}=(1-\unicode[STIX]{x1D6FC}^{k})\mathbf{G}(\mathbf{x}_{k})+\unicode[STIX]{x1D6FC}^{k}\mathbf{G}(\mathbf{x}_{k-1}),\end{eqnarray}$$

and note that

(5.10) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}^{k}=\displaystyle \frac{\mathbf{F}(\mathbf{x}_{k})^{T}(\mathbf{F}(\mathbf{x}_{k})-\mathbf{F}(\mathbf{x}_{k-1}))}{\Vert \mathbf{F}(\mathbf{x}_{k})-\mathbf{F}(\mathbf{x}_{k-1})\Vert ^{2}}.\end{eqnarray}$$

Now define

$$\begin{eqnarray}\mathbf{a}_{k}=\mathbf{G}(\mathbf{x}_{k+1})-\mathbf{G}((1-\unicode[STIX]{x1D6FC}^{k})\mathbf{x}_{k}+\unicode[STIX]{x1D6FC}^{k}\mathbf{x}_{k-1})\end{eqnarray}$$

and

$$\begin{eqnarray}\mathbf{b}_{k}=\mathbf{G}((1-\unicode[STIX]{x1D6FC}^{k})\mathbf{x}_{k}+\unicode[STIX]{x1D6FC}^{k}\mathbf{x}_{k-1})-\mathbf{x}_{k+1}.\end{eqnarray}$$

Clearly

(5.11) $$\begin{eqnarray}\mathbf{F}(\mathbf{x}_{k+1})=\mathbf{G}(\mathbf{x}_{k+1})-\mathbf{x}_{k+1}=\mathbf{a}_{k}+\mathbf{b}_{k}.\end{eqnarray}$$

We will obtain an estimate of $\mathbf{F}(\mathbf{x}_{k+1})$ by estimating $\mathbf{a}_{k}$ and $\mathbf{b}_{k}$ separately.

By definition of the Anderson iteration (5.9) and contractivity of $G$ ,

(5.12) $$\begin{eqnarray}\displaystyle \Vert \mathbf{a}_{k}\Vert & = & \displaystyle \Vert \mathbf{G}(\mathbf{x}_{k+1})-\mathbf{G}((1-\unicode[STIX]{x1D6FC}^{k})\mathbf{x}_{k}+\unicode[STIX]{x1D6FC}^{k}\mathbf{x}_{k-1})\Vert \nonumber\\ \displaystyle & \leq & \displaystyle c\Vert \mathbf{x}_{k+1}-(1-\unicode[STIX]{x1D6FC}^{k})\mathbf{x}_{k}-\unicode[STIX]{x1D6FC}^{k}\mathbf{x}_{k-1}\Vert \nonumber\\ \displaystyle & = & \displaystyle c\Vert (1-\unicode[STIX]{x1D6FC}^{K})(\mathbf{G}(\mathbf{x}_{K})-\mathbf{x}_{K})-\unicode[STIX]{x1D6FC}^{K}(\mathbf{G}(\mathbf{x}_{K-1})-\mathbf{x}_{K-1})\Vert \nonumber\\ \displaystyle & = & \displaystyle c\Vert (1-\unicode[STIX]{x1D6FC}^{k})\mathbf{F}(\mathbf{x}_{k})-\unicode[STIX]{x1D6FC}^{k}\mathbf{F}(\mathbf{x}_{k-1})\Vert \leq c\Vert \mathbf{F}(\mathbf{x}_{k})\Vert ,\end{eqnarray}$$

where the last inequality follows from the optimality property of the coefficients.

We now estimate $\mathbf{b}_{k}$ . Now let

$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{k}=\mathbf{x}_{k-1}-\mathbf{x}_{k}.\end{eqnarray}$$

To estimate $\mathbf{b}_{k}$ we note that

(5.13) $$\begin{eqnarray}\displaystyle \mathbf{b}_{k} & = & \displaystyle \mathbf{G}((1-\unicode[STIX]{x1D6FC}^{k})\mathbf{x}_{k}+\unicode[STIX]{x1D6FC}^{k}\mathbf{x}_{k-1})-(1-\unicode[STIX]{x1D6FC}^{k})\mathbf{G}(\mathbf{x}_{k})-\unicode[STIX]{x1D6FC}^{k}\mathbf{G}(\mathbf{x}_{k-1})\nonumber\\ \displaystyle & = & \displaystyle \mathbf{G}(\mathbf{x}_{k}+\unicode[STIX]{x1D6FC}^{k}\unicode[STIX]{x1D6FF}_{k})-\mathbf{G}(\mathbf{x}_{k})+\unicode[STIX]{x1D6FC}^{k}(\mathbf{G}(\mathbf{x}_{k})-\mathbf{G}(\mathbf{x}_{k-1}))\nonumber\\ \displaystyle & = & \displaystyle \int _{0}^{1}\mathbf{G}^{\prime }(\mathbf{x}_{k}+t\unicode[STIX]{x1D6FC}^{k}\unicode[STIX]{x1D6FF}_{K})\unicode[STIX]{x1D6FC}^{k}\unicode[STIX]{x1D6FF}_{k}\,\text{d}t-\unicode[STIX]{x1D6FC}^{K}\int _{0}^{1}\mathbf{G}^{\prime }(\mathbf{x}_{k}+t\unicode[STIX]{x1D6FF}_{k})\unicode[STIX]{x1D6FF}_{k}\,\text{d}t\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D6FC}^{k}\int _{0}^{1}[\mathbf{G}^{\prime }(\mathbf{x}_{k}+t\unicode[STIX]{x1D6FC}^{k}\unicode[STIX]{x1D6FF}_{k})-\mathbf{G}^{\prime }(\mathbf{x}_{k}+t\unicode[STIX]{x1D6FF}_{k})]\unicode[STIX]{x1D6FF}_{k}\,\text{d}t.\end{eqnarray}$$

Since $\mathbf{G}^{\prime }$ is continuous in ${\mathcal{B}}(\mathbf{x}^{\ast },\unicode[STIX]{x1D70C})$ , there is a non-decreasing function $\unicode[STIX]{x1D702}\in C[0,\infty )$ with $\unicode[STIX]{x1D702}(0)=0$ so that

(5.14) $$\begin{eqnarray}\Vert \mathbf{G}^{\prime }(\mathbf{x})-\mathbf{G}^{\prime }(\mathbf{x}^{\ast })\Vert \leq \unicode[STIX]{x1D702}(\Vert \mathbf{e}\Vert )\end{eqnarray}$$

for all $\mathbf{x}\in {\mathcal{B}}(\mathbf{x}^{\ast },\unicode[STIX]{x1D70C})$ . Hence, if $\mathbf{x}_{k}$ and $\mathbf{x}_{k-1}$ are both in ${\mathcal{B}}(\mathbf{x}^{\ast },\unicode[STIX]{x1D70C})$ (which is certainly true if $k=1$ ), then

(5.15) $$\begin{eqnarray}\Vert \mathbf{b}_{k}\Vert \leq 2\unicode[STIX]{x1D702}(\Vert \mathbf{e}_{k}\Vert +\Vert \unicode[STIX]{x1D6FF}_{k}\Vert )\Vert |\unicode[STIX]{x1D6FC}^{k}|\Vert \unicode[STIX]{x1D6FF}_{k}\Vert .\end{eqnarray}$$

Because $m=1$ , it is not difficult to estimate $\unicode[STIX]{x1D6FC}^{k}$ . Clearly,

$$\begin{eqnarray}\displaystyle \mathbf{F}(\mathbf{x}_{k})-\mathbf{F}(\mathbf{x}_{k-1}) & = & \displaystyle \mathbf{G}(\mathbf{x}_{k})-\mathbf{G}(\mathbf{x}_{k-1})+\unicode[STIX]{x1D6FF}_{k}=\unicode[STIX]{x1D6FF}_{k}-\int _{0}^{1}\mathbf{G}^{\prime }(\mathbf{x}_{k-1}-t\unicode[STIX]{x1D6FF}_{k})\unicode[STIX]{x1D6FF}_{k}\,\text{d}t\nonumber\\ \displaystyle & = & \displaystyle \biggl(I-\int _{0}^{1}\mathbf{G}^{\prime }(\mathbf{x}_{k-1}-t\unicode[STIX]{x1D6FF}_{k})\,\text{d}t\biggr)\unicode[STIX]{x1D6FF}_{k}.\nonumber\end{eqnarray}$$

Since $\Vert \mathbf{G}^{\prime }(u)\Vert \leq c$ for all $u\in {\mathcal{B}}(\mathbf{x}^{\ast },\unicode[STIX]{x1D70C})$ , we have

(5.16) $$\begin{eqnarray}\Vert \unicode[STIX]{x1D6FF}_{k}\Vert \leq \Vert \mathbf{F}(\mathbf{x}_{k})-\mathbf{F}(\mathbf{x}_{k-1})\Vert /(1-c).\end{eqnarray}$$

Combine (5.16) and (5.10) to obtain

(5.17) $$\begin{eqnarray}|\unicode[STIX]{x1D6FC}^{k}|\Vert \unicode[STIX]{x1D6FF}_{k}\Vert \leq \displaystyle \frac{\Vert \mathbf{F}(\mathbf{x}_{k})\Vert }{\Vert \mathbf{F}(\mathbf{x}_{k})-\mathbf{F}(\mathbf{x}_{k-1})\Vert }\Vert \unicode[STIX]{x1D6FF}_{k}\Vert \leq \displaystyle \frac{\Vert \mathbf{F}(\mathbf{x}_{k})\Vert }{1-c}.\end{eqnarray}$$

Hence

(5.18) $$\begin{eqnarray}\Vert \mathbf{F}(\mathbf{x}_{k+1})\Vert \leq \Vert \mathbf{F}(\mathbf{x}_{k})\Vert (c+2\unicode[STIX]{x1D702}(\Vert \mathbf{e}_{k}\Vert +\Vert \unicode[STIX]{x1D6FF}_{k}\Vert )\Vert ).\end{eqnarray}$$

So, given $c<{\hat{c}}<1$ , we may reduce $\unicode[STIX]{x1D70C}$ if needed so that

(5.19) $$\begin{eqnarray}\Vert \mathbf{F}(\mathbf{x}_{k+1})\Vert \leq {\hat{c}}\Vert \mathbf{F}(\mathbf{x}_{k})\Vert .\end{eqnarray}$$

We complete the proof by combining (5.18) and (5.19) to show that

$$\begin{eqnarray}\Vert \mathbf{F}(\mathbf{x}_{k+1})\Vert \leq \Vert \mathbf{F}(\mathbf{x}_{k})\Vert (c+o(1))\end{eqnarray}$$

as $k\rightarrow \infty$ .☐

5.2.3 Local convergence for general $m$

We must either assume or arrange that the $\ell ^{1}$ -norm of the coefficients be uniformly bounded to prove local convergence for $m>1$ or for any norm other than the $\ell ^{2}$ -norm. In addition, the convergence is r-linear rather than q-linear. We will state Chen and Kelley’s (Reference Chen and Kelley2017) extension of the local convergence result from Toth and Kelley (Reference Toth and Kelley2015). We refer to those papers for the proof.

Theorem 5.5. Assume that Assumption 5.3 holds. Assume that there is $M_{\unicode[STIX]{x1D6FC}}$ such that, for all $k\geq 0$ ,

(5.20) $$\begin{eqnarray}\mathop{\sum }_{j=1}^{m_{k}}|\unicode[STIX]{x1D6FC}_{j}|\leq M_{\unicode[STIX]{x1D6FC}}.\end{eqnarray}$$

Then if $\mathbf{x}_{0}$ is sufficiently near to $\mathbf{x}^{\ast }$ , the Anderson iterations converge, and

(5.21) $$\begin{eqnarray}\limsup _{k\rightarrow \infty }\biggl(\displaystyle \frac{\Vert \mathbf{F}(\mathbf{x}_{k})\Vert }{\Vert \mathbf{F}(\mathbf{x}_{0})\Vert }\biggr)^{1/k}\leq c.\end{eqnarray}$$

The assumption (5.20) that the $\ell ^{1}$ -norm of the coefficients is bounded can be enforced within the iteration by controlling $m_{k}$ . One way to do this is to reduce $m_{k}$ if the $\ell ^{1}$ -norm of the coefficients or (in the $\ell ^{2}$ -norm case) the conditioning of the least-squares problem (5.6) exceeds a predetermined bound. Walker and Ni (Reference Walker and Ni2011) and An et al. (Reference An, Jia and Walker2017) advocate limiting the condition number of the least-squares problem (5.6). Toth (Reference Toth2016) shows that if one does this the r-linear convergence improves to q-linear. However, the performance in practice of methods that enforce (5.20) is not always better than the original version. In Section 5.3 we discuss another approach which requires that the coefficients be non-negative (so $M_{\unicode[STIX]{x1D6FC}}\equiv 1$ ) as a way to globalize convergence.

5.2.4 $H$ -equation example

We return to the $H$ -equation with some of the results from Toth and Kelley (Reference Toth and Kelley2015). This computation exposes one of the many mysteries in Anderson acceleration.

We begin with a look at Newton–GMRES and Picard iteration, which are well known to converge even in the singular $c=1$ case. If $c=1$ the fixed-point map is not a contraction and, as we have seen, the Jacobian is singular at the solution. However, if the initial iterate is well chosen ( $H_{0}(\unicode[STIX]{x1D707})\equiv 1$ is a good choice), both the Picard and Newton iterations will converge.

The computations in this section are from an $N=500$ point mesh. We terminated the iterations when $\Vert \mathbf{F}(\mathbf{h}_{n})\Vert \leq 10^{-8}\Vert \mathbf{F}(\mathbf{h}_{0})\Vert$ . Table 5.1 reports cost in terms of function evaluations for both Newton–GMRES and Picard iterations for three values if $\unicode[STIX]{x1D714}$ . One can see the effect of the singular Jacobian for both Newton–GMRES and Picard iteration. In the statistics for Newton–GMRES we count both the function evaluations in the nonlinear iterations and those used in a finite-difference Jacobian-vector product.

Table 5.1.

Function evaluations for Newton–GMRES and fixed-point iteration.

In Table 5.2 we report the cost of Anderson( $m$ ) with $m=1,2,5$ for the same problems. We tabulate the number of iterations ITS needed to terminate, the maximum condition number $\unicode[STIX]{x1D705}$ of the least-squares problem, and the maximum $\ell ^{1}$ -norm $S_{\max }$ of the coefficients. $S_{\max }$ is not large and, at least for this problem, the assumption that (5.20) holds is reasonable.

Anderson( $m$ ) does far better for these problems than the theory predicts, costing less than even Newton–GMRES. Note that each iteration of Anderson costs a single function evaluation. For $m=5$ the least-squares problems become very ill-conditioned, which is unsurprising given that $\mathbf{G}^{\prime }$ is an integral operator. This ill-conditioning does not cause the iteration to fail, but does, especially when $\unicode[STIX]{x1D714}=1$ , increase the cost.

Table 5.2.

Iteration statistics for Anderson( $m$ ).

In Figure 5.1 we plot the residual histories for Anderson(1) and the three values of $\unicode[STIX]{x1D714}$ . The reader should compare Figure 5.1 to Figure 2.5.

Figure 5.1.

Anderson(1) for the $H$ -equation example.

Toth et al. (Reference Toth, Ellis, Evans, Hamilton, Kelley, Pawlowski and Slattery2017) contains mesh independence results for Anderson acceleration similar to those in Section 3.3.

5.3 The EDIIS globalization

Anderson acceleration performs poorly for some applications. One example is electronic structure computations for metallic systems where the HOMO–LUMO gap is small (Kudin et al. Reference Kudin, Scuseria and Cancès2002). In such cases one must use a small mixing parameter to ensure convergence. However, a small mixing parameter degrades the performance of the iteration. In addition, an accurate initial iterate is often necessary, and finding an acceptable initial iterate is often problematic.

One attempt to solve these problems for small systems is the EDIIS (energy DIIS) algorithm from Kudin et al. (Reference Kudin, Scuseria and Cancès2002). This is the form of Anderson acceleration in the Gaussian computational chemistry code (Frisch et al. Reference Frisch, Trucks, Schlegel, Scuseria, Robb, Cheeseman, Scalmani, Barone, Mennucci, Petersson, Nakatsuji, Caricato, Li, Hratchian, Izmaylov, Bloino, Zheng, Sonnenberg, Hada, Ehara, Toyota, Fukuda, Hasegawa, Ishida, Nakajima, Honda, Kitao, Nakai, Vreven, Montgomery, Peralta, Ogliaro, Bearpark, Heyd, Brothers, Kudin, Staroverov, Kobayashi, Normand, Raghavachari, Rendell, Burant, Iyengar, Tomasi, Cossi, Rega, Millam, Klene, Knox, Cross, Bakken, Adamo, Jaramillo, Gomperts, Stratmann, Yazyev, Austin, Cammi, Pomelli, Ochterski, Martin, Morokuma, Zakrzewski, Voth, Salvador, Dannenberg, Dapprich, Daniels, Farkas, Foresman, Ortiz, Cioslowski and Fox2009). EDIIS differs from Anderson acceleration in that the fixed-point problem is designed differently to minimize energy and a non-negativity constraint on the coefficients. So, the optimization problem becomes

(5.22) $$\begin{eqnarray}\min \biggl\|\mathop{\sum }_{j=0}^{m_{k}}\unicode[STIX]{x1D6FC}_{j}^{k}\mathbf{F}_{k-m_{k}+j}\biggr\|\quad \text{subject to }\mathop{\sum }_{j=0}^{m_{k}}\unicode[STIX]{x1D6FC}_{j}^{k}=1,~\unicode[STIX]{x1D6FC}_{j}^{k}\geq 0.\end{eqnarray}$$

The optimization problem (5.22) for EDIIS is much harder than the linear least-squares problem (5.6) in algorithm anderson. The approach of Kudin et al. (Reference Kudin, Scuseria and Cancès2002) is a direct examination of the boundary of the feasible simplex, which is not practical for a depth much greater than $m=3$ . Since $m$ is small in practice, expressing the optimization problem as a bound-constrained quadratic program is an efficient alternative. Moré and Toraldo (Reference Moré and Toraldo1991) survey the literature on this topic. For example, a bound-constrained quadratic programming code such as the MINQ code (Neumaier Reference Neumaier1998) is a reasonable choice. However, this approach squares the condition number. The classic method from Golub and Saunders (Reference Golub and Saunders1969) uses an active set method and QR factorization to avoid this.

In Chen and Kelley (Reference Chen and Kelley2017) we prove that adding the non-negativity constraint condition to the optimization makes the algorithm globally convergent. Theorem 5.6 and its proof make that precise.

Theorem 5.6. Let $\mathbf{G}$ be a contraction on a convex $D\subset \mathbb{R}^{N}$ with contractivity constant $c$ . Let $\mathbf{x}^{\ast }$ be the unique fixed point of $\mathbf{g}$ in $D$ . Then, for any $\mathbf{x}_{0}\in D$ , EDIIS $(m)$ converges to $\mathbf{x}^{\ast }$ r-linearly with r-factor

$$\begin{eqnarray}{\hat{c}}=c^{1/(m+1)}.\end{eqnarray}$$

In fact,

(5.23) $$\begin{eqnarray}\Vert \mathbf{e}_{k}\Vert \leq {\hat{c}}^{k}\Vert \mathbf{e}_{0}\Vert .\end{eqnarray}$$

Proof. The proof does not use the optimality condition and only requires that the iteration $\{\mathbf{x}_{k}\}$ have the form

(5.24) $$\begin{eqnarray}\mathbf{x}_{k+1}=\mathop{\sum }_{j=0}^{m_{k}}\unicode[STIX]{x1D6FC}_{j}^{k}\mathbf{G}(\mathbf{x}_{k-m_{k}+j}),\end{eqnarray}$$

where $m_{k}\leq m$ , $\unicode[STIX]{x1D6FC}_{j}^{k}\geq 0$ , and $\sum _{j=0}^{m_{k}}\unicode[STIX]{x1D6FC}_{j}^{k}=1$ .

We induct on $k$ . Clearly (5.23) holds for both $m_{k}=0$ , by definition, and $k=m_{k}=1$ because the iteration in that case is a single Picard iteration. Assume that the result holds for $k\leq K$ . Then (5.24) and $\sum _{j=0}^{m_{k}}\unicode[STIX]{x1D6FC}_{j}^{k}=1$ imply that

$$\begin{eqnarray}\mathbf{e}_{K+1}=\mathop{\sum }_{j=0}^{m_{k}}\unicode[STIX]{x1D6FC}_{j}^{k}(\mathbf{G}(\mathbf{x}_{k-m_{k}+j})-\mathbf{x}^{\ast }).\end{eqnarray}$$

Non-negativity of the $\unicode[STIX]{x1D6FC}^{k}$ then implies that

$$\begin{eqnarray}\displaystyle \hspace{35.39996pt}\Vert \mathbf{e}_{K+1}\Vert & \leq & \displaystyle \mathop{\sum }_{j=0}^{m_{K}}\unicode[STIX]{x1D6FC}_{j}^{K}\Vert \mathbf{G}(\mathbf{x}_{K-m_{K}+j})-\mathbf{x}^{\ast }\Vert \nonumber\\ \displaystyle & \leq & \displaystyle \mathop{\sum }_{j=0}^{m_{K}}\unicode[STIX]{x1D6FC}_{j}^{K}c\Vert \mathbf{x}_{K-m_{K}+j}-\mathbf{x}^{\ast }\Vert \nonumber\\ \displaystyle & \leq & \displaystyle c\mathop{\sum }_{j=0}^{m_{K}}\unicode[STIX]{x1D6FC}_{j}^{K}{\hat{c}}^{K-m_{K}+j}\Vert \mathbf{e}_{0}\Vert \leq {\hat{c}}^{K+1}(c{\hat{c}}^{-m-1})\leq {\hat{c}}^{K+1}.\hspace{35.39996pt}\square\nonumber\end{eqnarray}$$

The theorem says the iteration history of $m$ vectors will eventually be arbitrarily close to $\mathbf{x}^{\ast }$ . Hence restarting the iteration after sufficiently many EDIIS iterations will result in local convergence at the rate predicted by Theorem 5.5, which is better than (5.23). However, it is not clear how to decide when to restart. Theorem 5.7, the local convergence result from Chen and Kelley (Reference Chen and Kelley2017), says that one can simply continue with the EDIIS iteration and the local convergence estimate will hold.

The proof of Theorem 5.7 depends strongly on Theorem 5.6 to generate good initial data and the constrained optimization problem to guarantee that the $\ell ^{1}$ -norm of the coefficients is bounded.

6.1 Generalized derivatives and semi-smooth functions

We begin with a idealized version of Newton’s method for Lipschitz continuous functions. Let $\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{N}$ be open. Suppose $\mathbf{F}:\unicode[STIX]{x1D6FA}\rightarrow \mathbb{R}^{N}$ is locally Lipschitz continuous. Rademacher’s theorem says that $\mathbf{F}$ is Fréchet differentiable almost everywhere. The proof of this remarkable result can be found in Federer (Reference Federer1969, Theorem 3.1.6).

We let $D_{\mathbf{F}}$ denote the set of points where $\mathbf{F}$ is Fréchet differentiable. The generalized Jacobian (Clarke Reference Clarke1990) of $\mathbf{F}$ at $\mathbf{u}\in \mathbb{R}^{N}$ is the set

(6.1) $$\begin{eqnarray}\unicode[STIX]{x2202}\mathbf{F}(\mathbf{u})=\text{co}\Bigl\{\lim _{\mathbf{u}_{j}\rightarrow \mathbf{u};\mathbf{u}_{j}\in D_{\mathbf{F}}}\mathbf{F}^{\prime }(\mathbf{u}_{j})\Bigr\},\end{eqnarray}$$

where $\text{co}$ denotes the closed convex hull.

Consider the scalar function $f(x)=|x|$ . The function is differentiable except at $x=0$ . At $x=0$ , where $f$ is not differentiable, the generalized derivative is the interval $[-1,1]$ .

The chain rule is not a trivial matter for non-smooth problems, and there are several variations. As an example we will state Theorem 2.6.6 from Clarke (Reference Clarke1990) and one of its corollaries. We state the results using the formulation from Hintermüller (Reference Hintermüller2010).

In the theorem $\mathbf{F}\circ \mathbf{G}$ will denote the composition $\mathbf{F}\circ \mathbf{G}(\mathbf{x})\equiv \mathbf{F}(\mathbf{G}(\mathbf{x}))$ .

Theorem 6.1. Let $\mathbf{Q}=\mathbf{G}\circ \mathbf{F}$ , where $\mathbf{F}:\mathbb{R}^{N}\rightarrow \mathbb{R}^{N}$ is Lipschitz continuous in a neighbourhood of $\mathbf{x}$ and $\mathbf{G}:\mathbb{R}^{N}\rightarrow \mathbb{R}^{P}$ is Lipschitz continuous in a neighbourhood of $\mathbf{F}(\mathbf{x})$ . Then $\mathbf{Q}$ is Lipschitz continuous in a neighbourhood of $\mathbf{x}$ and, for all $\mathbf{v}\in \mathbb{R}^{N}$ ,

(6.2) $$\begin{eqnarray}\unicode[STIX]{x2202}\mathbf{Q}(\mathbf{x})\mathbf{v}\subset \text{co}\{\unicode[STIX]{x2202}\mathbf{G}(\mathbf{F}(\mathbf{x}))\unicode[STIX]{x2202}\mathbf{F}(\mathbf{x})\mathbf{v}\}.\end{eqnarray}$$

Moreover, if $\mathbf{G}$ is continuously differentiable, then the inclusion is an equality, that is, for all $\mathbf{v}\in \mathbb{R}^{N}$ ,

(6.3) $$\begin{eqnarray}\unicode[STIX]{x2202}\mathbf{Q}(\mathbf{x})\mathbf{v}=\mathbf{G}^{\prime }(\mathbf{F}(\mathbf{x}))\unicode[STIX]{x2202}\mathbf{F}(\mathbf{x})\mathbf{v}.\end{eqnarray}$$

If $\mathbf{G}$ is real-valued ( $P=1$ ), then the vector $\mathbf{v}$ can be omitted from (6.2) and (6.3).