## Subsection8.3.6Final touches for the Conjugate Gradient Method

We finish our discussion of the Conjugate Gradient Method by revisiting the stopping criteria and preconditioning.

### Subsubsection8.3.6.1Stopping criteria

In theory, the Conjugate Gradient Method requires at most $n$ iterations to achieve the condition where the residual is zero so that $x^{(k)}$ equals the exact solution. In practice, it is an iterative method due to the error introduced by floating point arithmetic. For this reason, the iteration proceeds while $\| r^{(k)} \|_2 \geq \meps \| b \|_2$ and some maximum number of iterations is not yet performed.

### Subsubsection8.3.6.2Preconditioning

In Subsection 8.2.5 we noted that the method of steepest Descent can be greatly accelerated by employing a preconditioner. The Conjugate Gradient Method can be greatly accelerated. While in theory the method requires at most $n$ iterations when $A$ is $n \times n \text{,}$ in practice a preconditioned Conjugate Gradient Method requires very few iterations.

###### Homework8.3.6.1.

Add preconditioning to the algorithm in Figure 8.3.5.2 (right).

Solution

\begin{equation*} A x = b \end{equation*}

we pick a SPD preconditioner $M = \tilde L \tilde L^T$ and instead solve the equivalent problem

\begin{equation*} \begin{array}[t]{c} \underbrace{ \tilde L^{-1} A \tilde L^{-T} } \\ \tilde A \end{array} \begin{array}[t]{c} \underbrace{ \tilde L^T x } \\ \tilde x \end{array} = \begin{array}[t]{c} \underbrace{ \tilde L^{-1} b . } \\ \tilde b \end{array} \end{equation*}

This changes the algorithm in Figure 8.3.5.2 (right) to

\begin{equation*} \begin{array}{l} {\bf Given:~} A, b, M = \tilde L \tilde L^T \\ \tilde x^{(0)} := 0 \\ \tilde A = \tilde L^{-1} A \tilde L^{-T} \\ \tilde r^{(0)} := \tilde L^{-1} b \\ k := 0 \\ {\bf while~} \tilde r^{(k)} \neq 0 \\ ~~~ {\bf if~} k = 0 \\ ~~~ ~~~ \tilde p^{(k)} = \tilde r^{(0)} \\ ~~~ {\bf else} \\ ~~~ ~~~ \tilde \gamma_k := (\tilde r^{(k)\,T} \tilde r^{(k)}) / (\tilde r^{(k-1)\,T} \tilde r^{(k-1)}) \\ ~~~ ~~~ \tilde p^{(k)} := \tilde r^{(k)} + \tilde \gamma_k \tilde p^{(k-1)} \\ ~~~ {\bf endif} \\ ~~~ \tilde \alpha_k := \frac{\tilde r^{(k)\,T} \tilde r^{(k)}}{\tilde p^{(k)\,T} \tilde A \tilde p^{(k)}} \\ ~~~ \tilde x^{(k+1)} := \tilde x^{(k)} + \tilde \alpha_k \tilde p^{(k)} \\ ~~~ \tilde r^{(k+1)} := \tilde r^{(k)} - \tilde \alpha_k \tilde A \tilde p^{(k)} \\ ~~~ k := k + 1 \\ {\bf endwhile} \end{array} \end{equation*}

Now, much like we did in the constructive solution to Homework 8.2.5.1 we now morph this into an algorithm that more directly computes $x^{(k+1)} \text{.}$ We start by substituting

\begin{equation*} \tilde A =\tilde L^{-1} A \tilde L^{-T}, \tilde x^{(k)} = \tilde L^T x^{(k)}, \tilde r^{(k)} = \tilde L^{-1} r^{(k)}, \tilde p^{(k)} = \tilde L^{T} p^{(k)}, \end{equation*}

which yields

