%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%  LaTeX file for demonstration of various ``texniques''.
%%%%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentstyle{article}

%--- I prefer to use more of the page than is allowed by defaut. Hence, ...
\textwidth=6.5in
\textheight=9.0in
\voffset=-0.75in
\hoffset=-0.75in
%-------------------

\begin{document}

%--------- Some of my favorite macros ------------------
\newtheorem{thm}{Theorem}
\newtheorem{cor}{Corollary}
\newtheorem{lm}{Lemma}
%
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bear}{\begin{eqnarray}}
\newcommand{\bears}{\begin{eqnarray*}}
\newcommand{\eear}{\end{eqnarray}}
\newcommand{\eears}{\end{eqnarray*}}
\newcommand{\eqdef}{\stackrel{\rm def}{=}}
\newcommand{\QED}{\hfill\rule{2mm}{2mm}\vspace{3mm}}
%
\newcommand{\sP}{{\cal P}}
\newcommand{\sR}{{\cal R}}
% ------------------------------------------------------

\title{Convergence properties of the \\
       Conjugate Gradient Method}

\author{J. G. Wade
       \thanks{ \today  }
}

% ---------------------------------------------------

\maketitle

\begin{abstract}
  The Conjugate Gradient method for quadratic minimization enjoys widespread 
popularity because it is relatively easy to implement and, at the same time, 
posseses good convergence properties. One of the best-known properties 
is that if the problem has $n$ degrees of freedom then 
convergence is obtained in at most $n$ steps (under exact arithmetic). 
This result certainly has some practical value. However, for very large $n$,
it may be infeasible or undesirable to carry out $n$ iterations. The good 
news is that it is often unnecessary also. 

In this talk we carry out a rather detailed convergence analysis, with two
goals. First, we want to study the {\em rate} of convergence, and 
second, we want to see why, under some circumstances, we can guanantee
convergence in fewer than $n$ steps.
\end{abstract}

% --------------------------------------------------
\section{Review of the method}
\setcounter{equation}{0}
 
 The notation and analysis here comes largely from \cite[\S 8.7]{SB}. 
See also \cite{L_blue} for
a similar analysis. For the c.g. method in infinite dimensions, see
\cite{L_red}. For excellent, ready-to-code algorithms and some 
additional enlightening analysis, see \cite{GV}. 

 The Conjugate gradient method is a non-linear iterative method for 
\[
   Ax=b,
\]
where $A$ is an $n\times n$ matrix, $A^T=A>0$, and $b\in {\bf R}^n$ 
is given. An equivalent formulation of the problem is that of 
minimizing the function
\bear \label{Fdef}
   F(z) &\eqdef& {1\over2}(b-Az)^TA^{-1}(b-Az) \\
        &=& {1\over2}z^TAz - b^Tz + {1\over2}b^TA^{-1}b \label{Fmin} 
\eear
Since $A^T=A>0$, from (\ref{Fmin}) we easily see that $x=\arg\min F(z)$ 
iff $Ax=b$.

\noindent\underline{Algorithm:}
\begin{tabbing}
Given \= $x_0$, generate $x_k, k=1,2,\dots$ by: \\
  \> $r_0 = b - Ax_0$,  \\
  \> $p_0 = r_0$.       \\
  \> For \= $k=0,1,2,\dots:$ \\
  \> \> If \= $p_0\neq 0$, then: \\
  \> \> \> (i) $a_k = (r_k^Tr_k)/(p_k^TAp_k)$ \\
  \> \> \> (ii) $x_{k+1} = x_k + a_kp_k$ \\
  \> \> \> (iii) $r_{k+1} = r_k - a_kAp_k$ \\
  \> \> \> (iv) $b_k = (r_{k+1}^Tr_{k+1})/(r_k^Tr_k)$ \\
  \> \> \> (v) $p_{k+1} = r_{k+1} + b_kp_k.$ \\
  \> \> End If. \\
  \> End While. \\
End.
\end{tabbing}

The cg method follows the pattern of a wide class 
of ``descent methods'' for minimization. At the $k^{th}$ step, we
(a) compute a ``descent direction'' $p_k$, and then (b) minimize
$F(x_k+sp_k)$ over $s>0$ to obtain $x_{k+1}$. 
One of the main features of the cg method is that the directions $p_k$ 
turn out to be mutually ``A-orthogonal'', that is, orthogonal in the
inner product
\[
  \langle x,y \rangle_A \eqdef x^TAy.
\]
Another main feature is this:
\begin{thm}
  If all arithmetic is carried out in infinite precision, then
