Skip to main content

Unit 8.3.6 Final touches for the Conjugate Gradient Method

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

Subsubsection 8.3.6.1 Stopping 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.

Subsubsection 8.3.6.2 Preconditioning

In Unit 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.

Homework 8.3.6.1.

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

Solution

To add preconditioning to

\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{.}\))