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

\newcommand{\x}{\boldsymbol{x}}
\newcommand{\y}{\boldsymbol{y}}
\renewcommand{\v}{\boldsymbol{v}}
\newcommand{\w}{\boldsymbol{w}}
\renewcommand{\d}{\boldsymbol{d}}
\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}}


\begin{document}

\title{Numerical Methods I}
\author{The Gradient and Conjugate Gradient Method}
\date{Marcel Oliver}
\maketitle

\section{Reformulation as an Optimization Problem}

Let $A \in \R^{n\times n}$ be a symmetric, positive definite matrix.
Then $A\x = \b$ if and only if $\x$ minimizes the function
\begin{equation}
  \label{e.phi}
  \Phi(\x) = \tfrac12 \, \x^T A \x - \x^T \b \,.
\end{equation}

\begin{proof}
Assume first that $\x$ is a minimizer of $\Phi$.  Thus,
$\Phi(\x+\lambda \, \v)$ must have a minimum at $\lambda=0$ for any
fixed vector $\v \in \R^n$.  In other words,
\begin{equation}
  \frac{d}{d\lambda} \Phi(\x+\lambda \, \v) \Bigl|_{\lambda=0} = 0 \,.
\end{equation}
We compute
\begin{align}
  \frac{d}{d\lambda} \Phi(\x+\lambda \, \v)
  & = \tfrac12 \, (\x + \lambda \, \v)^T A \v
      + \tfrac12 \, \v^T A (\x + \lambda \, \v)
      - \v^T \b \notag \\
  & = \v^T \, \bigl( A (\x + \lambda \, \v) -\b \bigr) \,.
  \label{e.1dmin}
\end{align}
Setting $\lambda=0$, we see that we must require
\begin{equation}
  \v^T \, ( A \x  - \b) = 0 \,.
  \label{e.optimality}
\end{equation}
Since $\v$ is arbitrary, this implies that $A\x = \b$.

Vice versa, assume that $A\x = \b$.  Then, for any $\y \in \R^n$,
\begin{align}
  \Phi(\y) - \Phi(\x)
  & = \tfrac12 \, \y^T A \y - \y^T \b -
      \tfrac12 \, \x^T A \x + \x^T \b  \notag \\
  & = \tfrac12 \, \y^T A \y - \y^T A \x -
      \tfrac12 \, \x^T A \x + \x^T A \x \notag \\
  & = \tfrac12 \, (\y-\x)^T A (\y-\x) \notag \\
  & \geq 0
\end{align}
with equality if and only if $\x=\y$ since $A$ is positive definite.
We conclude that $\x$ is the unique minimizer of $\Phi$.
\end{proof}

\section{The Gradient Method}

The algorithm works works as follows:
\begin{itemize}
\item Choose a descent direction $\d_k$;
\item Walk in the descent direction until you reach a minimum along the
line;
\item Repeat.
\end{itemize}
There are many ways of choosing descent directions.  The simplest is to
take the direction of steepest descent,
\begin{equation}
  \d_k = -\grad \Phi (\x_k) = \b - A \x_k = \r_k \,.
\end{equation}
Given $\x_k$, compute the next iterate via
\begin{equation}
  \x_{k+1} = \x_k + \alpha_k \, \d_k
\end{equation}
where $\alpha_k$ is chosen such that
\begin{equation}
  \frac{d}{d \alpha_k} \Phi( \x_k + \alpha_k \, \d_k ) = 0 \,.
\end{equation}
Following the computation leading up to \eqref{e.1dmin}, we find that
\begin{equation}
  0 = \d_k^T \, \bigl( A (\x_k + \alpha_k \, \d_k) -\b \bigr) \,,
\end{equation}
and therefore
\begin{equation}
  \label{e.alpha}
  \alpha_k = \frac{\d_k^T \r_k}{\d_k^T A \d_k} \,.
\end{equation}
In summary, one iteration of the gradient method consists of the steps
\begin{align}
  \r_k     & = \b - A \x_k \,, \\
  \alpha_k & = \frac{\r_k^T \r_k}{\r_k^T A \r_k} \,, \\
  \x_{k+1} & = \x_k + \alpha_k \, \r_k \,.
\end{align}


\section{The Conjugate Gradient Method}

The conjugate gradient method is based on the concept of optimality
with respect to a set of search directions.  Once the algorithm has
reached optimality in some direction, we allow only changes that are
in a certain sense orthogonal, thereby preserving optimality under
iteration.

Specifically, we say that a point $\x \in \R^n$ is optimal with
respect to a subspace $V \subset \R^n$ if $\Phi$ has a minimum at $\x$
along each line passing through $\x$ in a direction $\v \in V$.

Repeating the calculation that lead up to \eqref{e.optimality}, we
find that $\x$ is optimal with respect to $V$ if
\begin{equation}
  \r^T \v = 0 \quad \text{for all } \v \in V \,,
\end{equation}
where $\r = \b - A\x$ is the corresponding residual.