$x_k=x$ for some $k\leq n$. 
\end{thm}
This result certainly has some practical value. However, for very large $n$,
it may be infeasible or undesirable to carry out $n$ iterations. The good 
news is that it is often unnecessary also. Pursuant to this comment, 
we want to look at two refinements of the analysis. First, we want to 
study the {\em rate} at which $\Vert x_k-x\Vert_A$ converges to zero. 
Second, we want to see why, under some circumstances, we can guanantee
$x_k=x$ for some $k\leq m$ where $m<n$.
% --------------------------------------------------
\section{Estimating the rate of convergence}
 
 Throughout the following, we shall use the notation
\bears
 e_k &\eqdef& x-x_k, \\
 r_k &\eqdef& b-Ax_k.
\eears
Clearly, the ``error'' $e_k$ and the ``residual'' $r_k$ are related by
\[
   Ae_k=r_k. 
\]
Also, denote
\[
  \sR_k \eqdef span\{ r_0, r_1, \ldots, r_k \}.
  \sP_k \eqdef span\{ p_0, p_1, \ldots, p_k \}
\]
We introduce the concept of a ``Krylov subspace'', which for a given 
vector $z$ is defined by
\[
  K_i(z,A) \eqdef span\{z, Az, A^2z, \ldots, A^{i-1}z\}
\]
The cg method is referred to as a ``Krylov method'', because of
the following lemma.
\begin{lm} \label{lm1}
  \begin{enumerate}\begin{enumerate}
   \item $\sP_k = \sR_k.$
   \item $\sR_k = K_{k+1}(r_0,A).$
   \item $e_k\in \bar{K}_{k+1}(e_0,A)$, where 
         \[
           \bar{K}_i(z,A) \eqdef \{ v\in K_i(z,A) \;\mid\; 
                                 \mbox{coef. of $z$ is $1$}\}.
         \]
  \end{enumerate}\end{enumerate}
\end{lm}
{\sc Proof} of (a): The result is clearly true for $k=0$, so we assume that 
it holds for some $k\geq 0$ and proceed by induction to show that
$\sP_{k+1} = \sR_{k+1}.$ Fix $w\in\sP_{k+1}$. Then $w$ may be 
written as $w=w_k+\alpha p_{k+1}$ for some 
$w_k\in\sP_k=\sR_k\subset\sR_{k+1}$, so to show 
that $w\in\sR_{k+1}$, it suffices to show that $p_{k+1}\in \sR_{k+1}.$
But step (v) of the algorithm and the inductive assumption together imply
that this is true. Hence we have $\sP_{k+1}\subset\sR_{k+1}$. In 
a completely analogous manner we can also show $\sP_{k+1}\subset\sR_{k+1}$,
and hence $\sP_{k+1}=\sR_{k+1}$. 
\newline Proof of (b): Again the result clearly holds for $k=0$. So we
inductively assume that $\sR_k = K_{k+1}(r_0,A)$ and show shat 
$\sR_{k+1}=K_{k+2}(r_0,A).$ Since
\bears
   \sR_{k+1}&=&\sR_k\oplus span\{r_{k+1}\} \\
            &=& K_{k+1}(r_0,A)\oplus span\{r_{k+1}\},
\eears
it suffices to show shat $r_{k+1}\in K_{k+2}(r_0,A)$ or, equivalently, 
that both terms on the right side of step (iii) lie in 
$K_{k+2}(r_0,A)$. The first term lies there by the inductive assumption.
As for the second term, by (a) and the inductive assumption we have
\[
   p_k \in \sP_k = \sR_k - K_{k+1}(r_0,A) 
     = span\{ r_0, Ar_0, \ldots, A^kr_0\},
\]
so that 
$Ap_k\in span\{ Ar_0, A^2r_0, \ldots, A^{k+1}r_0\} \subset K_{k+2}(r_0,A).$
\newline Proof of (c): From step (ii), since 
$p_{k-1}\in\sP_{k-1}=K_k(r_0,A)$, 
we have 
\bears
   x_k &\in& x_0 + span\{ r_0, Ar_0, \ldots, A^{k-1}r_0\} \\
       &=&   x_0 + span\{ Ae_0, A^2e_0, \ldots, A^ke_0\},
