Newton’s Method for Nonlinear Equations and Unconstrained Minimization

Exercise 1

[1]:
import numpy

def F(x):
    return numpy.asfarray([2 * (x[0] + x[1])**2 + (x[0] - x[1])**2 - 8,
                           5 * x[0]**2 + (x[1] - 3)**2 - 9])

def J(x):
    return numpy.asfarray([[6 * x[0] + 2 * x[1], 2 * x[0] + 6 * x[1]],
                           [10 * x[0], 2 * x[1] - 6]])

x_c = numpy.asfarray([2, 0])
for i in range(2):
    s_k = numpy.dot(numpy.linalg.inv(J(x_c)), -F(x_c))
    x_k = x_c + s_k
    print('x_{0} = {1}'.format(i, x_c))
    x_c = x_k
x_0 = [2. 0.]
x_1 = [1.31578947 1.05263158]

Exercise 2

The Jacobian becomes ill-conditioned in the second iteration i.e. it becomes the zero matrix.

[2]:
import numpy

def F(x):
    return numpy.exp(x) - 1

def J(x):
    return numpy.asfarray([[numpy.exp(x[0]), 0],
                           [0, numpy.exp(x[1])]])

x_c = numpy.asfarray([-10, -10])
for i in range(2):
    _ = J(x_c)
    s_k = numpy.dot(numpy.linalg.inv(_), -F(x_c))
    x_k = x_c + s_k
    print('Iteration {0}'.format(i))
    print('x_{0} = {1}'.format(i, x_c))
    print('s_{0} = {1}'.format(i, s_k))
    print('J_{0} =\n{1}'.format(i, _))
    print('x_{0} = {1}'.format(i + 1, x_k))
    x_c = x_k
Iteration 0
x_0 = [-10. -10.]
s_0 = [22025.46579481 22025.46579481]
J_0 =
[[4.53999298e-05 0.00000000e+00]
 [0.00000000e+00 4.53999298e-05]]
x_1 = [22015.46579481 22015.46579481]
Iteration 1
x_1 = [22015.46579481 22015.46579481]
s_1 = [nan nan]
J_1 =
[[inf  0.]
 [ 0. inf]]
x_2 = [nan nan]

Exercise 3

Define \(f_1(x) = x, f_2(x) = x^2 + x,\) and \(f_3(x) = e^x - 1\).

(a)

\[f'_1(x) = 1 \Rightarrow f'_1(x_* = 0) = 1\]
\[f'_2(x) = 2x + 1 \Rightarrow f'_2(x_* = 0) = 1\]
\[f'_3(x) = e^x \Rightarrow f'_3(x_* = 0) = 1\]

(b)

The Lipschitz constant for each derivative at the root \(x_* = 0\) is

\[\begin{split}\left| f'_1(x) - f'_1(0) \right| &\leq \gamma_1 \left| x - 0 \right|\\ \left| \frac{1 - 1}{x} \right| &\leq \gamma_1\\ 0 &\leq \gamma_1\end{split}\]
\[\begin{split}\left| f'_2(x) - f'_2(0) \right| &\leq \gamma_2 \left| x - 0 \right|\\ \left| \frac{2x + 1 - 1}{x} \right| &\leq \gamma_2\\ 2 &\leq \gamma_2\end{split}\]
\[\begin{split}\left| f'_3(x) - f'_3(0) \right| &\leq \gamma_3 \left| x - 0 \right|\\ \left| \frac{e^x - 1}{x} \right| &\leq \gamma_3\\ &\leq e^x \leq \gamma_3\end{split}\]

The last inequality holds because the Bernoulli inequality gives

\[1 + x \leq \left( 1 + \frac{x}{n} \right)^n \xrightarrow{n \rightarrow \infty} \exp x.\]

Note that \(e^x\) is an upper bound when \(x > 0\) and becomes a lower bound when \(x < 0\).

(c)

The derivative of each function is bounded on the interval \([-a, a]\), so each function is Lipschitz continuous. The bound for each derivative is

