12) Newton methods

12) Newton methods#

Last time#

Convergence classes
Intro to Newton methods

Today#

Newton-Raphson Method
1.1 An implementation
Convergence of fixed-point (by mean value theorem)
Convergence of fixed-point (by Taylor series)
A fresh derivation of Newton’s method via convergence theory

using Plots
default(linewidth=4, legendfontsize=12)

f(x) = cos(x) - x
hasroot(f, a, b) = f(a) * f(b) < 0
function bisect_iter(f, a, b, tol)
    hist = Float64[]
    while abs(b - a) > tol
        mid = (a + b) / 2
        push!(hist, mid)
        if hasroot(f, a, mid)
            b = mid
        else
            a = mid
        end
    end
    hist
end

Precompiling packages...

    Plots Being precompiled by another process (pid: 21767, pidfile: /home/vbarra/.julia/compiled/v1.10/Plots/ld3vC_H1Fjh.ji.pidfile)

  23314.2 ms  ✓ Plots
    UnitfulExt Being precompiled by another process (pid: 21767, pidfile: /home/vbarra/.julia/compiled/v1.10/UnitfulExt/4ishe_H1Fjh.ji.pidfile)

Convergence classes#

Tip

Check this Reading

A convergent rootfinding algorithm produces a sequence of approximations $x_{k}$ such that

lim_{k \to \infty} x_{k} \to x_{*}

where

f (x_{*}) = 0

. For analysis, it is convenient to define the errors

e_{k} = x_{k} - x_{*}

. We say that an iterative algorithm is $q$ -linearly convergent if

lim_{k \to \infty} | e_{k + 1} | / | e_{k} | = ρ < 1.

(The

q

is for “quotient”.) A smaller convergence factor

ρ

represents faster convergence. A slightly weaker condition (

r

-linear convergence or just linear convergence - the “r” here is for “root”) is that

| e_{k} | \leq ϵ_{k}

for all sufficiently large

k

when the sequence

{ϵ_{k}}

converges

q

-linearly to 0.

We saw that:

Note

Why is $r$ -convergence considered a weaker form of convergence relative to $q$ -convergence? Notice that root convergence is concerned only with the overall rate of decrease of the error while quotient convergence requires the error to decrease at each iteration of the algorithm. Thus $q$ -convergence is a stronger form of convergence than $r$ -convergence
$q$ -convergence implies ( $⟹$ ) $r$ -convergence.

hist = bisect_iter(f, -1, 3, 1e-10)
r = hist[end] # What are we trusting?
hist = hist[1:end-1]
scatter( abs.(hist .- r), yscale=:log10)
ks = 1:length(hist)
ρ = 0.5
plot!(ks, 4 * (ρ .^ ks))

1. Newton-Raphson Method#

Much of numerical analysis reduces to Taylor series, the approximation