\eears
so that
\bears
   x_k-x = e_k &\in& e_0 + span\{ Ae_0, A^2e_0, \ldots, A^ke_0\} \\
               &=&   \bar{K}_{k+1}(e_0,A).
\eears
\QED

{\it Remark:} ``Krylov methods'' stand at the forefront of 
current research into iterative methods for problems in which $A$ is 
{\em nonsymetric}.

The following theorem is proved in \cite{SB}. In fact, in those authors' 
discussion of the cg method this result is the conceptual 
``point of departure''. In other words, they say ``wouldn't it be nice if
we could choose the $x_{k+1}$ in the following way'', and then they proceed
with the description and analysis of the method.
\begin{thm} \label{Fmin_thm}
  Given $x_k$ at the $k^th$ step of the algorithm, the $x_{k+1}$ satisfies
  \be \label{xkp1_min}
     F(x_{k+1}) = \min_{v\in\sR_k} F(x_k+v).
  \ee
\end{thm}
On the basis of this result and the Lemma, we can show:
\begin{thm} \label{polymin}
 \[
   \Vert e_{k} \Vert_A = \min_{v\in\bar{K}_k(e_0,A)} \Vert v \Vert_A.
 \]
\end{thm}

{\sc Proof:}
From (\ref{Fmin}) and the definition of the $A$-norm 
$\Vert v \Vert_A\eqdef\sqrt{\langle v,v\rangle_A}$, we have
$F(z) = \Vert x-z \Vert_A^2/2$. In light of this and Lemma \ref{lm1}(b),
(\ref{xkp1_min}) can be restated as
\bears
  \Vert e_{k+1}\Vert_A &=& \min_{v\in\sR_k}\Vert e_k+v \Vert_A \\
          &=& \min_{v\in\sR_k}\Vert e_k+v \Vert_A.
\eears
But 
\bears
   \{e_k+v\;\mid\; v\in K_{k+1}(r_0,A)\} 
   &=& span\{ e_k, r_0, Ar_0, \ldots, A^kr_0\} \\
   &=& span\{ e_k, Ae_0, A^2e_0, \ldots, A^{k+1}e_0\} \\
   &=& \bar{K}_{k+1}(e_0,A).
\eears
\QED

Now, denote by $\bar{\Pi}_k$ the set of polynomials $g(t)$ with degree
$\leq n$ and with the coefficient of $t^0$ equal to one. Then
Theorem \ref{polymin} implies that 
\[
   \Vert e_k \Vert_A = \min_{g\in\bar{\Pi}_k} Vert g(A)e_0 \Vert.
\]
Also, let $\sigma(A)$ denote the spectrum of $A$, let
$\lambda_i\in\sigma(A),i=1,/ldots,n$ be the (not necessarily distinct)
eigenvalues of $A$ arranged in nondescending order, 
and let $z_i$ be the corresponding eigenvectors 
(normalized in the Euclidean norm).
Then $e_0$ has representation
\be \label{e0}
   e_0 = \sum_{i=1}^n \hat{e}_iz_i, 
\ee
and
\[
   g(A)e_0 = \sum_{i=1}^n g(\lambda_i)\hat{e}_iz_i
\]
(why?), and 
\bears
   \Vert g(A)e_0 \Vert_A^2 
      &=& \sum_{i=1}^n g(\lambda_i)^2\hat{e}_i^2\lambda_i\\
      &\leq& \max_{\lambda\in\sigma(A)}g(\lambda)^2
             \sum_{i=1}^n \hat{e}_i^2\lambda_i \\
      &=& \max_{\lambda\in\sigma(A)}g(\lambda)^2
             \Vert e_0 \Vert_A^2.
\eears
Hence we have the following ``min-max'' (or, if you prefer, ``inf-sup''!) 
condition:
\begin{cor}
  \be \label{minmax}
     \Vert e_k \Vert_A \leq 
          \min_{g\in\bar{\Pi}_k} \max_{\lambda\in\sigma(A)}g(\lambda)^2
                   \Vert e_0 \Vert_A^2.
  \ee
\end{cor}

