#%config InlineBackend.figure_format = 'svg'
from pylab import *
3.5.1 Quasi-Newton methods ¶ Newton method requires the knowledge of the gradient. In many situations, it may not be possible to find the gradient. E.g., if f ( x ) f(x) f ( x ) [1] quasi-Newton methods . For a scalar function f ( x ) f(x) f ( x )
f ′ ( x k ) ≈ f ( x k + h ) − f ( x k ) h , 0 ≤ h ≪ 1 f'(x_k) \approx \frac{f(x_k+h) - f(x_k)}{h}, \qquad 0 \le h \ll 1 f ′ ( x k ) ≈ h f ( x k + h ) − f ( x k ) , 0 ≤ h ≪ 1 Since we want to generate a sequence of approximations to the root, we can use two successive values of x k x_k x k
f ′ ( x k ) ≈ f ( x k ) − f ( x k − 1 ) x k − x k − 1 f'(x_k) \approx \frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} f ′ ( x k ) ≈ x k − x k − 1 f ( x k ) − f ( x k − 1 ) Using this derivative in place of the exact derivative in Newton method leads to the secant method .
3.5.2 Secant method ¶ Given two starting values x 0 x_0 x 0 x 1 x_1 x 1 r r r f ( x ) f(x) f ( x ) ( x 0 , f ( x 0 ) ) (x_0,f(x_0)) ( x 0 , f ( x 0 )) ( x 1 , f ( x 1 ) ) (x_1, f(x_1)) ( x 1 , f ( x 1 )) x 2 x_2 x 2
x 2 = x 1 − f ( x 1 ) x 1 − x 0 f ( x 1 ) − f ( x 0 ) x_2 = x_1 - f(x_1) \frac{x_1 - x_0}{f(x_1) - f(x_0)} x 2 = x 1 − f ( x 1 ) f ( x 1 ) − f ( x 0 ) x 1 − x 0 Inductively, we obtain
x k + 1 = x k − f ( x k ) x k − x k − 1 f ( x k ) − f ( x k − 1 ) , k = 1 , 2 , … x_{k+1} = x_k - f(x_k) \frac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})}, \qquad k=1,2,\ldots x k + 1 = x k − f ( x k ) f ( x k ) − f ( x k − 1 ) x k − x k − 1 , k = 1 , 2 , … Note that this is just Newton iteration formulae where the derivative has been approximated by the finite difference formula.
3.5.3 Convergence analysis ¶ From the iteration formula, we get by subtracting r r r
r − x k + 1 = r − x k + f ( x k ) x k − x k − 1 f ( x k ) − f ( x k − 1 ) r - x_{k+1} = r - x_k + f(x_k) \frac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})} r − x k + 1 = r − x k + f ( x k ) f ( x k ) − f ( x k − 1 ) x k − x k − 1 Define the Newton divided differences
f [ x k − 1 , x k ] : = f ( x k ) − f ( x k − 1 ) x k − x k − 1 , f [ x k − 1 , x k , r ] : = f [ x k , r ] − f [ x k − 1 , x k ] r − x k − 1 f[x_{k-1},x_k] := \frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}}, \qquad f[x_{k-1},x_k,
r]:= \frac{f[x_k,r] - f[x_{k-1},x_k]}{r - x_{k-1}} f [ x k − 1 , x k ] := x k − x k − 1 f ( x k ) − f ( x k − 1 ) , f [ x k − 1 , x k , r ] := r − x k − 1 f [ x k , r ] − f [ x k − 1 , x k ] then it is easy to show that
f [ x k − 1 , x k ] = f ′ ( y k ) , for some y k between x k − 1 , x k f[x_{k-1},x_k] = f'(y_k), \qquad \textrm{for some $y_k$ between $x_{k-1}, x_k$} f [ x k − 1 , x k ] = f ′ ( y k ) , for some y k between x k − 1 , x k and
f [ x k − 1 , x k , r ] = 1 2 f ′ ′ ( z k ) for some z k between x k − 1 , x k , r f[x_{k-1},x_k,r] = \half f''(z_k) \qquad \textrm{for some $z_k$ between
$x_{k-1},x_k,r$} f [ x k − 1 , x k , r ] = 2 1 f ′′ ( z k ) for some z k between x k − 1 , x k , r Hence we get the error formula
r − x k + 1 = − ( r − x k − 1 ) ( r − x k ) f [ x k − 1 , x k , r ] 2 f [ x k − 1 , x k ] = − ( r − x k − 1 ) ( r − x k ) f ′ ′ ( z k ) 2 f ′ ( y k ) \begin{align}
r - x_{k+1} &= -(r - x_{k-1})(r - x_k) \frac{f[x_{k-1},x_k,r]}
{2f[x_{k-1},x_k]} \\
&= -(r - x_{k-1})(r - x_k)\frac{f''(z_k)}{2f'(y_k)}
\end{align} r − x k + 1 = − ( r − x k − 1 ) ( r − x k ) 2 f [ x k − 1 , x k ] f [ x k − 1 , x k , r ] = − ( r − x k − 1 ) ( r − x k ) 2 f ′ ( y k ) f ′′ ( z k ) If we can bound the last factor
∣ f ′ ′ ( z k ) 2 f ′ ( y k ) ∣ ≤ M < ∞ \left| \frac{f''(z_k)}{2f'(y_k)} \right| \le M < \infty ∣ ∣ 2 f ′ ( y k ) f ′′ ( z k ) ∣ ∣ ≤ M < ∞ then we get a bound on the error e k = ∣ r − x k ∣ e_k = |r - x_k| e k = ∣ r − x k ∣
e k + 1 ≤ M e k − 1 e k e_{k+1} \le M e_{k-1} e_k e k + 1 ≤ M e k − 1 e k Assume that f ( x ) f(x) f ( x ) f ′ ( x ) f'(x) f ′ ( x ) f ′ ′ ( x ) f''(x) f ′′ ( x ) x x x r r r f ( r ) = 0 f(r)=0 f ( r ) = 0 f ′ ( r ) ≠ 0 f'(r) \ne 0 f ′ ( r ) = 0 x 0 x_0 x 0 x 1 x_1 x 1 r r r r r r
p = 1 2 ( 1 + 5 ) ≈ 1.62 p = \half (1 + \sqrt{5}) \approx 1.62 p = 2 1 ( 1 + 5 ) ≈ 1.62 Recall that Newton method converges with order p = 2 p=2 p = 2
(1) Pick a sufficiently small interval I = [ r − δ , r + δ ] I = [r-\delta, r+\delta] I = [ r − δ , r + δ ] f ′ ( x ) ≠ 0 f'(x) \ne 0 f ′ ( x ) = 0 f ′ ( r ) ≠ 0 f'(r) \ne 0 f ′ ( r ) = 0 f ′ f' f ′
M = max x ∈ I ∣ f ′ ′ ( x ) ∣ 2 min x ∈ I ∣ f ′ ( x ) ∣ < ∞ M = \frac{\max_{x \in I} |f''(x)|}{2\min_{x \in I} |f'(x)|} < \infty M = 2 min x ∈ I ∣ f ′ ( x ) ∣ max x ∈ I ∣ f ′′ ( x ) ∣ < ∞ For any x 0 , x 1 ∈ I x_0, x_1 \in I x 0 , x 1 ∈ I
e 2 ≤ e 1 ⋅ M e 0 M e 2 ≤ M e 1 ⋅ M e 0 \begin{aligned}
e_2 &\le e_1 \cdot M e_0 \\
M e_2 &\le Me_1 \cdot Me_0
\end{aligned} e 2 M e 2 ≤ e 1 ⋅ M e 0 ≤ M e 1 ⋅ M e 0 Further assume that x 0 x_0 x 0 x 1 x_1 x 1
Δ = max { M e 0 , M e 1 } < 1 \Delta = \max\{ Me_0 , Me_1\} < 1 Δ = max { M e 0 , M e 1 } < 1 Then
M e 2 ≤ Δ 2 < 1 Me_2 \le \Delta^2 < 1 M e 2 ≤ Δ 2 < 1 Also
M e 2 ≤ Δ 2 < Δ M e_2 \le \Delta^2 < \Delta M e 2 ≤ Δ 2 < Δ ⟹ e 2 < Δ M = max { e 0 , e 1 } ≤ δ ⟹ x 2 ∈ I \begin{aligned}
& \implies e_2 < \frac{\Delta}{M} = \max\{e_0, e_1\} \le \delta \\
& \implies x_2 \in I
\end{aligned} ⟹ e 2 < M Δ = max { e 0 , e 1 } ≤ δ ⟹ x 2 ∈ I Applying the above argument inductively, we get
x k ∈ I M e k ≤ Δ < 1 } for k = 0 , 1 , 2 , 3 , … \left.
\begin{aligned}
x_k & \in I \\
M e_k & \le \Delta < 1
\end{aligned}
\right\} \quad \textrm{for } k = 0,1,2,3,\ldots x k M e k ∈ I ≤ Δ < 1 } for k = 0 , 1 , 2 , 3 , … To prove convergence, we apply the error inequality recursively.
M e 0 ≤ Δ M e 1 ≤ Δ M e 2 ≤ M e 1 ⋅ M e 0 ≤ Δ 2 M e 3 ≤ M e 2 ⋅ M e 1 ≤ Δ 3 M e 4 ≤ M e 3 ⋅ M e 2 ≤ Δ 5 M e 5 ≤ M e 4 ⋅ M e 3 ≤ Δ 8 \begin{aligned}
M e_0 & \le \Delta \\
M e_1 & \le \Delta \\
M e_2 & \le Me_1 \cdot Me_0 \le \Delta^2 \\
M e_3 & \le M e_2 \cdot M e_1 \le \Delta^3 \\
M e_4 & \le M e_3 \cdot M e_2 \le \Delta^5 \\
M e_5 & \le M e_4 \cdot M e_3 \le \Delta^8
\end{aligned} M e 0 M e 1 M e 2 M e 3 M e 4 M e 5 ≤ Δ ≤ Δ ≤ M e 1 ⋅ M e 0 ≤ Δ 2 ≤ M e 2 ⋅ M e 1 ≤ Δ 3 ≤ M e 3 ⋅ M e 2 ≤ Δ 5 ≤ M e 4 ⋅ M e 3 ≤ Δ 8 In general
M e k ≤ Δ q k M e_k \le \Delta^{q_k} M e k ≤ Δ q k with
q 0 = q 1 = 1 , q k = q k − 1 + q k − 2 q_0 = q_1 = 1, \qquad q_k = q_{k-1} + q_{k-2} q 0 = q 1 = 1 , q k = q k − 1 + q k − 2 We see that q k q_k q k Δ < 1 \Delta < 1 Δ < 1 q k → ∞ q_k \to \infty q k → ∞ k → ∞ k \to \infty k → ∞
lim k → ∞ e k ≤ 1 M Δ q k = 0 \lim_{k \to \infty} e_k \le \frac{1}{M} \Delta^{q_k} = 0 k → ∞ lim e k ≤ M 1 Δ q k = 0 Hence the secant method converges to some root.
(2) Now let us estimate the convergence rate p p p
e k ≈ A e k − 1 p e_k \approx A e_{k-1}^p e k ≈ A e k − 1 p and we will try to estimate p p p
e k − 1 ≈ ( e k A ) 1 / p e_{k-1} \approx \left( \frac{e_k}{A} \right)^{1/p} e k − 1 ≈ ( A e k ) 1/ p Using the error relation for secant method
e k + 1 ≈ M e k ⋅ e k − 1 ≈ M e k ⋅ 1 A 1 / p e k 1 / p A e k p ≈ M A − 1 / p e k 1 + 1 / p \begin{align}
e_{k+1} & \approx M e_k \cdot e_{k-1} \approx M e_k \cdot \frac{1}{A^{1/p}} e_k^{1/p} \\
A e_k^p & \approx M A^{-1/p} e_k^{1 + 1/p}
\end{align} e k + 1 A e k p ≈ M e k ⋅ e k − 1 ≈ M e k ⋅ A 1/ p 1 e k 1/ p ≈ M A − 1/ p e k 1 + 1/ p This implies that
p = 1 + 1 p and A ≈ M A − 1 / p p = 1 + \frac{1}{p} \qquad \textrm{and} \qquad A \approx M A^{-1/p} p = 1 + p 1 and A ≈ M A − 1/ p The first equation p 2 − p − 1 = 0 p^2 -p - 1 = 0 p 2 − p − 1 = 0
p = 1 ± 5 2 p = \frac{1 \pm \sqrt{5}}{2} p = 2 1 ± 5 We choose the positive root since we need p > 1 p > 1 p > 1
p = 1 + 5 2 ≈ 1.618 p = \frac{1 + \sqrt{5}}{2} \approx 1.618 p = 2 1 + 5 ≈ 1.618 This is less than quadratic convergence, but still better than linear convergence. From the second relation we get
A ≈ M p / ( 1 + p ) ≈ ∣ f ′ ′ ( r ) 2 f ′ ( r ) ∣ 5 − 1 2 A \approx M^{p/(1+p)} \approx \left| \frac{ f''(r)}{2f'(r)} \right|
^{\frac{\sqrt{5} - 1}{2}} A ≈ M p / ( 1 + p ) ≈ ∣ ∣ 2 f ′ ( r ) f ′′ ( r ) ∣ ∣ 2 5 − 1 3.5.4 Error estimate ¶ We can use the function value to decide on convergence but this may not be a good criterion if the function is very flat around the root. Additionally we can estimate the error in the root also. From mean value theorem
f ( x k ) = f ( x k ) − f ( r ) = f ′ ( ξ k ) ( x k − r ) f(x_k) = f(x_k) - f(r) = f'(\xi_k) (x_k - r) f ( x k ) = f ( x k ) − f ( r ) = f ′ ( ξ k ) ( x k − r ) where ξ k \xi_k ξ k x k x_k x k r r r
r − x k = − f ( x k ) f ′ ( ξ k ) r - x_k = - \frac{f(x_k)}{f'(\xi_k)} r − x k = − f ′ ( ξ k ) f ( x k ) If f ′ f' f ′ x k x_k x k r r r
f ′ ( ξ k ) ≈ f ( x k + 1 ) − f ( x k ) x k + 1 − x k f'(\xi_k) \approx \frac{f(x_{k+1}) - f(x_k)}{x_{k+1} - x_k} f ′ ( ξ k ) ≈ x k + 1 − x k f ( x k + 1 ) − f ( x k ) We get the error estimate
r − x k ≈ − f ( x k ) f ( x k + 1 ) − f ( x k ) x k + 1 − x k = x k + 1 − x k r - x_k \approx - \frac{f(x_k)}{\frac{f(x_{k+1}) - f(x_k)}{x_{k+1} - x_k}} = x_{k+1}
- x_k r − x k ≈ − x k + 1 − x k f ( x k + 1 ) − f ( x k ) f ( x k ) = x k + 1 − x k and the relative error estimate
r − x k r ≈ x k + 1 − x k x k + 1 \frac{r - x_k}{r} \approx \frac{x_{k+1} - x_k}{x_{k+1}} r r − x k ≈ x k + 1 x k + 1 − x k Hence we can use the test
∣ x k + 1 − x k ∣ < ϵ ∣ x k + 1 ∣ , 0 < ϵ ≪ 1 |x_{k+1} - x_k| < \epsilon |x_{k+1}|, \qquad 0 < \epsilon \ll 1 ∣ x k + 1 − x k ∣ < ϵ ∣ x k + 1 ∣ , 0 < ϵ ≪ 1 to decide on convergence. This is used in the code examples to stop the
iterations.
3.5.5 Accuracy of FD approximation ¶ Compute derivative of
f ( x ) = sin ( x ) f(x) = \sin(x) f ( x ) = sin ( x ) at x = 2 π x=2\pi x = 2 π
f ( x + h ) − f ( x ) h \frac{ f(x+h) - f(x) }{h} h f ( x + h ) − f ( x ) for h = 1 0 − 1 , 1 0 − 2 , … , 1 0 − 14 h = 10^{-1}, 10^{-2}, \ldots, 10^{-14} h = 1 0 − 1 , 1 0 − 2 , … , 1 0 − 14
f = lambda x: sin(x)
h = 10.0**arange(-1,-15,-1)
df= zeros(len(h))
x = 2.0*pi
f0= f(x)
for i in range(len(h)):
f1 = f(x+h[i])
df[i] = (f1 - f0)/h[i]
loglog(h,abs(df-1.0),'o-')
grid(True), xlabel('h'), ylabel('Error in derivative');Initially, as h h h h h h h h h
The secant method can be written as
x k + 1 = x k − f ( x k ) g k , g k = f ( x k ) − f ( x k − 1 ) x k − x k − 1 x_{k+1} = x_k - \frac{f(x_k)}{g_k}, \qquad g_k = \frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} x k + 1 = x k − g k f ( x k ) , g k = x k − x k − 1 f ( x k ) − f ( x k − 1 ) As x k → r x_k \to r x k → r x k − x k − 1 x_k - x_{k-1} x k − x k − 1 f ( x k ) − f ( x k − 1 ) f(x_k) - f(x_{k-1}) f ( x k ) − f ( x k − 1 ) g k g_k g k
f ( x + h ) − f ( x ) h = f ′ ( x ) + 1 2 h f ′ ′ ( ξ ) , ξ between x and x + h \frac{f(x+h) - f(x)}{h} = f'(x) + \half h f''(\xi), \qquad \textrm{$\xi$ between $x$ and $x+h$} h f ( x + h ) − f ( x ) = f ′ ( x ) + 2 1 h f ′′ ( ξ ) , ξ between x and x + h We expect that as h → 0 h \to 0 h → 0 f ( x ) f(x) f ( x )
f ^ ( x + h ) − f ^ ( x ) h \frac{\hat{f}(x+h) - \hat{f}(x)}{h} h f ^ ( x + h ) − f ^ ( x ) where, assuming f ( x ) = O ( 1 ) f(x) = O(1) f ( x ) = O ( 1 )
f ^ ( x + h ) = f ( x + h ) + ϵ 2 , f ^ ( x ) = f ( x ) + ϵ 1 \hat{f}(x+h) = f(x+h) + \epsilon_2, \quad \hat{f}(x) = f(x) + \epsilon_1 f ^ ( x + h ) = f ( x + h ) + ϵ 2 , f ^ ( x ) = f ( x ) + ϵ 1 and
ϵ 1 , ϵ 2 ≈ C u , the unit round of the computer \epsilon_1, \epsilon_2 \approx C \uround, \quad \textrm{the unit round of the computer} ϵ 1 , ϵ 2 ≈ C u , the unit round of the computer Hence
The error consists of
discretization error: 1 2 h f ′ ′ ( ξ ) \half h f''(\xi) 2 1 h f ′′ ( ξ ) round-off error: C u h \frac{C\uround}{h} h C u As h → 0 h \to 0 h → 0 h h h
h o p t = ( 2 C u f ′ ′ ) 1 2 h_{opt} = \left( \frac{2C \uround}{f''} \right)^\half h o pt = ( f ′′ 2 C u ) 2 1 For the above finite difference approximation, the optimum step size h = O ( u ) = O ( 1 0 − 7 ) h = O(\sqrt{\uround}) = \order{10^{-7}} h = O ( u ) = O ( 1 0 − 7 )
∣ x k − x k − 1 ∣ > ∣ x k − 1 + x k 2 ∣ u |x_k - x_{k-1}| > \left| \frac{x_{k-1} + x_k}{2} \right| \sqrt{\uround} ∣ x k − x k − 1 ∣ > ∣ ∣ 2 x k − 1 + x k ∣ ∣ u use the g k g_k g k
g k = f ( x k + h ) − f ( x k ) h , h = ∣ x k ∣ u g_k = \frac{f(x_k + h) - f(x_k)}{h}, \qquad h = |x_k| \sqrt{\uround} g k = h f ( x k + h ) − f ( x k ) , h = ∣ x k ∣ u 3.5.6 Complex step method ¶ We can take a complex step instead of a real step. Taylor expansion with a small complex step yields
f ( x + i h ) = f ( x ) + i h f ′ ( x ) − 1 2 h 2 f ′ ′ ( x ) + O ( i h 3 ) f(x+\ii h) = f(x) + \ii h f'(x) - \half h^2 f''(x) + O(\ii h^3) f ( x + i h ) = f ( x ) + i h f ′ ( x ) − 2 1 h 2 f ′′ ( x ) + O ( i h 3 ) Hence
imag f ( x + i h ) = h f ′ ( x ) + O ( h 3 ) \imag f(x+\ii h) = h f'(x) + O(h^3) imag f ( x + i h ) = h f ′ ( x ) + O ( h 3 ) and
imag f ( x + i h ) h = f ′ ( x ) + O ( h 2 ) \frac{\imag f(x+\ii h)}{h} = f'(x) + O(h^2) h imag f ( x + i h ) = f ′ ( x ) + O ( h 2 ) is a second order accurate approximation to the derivative. This formula does not involve subtraction of near equal quantities and hence it does not suffer from round-off error even if we take h h h h h h h = 1 0 − 20 h = 10^{-20} h = 1 0 − 20