Round off error in modified Newton-Raphson

Question

Newton's method for approximating roots of $f(x)=0$ is given by $$x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}, \quad n=0,1,2,\ldots,$$ where $x_n$ denotes the $n$th approximation of the root, and some sufficiently good initial approximation $x_0$ is required for convergence.

In general, this converges quadratically, but if the root, say $x^*$, satisfies $f'(x^*)=0$, the convergence is only linear. I'm aware that quadratic convergence can be restored by using the scheme $$x_{n+1}=x_n-\frac{F(x_n)}{F'(x_n)},$$ where $F(x_n)=f(x_n)/f'(x_n)$.

After substituting the definition of $F(x_n)$, it can be seen that $$x_{n+1}=x_n-\frac{f(x_n)f'(x_n)}{f'(x_n)f'(x_n)-f(x_n)f''(x_n)}.$$

I'm told that this raises another issue, namely that the denominator involves computing the difference of two nearly equal quantities and can lead to large rounding errors.

Questions

• Why are the quantities on the denominator nearly equal?
• How do we get around this issue?

Answer 1

1. The quantities in the denominator are not nearly equal.
2. Subtractive cancellation is not an issue here.

As @SimplyBeatifulArt has already pointed out we have $$ f'(x)f'(x) \approx 2f(x)f''(x)$$ near a double root of $f$. This eliminates the possibility of subtractive cancellation in the calculation of $$ f'(x)f'(x) - f(x) f''(x).$$ A general analysis of subtractive cancellation is given in this answer.

In general, the componentwise relative condition number of a subtraction $d = a - b \not = 0$ is given by $$\kappa = \frac{|a|+ |b|}{|a-b|},$$ while the componentwise relative condition number of a division $d = a/b$, $(b \not = 0)$ is given by $$\kappa = 2.$$ In particular, divisions are always well-conditioned in the relative sense.

Answer 2

Since both $f(x)\to0$ and $f'(x)\to0$, it is clear both terms on the denominator tend to $0$ quickly. Furthermore, if $f''(x)\to c\ne0$, then by Taylor's theorem one may use the estimations:

$$f(x)\sim\frac c2(x-x^\star)^2,\quad f'(x)\sim c(x-x^\star)$$

and the terms in the denominator are

$$f'(x)f'(x)\sim c^2(x-x^\star)^2,\quad f(x)f''(x)\sim\frac{c^2}2(x-x^\star)^2$$

and as you can see, they only differ by a factor of $2$ asymptotically. To try to avoid cancellation error, you can rewrite the expression as

$$x-\frac{f(x)}{f'(x)}\left(1-\frac{f(x)}{f'(x)}\frac{f''(x)}{f'(x)}\right)^{-1}$$

where one can see that

$$\frac{f(x)}{f'(x)}\frac{f''(x)}{f'(x)}\sim\frac12$$

$$\left(1-\frac{f(x)}{f'(x)}\frac{f''(x)}{f'(x)}\right)^{-1}\sim2$$

and that evaluations are more numerically stable since you get tricky divisions out of the way before other operations.

Some notable remarks:

• If it is known ahead of time that $f^{(n)}(x^\star)$ is $0$ for all $n<N$ and non-zero for $n=N$, then you can replace the term multiplying $f(x)/f'(x)$ by $N$. The proof is the same Taylor argument as above.

• Since it can be expected that the term multiplying $f(x)/f'(x)$ tends to a constant except in extreme examples (such as $\exp(-x^{-2})$, which has $f^{(n)}(0)=0$ for all $n$), it suffices to replace the term by the last numerically stable value you were able to compute, if you were to have issues evaluating it.

• If $f$ is known to be smooth, you can round the term multiplying $f(x)/f'(x)$ to the nearest integer.

• The above form is also advantageous if you can evaluate $f(x)/f'(x)$ and $f'(x)/f''(x)$ faster, or as fast as, or more numerically stable than $f(x),f'(x),f''(x)$ individually.