0

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?
MHW
  • 2,940
  • The idea to investigate a simple model equation and to look what exactly is happening there, both analytically and numerically, didn't occur to you? –  Dec 24 '20 at 21:37
  • @ProfessorVector It didn't occur to me, but I will try it now. – MHW Dec 24 '20 at 21:56

2 Answers2

1

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.

Simply Beautiful Art
  • 74,685
  • Note that a subtraction $d=a-b \not =0$ does not really experience cancellation when $|a| \ge 2|b|$ or $|b| \ge 2|a|$. Hence, OP is actually safe and OP's instructor made a mistake. – Carl Christian Dec 29 '20 at 23:36
  • Ah yes you are correct. I suppose the main issue then would stem from other factors, such as noise, which may cause the denominator to experience cancellation. In any case, I think most of my answer is still pretty useful. – Simply Beautiful Art Dec 30 '20 at 02:29
  • I do not understand why cancellation is an issue here at all. As long as we can compute $x_{n+1}$ by adding a small correction to $x_n$, cancellation is a very minor problem. In fact, we can afford to lose 8 of 16 digits in each correction and still experience quadratic convergence. Why is that? Suppose $x_0$ has $1$ correct digit. Then we need first 1, then 2, then 4, and finally 8 correct digits in the corrections to experience quadratic convergence. This gives us a total of 16 correct digits. – Carl Christian Dec 30 '20 at 20:40
  • I already agreed with you on that point in my last comment. As I stated, the only issues would come from other factors e.g. the terms on the denominator cannot be evaluated accurately on their own, and since they are in fact close in the absolute sense, this could raise issues. – Simply Beautiful Art Dec 30 '20 at 23:10
0
  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.

Carl Christian
  • 12,583
  • 1
  • 14
  • 37