\section{Newton's Method}\label{S:NewtonsMethod}

\subsection{Shortcomings of the Steepest Descent Method}

\noi A major drawback of the Steepest Descent Method is
its slow rate of convergence near the minimum. The simple positive definite
quadratic function, $\ds f(x_1,x_2)=  x_1^2+x_1x_2+\frac{3}{2}x_2^2+x_1-4x_2+6$, from
(\ref{E:SteepestDescentBasicPositiveDefiniteQuadratic})
illustrates this phenomenon, both
in Table \ref{Ta:SteepestDescent_steepest_descent_example}, where the relative change in
$\ds \left\|  \nabla f\left(\vecx_k\right)\right\|$ is largest for the first few
iterations, but also in
Figure \ref{Fi:unconstrained_numeric_quad_form_steepest_descent_vectors}. Superimposed on this diagram are four contours of $f$, the first three Steepest Descent Method iterates, and scaled versions of the corresponding optimal descent directions,
$-\nabla f(\vecx_k)$, where $k=0,1,2,3$. The descent direction vectors demonstrate a drawback with the Steepest Descent Method known
as ``zig-zagging'' or ``hemstitching\index{hemstiching}.'' This phenomenon, which frequently results for functions having mildly ``elliptic'' or
``elongated'' valleys, slows the rate of convergence.

\begin{figure}[h]
\centering
\psfrag{x1_axis_label}[Bc][Bc]{$x_1$}
\psfrag{x2_axis_label}[Bc][Bc]{$x_2$}
\psfrag{k0}[Bc][Bc]{$\vecx_0$}
\psfrag{k1}[Bc][Bc]{$\vecx_1$}
\psfrag{k2}[Bc][Bc]{$\vecx_2$}
\psfrag{k3}[Bc][Bc]{$\vecx_3$}
\psfrag{z0}[Bc][Bc]{$f(\vecx)=12.5$}
\psfrag{z1}[Bc][Bc]{$f(\vecx)=8$}
\psfrag{z2}[Bc][Bc]{$f(\vecx)=4$}
\psfrag{z3}[Bc][Bc]{$f(\vecx)=2$}
\includegraphics[width=.7\linewidth]
{images/f08x02x01.eps}
\caption{Steepest descent approximations $x_k$, $k=0,1,2,3$ for
\newline $\ds f(x_1,x_2)=x_1^2+x_1x_2+\frac{3}{2}x_2^2+x_1-4x_2+6$, along
with scaled optimal descent directions.}
\label{Fi:unconstrained_numeric_quad_form_steepest_descent_vectors}
\end{figure}



\subsection{Method Derivation}

\noi \emph{Newton's Method}\index{Newton's Method!derivation} is well-suited for addressing this problem,
although it does require that the objective function, $f$, be twice-differentiable, and it has its
own computational drawbacks. Suppose $S \subseteq \R^n$, $f:S \to \R$ is twice-differentiable, and $\vecx_0$ belongs
to the interior of $S$. For all $\| \vecd \|$ sufficiently small,
$\vecx_0+\vecd$ belongs to $S$ and second-order differentiability implies

\begin{equation}\label{E:unconstrained_numeric_Newton_quad_approx}
f(\vecx_0) - f(\vecx_0+\vecd) \approx -\nabla f(\vecx_0)^t \vecd - \frac{1}{2}\vecd^t H_f(\vecx_0)\vecd.
\end{equation}
If we intend to express the first iterate,
$\vecx_1$, in the form $\vecx_1=\vecx_0+\vecd$, then
the decrease in $f$ is maximized by differentiating
the right-hand side of (\ref{E:unconstrained_numeric_Newton_quad_approx})
with respect to the {\underbar{vector}} $\vecd$.  Doing so, setting the result equal to $\zerovec$ and solving for $\vecd$, we obtain
\begin{equation}\label{E:unconstrained_numeric_Newton_direction}
\vecd=-H_f(\vecx_0)^{-1}\nabla f(\vecx_0),
\end{equation}
provided $\ds H_f(\vecx_0)^{-1}$ exists.

\newpar
We call the vector $\vecd$ from (\ref{E:unconstrained_numeric_Newton_direction}) the \emph{Newton direction}\index{Newton's Method!Newton direction}
of $f$ at $\vecx_0$. It is the second-order analog of the optimal descent direction associated
with the Steepest Descent Method. Like the optimal descent direction, it produces a curve, ``embedded''
in the graph of $f$. Figure \ref{Fi:f08x02x02} depicts both the optimal descent and Newton directions starting from
$\ds \vecx_0=\begin{bmatrix} 1 \\ -1 \end{bmatrix}$ for the quadratic function,
$\ds f(x_1,x_2)=  x_1^2+x_1x_2+\frac{3}{2}x_2^2+x_1-4x_2+6$.


\begin{figure}[h]
\centering
\psfrag{x0}[Bc][Bc]{$\ds \vecx_0=\begin{bmatrix} 1 \\ -1 \end{bmatrix}$}
\psfrag{x1axislabel}[Bc][Bc]{$x_1$}
\psfrag{x2axislabel}[Bc][Bc]{$x_2$}
\psfrag{dSDM}[Bc][Bc]{Optimal Descent Direction}
\psfrag{dNM}[Bc][Bc]{Newton Direction}
\includegraphics[width=.7\linewidth]{images/f08x02x02.eps}
\caption{Illustration of optimal descent and Newton directions for
$\ds f(x_1,x_2)=x_1^2+x_1x_2+\frac{3}{2}x_2^2+x_1-4x_2+6$ starting at $\vecx_0$.}
\label{Fi:f08x02x02}
\end{figure}

\newpar
Thus, starting from $\vecx_0$
and moving in the Newton direction leads to
\newline $\ds \vecx_1=\vecx_0+\vecd=\vecx_0-H_f(\vecx_0)^{-1}\nabla f(\vecx_0)$.
Of course, assuming invertibility of the Hessian at each stage, this process can be repeated, thereby
leading to Newton's Method.
We now summarize the steps of the algorithm.

\newpar
{\underbar{Newton's Method}}\index{Newton's Method!summary of steps} \newline
To obtain an approximate solution of NLP (\ref{SteepestDescentStandardNLP}) under
the assumption $f$ is twice-differentiable on $S$:
\be
\item From the initial value $\vecx_0$ in $S$, calculate the Newton direction, $\vecd=-H_f(\vecx_0)^{-1}\nabla f(\vecx_0)$,
provided $\ds H_f(\vecx_0)^{-1}$ exists.
\item Define $\ds \vecx_1=\vecx_0+\vecd=\vecx_0-H_f(\vecx_0)^{-1}\nabla f(\vecx_0)$.
\item Proceed to (1) and replace $\vecx_0$ with
$\vecx_1$. Compute the Newton direction of $f$ at $\vecx_1$, $\vecx_2$, and so on.
\ee

\newpar
For example, recall the positive definite quadratic function from
(\ref{E:SteepestDescentBasicPositiveDefiniteQuadratic}). If
\newline $\ds A=\begin{bmatrix} 2 & 1 \\ 1 & 3 \end{bmatrix}$,
$\ds \vecb=\begin{bmatrix} 1 \\ -4 \end{bmatrix}$, and $c=6$,
this function is given by
\begin{align}\label{E:unconstrained_numeric_steepest_descent_quadratic_form_example}
f(\vecx)&=\frac{1}{2}\vecx^t A \vecx -\vecb^t\vecx +c  \notag \\
&= x_1^2+x_1x_2+\frac{3}{2}x_2^2+x_1-4x_2+6 \notag
\end{align}
If $\ds \vecx_0=\begin{bmatrix} 1 \\ -1 \end{bmatrix}$, then
$\ds \nabla f(\vecx_0)=\begin{bmatrix} 2 \\ -6 \end{bmatrix}$ and
$\ds H_f(\vecx_0)= \begin{bmatrix} 2 & 1 \\ 1 & 3 \end{bmatrix}$. Thus,
the Newton direction at $\vecx_0$ becomes
$\vecd=\begin{bmatrix} -2.4 \\ 2.8 \end{bmatrix}$
and $\ds \vecx_1=\vecx_0+\vecd= \begin{bmatrix} -1.4 \\ 1.8 \end{bmatrix}$.
Recall from Section \ref{S:SteepestDescent} that $\vecx_1$ is precisely the global minimum of $f$!
That the first iterate of Newton's Method leads to a global minimum is true in
general for all positive definite quadratic mappings. Verification of
this fact is fairly straightforward.

\begin{progress_ck}
Suppose that $f:\R^n \to \R$ is the positive definite quadratic function defined by
\begin{equation*}
f(\vecx)=\frac{1}{2}\vecx^tA\vecx-\vecb^t\vecx+c,
\end{equation*}
where $A$ is an $n$-by-$n$ positive definite matrix, $\vecb$ is in $\R^n$, and $c$ is a real scalar.
\be
\item Compute $\nabla f(\vecx)$, and use your result to determine the global minimum $\vecx_{\bigstar}$ in terms
of $A$ and $\vecb$.
\item Given $\vecx_0$ in $\R^n$, determine the Newton direction in terms of $\vecx_0$, $A$, and $\vecb$.
\item Show that the first iterate, $\vecx_1$, obtained using Newton's method is the global minimum obtained in (1).
\ee
\end{progress_ck}

\subsection{A Maple Implementation of Newton's Method}\label{SS:NewtonsMethodMapleImplementation}\index{Newton's Method!Maple implementation}

\noi The \mf{SteepestDescentMethod} procedure outlined in Section
\ref{S:SteepestDescent} is easily modified to create a new procedure,
\mf{NewtonsMethod}, which may be found in the worksheet \textbf{Newtons Method.mw}. The
procedure call takes the form
\begin{equation*}
\mbox{\mf{NewtonsMethod(\ma{function,initial,N,tolerance})},}
\end{equation*}
where \ma{function} is a twice-differentiable function to be minimized, \ma{initial} is the initial value, $\ds \vecx_0$, (written as a
list), \ma{N} is the maximum number of iterations, and \ma{tolerance} is the specified tolerance. The algorithm terminates when either the number of iterations is reached or the gradient norm is less than the desired tolerance. The returned value is the approximate minimum value. Intermediate
iterates leading to the approximate minimum are stored in a list, \mf{iterates}, a global variable accessible outside of the procedure.

\newpar
The following worksheet constructs the procedure and approximates the minimum of the \CM objective,
\newline
$\ds f(x_1,x_2 )= -1400 x_1^{\frac{1}{2}}x_2^{\frac{1}{3}} + 350x_1 +200x_2$
using $\ds \vecx_0=\begin{bmatrix} 10 \\ 10 \end{bmatrix}$, $N=5$, and a tolerance
of $\epsilon=.01$.

\begin{maplelist}
\item \mf{restart: with(VectorCalculus):with(LinearAlgebra):}
\item \mf{NewtonsMethod := proc (function, initial, N, tolerance) \newline
local x,f,Delf,Hf,d,j,epsilon; global iterates: \newline
\# Create local and global variables. \newline
x:=array(0 .. N):
\newline \# Define array of iterates. \newline
f:=unapply(function, [x1, x2]): \newline
 \# Create function. \newline
Delf:=unapply(Gradient(function, [x1, x2]), [x1, x2]): \newline
\# Create corresponding gradient function. \newline
Hf:=unapply(Hessian(function, [x1, x2]), [x1, x2]): \newline
\# Create corresponding Hessian matrix function. \newline
x[0]:=evalf(initial):
\newline \# Set initial array value.\newline
epsilon := 10:
\newline \# Set initial epsilon to a large value. \newline
j:=0: while (tolerance <= epsilon and j<=N-1) do \newline
\# Create loop structure to perform algorithm. \newline
d:=-convert(
\newline MatrixInverse(Hf(x[j][1],x[j][2])).Delf(x[j][1],x[j][2]),
\newline list): \newline
 \# Compute Newton direction. \newline
x[j+1]:=[x[j][1]+d[1],x[j][2]+d[2]]: \newline
\# Compute next iterate. \newline
epsilon:=evalf(Norm(Delf(x[j+1][1],x[j+1][2]),Euclidean)):\newline
\# Update epsilon using gradient.\newline
j:=j+1: \newline
\# Increase loop index. \newline
end do: \newline
iterates:=[seq(x[i],i=0..j)]:RETURN(x[j]):end: \newline
\# Construct iterates and return last array value.}
\item \mf{f:=(x1,x2)->-1400x1\^{}(1/2)x2\^{}(1/3)+350x1+200x2;
\newline \# Enter function}
\mo{f:=(x_1,x_2) \rightarrow -1400x_1^{\frac{1}{2}}x_2^{\frac{1}{3}}
+350x_1+200x_2}
\item \mf{NewtonsMethod(f(x1, x2), [10,10], 5, .01);\newline
\# Procedure call using initial value [10,10], bound of 
\newline 5 iterations, and tolerance, .01.}
\mo{ [87.10975551,101.6262013]}
\item \mf{iterates; \newline \# Print iterates.}
\begin{align*}
&[[10., 10.], [28.06287805, 29.11130610], [56.81932208, 61.72021291],\\
 &[80.07576281, 91.03969995], [86.73692387, 100.9216310], \\
 [87.10975551, 101.6262013]]
\end{align*}
\end{maplelist}

\subsection{Convergence Issues and Comparison with the Steepest Descent Method}\label{SS:NewtonsMethodConvergenceIssues}\index{Newton's Method!convergence issues}

\noi The preceding examples suggest that Newton's Method is more effective for approximating a local or
global minimum of a function than is the Steepest Descent Method. However, this result
is not true in general, and Newton's Method has its own computation drawbacks. Among these
are the following:
\be
\item Newton's Method assumes second-order differentiability of the objective function, $f$, and requires
evaluating the Hessian matrix, $H_f(\vecx_k)$, at each iteration.
\item Whereas the Steepest Descent Method requires solving an equation in a single variable at each iteration, Newton's Method requires solving a system of equations. This task frequently requires
    a much greater number of computations.
\item To obtain a unique solution to this system, we require that $H_f(\vecx_k)$ be invertible, which may not be the case. As a result, the method fails altogether. However,
even if $H_f(\vecx_0)$ is invertible, it may be ``almost non-invertible,'' as measured by having
eigenvalues extremely small in magnitude. Numeric results then become very sensitive
to computational roundoff error, a general phenomenon known as \emph{ill-conditioning}\index{ill-conditioning}.
\ee

\newpar
While Newton's Method has its drawbacks, the advantage of its use stems from the fact that
an initial value sufficiently close to a known local minimum yields a sequence of iterates converging to the minimum more rapidly to the minimum than a sequence obtained using the Steepest Descent Method. Theorem \ref{Th:NewtonsMethodConvergenceRateTheorem} summarizes this result.

\newpar
The statement of this theorem requires
us to recall the notion of a \emph{matrix norm}\index{matrix norm} as
summarized in Appendix \ref{C:Linear_Algebra_Appendix}. Any matrix norm, $\|\cdot\|$, on $\R^n$ satisfies
$\|A\vecx\| \leq \|A\| \|\vecx\|$ for every $\vecx$ in $\R^n$. As our discussion focuses almost exclusively on the
case when $A$ is symmetric, we will use, $\ds \|A\|=\rho(A)$, where $\rho(A)$ is the spectral radius, or maximum of the absolute values of the eigenvalues of $A$.

\newpar
Because the Hessian matrix, $H_f(\vecx)$, is a function of $\vecx$, so too will
be its matrix norm. We say that $H_f(\vecx)$ is Lipschitz continuous\index{Lipschitz continuous} in the matrix norm if
 there exists a constant, $M$, satisfying
 \begin{equation}\label{E:NewtonsMethodHessianLipschitz1}
\left\| H_f(\vecx_1)-H_f(\vecx_2) \right\| \leq M\|\vecx_1-\vecx_2 \|.
\end{equation}
Phrased another way, the norm of the matrix difference, $\ds H_f(\vecx_1)-H_f(\vecx_2)$,
is bounded by a constant times the distance from $\vecx_2$ to $\vecx_1$.
A consequence of (\ref{E:NewtonsMethodHessianLipschitz1}), whose proof we omit, is that
if $H_f(\vecx_1)$ is positive
definite and if $H_f(\vecx)$ is Lipschitz continuous in some neighborhood of $\vecx_1$, then $H_f(\vecx_2)$ is also positive definite and bounded below in norm
by some constant, $m$, independent of $\vecx_2$, provided $\ds \|\vecx_1-\vecx_2 \|$
is sufficiently small.

\newpar
With this background in mind, we now state an important convergence result.

\begin{theorem}\label{Th:NewtonsMethodConvergenceRateTheorem}
Assume $\ds S \subseteq \R^n$ is nonempty, open set. Suppose
 $f: S \to \R^n$ is continuously twice-differentiable, having a local
 minimum at $\vecx_{\bigstar}$, where $H_f(\vecx_{\bigstar})$
 is positive definite. Suppose that for $\vecx$ sufficiently
 close to $\vecx_{\bigstar}$,
 $H_f(\vecx)$ is Lipschitz continuous in the matrix norm, meaning
 there exists a constant, $M$, satisfying
 \begin{equation}\label{E:NewtonsMethodHessianLipschitz1}
\left\| H_f(\vecx_1)-H_f(\vecx_2) \right\| \leq M\|\vecx_1-\vecx_2 \|,
\end{equation}
for all $\vecx_1$ and $\vecx_2$ sufficiently close to $\vecx_{\bigstar}$.
Then, if $\vecx_0$ is sufficiently close to $\vecx_{\bigstar}$, the
sequence of Newton Method iterates satisfies
\begin{equation}\label{E:NewtonsMethodConvergenceRateInequality}
\left\| \vecx_{k+1}-\vecx_{\bigstar}\right\|
\leq C \left\|\vecx_k-\vecx_{\bigstar}\right\|^2 \mbox{ for $k=0,1,2,\ldots$},
\end{equation}
for some constant, $C$.
\end{theorem}
\begin{proof}
Choose $\delta$ in $(0,1)$ small enough to satisfy the following three conditions:
\be
\item  Inequality (\ref{E:NewtonsMethodHessianLipschitz1})
holds whenever $\ds \left\| \vecx_1 -\vecx_{\bigstar}\right\|$ and
$\ds \left\| \vecx_2 -\vecx_{\bigstar}\right\|$ are less than $\delta$.
\item $\ds H_f(\vecx)$ is positive definite if
$\ds \|\vecx-\vecx_{\bigstar} \|<\delta$
\item $\ds 0<m \leq \| H_f(\vecx)\|$ if
$\ds \|\vecx-\vecx_{\bigstar} \| < \delta $
\ee

\noi
Suppose $\vecx_0$ is chosen close enough to $\vecx_{\bigstar}$ so that
\begin{equation}\label{E:NewtonsMethodConvergenceRateEQ0}
\|\vecx_0-\vecx_{\bigstar}\| < \min \left\{\delta ,\sqrt{\frac{2m}{M}}\cdot \delta \right\}.
\end{equation}.

\newpar
By (2), $\ds H_f(\vecx) $ is invertible
so that
$\ds \vecx_1 = \vecx_0-H_f (\vecx_0)^{-1} \nabla f(\vecx_0)$
is defined. Since
 $\ds \nabla f(\vecx_{\bigstar})=\zerovec$, we have

\begin{align}\notag
\vecx_{1}-\vecx_{\bigstar}
&=\vecx_{0}- H_f(\vecx_0)^{-1}\nabla f(\vecx_0) -\vecx_{\bigstar}  \\
&=\vecx_0-\vecx_{\bigstar}+H_f(\vecx_0)^{-1}\left(\nabla f(\vecx_{\bigstar})-\nabla f(\vecx_0)\right).
\label{E:NewtonsMethodConvergenceRateProofEQ1}
\end{align}

\newpar
We now focus on the difference $\ds \left(\nabla f(\vecx_{\bigstar})-\nabla f(\vecx_0 )\right)$
in (\ref{E:NewtonsMethodConvergenceRateProofEQ1}) and introduce a vector-valued function
$\ds \phi: [0,1] \to \R^n$ given by
\begin{equation*}
\phi(t)=\nabla f \left(t\vecx_{\bigstar}+(1-t)\vecx_0 \right).
\end{equation*}
Note that $\ds \phi(0)=\nabla f(\vecx_0)$ and $\ds \phi(1)=\nabla f(\vecx_{\bigstar})$. An application
of the chain rule yields
\begin{equation*}
\phi'(t)= H_f\left(t\vecx_{\bigstar}+(1-t)\vecx_0 \right)\left(\vecx_{\bigstar}-\vecx_0\right),
\end{equation*}
which is continuously differentiable on $[0,1]$ since $f$ is continuously twice-differentiable on $S$.
By the Fundamental Theorem of Calculus,
\begin{equation}\label{E:NewtonsMethodConvergenceRateProofEQ2}
\nabla f(\vecx_{\bigstar})-\nabla f(\vecx_0)
=\int_{t=0}^{t=1} H_f\left(t\vecx_{\bigstar}+(1-t)\vecx_0 \right)\left(\vecx_{\bigstar}-\vecx_0\right) \,dt.
\end{equation}

\newpar
If we substitute
(\ref{E:NewtonsMethodConvergenceRateProofEQ2}) into (\ref{E:NewtonsMethodConvergenceRateProofEQ1}),
we obtain
\begin{align}\label{E:NewtonsMethodConvergenceRateProofEQ3}
\vecx_{1}-\vecx_{\bigstar} &=
\vecx_0-\vecx_{\bigstar}+H_f(\vecx_0)^{-1}\left(\nabla f(\vecx_{\bigstar})-\nabla f(\vecx_0)\right) \\
 &=H_f(\vecx_0)^{-1}
\left(
 H_f (\vecx_0 )\left(\vecx_0-\vecx_{\bigstar}\right)+
 \left(
 \int_{t=0}^{t=1} H_f\left(t\vecx_{\bigstar}+(1-t)\vecx_0 \right)\left(\vecx_{\bigstar}-\vecx_0\right) \,dt
 \right)
 \right) \notag \\
 &=H_f(\vecx_0)^{-1} \int_{t=0}^{t=1}
 \left( H_f\left(t\vecx_{\bigstar}+(1-t)\vecx_0 \right)-H_f (\vecx_0 )\right)
 \left(\vecx_{\bigstar}-\vecx_0\right) \,dt. \notag
 \end{align}

\newpar
We now apply the result from (\ref{E:NewtonsMethodConvergenceRateProofEQ3})
to bound the norm of the vector difference,
$\vecx_{1}-\vecx_{\bigstar}$. In this process we utilize the fact that the eigenvalues
of $H_f(\vecx_0)^{-1}$ are the reciprocals of the eigenvalues of $H_f(\vecx_0)$.
Hence, $\ds \left\| H_f(\vecx_0)^{-1} \right\| \leq m$.
Using this bound, along with properties of the definite integral, and the Lipschitz continuity of $H_f(\vecx)$, we obtain the following sequence of inequalities:

\begin{align}\label{E:NewtonsMethodConvergenceRateProofEQ4}
\left\|
\vecx_{1}-\vecx_{\bigstar}
\right\|
&=\left\|
H_f(\vecx_0)^{-1} \int_{t=0}^{t=1}
 \left(
 H_f\left(t\vecx_{\bigstar}+(1-t)\vecx_0 \right)
 -H_f (\vecx_0 )
 \right)
 \left(\vecx_{\bigstar}-\vecx_0 \right) \,dt
\right\| \notag \\
&\leq  \left\|
 H_f(\vecx_0)^{-1}
 \right\|
 \int_{t=0}^{t=1}
 \left\|
 H_f\left(t\vecx_{\bigstar}+(1-t)\vecx_0 \right)-H_f (\vecx_0 )
 \right\|
 \left\| \vecx_{\bigstar}-\vecx_0\right\|
 \,dt \notag \\
 &\leq  \left\|
 H_f(\vecx_0)^{-1}
 \right\|
 \int_{t=0}^{t=1}
 M \left\| t\vecx_{\bigstar} + (1-t)\vecx_0 - \vecx_0 \right\|
 \left\| \vecx_{\bigstar}-\vecx_0\right\|
 \,dt \notag \\
 &= \left\|
 H_f(\vecx_0)^{-1}
 \right\|
\int_{t=0}^{t=1} M  \left\| \vecx_{\bigstar}-\vecx_0 \right\|^2 t \, dt \notag \\
&= \frac{M}{2m} \left\| \vecx_{\bigstar}-\vecx_0 \right\|^2. \notag
\end{align}
Hence,
\begin{equation*}
\left\| \vecx_{1}-\vecx_{\bigstar}\right\|
\leq \frac{M}{2m}\left\|\vecx_0-\vecx_{\bigstar}\right\|^2.
\end{equation*}
Since $\ds \|\vecx_0-\vecx_{\bigstar}\| < \sqrt{\frac{2m}{M}}\cdot \delta$, we see that
$\ds \left\| \vecx_{1}-\vecx_{\bigstar}\right\| \leq \delta^2 < \delta$, so the preceding argument
can be applied again, this time starting at $\vecx_1$. Repeating the results,
we eventually obtain inequality (\ref{E:NewtonsMethodConvergenceRateInequality})
with $\ds C=\frac{M}{2m}$, which completes the proof.
\end{proof}

\newpar
Because
inequality (\ref{E:NewtonsMethodConvergenceRateInequality}) implies that
\begin{equation*}
\frac{\left\| \vecx_{k+1}-\vecx_{\bigstar}\right\|}{\left\|\vecx_k-\vecx_{\bigstar}\right\|^2}
\leq C,
\end{equation*}
Theorem \ref{Th:NewtonsMethodConvergenceRateTheorem} can be rephrased as saying the sequence of
iterates, starting at $\vecx_0$, converges to $\vecx_{\bigstar}$ at a
\emph{quadratic rate}\index{Newton's Method!quadratic rate of convergence}. Because
the denominator in the ratio is squared, this type of convergence is more rapid
than linear convergence, as defined following the statement of Theorem \ref{Th:SteepestDescent_ConvergenceRate}.

\newpar
Unfortunately, while Theorem \ref{Th:NewtonsMethodConvergenceRateTheorem} guarantees quadratic
convergence provided $\vecx_0$ is chosen sufficiently close $\vecx_{\bigstar}$, the proof given here does not
explicitly indicate ``how close'' is good enough. Part of the difficulty stems from the facts
we must establish that $H_f(\vecx)$ is Lipschitz continuous and also determine the value of
the Lipschitz constant, $M$. In most cases, Lipschitz continuity itself follows from ``smoothness'' conditions on $f$. For example, if all third-order partial derivatives of $f$ are continuous
in a neighborhood of $\vecx_{\bigstar}$, then Lipschitz continuity holds in that neighborhood. This follows from a higher-dimensional version of the Mean Value Theorem. Bounds on the third-order
partial derivatives can then be used to estimate $M$. In practice, however, this work is not necessary
as convergence of the iterates can occur even if $\vecx_0$ is relatively far from
$\vecx_{\bigstar}$.

\newpar
In the general case of a twice-differentiable function, Newton's Method and the Steepest Descent Method
frequently complement one another. The former method converges more rapidly when iterates are near the global
minimum and less rapidly when they are far away. For the Steepest Descent method, the opposite is true.

\newpar
Minimization of the \CM objective function,
\begin{equation*}
f(x_1,x_2 )= -1400 x_1^{\frac{1}{2}}x_2^{\frac{1}{3}} + 350x_1 +200x_2,
 \end{equation*}
 illustrates this principle. Table \ref{Ta:unconstrained_numeric_conpro_steepest_descent_newton} lists the first five
 iterations of both methods using an initial value of
 $\ds \vecx_0=\begin{bmatrix} 10 \\ 10 \end{bmatrix}$.


\begin{table}[h]
\caption{Results of the Steepest Descent and Newton's Methods applied to
the \CM objective function }
\label{Ta:unconstrained_numeric_conpro_steepest_descent_newton}
\begin{tabular}{|c|c|c|c|c|}
\hline
 & \multicolumn{2}{c|}{Steepest Descent Method} & \multicolumn{2}{c|}{Newton's Method}\\ \hline
$k$ & $\vecx_k^t$ &  $\ds \left\|  \vecx_k-\vecx_{\bigstar} \right\|$
& $\vecx_k^t$ &  $\ds \left\|  \vecx_k-\vecx_{\bigstar} \right\|$ \\
\hline
0 & $[10,10]$ & 119.76 & $[10,10]$ & 119.76 \\ \hline
                1  & $[91.680, 85.907]$ & 16.373 & $[28.063, 29.111]$ & 93.518 \\ \hline
                 2 & $[82.493, 95.792]$ & 7.444 & $[56.819, 61.720]$ & 50.103\\ \hline
                 3 & $[87.236, 100.200]$ & 1.435 & $[80.0760, 91.040]$ & 12.714\\ \hline
                4 & $[86.558, 100.930]$ & 0.892 & $[86.7370, 100.922]$ & 0.801\\ \hline
                5 & $[87.125, 101.455]$ & 0.175 & $[87.110, 101.626]$ & 0.004\\ \hline
\end{tabular}
\end{table}

\noi
Both methods demonstrate convergence of iterates to the global minimum, but they do so in a different manner. For small $k$,
the error $\ds \left\|  \vecx_k-\vecx_{\bigstar} \right\|$ using the Steepest Descent Method is much less than obtained
using Newton's Method.
Starting at $k=4$, this trend reverses. In fact, if $\epsilon=.01$, the Steepest Descent
Method requires nine iterations before $\nabla f(\vecx_k) < \epsilon$, whereas, Newton's Method requires only five.
That one of these methods
is more effective when iterates are far from the global
minimum and the other is more effective when iterates are near prompts
the question of whether it is possible to construct a single procedure that blends Steepest Descent
and Newton Methods
so as to be effective in either case. In the next section, we see that this
is the case and combine the
best of both techniques, creating a new procedure known as the
\emph{Levenberg-Marquardt Algorithm}, which
is particularly well-suited for approximating solutions of unconstrained NLPs, especially those arising in the area of nonlinear regression.

\section*{Exercises Section \protect \ref{S:NewtonsMethod}}

\addcontentsline{toc}{subsection}{Exercises}

\input{Exercises/ExercisesNewtonsMethod}
\endinput

