\documentclass[10pt]{article}
\usepackage{amsmath,amsfonts,amsthm}

\renewcommand{\d}{\boldsymbol{d}}
\newcommand{\x}{\boldsymbol{x}}
\newcommand{\y}{\boldsymbol{y}}
\renewcommand{\v}{\boldsymbol{v}}
\newcommand{\w}{\boldsymbol{w}}
\newcommand{\p}{\boldsymbol{p}}
\renewcommand{\r}{\boldsymbol{r}}
\renewcommand{\a}{\boldsymbol{a}}
\renewcommand{\b}{\boldsymbol{b}}
\renewcommand{\c}{\boldsymbol{c}}
\newcommand{\e}{\boldsymbol{e}}
\newcommand{\grad}{\boldsymbol{\nabla}}
\renewcommand{\Re}{\operatorname{Re}}
\renewcommand{\Im}{\operatorname{Im}}
\newcommand{\Range}{\operatorname{Range}}
\newcommand{\Ker}{\operatorname{Ker}}
\newcommand{\Span}{\operatorname{Span}}
\renewcommand{\dim}{\operatorname{dim}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}

\newtheorem{lemma}{Lemma}


\begin{document}

\title{Numerical Methods II}
\author{Convergence of Gradient and Conjugate Gradient Methods}
\date{Marcel Oliver}
\maketitle

\section{Notation}

We apply the Gradient and the Conjugate Gradient methods to the
problem of finding the minimum of the quadratic functional
\begin{equation}
  \label{e.phi}
  \Phi(\x) = \tfrac12 \, \x^T A \x - \x^T \b \,,
\end{equation}
where $A \in \R^{n\times n}$ is symmetric and positive definite.
Recall that both methods are \emph{descent methods}.  They construct a
minimizing sequence
\begin{equation}
  \label{e.descent}
  \x_{k+1} = \x_k + \alpha_k \, \d_k
\end{equation}
where $\d_k$ is a descent direction, i.e.\ $\d_k^T \grad \Phi(x_k) <
0$, so that
\begin{equation}
  \Phi(\x_k + \alpha \, \d_k) < \Phi (\x_k)
\end{equation}
for small positive values of $\alpha$.  For quadratic functionals
\eqref{e.phi}, the value of $\alpha$ that minimizes $\Phi$ along the
line through $\x_k$ in the direction $\d_k$ is easily found to be
\begin{equation}
  \alpha = \frac{\d_k^T \r_k}{\d_k^T A \d_k} \,.
\end{equation}

We denote the error in the $k$-th step by
\begin{equation}
  \e_k = \x - \x_k \,,
\end{equation}
$\x$ being the location of the true minimum, and define the
\emph{residual}
\begin{equation}
  \r_k = A \e_k = b - A \x_k \,.
  \label{e.re}
\end{equation}
The last equality is true since the location of the minimum is the
solution of $A\x=\b$.  See the notes from Numerical Methods~I for
details.

Finally, we introduce the $A$-norm
\begin{equation}
  \| \x \|_A^2 \equiv \x^T A \x \,.
\end{equation}


\section{Convergence of the Gradient Method}

In the gradient method we always walk down the direction of steepest
descent, i.e.
\begin{equation}
  \d_k = - \grad \Phi(\x_k) = \r_k \,.
\end{equation}
A direct computation of the error norm shows that
\begin{align}
  \| \e_{k+1} \|_A^2
  & = \| \x - (\x_k + \alpha_k \, \d_k) \|_A^2
      \notag \\
  & = \| \e_k - \alpha_k \, \d_k \|_A^2
      \notag \\
  & = \| \e_k \|_A^2 - 2 \, \alpha_k \, \e_k^T A \d_k +
      \alpha_k^2 \, \| \d_k \|_A^2
      \notag \\
  & = \| \e_k \|_A^2 - 2 \,
      \frac{\r_k^T \r_k}{\r_k^T A \r_k} \, \e_k^T A \r_k
    + \frac{(\r_k^T \r_k)^2}{(\r_k^T A \r_k)^2} \,
      \r_k^T A \r_k
      \notag \\
  & = \| \e_k \|_A^2 \,
      \biggl(
        1 - \frac{(\r_k^T \r_k)^2}{\r_k^T A \r_k \, \r_k^T A^{-1} \r_k}
      \biggr) \,.
      \label{e.gradient-estimate}
\end{align}
To proceed further, we need the following result.

\begin{lemma}[Kantorovich inequality]
  Let $A \in \R^{n\times n}$ be symmetric and positive definite, and
  let $0<\lambda_{\operatorname{min}} < \lambda_{\operatorname{max}}$
  denote its smallest and largest eigenvalue, respectively.  Then
