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 $$


  • 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)}\| $$


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?

  • $\begingroup$ Have you considered enforcing multiple criteria and stopping only when all of them are met? $\endgroup$ – Brian Borchers Mar 2 at 17:33
  • $\begingroup$ @BrianBorchers Yes. What I currently do is both ‖𝐱(π‘˜)βˆ’π±(π‘˜βˆ’1)‖≀𝛿‖𝐱(π‘˜)β€–+πœ– and ‖𝑓(𝐱(π‘˜))β€–β‰€πœ–, following advises from my (dear) supervisors. But I'm still wondering about the validity of the other test β€–(𝐱(π‘˜)βˆ’π±(π‘˜βˆ’1))⊘(∣∣𝐱(π‘˜)∣∣+1)‖≀𝛿, because it seems to me that it would better consider elements of the residues or of the gradients which are different by several orders of magnitude. Till now both approaches give me good convergence rates, but since the last one doesn't seem to be much used by others, I'm wondering if it's really robust. $\endgroup$ – Camille C Mar 2 at 19:30
  • 2
    $\begingroup$ You've brought up the issue of scaling which is also very important. You'll generally want to have the elements of the $x$ vector scaled so that they're of comparable magnitude and you'll want the elements of the vector $f(x)$ to be of similar magnitude. If there are very large changes in $x$ from one iteration to the next, but $f(x)$ isn't changing much, then this is an indication of very badly scaled or conditioned problem. $\endgroup$ – Brian Borchers Mar 2 at 20:13
  • 1
    $\begingroup$ You could use use individual relative errors when the system consists of clear components that might have different scales (for a solar system simulation in km-s units one would probably want to treat the planets separately, and then also positions and velocities). $\endgroup$ – Lutz Lehmann Mar 3 at 19:57
  • 1
    $\begingroup$ As a practical matter you might want to play it safe and include a criterion on relative change in $x$ in case your problem isn't as well scaled as you thought. $\endgroup$ – Brian Borchers Mar 3 at 21:56

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.

  • $\begingroup$ Thanks. I tested the different criteria. Absolute residuals-based combined or not with absolute or relative Newton-step-based criteria lead to convergence. But relative residuals-based criteria either never get satisfied or lead to bad results. I obtain the best convergence rates using absolute Newton-step-based criteria, norm of the relative differences between Newton steps, or the intersection of both, i.e. ‖𝐱(π‘˜)βˆ’π±(π‘˜βˆ’1)β€–β‰€πœ– or β€–(𝐱(π‘˜)βˆ’π±(π‘˜βˆ’1))⊘(∣∣𝐱(π‘˜)∣∣+1)‖≀𝛿 or ‖𝐱(π‘˜)βˆ’π±(π‘˜βˆ’1)β€–β‰€πœ– ∩ β€–(𝐱(π‘˜)βˆ’π±(π‘˜βˆ’1))⊘(∣∣𝐱(π‘˜)∣∣+1)‖≀𝛿. $\endgroup$ – Camille C Mar 3 at 11:55
  • $\begingroup$ I also compared iteration numbers, and found finally that the most effective criterion for my problem is ‖𝐱(π‘˜)βˆ’π±(π‘˜βˆ’1)‖≀𝛿‖𝐱(π‘˜)β€–+πœ–. So I decided, as you suggested, to use this criterion and to also use criterion ‖𝑓(𝐱(π‘˜))β€–β‰€πœ– for safety (β€–.β€– = infinity-norm). $\endgroup$ – Camille C Mar 3 at 13:36
  • 1
    $\begingroup$ @camille-c For a propre criterion on the size of the Newton step, it might be better to use something similar to what is done for the relative error in embedded Runge-Kutta methods: $||\frac{\delta \vec{x}}{rtol \vec{x} + atol}||<1$, where you can set $atol=rtol=\epsilon$, and the division is done element by element. This is quite a robust way to handle different orders of magnitude in your variables. $\endgroup$ – Laurent90 Mar 3 at 18:03
  • $\begingroup$ Thanks. As I answered to Brian Borchers, using a criterion on relative change in π‘₯ for my problem currently increases the number of iterations needed ; it leads to smaller residues, but I don't need such a level of accuracy by now. So I will keep your criterion in mind for the case I may encounter convergence difficulties, either because of extreme parameters and/or if I add more complex processes or interdependences. $\endgroup$ – Camille C Mar 4 at 9:41
  • $\begingroup$ @CamilleC Yes, that depends indeed on your objective. The above criterion may be useful if you want to compare the solutions (obtained via a Newton method) of two optimisations with very similar inputs. If each Newton is not converged enough, the difference between the two solutions may be polluted by the poor convergence. I don't know if that applies to your case. $\endgroup$ – Laurent90 Mar 4 at 10:39

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.

  • $\begingroup$ Thanks. Indeed, I already saw your previous good answer to this post, and it actually helped me find the πœ€ you're talking about. Since the units in my problem are fixed, and all the elements of my solution-vector have the same magnitude order (even if they may represent different variables), then it's easy to choose my πœ€. $\endgroup$ – Camille C Mar 3 at 16:24

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.

  • $\begingroup$ Thanks. I'm not a specialist of generalized Poisson equation or minimal surface problem, but what is the difference between your Ξ”<πœ–β‹…π‘†(𝑒) and ‖𝐱(π‘˜)βˆ’π±(π‘˜βˆ’1)‖≀𝛿‖𝐱(π‘˜)β€–? $\endgroup$ – Camille C May 5 at 22:00
  • $\begingroup$ There are a few differences. Relative change in the solution norm requires you to at least compute one more iterate, but the ratio of expected decrease to problem scale can be evaluated purely from the initial guess. So for example if you call a solve routine and then call solve again using the answer from the first call as the initial guess, a relative decrease criterion will do more work whereas an expected decrease / scale criterion will not. $\endgroup$ – Daniel Shapero May 6 at 1:10
  • $\begingroup$ I should add that this isn't the be-all, end-all approach to stopping criteria and it's a good idea to combine a few different criteria. Just that if you know something special about the problem you're solving then you can very often use that to do better than a completely general approach. $\endgroup$ – Daniel Shapero May 6 at 1:12
  • $\begingroup$ So, if I well understood, the use of an expected decrease / scale criterion may prevent from computing the last "useless" iteration, i.e. the iteration computed just to test for the convergence. Am I right? $\endgroup$ – Camille C May 6 at 8:24
  • $\begingroup$ Yes that's correct! I like to think of this as a kind of idempotence. Again, your mileage may vary with this, just food for thought! $\endgroup$ – Daniel Shapero May 6 at 15:28

Your Answer

By clicking β€œPost Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.