\[\begin{split}\left| f'_1(x) \right| &= 1\\ &\geq \rho_1 = 1\end{split}\]
\[\begin{split}\left| f'_2(x) \right| &= \left| 2x + 1 \right|\\ &\geq \rho_2 = 0\end{split}\]
\[\begin{split}\left| f'_3(x) \right| &= \left| e^x \right|\\ &\geq \rho_3 = e^{-a}\end{split}\]

The smallest region of convergence is defined to be \(\eta = \min \left\{ \hat{\eta}, 2 \rho / \gamma \right\}\) where \(\hat{\eta}\) is defined to be the radius of the largest open interval around \(x_*\).

(d)

\(f_1\) will converge to \(x_* = 0\) in the interval \((-\infty, \infty)\) since it only has one root and \(\lim_{\gamma \rightarrow 0^+} \frac{2 \rho_1}{\gamma_1} = \infty\).

\(f_2\) violates the assumption that \(\rho_2 > 0\) at \(x = -0.5\) causing Newton’s method to fail because the derivative is zero. Splitting the interval into \((-\infty, -0.5)\) and \((-0.5, \infty)\) sidesteps that issue. Applying Theorem 2.4.3 to each region shows that any point in those regions will converge.

\(f_3\) will converge to \(x_* = 0\) in the interval \((-\infty, \infty)\) since it only has one root in theory and \(\lim_{a \rightarrow \infty} \frac{2 \rho_3}{\gamma_3} = \lim_{a \rightarrow \infty} 2 e^{-2a} = 0\). However, the floating-point precision restricts the region to \([-6.5, 709]\).

[3]:
import numpy

def Newtons_Method(f, f_prime, x_c):
    return x_c - f(x_c) / f_prime(x_c)

f = [lambda x: x,
     lambda x: x**2 + x,
     lambda x: numpy.exp(x) - 1
    ]
f_prime = [lambda x: 1,
           lambda x: 2 * x + 1,
           lambda x: numpy.exp(x)]
x_c = [-7.5, -0.6, 7.5]
for i in range(16):
    for j in range(len(x_c)):
        x_c[j] = Newtons_Method(f[j], f_prime[j], x_c[j])
    print('x_{0} = {1}'.format(i + 1, x_c))
x_1 = [0.0, -1.8000000000000003, 6.500553084370148]
x_2 = [0.0, -1.2461538461538462, 5.502055692264317]
x_3 = [0.0, -1.0406026962727994, 4.506134071187519]
x_4 = [0.0, -1.0015247601944897, 3.5171751329216505]
x_5 = [0.0, -1.000002317825395, 2.546858300770652]
x_6 = [0.0, -1.0000000000053721, 1.6251856616292897]
x_7 = [0.0, -1.0, 0.8220607812846256]
x_8 = [0.0, -1.0, 0.2615857370546373]
x_9 = [0.0, -1.0, 0.031415606703852794]
x_10 = [0.0, -1.0, 0.0004883429491289379]
x_11 = [0.0, -1.0, 1.192200103575106e-07]
x_12 = [0.0, -1.0, 7.157097849691172e-15]
x_13 = [0.0, -1.0, 5.16704920902205e-17]
x_14 = [0.0, -1.0, 5.16704920902205e-17]
x_15 = [0.0, -1.0, 5.16704920902205e-17]
x_16 = [0.0, -1.0, 5.16704920902205e-17]

Exercise 4

\[\begin{split}J(x) = \begin{bmatrix} 1 & 0 & 0\\ 0 & 2 x_2 + 1 & 0\\ 0 & 0 & e^{x_3} \end{bmatrix}\end{split}\]

(a)

\[\begin{split}J(x_*) = \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix}\end{split}\]

(b)

The Lipschitz constant on \(J(x) \in \text{Lip}_\gamma(N(x_*, a))\) is

\[\begin{split}\left\Vert J(x) - J(x_*) \right\Vert &\leq \gamma \left\Vert x - x_* \right\Vert & \quad & \text{(4.1.11)}\\ \left\Vert \begin{bmatrix} 0 & 0 & 0\\ 0 & 2 x_2 & 0\\ 0 & 0 & e^{x_3} - 1 \end{bmatrix} \right\Vert &< \gamma a & \quad & d(x_*, x) < a\\ \frac{ \max \left\{ \left\vert 2 x_2 \right\vert, \left\vert e^{x_3} - 1 \right\vert \right\} }{a} &< \gamma & \quad & \text{1-norm is an upper bound for the numerator}\end{split}\]

(c)

1-norm is the largest possible value for \(\left\Vert J(x_*)^{-1} \right\Vert\), so the region of convergence is defined as \(\epsilon = \min \left\{ a, \frac{1}{2 \beta \gamma} \right\}\) where \(0 < \left\Vert J(x_*)^{-1} \right\Vert \leq 1 \leq \beta\). Hence

\[\epsilon = \min \left\{ a, \frac{a}{ \max \left\{ \left\vert 4 x_2 \right\vert, \left\vert 2e^{x_3} - 2 \right\vert \right\} } \right\}\]

(d)

When \((x_0)_3 = 0, \epsilon = \min \left\{ a, \frac{a}{\left\vert 4 x_2 \right\vert} \right\}\). When \((x_0)_2 = (x_0)_3 = 0, \epsilon = a\) because as \((x_0)_2 \rightarrow 0\) or \((x_0)_3 \rightarrow 0\), the denominator grows infinitely large.

Exercise 5

Proof by induction starts at the base case of \(x_0\), so the goal is to show \(x_1 = x_0 - J(x_0)^{-1} F(x_0)\) is well-defined i.e. \(J(x_0)\) is nonsingular.

Since \(J(x_*)^{-1}\) exists and \(\left\Vert J(x_*)^{-1} \right\Vert \leq \beta\) where \(\beta > 0\) by assumption, \(J(x_*)\) is nonsingular. Theorem 3.14 states that \(J(x_0)\) is nonsingular when \(\left\Vert J(x_*)^{-1} \left( J(x_0) - J(x_*) \right) \right\Vert < 1\) holds. Recall that \(J\) is Holder continuous within \(N(x_*, r)\) where \(r > 0\) and \(\alpha \in (0, 1]\) such that

\[\begin{split}\left\Vert J(x_*)^{-1} \left( J(x_0) - J(x_*) \right) \right\Vert &\leq \left\Vert J(x_*)^{-1} \right\Vert \left\Vert J(x_0) - J(x_*) \right\Vert\\ &\leq \beta \gamma \left\Vert x_0 - x_* \right\Vert^\alpha\\ &\leq \beta \gamma r^\alpha\end{split}\]

Picking \(\epsilon = \min \left\{ r, \left( \frac{\alpha}{(1 + \alpha) \beta \gamma} \right)^{1 / \alpha} \right\}\) satisfies Theorem 3.14 in addition to revealing the elegant inequality \(\beta \gamma \epsilon^\alpha \leq \frac{\alpha}{1 + \alpha} < 1\) that simplifies later derivations. Thus \(J(x_0)\) is nonsingular and its induced matrix norm is bounded above by

\[\begin{split}\left\Vert J(x_0) \right\Vert &\leq \frac{ \left\Vert J(x_*)^{-1} \right\Vert }{ 1 - \left\Vert J(x_*)^{-1} \left( J(x_0) - J(x_*) \right) \right\Vert } & \quad & \text{perturbation relation 3.1.20}\\ &\leq \beta \frac{1}{1 - \beta \gamma \epsilon^\alpha}\\ &\leq \beta (1 + \alpha).\end{split}\]

Observe that the sequence \(x_1, x_2, \ldots\) converges to \(x_*\) because

\[\begin{split}x_1 - x_* &= x_0 - J(x_0)^{-1} F(x_0) - x_*\\ &= x_0 - J(x_0)^{-1} \lbrack F(x_0) - F(x_*) \rbrack - x_* & \quad & F(x_*) = 0\\ &= J(x_0)^{-1} \lbrack F(x_*) - F(x_0) - J(x_0) (x_* - x_0) \rbrack\\ \left\Vert x_1 - x_* \right\Vert &= \left\Vert J(x_0)^{-1} \lbrack F(x_*) - F(x_0) - J(x_0) (x_* - x_0) \rbrack \right\Vert\\ \left\Vert x_1 - x_* \right\Vert &\leq \left\Vert J(x_0)^{-1} \right\Vert \left\Vert F(x_*) - F(x_0) - J(x_0) (x_* - x_0) \right\Vert & \quad & \text{(3.1.10)}\\ &\leq \beta (1 + \alpha) \left\Vert x_* - x_0 \right\Vert^{1 + \alpha} \frac{\gamma}{1 + \alpha} & \quad & p = x_* - x_0\\ &= \beta \gamma \left\Vert x_0 - x_* \right\Vert^{1 + \alpha} & \quad & \text{(3.1.1b).}\end{split}\]

To understand why \(p = x_* - x_0\), see Exercise 4.8.

Exercise 6

By definition of differentiability and continuity,

\[f(t) = \frac{\gamma}{2} t^2 - \frac{t}{\beta} + \frac{\eta}{\beta}\]

is continuously differentiable \(\forall t \in \mathbb{R}\), \(\gamma, \eta \in \mathbb{R}^+_0\), and \(\beta \in \mathbb{R}^+\). Assuming \(t_0 = 0\), \(N(t_0, r)\) is satisfied \(\forall r \in (0, \infty)\).

Since \(f(t): \mathbb{R} \rightarrow \mathbb{R}\), the vector norm and induced operator norm reduce respectively to an absolute value and \(J(t) = f'(t) = \gamma t - \frac{1}{\beta}\).

Notice that \(J \in \text{Lip}_\gamma(N(t_0, r))\) because \(\left\vert f'(t) \right\vert < \left\vert \gamma r \right\vert + \left\vert \frac{1}{\beta} \right\vert\).

The initial conditions for the induced operator norms

\[\left\Vert J(t_0)^{-1} \right\Vert = \left\vert \left( \gamma t_0 - \frac{1}{\beta} \right)^{-1} \right\vert = \beta\]
\[\left\Vert J(t_0)^{-1} F(t_0) \right\Vert = \left\vert \left( \gamma t_0 - \frac{1}{\beta} \right)^{-1} \left( \frac{\gamma}{2} t_0^2 - \frac{t_0}{\beta} + \frac{\eta}{\beta} \right) \right\vert = \eta\]

are satisfied by definition of \(f(t)\).

Define \(\gamma_\text{Rel} = \beta \gamma, \alpha = \gamma_\text{Rel} \eta\) for ease of notation. Analytically solving for the roots of the univariate quadratic function \(f(t_*) = 0\) gives

\[\begin{split}t_* &= \frac{ \beta^{-1} \pm \sqrt{\beta^{-2} - 2 \gamma \eta \beta^{-1}} }{\gamma}\\ &= \frac{\beta^{-1} \pm \sqrt{\frac{1 - 2 \alpha}{\beta^2}}}{\gamma}\\ &= \frac{1 \pm \sqrt{1 - 2 \alpha}}{\beta \gamma}.\end{split}\]

Thus \(\alpha \leq \frac{1}{2}\) is a necessary condition for \(f(t)\) to have real roots.

The following results will be useful in proving \(\{t_k\}\) converges.

Proof of \(\left\vert f'(t_{k + 1}) \right\vert \geq \left\vert \frac{f'(t_k)}{2} \right\vert\) via induction. The following shows that the base case holds:

\[\begin{split}\left\vert f'(t_1) \right\vert &= \left\vert \gamma t_1 - \frac{1}{\beta}\right\vert\\ &= \left\vert \gamma \left( t_0 - J(t_0)^{-1} F(t_0) \right) - \frac{1}{\beta} \right\vert\\ &= \left\vert \gamma \eta - \frac{1}{\beta}\right\vert\\ &= \left\vert \frac{\alpha - 1}{\beta}\right\vert\\ &= \frac{1 - \alpha}{\beta}\\ &\geq \left\vert \frac{f'(t_0)}{2} \right\vert & \quad & \frac{1}{2 \beta} \leq \frac{1 - \alpha}{\beta} \leq \frac{1}{\beta}\end{split}\]

Assume \(\left\vert f'(t_k) \right\vert \geq \left\vert \frac{f'(t_{k - 1})}{2} \right\vert\) holds. The inductive step holds because

\[\begin{split}\left\vert f'(t_{k + 1}) \right\vert &= \left\vert \gamma t_{k + 1} - \frac{1}{\beta}\right\vert\\ &= \left\vert \gamma \left( t_k - J(t_k)^{-1} F(t_k) \right) - \frac{1}{\beta} \right\vert\\ &= \left\vert \gamma \left( t_k - \left( \frac{\beta \gamma t_k - 1}{\beta} \right)^{-1} \left( \frac{\beta \gamma t_k^2 - 2 t_k + 2 \eta}{2 \beta} \right) \right) - \frac{1}{\beta} \right\vert\\ &= \left\vert \gamma t_k - \frac{ \frac{\beta \gamma^2 t_k^2}{2} - \gamma t_k + \gamma \eta }{\beta \gamma t_k - 1} - \frac{1}{\beta} \right\vert\\ &= \left\vert \gamma t_k - \frac{ \frac{\gamma_\text{Rel}^2 t_k^2}{2} - \gamma_\text{Rel} t_k + \alpha + \gamma_\text{Rel} t_k - 1 }{\beta \gamma_\text{Rel} t_k - \beta} \right\vert\\ &= \left\vert \gamma t_k + \frac{ 1 - \alpha - \frac{\gamma_\text{Rel}^2 t_k^2}{2} }{\beta \gamma_\text{Rel} t_k - \beta} \right\vert\\ &= \left\vert \gamma t_k + \frac{ 1 - \alpha - \left( \beta^2 \gamma t_k - \beta \right) \frac{\gamma t_k}{2} - \frac{\beta \gamma t_k}{2} }{\beta^2 \gamma t_k - \beta} \right\vert\\ &= \left\vert \frac{\gamma t_k}{2} - \frac{ 1 - \alpha - \frac{\gamma_\text{Rel} t_k}{2} }{\beta \left( 1 - \gamma_\text{Rel} t_k \right)} \right\vert\\ &= \left\vert \frac{\gamma t_k}{2} - \frac{1}{2 \beta} \frac{ 2 - 2 \alpha - \gamma_\text{Rel} t_k }{1 - \gamma_\text{Rel} t_k} \right\vert\\ &\geq \left\vert \frac{\gamma t_k}{2} - \frac{1}{2 \beta} \right\vert\\ &= \left\vert \frac{f'(t_k)}{2} \right\vert.\end{split}\]