\begin{equation}
  \min_{\y \neq 0}
    \frac{(\y^T \y)^2}{\y^T A \y \, \y^T A^{-1} \y}
    = \frac{4 \, \lambda_{\operatorname{min}} \,
      \lambda_{\operatorname{max}}}{(\lambda_{\operatorname{min}} +
      \lambda_{\operatorname{max}})^2} \,.
  \label{e.kantorovich}
\end{equation}
\end{lemma}

Applying the Kantorovich inequality to the right side of
\eqref{e.gradient-estimate} and noting that
\begin{equation}
  1 - \frac{4 \, \lambda_{\operatorname{min}} \,
      \lambda_{\operatorname{max}}}{(\lambda_{\operatorname{min}} +
      \lambda_{\operatorname{max}})^2}
  = \frac{(\lambda_{\operatorname{min}} -
      \lambda_{\operatorname{max}})^2}{(\lambda_{\operatorname{min}} +
      \lambda_{\operatorname{max}})^2} \,,
\end{equation}
we find that
\begin{equation}
  \| \e_{k+1} \|_A
  \leq \frac{\lambda_{\operatorname{max}} -
       \lambda_{\operatorname{min}}}{\lambda_{\operatorname{max}} +
       \lambda_{\operatorname{min}}} \,
       \| \e_{k} \|_A
  \equiv \frac{\kappa - 1}{\kappa + 1} \, \| \e_{k} \|_A \,,
  \label{e.error-grad}
\end{equation}
where $\kappa = \lambda_{\operatorname{max}} /
\lambda_{\operatorname{min}}$ is the so-called spectral condition
number of $A$.

\begin{proof}[Proof of the Kantorovich inequality]
The Kantorovich inequality \eqref{e.kantorovich} is invariant under
rescaling of $\y$.  It is therefore sufficient to prove it for
arbitrary unit vectors.  Let $\v_1, \dots, \v_n$ be an orthonormal
basis of $\R^n$ consisting of eigenvectors of $A$ with corresponding
eigenvalues $\lambda_1, \dots, \lambda_n$.  We can write
\begin{equation}
  \y = \sum_{i=1}^n y_i \, \v_i
\end{equation}
so that
\begin{gather}
  \| \y \|^2_2 = \sum_{i=1}^n y_i^2 = 1 \,, \\
  \y^T A \y = \sum_{i=1}^n y_i^2 \, \lambda_i \,, \\
  \y^T A^{-1} \y = \sum_{i=1}^n y_i^2 \, \lambda_i^{-1} \,.
\end{gather}
Thus, proving the Kantorovich inequality reduces to solving the
constraint optimization problem \emph{find the maximum of}
\begin{equation}
  g(\y) = \y^T A \y \, \y^T A^{-1} \y
  \label{e.g}
\end{equation}
\emph{under the constraint}
\begin{equation}
  h(\y) = \y^T \y -1 = 0 \,.
\end{equation}
Such maximum must necessarily satisfy
\begin{equation}
  \grad g(\y) = \mu \, \grad h(\y) \,,
  \label{e.lagrange}
\end{equation}
where $\mu$ is the Lagrange multiplier.  By direct computation,
\begin{align}
  \frac{\partial g}{\partial y_i}
  & = 2 \, y_i \, \lambda_i \, \y^T A^{-1} \y
      + 2 \, y_i \, \lambda_i^{-1} \, \y^T A \y \,, \\
  \frac{\partial h}{\partial y_i}
  & = 2 \, y_i \,,
\end{align}
so that \eqref{e.lagrange} reads
\begin{gather}
  y_i \, \lambda_i \, \y^T A^{-1} \y
      + y_i \, \lambda_i^{-1} \, \y^T A \y 
  = \mu \, y_i \,,
  \label{e.lagrange3}
\intertext{or}
  \bigl(
    \lambda_i^2 \, \y^T A^{-1} \y
    - \lambda_i \, \mu + \y^T A \y
  \bigr) \,
  \frac{y_i}{\lambda_i} = 0 \,.
  \label{e.lagrange2}
