Which absolute and/or relative stopping criteria do use for Newton's method?

Question

I saw many stopping criteria for Newton's method all around Web and books.

Some are defined from the residuals:

• of either current iteration only:

$$ \|f(\mathbf{x}^{(k)})\| \leq \epsilon $$

(https://fr.wikipedia.org/wiki/M%C3%A9thode_de_Newton#Crit%C3%A8re_d'arr%C3%AAt)

• or of current and previous iterations, mixing both absolute and relative tolerances:

$$ \|f(\mathbf{x}^{(k)})\| \leq \delta \|f(\mathbf{x}^{(k-1)})\| + \epsilon $$

(https://archive.siam.org/books/textbooks/fr16_book.pdf, p.73).

Some other criteria are defined from the gradients:

• of absolute differences:

$$ \|\mathbf{x}^{(k)}- \mathbf{x}^{(k-1)}\| \leq \epsilon $$

(https://fr.wikipedia.org/wiki/M%C3%A9thode_de_Newton#Crit%C3%A8re_d'arr%C3%AAt ; https://www.hds.utc.fr/~tdenoeux/dokuwiki/_media/en/univariate_optimization.pdf, p.20),

• of relative differences:

$$ \|\mathbf{x}^{(k)} - \mathbf{x}^{(k-1)}\| \leq \delta \|\mathbf{x}^{(k)}\| $$

(https://www.hds.utc.fr/~tdenoeux/dokuwiki/_media/en/univariate_optimization.pdf, p.20 ; https://scicomp.stackexchange.com/a/30606/38451),

• or of both absolute and relative differences, using absolute and relative tolerances:

$$ \|\mathbf{x}^{(k)} - \mathbf{x}^{(k-1)}\| \leq \delta \|\mathbf{x}^{(k)}\| + \epsilon $$

(https://archive.siam.org/books/textbooks/fr18_book.pdf, p.16).

For the vector norm, some authors use the infinity norm $\|.\|_\infty$, while some others use the L-1 norm $\|.\|_1$.

Finally, I saw that it's also possible to compute the norm of the relative errors rather than the relative error of the norms, e.g. (using $\oslash$ for element-wise vector division):

$$ \| \left( \mathbf{x}^{(k)} - \mathbf{x}^{(k-1)} \right) \oslash \mathbf{x}^{(k)} \| \leq \delta $$

rather than:

$$ \|\mathbf{x}^{(k)} - \mathbf{x}^{(k-1)}\| \leq \delta \|\mathbf{x}^{(k)}\| $$

(https://www.rocq.inria.fr/modulef/Doc/GB/Guide6-10/node21.html),

or, to avoid division by zero:

$$ \| \left( \mathbf{x}^{(k)} - \mathbf{x}^{(k-1)} \right) \oslash \left( \left| \mathbf{x}^{(k)} \right| + 1 \right) \| \leq \delta $$

So my question is: which one of these criteria should I use? Intuitively I would say that checking both residuals and gradients, using both absolute and relative tolerances, and using infinity norm would be the best. But then, why "simpler" criteria are still so often used for many applications? If there is no universal termination criteria applicable to all situations, how could I choose the one which is the best for me?

Answer 1

This is not a definite answer, but I highly recommend the work of Deuflhard on Newton methods (e.g. this report, this book, and another report I cannot find just yet).

The basic idea is the following: a criterion based on the norm of the Newton step is scaling independent, as the sequence of Newton iterates is actually independent of the scalings used for your variables (a property called affine invariance). If the size of the Newton step relative to your solution vector becomes lower than a low tolerance (say $10^{-9}$), this means that the next Newton iterates will be very close to the current one, hence you can consider that the solution is sufficiently converged. On the opposite, a stopping criterion based on the norms of the residual may be more erroneous depending on the scaling of the residuals. In my codes, I usually rely on a stopping criterion based on the relative norm of the Newton step, and just for safety I also use a criterion on the residual norm, but the latter is never triggered in practice.

A practical example is the case of adaptive implicit time integration schemes for ODEs (BDF, Radau...). In these algorithms, the stopping criterion for the Newton loop is generally based on the relative size of the Newton step, with some tweaks such as estimating the convergence rate to see if the Newton iterates are indeed likely to converge. If the user wants an integration error lower than $10^{-n}$, then the Newton algorithm should itself converge with a tighter tolerance so as to not pollute the error estimate (say $<0.3 \times 10^{-n}$). You can find some info on that in this book (in the case of Radau). For me personally, that is the reason why I use a criterion on the relative norm of the Newton step. Otherwise, it would not always be easy to determine the "optimal" tolerance on the residual norm (low enough so as to ensure the error estimate is not polluted, but not too low so that convergence is possible in a reasonable number of iterations).

In your case, I guess a good solution is to monitor the different stopping criteria to see which one seems most relevant and reliable, and then only use this one.

Answer 2

Some of the confusion comes about because mathematicians often consider artificial functions that happen to have size $O(1)$ whereas practitioners use functions and variables that just happen to be whatever they are, often vastly different from one.

One of the considerations that is important in practice is to make sure that the choice of physical units does not affect when to stop an iteration. For example, if you are optimizing a situation where $x$ denotes lengths, then the condition $$ \|x_k-x_{k-1}\| \le 10^{-9} $$ does not make any sense because you will only stop if distances change from one iterate to another by no more than one nanometer if you happen to measure distances in meters, but the iteration will already stop if distances change by less than 10,000 km if you happen to measure distances in lightyears.

One can make that same kind of argument to argue that a criterion such as $$ \|F(x_k)-F(x_{k-1})\| \le 10^{-9} $$ does not work because when the iteration stops depends on what units you choose to measure $F$.

As a consequence of these sorts of considerations, it makes sense to conclude that the tolerance used on the right hand side has to somehow take into account what units we are using, and in particular that when the iteration is stopped should not depend on the choice of units. Some of the criteria you list achieve this by putting something proportional to $x$ or $F(x)$ on the right hand side.

In the end, the "right" criteria are those where you would stop if $$ \|x_k-x_{k-1}\| \le \varepsilon L $$ where $L$ is a typical length scale (or, more generally: a typical scale for whatever it is that $x$ measures), or if $$ \|F(x_k)-F(x_{k-1})\| \le \varepsilon \Phi $$ where $\Phi$ is a typical scale for whatever it is that $F$ measures.

What $\Phi$ or $L$ are depends on the problem, and should be chosen based on practitioner knowledge of the problem. For example, if you are trying to optimize the shape of a car, then maybe $L=5m$ is a good scale if you are considering the overall shape of the car, but if you are optimizing the shape of a side mirror, then $L=20cm$ would be the right choice.

In actual practice, people don't want to go through these sorts of considerations, and so software creators have come up with ways of determining values for $L$ or $\Phi$ automatically. The way you would do that is ask "Over what scales of $x$ does $F(x)$ change significantly?" in some sense. The methods due to Deuflhard mentioned in one of the other comments are just mathematical representations of this sort of idea.

Answer 3

To add to Wolfgang's answer, some problems have specific structure that you can use to design physically-based convergence criteria. Suppose you're solving the generalized Poisson equation or some other elliptic partial differential equation that can be derived through minimization of some convex action functional $J(u)$. Very often, the action can be split into two terms, one of which is strictly positive. For example, the action for the generalized Poisson equation is

$$J(u) = \int_\Omega\left(\frac{1}{2}k|\nabla u|^2 - fu\right)\mathrm dx.$$

This is a little silly because it's quadratic, so you could also consider the minimal surface problem

$$J(u) = \int_\Omega\left(\sqrt{1 + \ell^2|\nabla u|^2} - fu\right)\mathrm dx.$$

The first term in both objective functionals is strictly positive; you can call this a scale functional $S(u)$, although admittedly I made up this terminology.

Now suppose you've computed a search direction $v$ by solving the linear system

$$\mathrm d^2J(u)\cdot v = -\mathrm dJ(u)$$

and you want to know whether to keep going. The scalar quantity

$$\Delta \equiv \frac{1}{2}\langle \mathrm d^2J(u)\cdot v, v\rangle$$

is (because of convexity) strictly positive and tells you how much you'd expect the objective to decrease if you took one more iteration of Newton's method, assuming that the objective is approximated well by its second-order Taylor expansion at $u$. But without (as Wolfgang points out) some idea of the problem scale, you don't know a priori whether that expected decrease is large or small. You could try to use the objective functional itself, but because of the forcing term $J(u)$ could be positive or negative or zero depending on the boundary conditions. But the first term of the objective -- the scale functional $S(u)$ -- gives you exactly the positive scale you need. You can then define a stopping criterion by the size of the Newton decrement $\Delta$ relative to the scale functional:

$$\Delta < \epsilon\cdot S(u)$$

Crucially, the scale is diagnosed from the current solution guess and doesn't have to be picked a priori.

The key insight is that the objective is separable into $J(u) = S(u) + $ forcing. These kinds of problem-specific tricks can be easy to miss when there's a disconnect between the people developing optimization algorithms or nonlinear solvers and the domain specialists who are using them.