Now assume that $\x_k$ is optimal with respect to some subspace $V_k$.
We would like to find a condition on a new descent direction $\d_{k}$
so that the next iterate under
\begin{equation}
  \x_{k+1} = \x_k + \alpha_k \, \d_k
  \label{e.iterate}
\end{equation}
is optimal not only with respect to the new descent direction $\d_k$,
as in the gradient method, but also with respect to all old descent
directions $V_k$.  Optimality with respect to the new direction
implies, as for the gradient method,
\begin{equation}
  \alpha_k = \frac{\d_k^T \r_k}{\d_k^T A \d_k} \,.
\end{equation}
To get optimality conditions with respect to $V_k$, we multiply
\eqref{e.iterate} by $-A$ from the left and add $\b$, so that
\begin{equation}
  \r_{k+1} = \r_k - \alpha_k \, A \d_k \,.
  \label{e.rk1}
\end{equation}
Optimality with respect to $V_k$ for $\x_{k+1}$ and $\x_k$ means that
\begin{equation}
  \r_{k+1}^T \v = \r_k^T \v = 0 \quad \text{for all } \v \in V_k \,,
\end{equation}
and therefore
\begin{equation}
  \v^T A \d_k = 0 \quad \text{for all } \v \in V_k \,.
\end{equation}
In other words, $\d_k$ must be $A$-orthogonal to all directions in
$V_k$.

The conjugate gradient method now works as follows. The first descent
direction is chosen as for the gradient method, namely $\d_1=\r_1$.
Each subsequent descent direction is the $A$-orthogonalization of
$\r_k$ with respect to the space of old descent directions
\begin{equation}
  V_k = \Span \{ \d_1, \dots, \d_{k-1} \} \,.
\end{equation}
Following this construction, each new descent direction $\d_k$ is
ultimately a linear combination of all previous residuals
$\r_1,\dots,\r_k$.  In particular, we see that
\begin{equation}
  V_k = \Span \{ \r_1, \dots, \r_{k-1} \} \,.
\end{equation}
Thus, equation \eqref{e.rk1} gives a recursion relation for the spaces
$V_k$, namely\footnote{In fact, it is easy to see that equality holds
unless the algorithm terminates with the exact solution, and that 
\begin{align*}
  V_{k+1} &= \Span \{ \r_1, \dots, \r_{k-1}, A\d_{k-1} \} \\
      &= \Span \{ A^0 \r_1, A^1 \r_1, \dots, A^{k-1} \r_1 \} \\
      &= \Span \{ A^0 \d_1, A^1 \d_1, \dots, A^{k-1} \d_1 \}  \,.
\end{align*}}
\begin{equation}
  V_{k+1} \subset V_k \oplus A V_k \,.
  \label{e.vrecursion}
\end{equation}

We now analyze the $A$-orthogonality condition when stepping from
$\x_k$ to $\x_{k+1}$ in detail.  First, $\x_k$ is already optimal with
respect to $V_k$---this has been achieved in the previous step of the
iteration---so that
\begin{equation}
  \r_k^T \v = 0 \quad \text{for all } \v \in V_k \,.
\end{equation}
Due to \eqref{e.vrecursion}, this implies, in particular, that
\begin{equation}
  \r_k^T A \v = 0 \quad \text{for all } \v \in V_{k-1} \,,
\end{equation}
which is to say that
\begin{equation}
  \r_k^T A \d_j =0 \quad \text{for all } j = 1, \dots, k-2 \,.
\end{equation}
In other words, $\r_k$ is already $A$-orthogonal to all but the
$(k-1)$-st previous descent direction.  This leads to a tremendous
simplification of the orthogonalization step; the
Gram--Schmidt-procedure needs only one projection, so that
\begin{equation}
  \d_k = \r_k -
         \frac{\r_k^T A \d_{k-1}}{\d_{k-1}^T A \d_{k-1}} \, \d_{k-1} 
\end{equation}
for $k \geq 2$.

We summarize the conjugate gradient method:
\begin{align}
  \r_k     & = \b - A \x_k \,, \\
  \d_k & = 
  \begin{cases}
    \r_1 & \text{for } k = 1 \\
    \r_k - \dfrac{\r_k^T A \d_{k-1}}{\d_{k-1}^T A \d_{k-1}} \,
         \d_{k-1} & \text{for } k \geq 2 \,,
  \end{cases} \\
  \alpha_k & = \frac{\d_k^T \r_k}{\d_k^T A \d_k} \,, \\
  \x_{k+1} & = \x_k + \alpha_k \, \d_k \,.
\end{align}
In exact arithmetic, one can show that either the dimension of $V_k$
increases by one each iteration, or the exact solution is reached.
This implies that the search space is exhausted after at most $n$
iteration and the algorithm must terminate with the exact answer.





\end{document}


Thus, we not only have $\r_k^T\d_j=0$ for $j=1,\dots,k-1$ as assumed
in the construction, but also
\begin{equation}
  \r_k^T\r_j=0 \quad \text{for } j=1,\dots,k-1 \,.
\end{equation}
Moreover, 