The last step can be proven via contradiction. Assume \(0 \leq \alpha \leq \frac{1}{2}\) and

\[\left\vert \gamma_\text{Rel} t_k - 1 \right\vert > \left\vert \gamma_\text{Rel} t_k - \frac{2 - 2 \alpha - \gamma_\text{Rel} t_k}{1 - \gamma_\text{Rel} t_k} \right\vert.\]

Given \(x = \gamma_\text{Rel} t_k\),

\[\begin{split}\left\vert x - 1 \right\vert &> \left\vert x - \frac{2 - 2 \alpha - x}{1 - x} \right\vert\\ \left\vert x - 1 \right\vert \left\vert 1 - x \right\vert &> \left\vert x - x^2 - 2 + 2 \alpha + x \right\vert\\ \left\vert x - 1 \right\vert^4 &> \left\vert x - x^2 - 2 + 2 \alpha + x \right\vert^2\\ (4 \alpha - 2) x^2 + (4 - 8 \alpha) x - 4 \alpha^2 + 8 \alpha - 3 &> 0\\ (2 \alpha - 1) (2x^2 - 4x + 3 - 2 \alpha) &> 0.\end{split}\]

Solving for roots of the quadratic equation yields

\[x = \frac{-(-4) \pm \sqrt{(-4)^2 - 4 (2) (3 - 2 \alpha)}}{2 (2)} = 1 \pm \sqrt{\alpha - \frac{1}{2}}\]

which means \(\alpha > \frac{1}{2}\) must hold to yield a real solution, but this contradicts the assumption that \(\alpha \leq \frac{1}{2}\). Therefore