\end{gather}
The expression in parenthesis is a quadratic equation in $\lambda_i$,
and can only be zero for at most two distinct eigenvalues.  Therefore,
there can be at most two non-zero components $y_i$ and $y_j$.  (If an
eigenvalue occurs with multiplicity larger than one, $g(\y)$ depends
only on the norm of the projection of $\y$ onto the corresponding
eigenspace, so the statement remains valid.)  Dividing
\eqref{e.lagrange3} by $y_i$ and equating the expressions for indices
$i$ and $j$, we obtain
\begin{gather}
    \lambda_i \,
    \biggl( \frac{y_i^2}{\lambda_i}  + \frac{y_j^2}{\lambda_j} \biggr) 
  + \frac{y_i^2 \, \lambda_i  + y_j^2 \, \lambda_j}{\lambda_i}
  = \lambda_j \,
    \biggl( \frac{y_i^2}{\lambda_i}  + \frac{y_j^2}{\lambda_j} \biggr) 
  + \frac{y_i^2 \, \lambda_i  + y_j^2 \, \lambda_j}{\lambda_j} \,,
\intertext{or}
  (y_j^2 - y_i^2) \,
  \biggl(
   \frac{\lambda_i}{\lambda_j} + \frac{\lambda_j}{\lambda_i} - 2
  \biggr)
  = (y_j^2 - y_i^2) \,
    \frac{(\lambda_i-\lambda_j)^2}{\lambda_i \lambda_j} \,.
\end{gather}
Therefore, $y_i^2=y_j^2=\tfrac12$ unless $\lambda_i=\lambda_j$ and, by
direct substitution into \eqref{e.g}, candidates for the maximum value
are
\begin{equation}
  \frac14 \, (\lambda_i + \lambda_j) \,
  \biggl( \frac1{\lambda_i} + \frac1{\lambda_j} \biggr)
  = \frac14 \,
    \biggl( \frac{\lambda_i}{\lambda_j} + 1 \biggr) \,
    \biggl( \frac{\lambda_j}{\lambda_i} + 1 \biggr)
  = \frac{(\lambda_i + \lambda_j)^2}{4 \, \lambda_i \, \lambda_j} \,.
\end{equation}
The expression in the middle shows that it is increasing in the ratio
$\lambda_i/\lambda_j$ when $\lambda_i>\lambda_j$; we thus take
$\lambda_i = \lambda_{\operatorname{max}}$ and
$\lambda_j=\lambda_{\operatorname{min}}$ to complete the proof.
\end{proof}



\section{Convergence of the CG Method}

Recall that one of the key properties of the Conjugate Gradient method
is that each new iterate is optimal with respect to \emph{all} descent
directions from the so-called Krylov subspace
\begin{align}
  V_k & = \Span \{ \d_1, \dots, \d_{k-1} \} \notag \\
      & = \Span \{ A^0 \r_1, \dots, A^{k-2} \r_1 \} \,.
\end{align}
For details, see the handout on the derivation of the CG method.
According to \eqref{e.re}, we can also write
\begin{equation}
  V_k  = \Span \{ A^1 \e_1, \dots, A^{k-1} \e_1 \} \,.
\end{equation}
Recall also the recursion for the residual,
\begin{equation}
  \r_{k+1} = \r_k - \alpha_k \, A \d_k
           \equiv \r_k + A \w_{k+1}
 \end{equation}
where $\w_{k+1}$ is some vector from $V_{k+1}$.  By definition of
$\e_k$, we can also write
\begin{equation}
  \e_{k+1} = \e_k + \w_{k+1} \,.
\end{equation}
We conclude inductively that $\e_{k} - \e_1 \in
V_{k}$ for $k=2,3,\dots$, and we can write
\begin{align}
  \e_k & = \e_1 + \sum_{i=1}^{k-1} \gamma_i \, A^i \, \e_1
           \notag \\
       & \equiv \phi_k (A) \, \e_1
  \label{e.ephi}
\end{align}
where
\begin{equation}
  \phi_k(x) = \sum_{i=0}^{k-1} \gamma_i \, x^i
\end{equation}
is some polynomial of degree $k-1$ with $\phi_k(0)=1$.

Since
\begin{align}
  \| \e_k \|_A^2
  & = (\x - \x_k)^T A (\x - \x_k)
      \notag \\
  & = \x^T A \x - 2 \, \x_k^T A \x + \x_k^T A \x_k
      \notag \\
  & = \x^T \b - 2 \, \x_k^T \b + \x_k^T A \x_k
      \notag \\
  & = 2 \, \Phi(\x_k) + \x^T \b \,,
\end{align}
last term on the right being independent of $\x_k$, optimality of
$\x_k$ with respect to a certain subspace is equivalent to optimality
of the $A$-norm of $\e_k$ with respect to the same subspace.  In other
words, CG is constructed in such a way that the polynomial $\phi_k$
which appears in \eqref{e.ephi} is the polynomial that minimizes the
$A$ norm of $\e_k$ among all polynomials or the same or lesser degree.

