Bossy Lobster

A blog by Danny Hermes; musing on tech, mathematics, etc.

Newton's (Method's) Failure

Finding zeros of any old function is a common task, and using Newton's method is one of the best tools for carrying out this task. I've even written an old post that used this method.

However, Newton's method loses some of it's luster around repeated roots. Consider the iteration

xn+1=xnf(xn)f(xn)x_{n + 1} = x_n - \frac{f(x_n)}{f'(x_n)}

On inspection, points where f(xn)=0f'(x_n) = 0 seem to be problematic. However, when these points are also roots, the division 0/00 / 0 may no longer be a problem. The case where

f(xn)=f(xn)=0f(x_n) = f'(x_n) = 0

is exactly the case of a repeated root.

For a practical example, consider f(x)=(x1)2f(x) = (x - 1)^2 which has corresponding iteration function

g(x)=x(x1)22(x1)=x+12.g(x) = x - \frac{(x - 1)^2}{2(x - 1)} = \frac{x + 1}{2}.

Starting with x0=1+20=2x_0 = 1 + 2^{-0} = 2, in exact arithmetic (i.e. no rounding) we'll always have xn=1+2nx_n = 1 + 2^{-n} since

1+xn2=2+2n2=1+2n1.\frac{1 + x_n}{2} = \frac{2 + 2^{-n}}{2} = 1 + 2^{-n-1}.

This sequence converges to 11 and the error halves every term. However, in double precision, the sequence stops once n=53n = 53: at that point 1253=01 \oplus 2^{-53} = 0.

However, if we told the computer that f(x)=(x1)2f(x) = (x - 1)^2, there would be no reason to use Newton's method, the roots are clearly 11 and 11. Instead, the "data" the computer is given are the coefficients in x22x+1x^2 - 2x + 1 and each term is computed in floating point

xn+1=xn(NnDn).x_{n + 1} = x_n \ominus \left(N_n \oslash D_n\right).

where the numerator and denominator are given (with the help of Horner's method) by

Nn=((xn2)xn)1,Dn=(2xn)2.N_n = \left((x_n \ominus 2) \otimes x_n\right) \oplus 1, \quad D_n = \left(2 \otimes x_n\right) \ominus 2.

Each floating point operation (,,,\otimes, \oslash, \oplus, \ominus) rounds the "exact" result to the nearest floating point number. For x1,,x27x_1, \, \ldots, \, x_{27}, there is no rounding needed and we get the same values that we got in exact arithmetic.

The real problem happens when n=27n = 27. At that point

(xn2)xn=((1227))(1+227)(x_n \ominus 2) \otimes x_n = \left(-\left(1 - 2^{-27}\right)\right) \otimes \left(1 + 2^{-27}\right)

and the exact result 12541 - 2^{-54} needs to be rounded into floating point:

N27=((1254))1=(1)1=0N_{27} = \left(-(1 \ominus 2^{-54})\right) \oplus 1 = \left(-1\right) \oplus 1 = 0

which forces x28=x270x_{28} = x_{27} \ominus 0. After this point, the iteration stalls and every term after is equal to x27x_{27}.

Thus when the algorithm terminates on a computer (since two consecutive terms are identical), the actual error is 2272^{-27}. This error is much too large! Using the standard floating point type (double), the machine precision is ε=252\varepsilon = 2^{-52} and our error is approximately ε\sqrt{\varepsilon}. In other words, we expected 52 bits of precision and only get 27!


In general, a root repeated mm times (i.e. with multiplicity mm) will only be accurate to 52m\frac{52}{m} bits, a very large penalty!

We can see this phenomenon in action, we define the following function to do one Newton step using Horner's method:

def horner_newton(coeffs, val):
    degree = len(coeffs) - 1
    deriv_factor = float(degree)
    fx = fpx = 0.0
    for i in range(degree):
        fx = val * fx + coeffs[i]
        fpx = val * fpx + deriv_factor * coeffs[i]
        deriv_factor -= 1.0
    # Add the last term to fx.
    fx = val * fx + coeffs[-1]
    return val - fx / fpx

Applying it to (x22)2=x44x2+4(x^2 - 2)^2 = x^4 - 4x^2 + 4 — which has two double roots — we get

>>> values = [1.0]
>>> coeffs = [1, 0, -4, 0, 4]
>>> next_val = horner_newton(coeffs, values[-1])
>>> while next_val not in values:
...     values.append(next_val)
...     next_val = horner_newton(coeffs, values[-1])
...
>>> len(values)
28
>>> values[-1]
1.4142135634821
>>> math.log(abs(values[-1] - math.sqrt(2)), 2)
-29.74808712981571

i.e. the error is around 229ε2^{-29} \approx \sqrt{\varepsilon}.

Trying out with a triple root (x1)3=x33x2+3x1(x - 1)^3 = x^3 - 3 x^2 + 3x - 1 we get

>>> values = [2.0]
>>> coeffs = [1, -3, 3, -1]
>>> next_val = horner_newton(coeffs, values[-1])
>>> while next_val not in values:
...     values.append(next_val)
...     next_val = horner_newton(coeffs, values[-1])
...
>>> len(values)
31
>>> values[-1]
1.0000046685597561
>>> math.log(abs(values[-1] - 1), 2)
-17.70859102007032

with an error around 217ε32^{-17} \approx \sqrt3 as predicted.

Comments