\begin{equation*} \begin{array}{l} {\bf Given:~} A, b, M = \tilde L \tilde L^T \\ \tilde L^T x^{(0)} := 0 \\ \tilde L^{-1} r^{(0)} := \tilde L^{-1} b \\ k := 0 \\ {\bf while~} \tilde L^{-1} r^{(k)} \neq 0 \\ ~~~ {\bf if~} k = 0 \\ ~~~ ~~~ \tilde L^{T} p^{(k)} = \tilde L^{-1} r^{(0)} \\ ~~~ {\bf else} \\ ~~~ ~~~ \tilde \gamma_k := ((\tilde L^{-1} r^{(k)})^T \tilde L^{-1} r^{(k)}) / (\tilde L^{-1} r^{(k-1)})^T \tilde L^{-1} r^{(k-1)}) \\ ~~~ ~~~ \tilde L^{T} p^{(k)} := \tilde L^{-1} r^{(k)} + \tilde \gamma_k \tilde L^{T} p^{(k-1)} \\ ~~~ {\bf endif} \\ ~~~ \tilde \alpha_k := \frac{(\tilde L^{-1} r^{(k)})^T \tilde L^{-1} r^{(k)}}{((\tilde L^{T} p^{(k)})^T \tilde L^{-1} A \tilde L^{-T} \tilde L^{T} p^{(k)}} \\ ~~~ \tilde L^T x^{(k+1)} := \tilde L^T x^{(k)} + \tilde \alpha_k \tilde L^{T} p^{(k)} \\ ~~~ \tilde L^{-1} r^{(k+1)} := \tilde L^{-1} r^{(k)} - \tilde \alpha_k \tilde L^{-1} A \tilde L^{-T} \tilde L^{T} p^{(k)} \\ ~~~ k := k + 1 \\ {\bf endwhile} \end{array} \end{equation*}

If we now simplify and manipulate various parts of this algorithm we get

\begin{equation*} \begin{array}{l} {\bf Given:~} A, b, M = \tilde L \tilde L^T \\ x^{(0)} := 0 \\ r^{(0)} := b \\ k := 0 \\ {\bf while~} r^{(k)} \neq 0 \\ ~~~ {\bf if~} k = 0 \\ ~~~ ~~~ p^{(k)} = M^{-1} r^{(0)} \\ ~~~ {\bf else} \\ ~~~ ~~~ \tilde \gamma_k := (r^{(k)\,T} M^{-1} r^{(k)}) / (r^{(k-1)\,T} M^{-1} r^{(k-1)}) \\ ~~~ ~~~ p^{(k)} := M^{-1} r^{(k)} + \tilde \gamma_k p^{(k-1)} \\ ~~~ {\bf endif} \\ ~~~ \tilde \alpha_k := \frac{r^{(k)\,T} M^{-1} r^{(k)}}{p^{(k)\,T} A p^{(k)}} \\ ~~~ x^{(k+1)} := x^{(k)} + \tilde \alpha_k p^{(k)} \\ ~~~ r^{(k+1)} := r^{(k)} - \tilde \alpha_k A p^{(k)} \\ ~~~ k := k + 1 \\ {\bf endwhile} \end{array} \end{equation*}

Finally, we avoid the recomputing of $M^{-1} r^{(k)}$ and $A p^{(k)}$ by introducing $z^{(k)}$ and $q^{(k)} \text{:}$

\begin{equation*} \begin{array}{l} {\bf Given:~} A, b, M = \tilde L \tilde L^T \\ x^{(0)} := 0 \\ r^{(0)} := b \\ k := 0 \\ {\bf while~} r^{(k)} \neq 0 \\ ~~~ z^{(k)} := M^{-1} r^{(k)} \\ ~~~ {\bf if~} k = 0 \\ ~~~ ~~~ p^{(k)} = z^{(0)} \\ ~~~ {\bf else} \\ ~~~ ~~~ \tilde \gamma_k := (r^{(k)\,T} z^{(k)}) / (r^{(k-1)\,T} z^{(k-1)}) \\ ~~~ ~~~ p^{(k)} := z^{(k)} + \tilde \gamma_k p^{(k-1)} \\ ~~~ {\bf endif} \\ ~~~ q^{(k)} := A p^{(k)} \\ ~~~ \tilde \alpha_k := \frac{r^{(k)\,T} z^{(k)}}{p^{(k)\,T} q^{(k)}} \\ ~~~ x^{(k+1)} := x^{(k)} + \tilde \alpha_k p^{(k)} \\ ~~~ r^{(k+1)} := r^{(k)} - \tilde \alpha_k q^{(k)} \\ ~~~ k := k + 1 \\ {\bf endwhile} \end{array} \end{equation*}

(Obviously, there are a few other things that can be done to avoid unnecessary recomputations of $r^{(k)\,T} z^{(k)} \text{.}$)