As in the proof of the Kantorovich inequality we express $\e_1$ in
terms of the orthonormal eigenvectors of $A$, i.e.
\begin{gather}
  \e_1 = \sum_{j=1}^n y_j \, \v_j  \,,
\intertext{so that}
  \e_k = \sum_{j=1}^n y_j \, \phi_k(\lambda_j) \, \v_j  \,,\\
  \| \e_k \|_A^2
  = \sum_{j=1}^n y_j^2 \, \phi_k^2(\lambda_j) \, \lambda_j \,.
\end{gather}
Let $P_k$ denote the vector space of polynomials of degree less or
equal to $k$.  The optimality condition of CG can therefore be
expressed as
\begin{align}
  \| \e_k \|_A^2
  & = \min_{\substack{\phi \in P_{k-1} \\ \phi(0)=1}} \,
      \sum_{j=1}^n y_j^2 \, \phi^2(\lambda_j) \, \lambda_j
      \notag \\
  & \leq \min_{\substack{\phi \in P_{k-1} \\ \phi(0)=1}} \,
      \max_{\lambda \in
        [\lambda_{\operatorname{min}},\lambda_{\operatorname{max}}]} 
      \, \phi^2(\lambda) \,
      \sum_{j=1}^n y_j^2 \, \lambda_j
      \notag \\
  & = \min_{\substack{\phi \in P_{k-1} \\ \phi(0)=1}} \,
      \max_{\lambda \in
        [\lambda_{\operatorname{min}},\lambda_{\operatorname{max}}]}
      \, \phi^2(\lambda) \,
      \| \e_1 \|_A^2 \,.
  \label{e.error-cg}
\end{align}
Let $T_k$ denote the Chebychev polynomial of degree $k$.  Then the
following is true.

\begin{lemma} \label{l.2}
If $0<\lambda_{\operatorname{min}}<\lambda_{\operatorname{max}}$, then
\begin{equation}
  \min_{\substack{\phi \in P_{k} \\ \phi(0)=1}} \,
  \max_{\lambda \in
    [\lambda_{\operatorname{min}},\lambda_{\operatorname{max}}]} \,
      |\phi(\lambda)|
  = T_k \biggl(
      \frac{\lambda_{\operatorname{max}}+\lambda_{\operatorname{min}}}% 
           {\lambda_{\operatorname{max}}-\lambda_{\operatorname{min}}}
    \biggr)^{-1}  
  \leq 2 \,
       \biggl(
         \frac{\sqrt \kappa - 1}{\sqrt \kappa + 1}
       \biggr)^k \,,
  \label{e.l2}
\end{equation}
where $\kappa = \lambda_{\operatorname{max}} /
\lambda_{\operatorname{min}}$.
\end{lemma}

Inserting this result into \eqref{e.error-cg}, we find that
\begin{equation}
  \| \e_k^{\operatorname{cg}} \|_A
  \leq 2 \,
       \biggl(
         \frac{\sqrt \kappa - 1}{\sqrt \kappa + 1}
       \biggr)^{k-1}
       \| \e_1 \|_A \,,
\end{equation}
whereas for the gradient method, from equation \eqref{e.error-grad},
\begin{equation}
  \| \e_k^{\operatorname{grad}} \|_A
  \leq \biggl(
         \frac{\kappa - 1}{\kappa + 1}
       \biggr)^{k-1}
       \| \e_1 \|_A \,.
\end{equation}
This shows that for poorly conditioned matrices the rate of
convergence of the conjugate gradient method is much better than that
of the gradient method.

\begin{proof}[Proof of Lemma~\ref{l.2}]
Recall that the Chebychev polynomials of order $k$ on the interval
$[-1,1]$ can be written in the form
\begin{equation}
  T_k(x) = \cos(k \, \arccos x) \,.
\end{equation}
It is clear that $|T_k| \leq 1$ on this interval.  Moreover,
$T_k(1)=1$, $T_k(-1) = (-1)^k$, and $T_k$ has $k-1$ distinct extrema
on $(-1,1)$, alternating between $1$ and $-1$.

To realize the minimum on the left of \eqref{e.l2}, we need a
polynomial that is uniformly small on the interval
$[\lambda_{\operatorname{min}}, \lambda_{\operatorname{max}}]$.  Since
$T_k$ is of uniform size on $[-1,1]$, we are tempted to use it as a
template for constructing a candidate minimizing polynomial by
remapping the interval $[\lambda_{\operatorname{min}},
\lambda_{\operatorname{max}}]$ onto $[-1,1]$ by a linear affine change
of variable $\ell$.  Setting
\begin{equation}
  \phi(x) = \frac{T_k(\ell(x))}{T_k(\ell(0))}
  \label{e.minimizer}
