Newton's Method for finding a root

Newton's method uses a clever insight to iteratively home in on the root of a function $f$. The central idea is to approximate $f$ by its tangent at some initial position $x_0$:

$$ y = f'(x_0) (x-x_0) + f(x_0) $$

The $x$-intercept of this line is then closer to the root than the starting position $x_0$. That is, we need to solve the linear relation

$$ f'(x_n)(x_1-x_0) + f(x_0) $$

for the updated position $x_1 = x_0 - f(x_0)/f'(x_0)$. Repeating this sequence

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

will yield a fixed point, which is the root of $f$ if one exists in the vicinity of $x_0$.



In [2]:

    
def newtons_method(f, df, x0, tol=1E-6):
    x_n = x0    
    while abs(f(x_n)) > tol:
        x_n = x_n - f(x_n)/df(x_n)
    return x_n

Minimizing a function

As the maximum and minimum of a function are defined by $f'(x) = 0$, we can use Newton's method to find extremal points by applying it to the first derivative. Let's try this with a simply function with known minimum:



In [10]:

    
# define a test function
def f(x):
    return (x-3)**2 - 9

def df(x):
    return 2*(x-3)

def df2(x):
    return 2.



In [14]:

    
root = newtons_method(f, df, x0=0.1)
print ("root {0}, f(root) = {1}".format(root, f(root)))









    



root -4.092847520435343e-14, f(root) = 2.4513724383723456e-13



In [16]:

    
minimum = newtons_method(df, df2, x0=0.1)
print ("minimum {0}, f'(minimum) = {1}".format(minimum, df(minimum)))









    



minimum 3.0, f'(minimum) = 0.0

There is an important qualifier in the statement about fixed points: a root needs to exist in the vicinity of $x_0$! Let's see what happens if that's not the case:



In [17]:

    
def f(x):
    return (x-3)**2 + 1
newtons_method(f, df, x0=0.1)









    



---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-17-cc1c41987083> in <module>()
      1 def f(x):
      2     return (x-3)**2 + 1
----> 3 newtons_method(f, df, x0=0.1)

<ipython-input-2-0945878377e7> in newtons_method(f, df, x0, tol)
      2     x_n = x0
      3     while abs(f(x_n)) > tol:
----> 4         x_n = x_n - f(x_n)/df(x_n)
      5     return x_n

<ipython-input-10-720123a4229c> in df(x)
      3     return (x-3)**2 - 9
      4 
----> 5 def df(x):
      6     return 2*(x-3)
      7 

KeyboardInterrupt:

With a little more defensive programming we can make sure that the function will terminate after a given number of iterations:



In [21]:

    
def newtons_method2(f, df, x0, tol=1E-6, maxiter=100000):
    x_n = x0    
    for _ in range(maxiter):
        x_n = x_n - f(x_n)/df(x_n)
        
        if abs(f(x_n)) < tol:
            return x_n
        
    raise RuntimeError("Failed to find a minimum within {} iterations ".format(maxiter))



In [22]:

    
newtons_method2(f, df, x0=0.1)









    



---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-22-8ff7ac5b1e25> in <module>()
----> 1 newtons_method2(f, df, x0=0.1)

<ipython-input-21-60eb0a1f76a6> in newtons_method2(f, df, x0, tol, maxiter)
      7             return x_n
      8 
----> 9     raise RuntimeError("Failed to find a minimum within {} iterations ".format(maxiter))

RuntimeError: Failed to find a minimum within 100000 iterations