\[\left\vert \gamma_\text{Rel} t_k - 1 \right\vert \leq \left\vert \gamma_\text{Rel} t_k - \frac{ 2 - 2 \alpha - \gamma_\text{Rel} t_k }{1 - \gamma_\text{Rel} t_k} \right\vert \quad \text{and} \quad \left\vert f'(t_{k + 1}) \right\vert \geq \left\vert 2^{-(k + 1)} \beta^{-1} \right\vert.\]

The upper bound \(\left\vert f'(t_{k + 1}) \right\vert \leq \left\vert f'(t_k) \right\vert\) can also be derived as

\[\begin{split}2 \left\vert x - 1 \right\vert &\geq \left\vert x - \frac{2 - 2 \alpha - x}{1 - x} \right\vert\\ 3x^4 - 12x^3 + (16 + 4 \alpha)x^2 - (8 + 8 \alpha)x + 8 \alpha - 4 \alpha^2 &\geq 0\\ (x^2 - 2x + 2 \alpha) (3x^2 - 6x + 4 - 2 \alpha) &\geq 0.\end{split}\]

Solving for roots of the quartic equation yields

\[x = \frac{-(-2) \pm \sqrt{(-2)^2 - 4 (1) (2 \alpha)}}{2 (1)} = 1 \pm \sqrt{1 - 2 \alpha}\]

and

\[x = \frac{-(-6) \pm \sqrt{(-6)^2 - 4 (3) (4 - 2 \alpha)}}{2 (3)} = 1 \pm \sqrt{\frac{2 \alpha - 1}{3}}.\]

Since \(0 \leq \alpha \leq \frac{1}{2}\), the second quadratic formula is inconsequential for real roots. Therefore the sequence \(\{ t_k \}\) is well-defined with \(N(t_0, r)\) because

\[\begin{split}t_{k + 1} - t_* &= t_k - J(t_k)^{-1} F(t_k) - t_*\\ &= J(t_k)^{-1} \left[ F(t_*) - F(t_k) - J(t_k) (t_* - t_k) \right]\\ \left\vert t_{k + 1} - t_* \right\vert &\leq 2^k \beta \frac{\gamma}{2} \left\vert t_* - t_k \right\vert^2 & \quad & \text{Lemma 4.1.12}\\ &= 2^{k - 1} \beta \gamma \left\vert t_k - t_* \right\vert^2 & \quad & \text{(3.1.1b)}\\ &\leq 2^{k - 1} r^2 \frac{\eta}{\alpha}.\end{split}\]

This inequality can be further simplified assuming \(r \geq r_0\)

\[\begin{split}\left\vert t_1 - t_* \right\vert &\leq 2^{-1} r^2 \frac{\eta}{\alpha}\\ &\leq 2^{-1} \left( 1 - \sqrt{1 - 2 \alpha} \right)^2 \frac{\eta}{\alpha}.\end{split}\]

Unrolling the recursion for \(t_k\) yields

\[\begin{split}\left\vert t_k - t_* \right\vert &\leq 2^{-k} \left( 1 - \sqrt{1 - 2 \alpha} \right)^{2^k} \frac{\eta}{\alpha}\\ &\leq (2 \alpha)^{2^k} \frac{\eta}{\alpha}.\end{split}\]

The last step can be proven via induction. The base case when \(k = 0\) gives \(1 - \sqrt{1 - 2 \alpha} \leq 2 \alpha\). Let us prove this via contradiction: assume \(0 \leq \alpha \leq \frac{1}{2}\) and \(1 - \sqrt{1 - 2 \alpha} > 2 \alpha\).

\[\begin{split}1 - \sqrt{1 - 2 \alpha} &> 2 \alpha\\ \left( 1 - 2 \alpha \right)^2 &> \sqrt{1 - 2 \alpha}^2\\ 4 \alpha^2 - 4 \alpha + 1 &> 1 - 2 \alpha\\ 4 \alpha^2 - 2 \alpha > 0\\ \alpha (2\alpha - 1) > 0\end{split}\]

shows that \(\alpha > \frac{1}{2}\) must hold to satisfy the inequality, but this contradicts the assumption that \(\alpha \leq \frac{1}{2}\). Thus the base case is true.

For the inductive step, assume \(2^{-k} \left( 1 - \sqrt{1 - 2 \alpha} \right)^{2^k} \leq \left( 2 \alpha \right)^{2^k}\) for \(0 \leq \alpha \leq \frac{1}{2}\).

\[\begin{split}2^{-(k + 1)} \left( 1 - \sqrt{1 - 2 \alpha} \right)^{2^{k + 1}} &\leq \left( 2 \alpha \right)^{2^{k + 1}}\\ -(k + 1) + 2^{k + 1} \log_2 \left( 1 - \sqrt{1 - 2 \alpha} \right) &\leq 2^{k + 1} \log_2 \left( 2 \alpha \right)\\ \frac{-(k + 1)}{2^{k + 1}} + \log_2 \left( 1 - \sqrt{1 - 2 \alpha} \right) &\leq \log_2 \left( 2 \alpha \right)\\ 2^{\frac{-(k + 1)}{2^{k + 1}}} \left( 1 - \sqrt{1 - 2 \alpha} \right) &\leq 2 \alpha\\ c \left( 1 - \sqrt{1 - 2 \alpha} \right) &\leq 2 \alpha & \quad & c = 2^{\frac{-(k + 1)}{2^{k + 1}}} > 0.\end{split}\]

By inspection, as \(k \rightarrow \infty\), \(c \rightarrow 1\), which completes the inductive step.

Exercise 7

\(\left\Vert J(x_n)^{-1} \right\Vert \leq \left\vert f'(t_n)^{-1} \right\vert\)

When \(n = 0\), Exercise 6 shows that \(\left\vert f'(t_0)^{-1} \right\vert = \beta\) and \(F\) is presumed to satisfy the assumptions of the Kantorovich theorem

\[\left\Vert J(x_0)^{-1} \right\Vert \leq \left\vert f'(t_0)^{-1} \right\vert.\]

For the inductive step \(n = k\), recall Theorem 3.14 states that \(J(x_k)\) is nonsingular when \(\left\Vert J(x_0)^{-1} \left( J(x_k) - J(x_0) \right) \right\Vert < 1\) holds.

\[\begin{split}\left\Vert J(x_0)^{-1} \left( J(x_k) - J(x_0) \right) \right\Vert &\leq \left\Vert J(x_0)^{-1} \right\Vert \left\Vert J(x_k) - J(x_0) \right\Vert\\ &\leq \beta \gamma r & \quad & J \in \text{Lip}_\gamma(N(x_0, r > 0))\end{split}\]

Picking \(r = \frac{1 - \sqrt{1 - 2 \alpha}}{\beta \gamma}\) where \(0 \leq \alpha \leq \frac{1}{2}\) satisfies Theorem 3.14. Thus \(J(x_k)\) is nonsingular and its induced matrix norm is bounded above by

\[\begin{split}\left\Vert J(x_k)^{-1} \right\Vert &\leq \frac{ \left\Vert J(x_0)^{-1} \right\Vert }{ 1 - \left\Vert J(x_0)^{-1} \left( J(x_k) - J(x_0) \right) \right\Vert } & \quad & \text{(3.1.20)}\\ &\leq \frac{\left\vert f'(t_0)^{-1} \right\vert}{ 1 - c} & \quad & 0 \leq c = 1 - \left\Vert J(x_0)^{-1} \left( J(x_k) - J(x_0) \right) \right\Vert < 1\\ &\leq \frac{\left\vert f'(t_k)^{-1} \right\vert}{1 - c}\\ (1 - c) \left\Vert J(x_k)^{-1} \right\Vert &\leq \left\vert f'(t_k)^{-1} \right\vert.\end{split}\]

The second to last inequality is explained in Exercise 6.

\(\left\Vert x_{n + 1} - x_n \right\Vert \leq \left\vert t_{n + 1} - t_n \right\vert\)

When \(n = 0\), Exercise 6 shows that \(\left\Vert J(t_0)^{-1} F(t_0) \right\Vert = \eta\) and \(F\) is presumed to satisfy the assumptions of the Kantorovich theorem

\[\left\Vert x_1 - x_0 \right\Vert \leq \left\vert t_1 - t_0 \right\vert.\]

The following relations will be useful for the induction.

\[\begin{split}\left\Vert F(x_k) \right\Vert &= \left\Vert F(x_k) - F(x_{k - 1}) - J(x_{k - 1}) (x_k - x_{k - 1}) \right\Vert & \quad & \text{(4.1.10)}\\ &= \left\Vert F(x_{k - 1} + p) - F(x_{k - 1}) - J(x_{k - 1}) p \right\Vert & \quad & p = x_k - x_{k - 1}\\ &\leq \frac{\gamma}{2} \left\Vert x_k - x_{k - 1} \right\Vert^2 & \quad & \text{Lemma 4.1.12}\end{split}\]

Applying the same analysis to \(f(t_k) = \frac{\gamma}{2} (t_k - t_{k - 1})^2\) gives

\[\begin{split}\left\vert f(t_k) \right\vert &= \left\vert f(t_k) - f(t_{k - 1}) - f'(t_{k - 1}) (t_k - t_{k - 1}) \right\vert & \quad & \text{(4.1.10)}\\ &= \left\vert \frac{\gamma}{2} t_k^2 - \frac{t_k}{\beta} + \frac{\eta}{\beta} - \frac{\gamma}{2} t_{k - 1}^2 + \frac{t_{k - 1}}{\beta} - \frac{\eta}{\beta} - \left( \gamma t_{k - 1} - \frac{1}{\beta} \right) (t_k - t_{k - 1}) \right\vert\\ &= \left\vert \frac{\gamma}{2} t_k^2 - \frac{t_k}{\beta} - \frac{\gamma}{2} t_{k - 1}^2 + \frac{t_{k - 1}}{\beta} - \gamma t_{k - 1} t_k + \gamma t_{k - 1}^2 + \frac{t_k}{\beta} - \frac{t_{k - 1}}{\beta} \right\vert\\ &= \left\vert \frac{\gamma}{2} t_k^2 - \gamma t_k t_{k - 1} + \frac{\gamma}{2} t_{k - 1}^2 \right\vert\\ &= \frac{\gamma}{2} \left( t_k - t_{k - 1} \right)^2.\end{split}\]

The inductive step for \(n = k\) is

\[\begin{split}x_{k + 1} &= x_k - J(x_k)^{-1} F(x_k)\\ \left\Vert x_{k + 1} - x_k \right\Vert &= \left\Vert -J(x_k)^{-1} F(x_k) \right\Vert\\ &\leq \left\vert f'(t_k)^{-1} \right\vert \left\Vert F(x_k) \right\Vert\\ &= \left\vert t_{k + 1} - t_k \right\vert \frac{\left\Vert F(x_k) \right\Vert}{\left\vert f(t_k) \right\vert} & \quad & \left\vert t_{k + 1} - t_k \right\vert = \left\vert f'(t_k)^{-1} f(t_k)\right\vert\\ &\leq \left\vert t_{k + 1} - t_k \right\vert & \quad & \text{ratio is less than one by the inductive step.}\end{split}\]

Exercise 8

As stated in the section after Theorem 5.3.1, Exercise 6 and Exercise 7 shows that \(\left\Vert x_* - x_k \right\Vert\) is bounded above by \(\left\vert t_* - t_k \right\vert\); hence it is a convergent sequence. Every convergent sequence is a Cauchy sequence, so \(\{ x_k \}\) is a Cauchy sequence that will converge to some \(x_*\). This technique of proof is known as majorization where \(\{ t_k \}\) is said to majorize \(\{ x_k \}\).

See [Ort68] for a very concise and formal proof.

Exercise 9

See [Bry68] for an elegant and concise proof that enables the evaluation of nonlinear integral equations using numerical quadrature.

Exercise 10

The contractive mapping theorem is also known as the Banach fixed-point theorem or the contraction mapping principle.

Given \(G: D \rightarrow D\) where \(D\) is a closed subset of \(\mathbb{R}^n\), for any starting point \(x_0 \in D\), the sequence \(\{ x_k \}\) generated by \(x_{k + 1} = G(x_k)\) remains in \(D\) for \(k = 0, 1, \ldots\).

(a)

Observe that recursively applying

\[\left\Vert x_{k + 1} - x_k \right\Vert = \left\Vert G(x_k) - G(x_{k - 1}) \right\Vert \leq \alpha \left\Vert x_k - x_{k - 1} \right\Vert\]

yields

\[\left\Vert x_{k + 1} - x_k \right\Vert \leq \alpha^k \left\Vert x_1 - x_0 \right\Vert\]

where \(\alpha \in [0,1)\).

For any \(j \geq 0\),

\[\begin{split}\sum_{i = 0}^j \left\Vert x_{i + 1} - x_i \right\Vert &\leq \left\Vert x_1 - x_0 \right\Vert \sum_{i = 0}^j \alpha^i\\ &= \left\Vert x_1 - x_0 \right\Vert \frac{1 - \alpha^{j + 1}}{1 - \alpha} & \quad & \text{geometric series formula}\\ &\leq \frac{\left\Vert x_1 - x_0 \right\Vert}{1 - \alpha} & \quad & \text{as } j \rightarrow \infty.\end{split}\]

(b)

Recall that a sequence \(\{ x_i \}\) is called a Cauchy sequence if \(\forall \epsilon > 0 \, \exists N \in \mathbb{N}\) such that \(\forall \, m, n \geq N\), then \(\left\Vert x_m - x_n \right\Vert < \epsilon\).

Without loss of generality, let \(m \leq n\) with \(k = n - m \geq 0\).

\[\begin{split}\left\Vert x_m - x_n \right\Vert &= \left\Vert x_m - x_{m + k} \right\Vert\\ &= \left\Vert (x_m - x_{m + 1}) + (x_{m + 1} - x_{m + 2}) + \ldots + (x_{m + k - 1} - x_{m + k}) \right\Vert\\ &\leq \left\Vert x_{m + 1} - x_m \right\Vert + \left\Vert x_{m + 2} - x_{m + 1} \right\Vert + \ldots + \left\Vert x_{m + k} - x_{m + k - 1} \right\Vert & \quad & \text{triangle inequality (3.1.1b)}\\ &\leq \left\Vert x_1 - x_0 \right\Vert \sum_{i = m}^{m + k - 1} \alpha^i & \quad & \text{(a)}\\ &= \left\Vert x_1 - x_0 \right\Vert \frac{\alpha^m - \alpha^{m + k}}{1 - \alpha} & \quad & \text{geometric series formula}\\ &\leq \left\Vert x_1 - x_0 \right\Vert \frac{\alpha^N}{1 - \alpha}.\end{split}\]

Hence \(\{ x_i \}\) is a Cauchy sequence as long as \(N\) satisfies

\[\begin{split}\left\Vert x_1 - x_0 \right\Vert \frac{\alpha^N}{1 - \alpha} &< \epsilon\\ \alpha^N &< \epsilon \frac{1 - \alpha}{\left\Vert x_1 - x_0 \right\Vert}\\ N &> \log_{\alpha}\left( \epsilon \frac{1 - \alpha}{\left\Vert x_1 - x_0 \right\Vert} \right).\end{split}\]

The following observations will be useful in proving \(G\) has a fixpoint \(x_*\) where a fixed point (a.k.a. fixpoint, invariant point) of a function is defined as \(G(x_*) = x_*\).

A function is said to be a contraction map if it satisfies the assumptions of the Contractive Mapping Theorem. Notice that the inequality (5.3.2) is a restricted definition of Lipschitz continuity, which means a contraction map is continuous on \(D\).

As shown in (a), \(\{ x_k \}\) converges to some \(x_*\) because \(\lim_{k \rightarrow \infty} \alpha^k = 0\). Therefore,

\[\begin{split}\lim_{k \rightarrow \infty} x_{k + 1} &= \lim_{k \rightarrow \infty} G(x_k)\\ &= G(\lim_{k \rightarrow \infty} x_k) & \quad & \text{G is continuous}\\ x_* &= G(x_*).\end{split}\]

(c)

Proof by contradiction is the easily method to show that \(x_* \in D\) is unique.

Assume that \(\exists \, y_* \in D\) such that \(y_* \neq x_*\).

\[\begin{split}\left\Vert G(x_*) - G(y_*) \right\Vert &\leq \alpha \left\Vert x_* - y_* \right\Vert & \quad & \text{(5.3.2)}\\ \left\Vert x_* - y_* \right\Vert &\leq \alpha \left\Vert x_* - y_* \right\Vert & \quad & \text{(b)}\\ (1 - \alpha) \left\Vert x_* - y_* \right\Vert &\leq 0.\end{split}\]

Since \(\alpha \in [0, 1)\), \(\left\Vert x_* - y_* \right\Vert = 0\), but this contradicts the assumption \(y_* \neq x_*\). Thus a contraction map has a unique fixed point.

By definition

\[\left\Vert x_{k + 1} - x_* \right\Vert = \left\Vert G(x_k) - G(x_*) \right\Vert \leq \alpha \left\Vert x_k - x_* \right\Vert,\]

\(\{ x_k \}\) converges q-linearly to \(x_*\).

By applying the results of (b) with \(m = k\) as \(n \rightarrow \infty\), the error bound can be derived as

\[\begin{split}\lim_{n \rightarrow \infty} \left\Vert x_k - x_n \right\Vert &\leq \lim_{n \rightarrow \infty} \left\Vert x_1 - x_0 \right\Vert \frac{\alpha^k - \alpha^{n}}{1 - \alpha}\\ \left\Vert x_k - x_* \right\Vert &\leq \frac{\eta \alpha^k}{1 - \alpha} & \quad & \left\Vert x_1 - x_0 \right\Vert \leq \eta.\end{split}\]

The last inequality holds because the norm’s continuity enables moving the limit inside.

(d)

As shown in [Pet], the other version of the Contractive Mapping Theorem no longer assumes a self-map.

Let \(G \colon D \rightarrow \mathbb{R}^n\) where \(D\) is a closed subset of \(\mathbb{R}^n\). If for some norm \(\left\Vert \cdot \right\Vert\), there exists \(\alpha \in [0, 1)\) such that

\[\left\Vert G(x) - G(y) \right\Vert \leq \alpha \left\Vert x - y \right\Vert \qquad \forall \, x, y \in D,\]

then the rest of the Contractive Mapping Theorem still holds. The goal is to show the sequence \(\{ x_k \}\) generated by \(x_{k + 1} = G(x_k)\) remains in \(D\) for \(k = 0, 1, \ldots\) with \(x_0 \in D\).

Since \(D\) is a closed subset and the results of (c) defines a closed neighborhood of radius \(\frac{\eta \alpha^k}{1 - \alpha}\) around \(x_*\), as long as each subsequent element in the sequence are within that iteration’s specified neighborhood, the original assumptions of the Contractive Mapping Theorem are satisfied and the sequence will converge to \(x_*\).

Exercise 11

Define \(x_{i + 1} = G(x_i) = x_i - \frac{f(x_i)}{a}\) where \(f(x) = x^2 - 1\) and \(a > 1\). By inspection, the roots of \(f\) are at \(x = \pm 1\) and \(G \colon \mathbb{R} \rightarrow \mathbb{R}\), but \(x_i\) will not converge to the roots \(\forall x_0 \in \mathbb{R}\). The goal is to find \(D \in \mathbb{R}\) such that \(x_0 \in D\) will converge to \(x_* = 1\) and

\[\begin{split}\left\Vert G(x_0) - G(x_*) \right\Vert &= \left\vert x_0 - \frac{f(x_0)}{a} - x_*\right\vert\\ &= \left\vert x_0 - x_* - \frac{x_0^2 - 1}{a} \right\vert\\ &= \left\vert x_0 - x_* \right\vert \left\vert 1 - \frac{x_0 + x_*}{a} \right\vert & \quad & x_* = 1\\ &< \left\vert x_0 - x_* \right\vert.\end{split}\]

In order to establish the last inequality, \(x_0\) needs to satisfy

\[\begin{split}\left\vert 1 - \frac{x_0 + x_*}{a} \right\vert &< 1\\ 1 - \frac{2 (x_0 + x_*)}{a} + \frac{(x_0 + x_*)^2}{a^2} &< 1\\ 0 &< -(x_0 + x_*)^2 + 2a (x_0 + x_*)\\ 0 &< -x_0^2 + (2a - 2x_*) x_0 + 2ax_* - x_*^2.\end{split}\]

Solving the quadratic equation for the \(x_0\) yields

\[x_0 = \frac{ -(2a - 2x_*) \pm \sqrt{(2a - 2x_*)^2 - 4 (-1) (2ax_* - x_*^2)} }{2 (-1)} = a - x_* \pm (-a).\]

Hence the inequality will be satisfied i.e. \(\exists \, \alpha \in [0, 1)\) such that \(\left\Vert G(x_0) - G(x_*) \right\Vert \leq \alpha \left\vert x_0 - x_* \right\vert\) when \(x_0 \in (-x_*, 2a - x_*)\).

To show \(G \colon D \rightarrow D\) holds by principle of induction, assume

\[\left\Vert x_i - x_* \right\Vert = \left\Vert G(x_{i - 1}) - G(x_*) \right\Vert \leq \alpha \left\vert x_{i - 1} - x_* \right\vert.\]

and

\[\begin{split}\left\Vert G(x_i) - G(x_*) \right\Vert &= \left\vert x_i - \frac{f(x_i)}{a} - x_* \right\vert\\ &= \left\vert x_i - x_* - \frac{(x_i + x_*) (x_i - x_*)}{a} \right\vert & \quad & x_* = 1\\ &= \left\vert x_i - x_* \right\vert \left\vert 1 - \frac{(x_i + x_*)}{a} \right\vert\\ &< \left\vert x_i - x_* \right\vert.\end{split}\]

Consequently,

\[\begin{split}\left\vert 1 - \frac{(x_i + x_*)}{a} \right\vert &< 1\\ 1 - \frac{2 (x_i + x_*)}{a} + \frac{(x_i + x_*)^2}{a^2} &< 1\\ 0 &< -(x_i + x_*)^2 + 2a (x_i + x_*)\\ 0 &< -x_i^2 + (2a - 2x_*) x_i + 2ax_* - x_*^2\\ 0 &< -\left( x_{i - 1} - \frac{x_{i - 1}^2 - 1}{a} \right)^2 + (2a - 2x_*) \left( x_{i - 1} - \frac{x_{i - 1}^2 - 1}{a} \right) + 2ax_* - x_*^2\\ 0 &< -\frac{\left( -ax_{i - 1} - ax_* + x_{i - 1}^2 - 1 \right) \left( 2a^2 - ax_{i - 1} - ax_* + x_{i - 1}^2 - 1 \right)}{a^2}\\ 0 &< \frac{ \left( x_{i - 1} + 1 \right) \left( a - x_{i - 1} + 1 \right) \left( 2a^2 - ax_{i - 1} - a + x_{i - 1}^2 - 1 \right)}{a^2} & \quad & x_* = 1.\end{split}\]

Since \(a > 1\), \(\left\Vert G(x_i) - G(x_*) \right\Vert \leq \alpha \left\vert x_i - x_* \right\vert\) when \(x_{i - 1} \in (-1, a + 1)\).

Thus defining \(\delta = \min\left\{ 2, a + 1, 2a - 1 \right\}\) and applying the results of Exercise 10d proves that the sequence \(\{ x_i \}\) converges q-linearly to \(x_* = 1\) with

\[\begin{split}\lim_{i \rightarrow \infty} \left\vert \frac{x_{i + 1} - x_*}{x_i - x_*} \right\vert &= \lim_{i \rightarrow \infty} \left\vert \frac{x_i - \frac{x_i^2 - 1}{a} - x_*}{x_i - x_*} \right\vert\\ &= \lim_{i \rightarrow \infty} \left\vert 1 - \frac{x_i^2 - 1}{a (x_i - x_*)} \right\vert\\ &= \lim_{i \rightarrow \infty} \left\vert 1 - \frac{x_i + 1}{a} \right\vert & \quad & x_* = 1\\ &= \left\vert 1 - \lim_{i \rightarrow \infty} \frac{x_i + 1}{a} \right\vert & \quad & \text{absolute value is continuous}\\ &= \left\vert 1 - \frac{2}{a} \right\vert\\ &> 0 & \quad & \text{unless } a = 2.\end{split}\]

Exercise 12

If \(\lim_{k \rightarrow \infty} h_k = 0\), then the convergence is q-superlinear.

The proof is similar to the one in the text up to applying Lemma 4.2.1.

\[\begin{split}\left\Vert x_k - x_* \right\Vert &\leq \left\Vert A^{-1}_k \right\Vert \left( \left\Vert F(x_*) - F(x_k) - J(x_k) (x_* - x_k) \right\Vert + \left\Vert A_k - J(x_k) \right\Vert \left\Vert x_* - x_k \right\Vert \right)\\ &\leq 2 \beta \left( \frac{\gamma}{2} \left\Vert x_* - x_k \right\Vert^2 + \frac{\gamma \left\vert h_k \right\vert}{2} \left\Vert x_* - x_k \right\Vert \right) & \quad & \text{Lemma 4.2.1}\\ &= c_k \left\Vert x_* - x_k \right\Vert & \quad & c_k = \beta \gamma \left( \left\Vert x_* - x_k \right\Vert + \left\vert h_k \right\vert \right).\end{split}\]

Since \(\lim_{k \rightarrow \infty} x_k = x_*\) and by assumption \(\lim_{k \rightarrow \infty} h_k = 0\), \(\lim_{k \rightarrow \infty} c_k = 0\) i.e. \(\{ x_k \}\) converges q-superlinearly.

Notice that the expressions (5.4.2) and (5.4.3) are equivalent via Lemma 4.1.16:

\[\begin{split}\left\vert h_k \right\vert &\leq c_2 \left\Vert F(x_k) \right\Vert & \quad & \text{(5.4.3)}\\ &= c_2 \left\Vert F(x_k) - F(x_*) \right\Vert & \quad & F(x_*) = 0\\ &\leq c_2 \beta \left\Vert x_k - x_* \right\Vert & \quad & 4.1.16\\ &= c_1 \left\Vert x_k - x_* \right\Vert & \quad & c_1 = c_2 \beta.\end{split}\]

To show q-quadratic convergence, the proof consists of applying the assumption that either (5.4.2) or (5.4.3) holds.

Exercise 13

Given \(x_k = \begin{pmatrix} 10^7 & 10^{-7} \end{pmatrix}^\top\), the analytical solution for \(J(x_k)\) is

\[\begin{split}J(x_k) = \nabla F(x_k)^\top = \begin{bmatrix} 2 x_1 & 0\\ 0 & 2 x_2 \end{bmatrix}.\end{split}\]

Approximating \(J(x_k)\) with (5.4.1) gives

\[\begin{split}A(x_k) = \frac{1}{h_k} \begin{bmatrix} \left( x_1 + h_k \right)^2 - \left( x_1 \right)^2 & \left( x_1 \right)^2 - \left( x_1 \right)^2\\ \left( x_2 \right)^2 - \left( x_2 \right)^2 & \left( x_2 + h_k \right)^2 - \left( x_2 \right)^2 \end{bmatrix} = \begin{bmatrix} 2 x_1 + h_k & 0\\ 0 & 2 x_2 + h_k \end{bmatrix}.\end{split}\]

When \(h_k = 1\), the approximation

\[\begin{split}A(x_k) = \begin{bmatrix} 2 x_1 + 1 & 0\\ 0 & 2 x_2 + 1 \end{bmatrix}\end{split}\]

has a relative error of

\[\frac{ \left\vert 2 x_1 - \left( 2 x_1 + 1 \right) \right\vert }{\left\vert 2 x_1 \right\vert} = \frac{1}{2} 10^{-7} \quad \text{and} \quad \frac{ \left\vert 2 x_2 - \left( 2 x_2 + 1 \right) \right\vert }{\left\vert 2 x_2 \right\vert} = \frac{1}{2} 10^7.\]

When \(h_k = 10^{-14}\), the approximation

\[\begin{split}A(x_k) = \begin{bmatrix} 2 x_1 + 10^{-14} & 0\\ 0 & 2 x_2 + 10^{-14} \end{bmatrix}\end{split}\]

has a relative error of

\[\frac{ \left\vert 2 x_1 - \left( 2 x_1 + 10^{-14} \right) \right\vert }{\left\vert 2 x_1 \right\vert} = \frac{1}{2} 10^{-21} \quad \text{and} \quad \frac{ \left\vert 2 x_2 - \left( 2 x_2 + 10^{-14} \right) \right\vert }{\left\vert 2 x_2 \right\vert} = \frac{1}{2} 10^{-7}.\]

Recall that a computer with \(14\) base-\(10\) digits has a machine epsilon \(\eta \approx 10^{-14}\). The first approximation incorrectly computed \(J(x_k)_{22}\) while the second approximation will experience underflow during the evaluation of \(A(x_k)_{22}\). Setting

\[h_k = \text{sign}(x_k) \sqrt{\eta} \max\{ \vert x_k \vert, \text{typ}x_k \}\]

is a good compromise.

Exercise 14

Let \(p \in \mathbb{R}, p \neq 0\). The step size proposed in (5.4.10) is

\[h_j = \sqrt{\eta} \max\{ |x_j|, \text{typ}x_j \} \text{sign}(x_j)\]

(a)

Let \(\hat{F} = f_1(x) = x^p\) and \(\gamma = f_1''(x) = p (p - 1) x^p\). The optimal step size given in (5.4.13) is

\[h_* = \left( \frac{4 (\eta + \epsilon)}{\gamma} \hat{F} \right)^{1/2} = \left( \frac{4 (\eta + \epsilon)}{p^2 - p} \right)^{1/2} = \frac{2}{\sqrt{p^2 - p}} \sqrt{\eta + \epsilon}.\]

(b)

Let \(\hat{F} = f_2(x) = e^{px}\) and \(\gamma = f_2''(x) = p^2 e^{px}\). The optimal step size given in (5.4.13) is

\[h_* = \left( \frac{4 (\eta + \epsilon)}{\gamma} \hat{F} \right)^{1/2} = \left( \frac{4 (\eta + \epsilon)}{p^2} \right)^{1/2} = \frac{2}{p} \sqrt{\eta + \epsilon}.\]

(c)

Both step sizes are scale independent, which is not ideal.

(d)

The forward finite-difference approximation to \(f_2'(x) = p e^{px}\) is

\[\frac{f_2(x + h) - f_2(x)}{h} = e^{px} \frac{e^{ph} - 1}{h}\]

with a relative error of

\[\frac{ \left\vert p e^{px} - \left( e^{px} \frac{e^{ph} - 1}{h} \right) \right\vert }{ \left\vert p e^{px} \right\vert } = \left\vert 1 - \frac{e^{ph} - 1}{ph} \right\vert.\]

Good step sizes will tend to make \(e^{ph} \rightarrow 1\), so the quantity will either be \(0\) or a very small constant compared to \(e^{px}\). This means when \(x\) is a very large positive number, either overflow occurs or the approximation will be off as described by the relative error.

Exercise 15

Lemma 4.3.1 asserts that a minimizer of

\[f(x) = f(x_1, x_2) = \frac{1}{2} (x_1^2 - x_2)^2 + \frac{1}{2} (1 - x_1)^2\]

needs to satisfy

\[\begin{split}\nabla f(x) = \begin{bmatrix} 2 x_1^3 - 2 x_1 x_2 + x_1 - 1\\ x_2 - x_1^2 \end{bmatrix} = 0.\end{split}\]

\(x = [1, 1]\) is one such minimizer. By definition 4.1.4, the Hessian is

\[\begin{split}\nabla^2 f(x) = \begin{bmatrix} 6 x_1^2 - 2 x_2 + 1 & -2 x_1\\ -2 x_1 & 1 \end{bmatrix} = 0.\end{split}\]

Based on the results of \(f(x_0)\) and \(f(x_1)\), it does seem like a good step.

[4]:
import numpy

def f(v):
    x = v.flat
    return 0.5 * (x[0]**2 - x[1])**2 + 0.5 * (1 - x[0])**2

def fp(v):
    x = v.flat
    return numpy.asmatrix([2 * x[0]**3 - 2 * x[0] * x[1] + x[0] - 1,
                           x[1] - x[0]**2]).T

def fpp(v):
    x = v.flat
    return numpy.asmatrix([[6 * x[0]**2 - 2 * x[1] + 1, -2 * x[0]],
                           [-2 * x[0], 1]])

x_k = numpy.asmatrix([2, 2]).T
for i in range(2):
    print('x_{0}: {1}'.format(i, f(x_k)))
    _ = x_k - numpy.linalg.inv(fpp(x_k)) * fp(x_k)
    x_k = _
x_0: 2.5
x_1: 0.3208000000000001

Exercise 16

Given \(f(x) = \sin(x)\), \(f'(x) = \cos(x)\) and \(f''(x) = -\sin(x)\).

Applying Newton’s method yields

\[x_1 = x_0 - \frac{\cos(x_0)}{-sin(x_0)} = x_0 + \frac{1}{\tan(x_0)}.\]

Define \(0 < \epsilon \ll 1\). When \(x_0 = -\epsilon\),

\[\begin{split}x_1 &= -\epsilon + \frac{1}{\tan(-\epsilon)}\\ &\cong -\epsilon - \frac{1}{\epsilon} & \quad & \text{small-angle approximation of tangent}\\ &= -\frac{\epsilon^2 + 1}{\epsilon}\\ &\cong -\frac{1}{\epsilon}.\end{split}\]

When \(x_0 = \epsilon\),

\[\begin{split}x_1 &= \epsilon - \frac{1}{\tan(\epsilon)}\\ &\cong \epsilon - \frac{1}{\epsilon} & \quad & \text{small-angle approximation of tangent}\\ &= \frac{\epsilon^2 - 1}{\epsilon}\\ &\cong -\frac{1}{\epsilon}.\end{split}\]

Exercise 17

(a)

The gradient and Hessian of \(f_1(x) = -x_1^2 - x_2^2\) are respectively

\[\begin{split}\nabla f_1(x) = \begin{bmatrix} -2 x_1\\ -2 x_2 \end{bmatrix} \quad \text{and} \quad \nabla^2 f_1(x) = \begin{bmatrix} -2 & 0\\ 0 & -2 \end{bmatrix}.\end{split}\]

Since \(\nabla^2 f_1(x)\) has an eigenvalue \(\lambda = -2\) of multiplicity \(2\), it is not positive definite. This means

\[H_k = \nabla^2 f_1(x) + \mu_k I = (\mu_k - 2) I \quad \Rightarrow \quad H_k^{-1} = \frac{1}{\mu_k - 2} I\]

and

\[\begin{split}x_{k + 1} = x_k - H_k^{-1}(x_k) \nabla f_1(x_k) = \begin{bmatrix} x_{k_1} + x_{k_1} \left( \frac{\mu_k}{2} - 1 \right)^{-1}\\ x_{k_2} + x_{k_2} \left( \frac{\mu_k}{2} - 1 \right)^{-1} \end{bmatrix}.\end{split}\]

From the assumption of algorithm A5.5.1, the denominator is positive. Thus the algorithm will grow towards to \(\pm \infty\) depending on \(\text{sgn}(x_0)\) unless \(x_0 = x_*\).

(b)

The gradient and Hessian of \(f_2(x) = x_1^2 - x_2^2\) are respectively

\[\begin{split}\nabla f_2(x) = \begin{bmatrix} 2 x_1\\ -2 x_2 \end{bmatrix} \quad \text{and} \quad \nabla^2 f_2(x) = \begin{bmatrix} 2 & 0\\ 0 & -2 \end{bmatrix}.\end{split}\]

Since the eigenvalues of \(\nabla^2 f_2(x)\) are \(\lambda_0 = 2\) and \(\lambda_1 = -2\), it is not positive definite. This means

\[\begin{split}H_k = \nabla^2 f_2(x) + \mu_k I \quad \Rightarrow \quad H_k^{-1} = \begin{bmatrix} \frac{1}{\mu_k + 2} & 0\\ 0 & \frac{1}{\mu_k -2} \end{bmatrix}\end{split}\]

and

\[\begin{split}x_{k + 1} = x_k - H_k^{-1}(x_k) \nabla f_1(x_k) = \begin{bmatrix} x_{k_1} - x_{k_1} \left( \frac{\mu_k}{2} + 1 \right)^{-1}\\ x_{k_2} + x_{k_2} \left( \frac{\mu_k}{2} - 1 \right)^{-1} \end{bmatrix}\end{split}\]

From the assumption of algorithm A5.5.1, the denominator is positive. Thus the algorithm will grow towards to \(\pm \infty\) depending on \(\text{sgn}(x_0)\) unless \(x_{0_2} = 0\).

Exercise 18

[5]:
import numpy

def f(x):
    return 0.5 * (x[0]**2 - x[1])**2 + 0.5 * (1 - x[0])**2

def fp(v):
    x = v.flat
    return numpy.asmatrix([2 * x[0]**3 - 2 * x[0] * x[1] + x[0] - 1,
                           x[1] - x[0]**2]).T

def fpp(v):
    x = v.flat
    return numpy.asmatrix([[6 * x[0]**2 - 2 * x[1] + 1, -2 * x[0]],
                           [-2 * x[0], 1]])

def afpp(v, root_eta):
    x_c = v.flat
    h = root_eta * x_c
    e = [numpy.asfarray([1, 0]), numpy.asfarray([0, 1])]

    A = numpy.asmatrix([[0.0, 0.0], [0.0, 0.0]])
    for i in range(2):
        for j in range(2):
            A[i,j] = f(x_c + h[i] * e[i] + h[j] * e[j])
            A[i,j] -= f(x_c + h[i] * e[i])
            A[i,j] -= f(x_c + h[j] * e[j])
            A[i,j] += f(x_c)
            A[i,j] /= h[i] * h[j]
    return A

x_k = numpy.asmatrix([2.0, 2.0]).T

solutions = [[], [], []]

for s in range(3):
    x_k = numpy.asmatrix([2.0, 2.0]).T
    for i in range(8):
        if 0 == s:
            A = fpp(x_k)
        elif 1 == s:
            A = afpp(x_k, numpy.finfo(float).eps**(1/2))
        else:
            A = afpp(x_k, numpy.finfo(float).eps**(1/3))
        _ = x_k - numpy.linalg.inv(A) * fp(x_k)
        x_k = _
        solutions[s].append(f(x_k).flat[0])

row_layout = '{0:<25} {1:<25} {2:<25}'
print(row_layout.format('Analytic', 'macheps^(1/2)', 'macheps^(1/3)'))
for i in range(len(solutions[0])):
    print(row_layout.format(solutions[0][i], solutions[1][i], solutions[2][i]))
Analytic                  macheps^(1/2)             macheps^(1/3)
0.3208000000000001        0.3208000000000001        0.32080752264304735
0.1522899437566909        0.1203483182966636        0.1522190630342569
0.00048098909417033226    0.012308840948738988      0.00048274254706565654
4.6037027257078434e-07    0.0005726840749779877     4.6389527426665865e-07
4.471929271276434e-15     6.465667087177092e-07     5.707975194043212e-15
3.944304526105059e-29     9.628815678752908e-13     1.369421165569959e-23
0.0                       2.77706417041588e-24      2.465190328815662e-32
0.0                       0.0                       0.0

Exercise 19

Recall that the Hessian of \(f \colon \mathbb{R}^n \rightarrow \mathbb{R}\) is equal to the Jacobian of \(\nabla f\).

\[\begin{split}x_+ - x_* &= x_c - A_c^{-1} g_c - x_*\\ &= A_c^{-1} [ A_c (x_c - x_*) - g_c ]\\ &= A_c^{-1} \left[ \left( A_c - \nabla^2 f(x_c) \right) (x_c - x_*) - g_c + \nabla^2 f(x_c) (x_c - x_*) \right]\\ &= A_c^{-1} \left[ \left( A_c - \nabla^2 f(x_c) \right) (x_c - x_*) - (g_c - \nabla f(x_c)) + \nabla f(x_*) - \nabla f(x_c) - \nabla^2 f(x_c) (x_* - x_c) \right]\\ \left\Vert x_+ - x_* \right\Vert &= \left\Vert A_c^{-1} \left[ \left( A_c - \nabla^2 f(x_c) \right) (x_c - x_*) - (g_c - \nabla f(x_c)) + \nabla f(x_*) - \nabla f(x_c) - \nabla^2 f(x_c) (x_* - x_c) \right] \right\Vert\\ \left\Vert x_+ - x_* \right\Vert &\leq \left\Vert A_c^{-1} \right\Vert \left[ \left\Vert A_c - \nabla^2 f(x_c) \right\Vert \left\Vert x_c - x_* \right\Vert + \left\Vert -(g_c - \nabla f(x_c)) \right\Vert + \left\Vert \nabla f(x_*) - \nabla f(x_c) - \nabla^2 f(x_c) (x_* - x_c) \right\Vert \right]\\ \left\Vert x_+ - x_* \right\Vert &\leq \left\Vert A_c^{-1} \right\Vert \left[ \left\Vert g_c - \nabla f(x_c) \right\Vert + \left\Vert A_c - \nabla^2 f(x_c) \right\Vert \left\Vert e_c \right\Vert + \frac{\gamma}{2} \left\Vert e_c^2 \right\Vert \right] & \quad & \text{Lemma 4.1.12.}\end{split}\]

Exercise 20

A5.6.3, A5.6.4, and A5.6.2 uses \(n + 1\), \(2n\), and \(1 + 2n + n^2 / 2\) function evaluations respectively. The overlap between each of these algorithms is at most \(n\) function evaluations. Therefore,

\[\text{A5.6.3} + \text{A5.6.2} = (n + 1) + (1 + 2n + n^2 / 2) - n = 2 + 2n + n^2 / 2\]

and

\[\text{A5.6.4} + \text{A5.6.2} = (2n) + (1 + 2n + n^2 / 2) - n = 1 + 3n + n^2 / 2.\]

The combined algorithm may yield answers with less accuracy and possibly lose cache locality due to more complex code.

References

Bry68

Charles A Bryan. Approximate solutions to nonlinear integral equations. SIAM Journal on Numerical Analysis, 5(1):151–155, 1968.

Ort68

James M Ortega. The newton-kantorovich theorem. The American Mathematical Monthly, 75(6):658–660, 1968.

Pet

Tobias von Petersdorf. Fixed point iteration and contraction mapping theorem. http://terpconnect.umd.edu/ petersd/466/fixedpoint.pdf. Accessed on 2017-04-10.