\end{equation}
will then also satisfy $\phi(0)=1$.  The ansatz $\ell(x)= mx+c$
together with the requirement that $\ell(\lambda_{\operatorname{min}})
= 1$ and $\ell(\lambda_{\operatorname{max}}) = -1$ gives
\begin{equation}
  \ell(x) = \frac{\lambda_{\operatorname{max}} +
  \lambda_{\operatorname{min}} - 2 \lambda}%
  {\lambda_{\operatorname{max}} -
  \lambda_{\operatorname{min}}} \,.
\end{equation}
Therefore,
\begin{equation}
  \max_{\lambda \in
    [\lambda_{\operatorname{min}},\lambda_{\operatorname{max}}]} \,
      |\phi(\lambda)|
  = \frac1{|T_k(\ell(0))|} \,
    \max_{x \in [-1,1]} \, |T_k(x)|
  = T_k \biggl(
      \frac{\lambda_{\operatorname{max}}+\lambda_{\operatorname{min}}}% 
           {\lambda_{\operatorname{max}}-\lambda_{\operatorname{min}}}
    \biggr)^{-1}
  \label{e.candidate}
\end{equation}
To complete the proof of the equality in \eqref{e.l2}, we need to show
that no other polynomial of degree $k$ can yield a bound lower than
\eqref{e.candidate}.  Assume the contrary, and let $\psi$ denote a
polynomial satisfying $\psi(0)=1$ and
\begin{equation}
  \max_{\lambda \in
    [\lambda_{\operatorname{min}},\lambda_{\operatorname{max}}]} \,
      |\psi(\lambda)|
  < \max_{\lambda \in
    [\lambda_{\operatorname{min}},\lambda_{\operatorname{max}}]} \,
      |\phi(\lambda)| \,.
\end{equation}
Then $\psi(0)-\phi(0)=0$, and the graph of $\psi$ must intersect the
graph of $\phi$ exactly $k$ times in the interval $
[\lambda_{\operatorname{min}}, \lambda_{\operatorname{max}}]$.  This
means that $\phi-\psi$ is a polynomial of degree $k$ with $k+1$ zeroes,
and must therefore be identically zero.

The proof of the inequality in \eqref{e.l2} is a homework exercise.
\end{proof}



\end{document}


First, multiplying this last equation by $y_i$ summing
over all indices $i$ and noting that $\y^T\y=1$, we find that
\begin{equation}
  \mu = \sum_{i=1}^n y_i^2 \, \lambda_i \, \y^T A^{-1} \y
      + \sum_{i=1}^n y_i^2 \, \lambda_i^{-1} \, \y^T A \y
      = 2 \, g(\y) \,.
\end{equation}
Moreover, for any $i$ such that $y_i \neq 0$ we can divide
\eqref{e.lagrange2} by $y_i$, so that
\begin{equation}
  g(\y) = \tfrac12 \, \lambda_i \, \y^T A^{-1} \y
      + \tfrac12 \, \lambda_i^{-1} \, \y^T A \y \,.
  \label{e.g2}
\end{equation}
This relation holds true at each local constraint extremum.  However,
we can also consider this expression as an arbitrary function $G(\y)$.
Clearly, the maximum of $G$ provides an upper bound for the maximum of
$g$. 

and maximizing $g$ reduces to maximizing the---now only
quadratic---expression on the right.  We thus play the Lagrange
multiplier game once more.  Differentiation of \eqref{e.g2} gives
\begin{equation}
  \frac{\partial g}{\partial y_j}
  = \frac{\lambda_i}{\lambda_j} \, y_j +
    \frac{\lambda_j}{\lambda_i} \, y_j \,,
\end{equation}
so that \eqref{e.lagrange} now reads
\begin{equation}
  \frac{\lambda_i}{\lambda_j} \, y_j +
  \frac{\lambda_j}{\lambda_i} \, y_j
  = 2 \mu \, y_j \,,
  \label{e.mu2}
\end{equation}
As before, we multiply with $y_j$, sum over all $j$, and compare with
equation \eqref{e.g2}, whence
\begin{equation}
  2 \mu = \lambda_i \, \y^T A^{-1} \y + \lambda_i^{-1} \, \y^T A \y
  = 2 \, g(\y) \,.
\end{equation}
Therefore, for every $j$ such that $y_j \neq 0$ we can divide
\eqref{e.mu2} by $y_j$ to find that
\begin{equation}
  g(\y) =  \frac{\lambda_i}{\lambda_j} +
  \frac{\lambda_j}{\lambda_i} \,.
\end{equation}