f (x) = f (x_{0}) + f^{'} (x_{0}) (x - x_{0}) + f^{″} (x_{0}) (x - x_{0})^{2} / 2 + \underset{O ((x - x_{0})^{3})}{\underset{⏟}{\dots}}

centered on some reference point

x_{0}

.

In numerical computation, it is exceedingly rare to look beyond the first-order approximation

{\tilde{f}}_{x_{0}} (x) = f (x_{0}) + f^{'} (x_{0}) (x - x_{0}) .

Since

{\tilde{f}}_{x_{0}} (x)

is a linear function, we can explicitly compute the unique solution of

{\tilde{f}}_{x_{0}} (x) = 0

as

x = x_{0} - \frac{f (x_{0})}{f^{'} (x_{0})} .

This is Newton’s Method (aka Newton-Raphson or Newton-Raphson-Simpson) for finding the roots of differentiable functions.

1.1 An implementation#

function newton(f, fp, x0; tol=1e-8, verbose=false)
    x = x0
    for k in 1:100 # max number of iterations
        fx = f(x)
        fpx = fp(x)
        if verbose
            println("iteration $k: x=$x,  f(x)=$fx,  f'(x)=$fpx")
        end
        if abs(fx) < tol
            return x, fx, k
        end
        x = x - fx / fpx
    end
end

f(x) = cos(x) - x
fp(x) = -sin(x) - 1
newton(f, fp, 1; tol=1e-15, verbose=true)

That’s really fast!#

10 digits of accuracy in 4 iterations.
How is this convergence test different from the one we used for bisection?
How can this break down?

x_{k + 1} = x_{k} - \frac{f (x_{k})}{f^{'} (x_{k})}

newton(f, fp, -pi/2+.1; verbose=true)

2. Convergence of fixed-point (by mean value theorem)#

Consider the iteration

x_{k + 1} = g (x_{k})

where

g

is a continuously differentiable function. Suppose that there exists a fixed point

x_{*} = g (x_{*})

. By the mean value theorem, we have that

x_{k + 1} - x_{*} = g (x_{k}) - g (x_{*}) = g^{'} (c_{k}) (x_{k} - x_{*})

for some

c_{i}

between

x_{k}

and

x_{*}

.

Taking absolute values,

| e_{k + 1} | = | g^{'} (c_{k}) | | e_{k} |,

which converges to zero if

| g^{'} (c_{k}) | < 1

.

3. Convergence of fixed-point (by Taylor series)#

Consider the iteration

x_{k + 1} = g (x_{k})

where

g

is a continuously differentiable function. Suppose that there exists a fixed point

x_{*} = g (x_{*})

. There exists a Taylor series at

x_{*}

,

g (x_{k}) = g (x_{*}) + g^{'} (x_{*}) (x_{k} - x_{*}) + O ((x_{k} - x_{*})^{2})

and thus

(1)#

\begin{aligned} x_{k + 1} - x_{*} & = g (x_{k}) - g (x_{*}) \\ = g^{'} (x_{*}) (x_{k} - x_{*}) + O ((x_{k} - x_{*})^{2}) . \end{aligned}

In terms of the error $e_{k} = x_{k} - x_{*}$ ,

| \frac{e_{k + 1}}{e_{k}} | = | g^{'} (x_{*}) | + O (e_{k}) .

Poll 12.1: Convergence class of the fixed-point method#

Is this convergence:

A=q-linear
B=r-linear
C=neither

Recall the definition of q-linear convergence

lim_{k \to \infty} | \frac{e_{k + 1}}{e_{k}} | = ρ < 1.

So, clearly, we need

| g^{'} (x_{*}) | < 1

to have convergence of the fixed-point iteration method. Otherwise, it diverges.

Aside:#

Big $O$ (“big oh”) notation

Limit $n \to \infty$ :#

We’d say an algorithm costs $O (n^{2})$ if its running time on input of size $n$ is less than $c n^{2}$ for some constant $c$ and sufficiently large $n$ .

Sometimes we write $cost (algorithm, n) = O (n^{2})$ or (preferably) $cost (algorithm) \in O (n^{2})$ .

Note that $O (\log n) \subset O (n) \subset O (n \log n) \subset O (n^{2}) \subset \dots$ .

We say the algorithm is in $Θ (n^{2})$ (“big theta”) if

c_{1} n^{2} < cost (algorithm) < c_{2} n^{2}

for some positive constants

c_{1}, c_{2}

and sufficiently large

n

.

Limit $h \to 0$ :#

In numerical analysis, we often have a small real number, and now the definitions take the limit as the small number goes to zero. So we say a term in an expression is in $O (h^{2})$ if

lim_{h \to 0} \frac{term (h)}{h^{2}} < \infty .

Big $O$ terms can be manipulated as

(2)#

\begin{aligned} h O (h^{k}) & = O (h^{k + 1}) \\ O (h^{k}) / h & = O (h^{k - 1}) \\ c O (h^{k}) & = O (h^{k}) \\ O (h^{k}) - O (h^{k}) & = ? \end{aligned}

Example of a fixed-point iteration#

We wanted to solve $\cos x - x = 0$ , which occurs when $g (x) = \cos x$ is a fixed point.

xstar, _ = newton(f, fp, 1.)
g(x) = cos(x)
gprime(x) = -sin(x)
@show xstar
@show gprime(xstar)
plot([x->x, g], xlims=(-2, 3), label = ["y=x" "y=g(x)"])
scatter!([xstar], [xstar],
    label="\$x_*\$")

function fixed_point(g, x, n)
    xs = [x]
    for k in 1:n
        x = g(x)
        append!(xs, x)
    end
    xs
end

xs = fixed_point(g, 2., 15)
plot!(xs, g.(xs), seriestype=:path, marker=:auto)

Verifying fixed-point method convergence theory#

| \frac{e_{k + 1}}{e_{k}} | \to | g^{'} (x_{*}) |

@show gprime(xstar)
es = xs .- xstar
es[2:end] ./ es[1:end-1]

scatter(abs.(es), yscale=:log10, label="fixed point")
plot!(k -> abs(gprime(xstar))^k, label="\$|g'|^k\$")

Reading: FNC Fixed-point iteration #

Plotting Newton convergence#

function newton_hist(f, fp, x0; tol=1e-12)
    x = x0
    hist = []
    for k in 1:100 # max number of iterations
        fx = f(x)
        fpx = fp(x)
        push!(hist, [x fx fpx])
        if abs(fx) < tol
            return vcat(hist...) # concatenate arrays
        end
        x = x - fx / fpx
    end
end

xs = newton_hist(f, fp, 1.97)
@show x_star = xs[end,1]
@show xs[1:end-1,1] .- x_star
plot(xs[1:end-1,1] .- x_star, yscale=:log10, marker=:auto, xlabel = "# of iterations", ylabel = "abs error")

Poll 12.2: Convergence of Newton’s method#

Is this convergence:

A=q-linear
B=r-linear
C=neither?

Formulations are not unique (constants):#

If $x = g (x)$ then

x = \underset{g_{2}}{\underset{⏟}{x + c (g (x) - x)}}

for any constant

c \neq 0

. Can we choose

c

to make

| g_{2}^{'} (x_{*}) |

small?

c = .5
g2(x) = x + c * (cos(x) - x)
g2p(x) = 1 + c * (-sin(x) - 1)
@show g2p(xstar)
plot([x->x, g, g2], ylims=(-5, 5), label=["y=x" "y=g(x)" "y=g2(x)"])

xs = fixed_point(g2, 1., 15)
xs .- xstar

Formulations are not unique (functions)#

If $x = g (x)$ then

x = \underset{g_{3} (x)}{\underset{⏟}{x + h (x) (g (x) - x)}}

for any smooth

h (x) \neq 0

. Can we choose

h (x)

to make

| g_{3}^{'} (x) |

small?

h(x) = -1 / (gprime(x) - 1)
g3(x) = x + h(x) * (g(x) - x)
plot([x-> x, cos, g2, g3], ylims=(-5, 5), label=["y=x" "y=g(x)" "y=g2(x)" "y=g3(x)"])

Exercise 12.1:#

Compute $| g_{3}^{'} (x_{*}) |$ .

Note:

We don’t know $g^{'} (x_{*})$ in advance because we don’t know $x_{*}$ yet.
This method converges very fast
We actually just derived Newton’s method.

4. A fresh derivation of Newton’s method via convergence theory#

A rootfinding problem $f (x) = 0$ can be converted to a fixed point problem
$x = x + f (x) =: g (x)$
but there is no guarantee that $g^{'} (x_{*}) = 1 + f^{'} (x_{*})$ will have magnitude less than 1.
Problem-specific algebraic manipulation can be used to make $| g^{'} (x_{*}) |$ small.
$x = x + h (x) f (x)$ is also a valid formulation for any $h (x)$ bounded away from $0$ .
Can we choose $h (x)$ such that
$g^{'} (x) = 1 + h^{'} (x) f (x) + h (x) f^{'} (x) = 0$
when $f (x) = 0$ ?

In other words,

x_{k + 1} = x_{k} + \underset{h (x_{k})}{\underset{⏟}{\frac{- 1}{f^{'} (x_{k})}}} f (x_{k}) .

Quadratic convergence!#

| \frac{e_{k + 1}}{e_{k}} | \to | g^{'} (x_{*}) |

What does it mean that $g^{'} (x_{*}) = 0$ ?
It turns out that Newton’s method has locally quadratic convergence to simple roots,
$lim_{k \to \infty} \frac{| e_{k + 1} |}{| e_{k} |^{2}} < \infty .$
“The number of correct digits doubles each iteration.”
Now that we know how to make a good guess accurate, the effort lies in getting a good starting guess.

Rootfinding summary and outlook#

Newton methods are immensely successful
- Convergence theory is local; we need good initial guesses
- Computing the derivative $f^{'} (x)$ is intrusive (and can be error-prone)
  - Avoided by secant methods (approximate the derivative)
  - Algorithmic or numerical differentiation can also be used instead of computing the derivative by hand
- Bisection is robust when conditions are met
- When does Newton diverge?
More topics
- Do we find all the roots?
- Times when Newton converges slowly

Exercise 12.2:#

Solve the fixed point iteration $g (x_{n}) = (2 + {(x_{n} - 2)}^{2})$ , up to a tolerance $10^{- 6}$ or maximum number of iterations $N = 50$ , with (a) $x_{0} = 2.8$ , (b) $x_{0} = 3.1$ . Which execution is better? And why?
Solve the fixed point iteration $g (x_{n}) = (2 + {(x_{n} - 2)}^{3})$ , up to a tolerance $10^{- 6}$ or maximum number of iterations $N = 50$ , with (a) $x_{0} = 2.8$ , (b) $x_{0} = 3.1$ . Which execution is better? And why?

12) Newton methods

Contents

12) Newton methods#

Last time#

Today#

Convergence classes#

1. Newton-Raphson Method#

1.1 An implementation#

That’s really fast!#

2. Convergence of fixed-point (by mean value theorem)#

3. Convergence of fixed-point (by Taylor series)#

Poll 12.1: Convergence class of the fixed-point method#

Aside:#

Limit n→∞:#

Limit h→0:#

Example of a fixed-point iteration#

Verifying fixed-point method convergence theory#

Reading: FNC Fixed-point iteration#

Plotting Newton convergence#

Poll 12.2: Convergence of Newton’s method#

Formulations are not unique (constants):#

Formulations are not unique (functions)#

Exercise 12.1:#

4. A fresh derivation of Newton’s method via convergence theory#

Quadratic convergence!#

Rootfinding summary and outlook#

Exercise 12.2:#

Limit $n \to \infty$ :#

Limit $h \to 0$ :#

Reading: FNC Fixed-point iteration #