This result shows that the method converges in at most $n$ steps, since $A$
has at most $n$ distinct eigenvalues and we can certainly choose a 
$g\in\bar{\Pi}_n$ which vanishes at those eigenvalues. Of course, 
Theorem \ref{Fmin_thm}, which we {\em used} to establish this the Corollary, 
itself can be used more directly to show n-step convergence). By the same 
reasoning, however, we get the following stronger statement:
\begin{cor} \label{mstep}
  If  $A$ has at most $m$ distinct eigenvalues, then the cg method in
exact arithmetic converges in at most $m$ steps.
\end{cor}

We can also use (\ref{minmax}) to obtain an estimate of the convergence 
rate of the method. Define $g_k\in\bar{\Pi}_k$ by
\[
   g_k(t) = \beta_k T_k
            \Bigl(2({t-\lambda_1\over\lambda_n-\lambda_1} - 1\Bigl),
\]
where $T_k$ is the $k^{th}$ Chebyshev polynomial and $\beta_k$ is chosen so 
that $g_k(0)=\beta_kT_k(-\kappa) = 1$ (so that $g_k\in\bar{\Pi}_k$). Here, 
\[
   \kappa\eqdef {\lambda_n+\lambda_1\over\lambda_n-\lambda_1}
\]
is the referred to as the ``spectral condition number'' of $A$. 
The reason for this choice of 
$g_k$ is that, amoung the many very useful properties of this family of 
polynomials, the  Chebyshev polynomials  satisfy
\[
   \max_{s\in[-1,1]} \vert T_k(s) \vert = 1.
\]
Another, less well-known property is 
\[
  T_k(s) \geq (|s| + \sqrt{s^2-1})^k/2.
\]
From this we get that $\beta_k\leq (\kappa + \sqrt{\kappa^2-1})^{-k}/2$, and 
hence
\[
   \Vert e_k \Vert_A \leq (\kappa + \sqrt{\kappa^2-1})^{-k}
                            \Vert e_0 \Vert_A/2
\]
or, 
\be \label{estimate}
   \Vert e_k \Vert_A \leq \kappa^{-k}\Vert e_0 \Vert_A.
\ee
 
% --------------------------------------------------
\section{Numerical Examples}

 The standard finite difference approximation of the boundary-value problem
\bears
 [Ax](\xi) \eqdef \Bigl(a(\xi)x'(\xi)\Bigr)' &=&= b(\xi), \\
                                        x(0) &=& 0, \\
                                        x(1) &=& 0.
\eears
on a partion of $[0,1]$ into $n+1$ subintervals leads to an algbraic 
problem of the form
\be \label{discrete}
   A^{(n)}x^{(n)} = b^{(n)}.
\ee
As $n\rightarrow \infty$, we have $\lambda^{(n)}_1\rightarrow \lambda_1$, 
and $\lambda^{(n)}_n = {\cal O}(n^2)$. Hence the spectral condition 
number growns like $n^2$, and if the cg method to to solve (\ref{discrete}),
we use we expect to see performance degradation for large $n$.

As a test we take $a(\xi) = 1 + \exp(-10\xi^2)$, $b(\xi)\equiv 0$, and
$\hat{e}_i^{(n)}=1/\sqrt{n\lambda_i}$ in (\ref{e0}). 
This choice of in initial guess gives us 
\[
  \Vert e_0^{(n)} \Vert_A^2 = \sum_{i=1}^n (\hat{e}_i^{(n)})^2\lambda_i
                            = {1\over n}\sum_{i=1}^n 1,
\]
so that $\Vert e_0^{(n)} \Vert_A=1$ and the initial error is 
``equally distributed'' over the spectrum.

The figure below shows semilog plots of $\Vert e^{(n)}_k \Vert_A$
versus iteration number $k$ for $k=1, 2, \ldots, 10$, and for various 
values of $n$. 

\begin{thebibliography}{99}

\bibitem{GV}
  G.H. Golub and C. VanLoan,
  {\sl Matrix Computations},
  Johns Hopkins Press.

\bibitem{L_blue}
  Luenberger, David G., {\it Introduction to Linear and Nonlinear 
  Programming}, Addison-Wesley, New York, 1973.
  Wiley, New York, 1969.

\bibitem{L_red}
  Luenberger, David G., {\it Optimization by Vector Space Methods}, 
  Wiley, New York, 1969.

\bibitem{SB}
   Stoer, J. and Bulirsch, R., {\it Introduction to Numerical Analysis, 
   $2^{nd}$ ed.},    Springer-Verlag, New York, 1993.


\end{thebibliography}
\end{document}