\[ \newcommand{\E}{\text{E}} \newcommand{\mbf}{\mathbf} \newcommand{\bo}{\mathbf} \newcommand{\bs}{\boldsymbol} \newcommand{\Var}{\text{Var}} \newcommand{\Cov}{\text{Cov}} \newcommand{\e}{\frac{1}{\sigma^2_e}} \newcommand{\f}{\frac{1}{\sigma^2_{\alpha}}} \]

11. Exercise: The Conjugate Gradient Method

In the conjugate gradient method, the solution at iteration \(n+1\) is:

\[ \mathbf{x}_{n+1} = \mathbf{x}_n + \alpha_n\mathbf{d}_n, \]

where

(11.1)\[ \alpha_n = \frac{-\mathbf{r}_n'\mathbf{r}_n}{\mathbf{d}_n'\mathbf{Ad}_n}, \]

(11.2)\[ \mathbf{r}_n = \mathbf{Ax}_n - \mathbf{b}, \]

(11.3)\[ \mathbf{d}_n = \mathbf{r}_n - \beta_{n-1}\mathbf{d}_{n-1}, \]

and

(11.4)\[ \beta_{n-1} = \frac{-\mathbf{r}_n'\mathbf{r}_n} {\mathbf{r}'_{n-1}\mathbf{r}_{n-1}}. \]

It can be shown that the residual can be computed as

(11.5)\[ \mathbf{r}_n = \mathbf{r}_{n-1} + \alpha_{n-1}\mathbf{Ad}_{n-1}, \]

and thus avoiding computation of \(\mathbf{Ax}_n\). However, using (11.5) leads to the accumulation of rounding errors. Thus, it is recommended that (11.2) is used every 50 iterations.

Unlike the Jacobi method, it is not intuitively obvious why the conjugate gradient method works. Following is an attempt to explain why the method works.

In the conjugate gradient method, the value of \(\bo{x}\) that minimizes

(11.6)\[\begin{split} \label{eq:10} \begin{split} f(\bo{x}) &= \frac{1}{2}\bo{x}'\bo{Ax} -\bo{b}'\bo{x}\\ \end{split} \end{split}\]

is obtained by line minimizations in \(n\) linearly independent directions, where \(n\) is the order of the symmetric positive definite matrix \(\bo{A}\). Note that after minimization of \(f(\bo{x})\) in any direction \(\bo{d}_i\), the gradient \(\bo{r}_{i+1}\) will be orthogonal to \(\bo{d}_i\).

In the conjugate gradient method, each direction \(\bo{d}_i\) is is chosen such that the gradient \(\bo{r}_{i+1}\) will also be orthogonal to all the directions \(\bo{d}_j\), for \(j<i\), that have already been used for line minimization. If the directions are also linearly independent, after \(n\) line minimizations, the function will be at its minimum in \(n\) linearly independent directions. Further, at this point, the gradient

\[ \bo{r} = \bo{Ax} - \bo{b} \]

is orthogonal to the \(n\) direction vectors. Thus, it must be the zero vector.

Following is a description of how the direction vectors are computed. Suppose the search is initiated at \(\bo{x}_0 = \bo{0}\). At this point, the gradient of the \(f(\bo{x})\) is

(11.7)\[\begin{split} \label{eq:11} \begin{split} \bo{r}_0 &= \bo{Ax}_o - \bo{b}\\ &= -b. \end{split} \end{split}\]

Let \(\bo{d}_0 = \bo{r}_0\) be the first direction for line minimization. After minimization of \(f(\bo{x})\) in the direction \(\bo{d}_0\), the value of \(\bo{x}\) is

(11.8)\[ \label{eq:12} \bo{x}_1 = \bo{x}_0 + \alpha_0\bo{d}_0. \]

At \(\bo{x}_1\), the gradient of the function is

(11.9)\[\begin{split} \label{eq:13} \begin{split} \bo{r}_1 &= \bo{Ax}_1 - \bo{b}\\ &= \bo{A}(\bo{x}_0 + \alpha_0\bo{d}_0) - \bo{b}\\ &= \bo{r}_0 + \alpha_0\bo{Ad}_0, \end{split} \end{split}\]

and \(\bo{r}_1\) is orthogonal to \(\bo{d}_0\). Writing the product \(\bo{d}_0\bo{r}_1=0\) as

(11.10)\[\begin{split} \label{eq:14} \begin{split} \bo{d}_0'\bo{r}_1 &= \bo{d}_0'\bo{r}_0 + \alpha_0\bo{d}_0'\bo{Ad}_0\\ &= 0 \end{split} \end{split}\]

shows that

(11.11)\[ \label{eq:15} \alpha_0 = \frac{- \bo{d}'_0\bo{r}_0}{\bo{d}_0'\bo{Ad}_0}, \]

and in general,

(11.12)\[ \label{eq:16} \alpha_i = \frac{- \bo{d}_i'\bo{r}_i}{\bo{d}_i'\bo{Ad}_i}. \]

At this point, line minimization proceeds in the direction \(\bo{d}_1\), and the value of \(\bo{x}\) at the minimum is

(11.14)\[ \label{eq:17} \bo{x}_2 = \bo{x}_1 + \alpha_1\bo{d}_1. \]

The gradient of the function at \(\bo{x}_2\) is

(11.14)\[\begin{split} \label{eq:18} \begin{split} \bo{r}_2 &= \bo{Ax}_2 - \bo{b}\\ &= \bo{r}_1 + \alpha_1\bo{Ad}_1, \end{split} \end{split}\]

and \(\bo{r}_2\) is orthogonal to \(\bo{d}_1\). In the conjugate gradient method, the direction \(\bo{d}_1\) is chosen such that \(\bo{r}_2\) will also be orthogonal to the direction \(\bo{d}_0\). Thus the product

(11.15)\[\begin{split} \label{eq:19} \begin{split} \bo{d}_0'\bo{r}_2 &= \bo{d}_0'\bo{r}_1 + \alpha_1\bo{d}_0'\bo{Ad}_1\\ &= \alpha_1\bo{d}_0'\bo{Ad}_1\\ &= 0. \end{split} \end{split}\]

So, the direction, \(\bo{d}_1\) must satisfy the condition

\[ \bo{d}_0'\bo{Ad}_1 = 0 \]

in order for \(\bo{r}_2\) to be orthogonal to \(\bo{d}_0\). To accomplish this, following the Grahm-Schmidt procedure, \(\bo{d}_1\) is written as

(11.16)\[ \label{eq:20} \bo{d}_1 = \bo{r}_1 - \beta_{10}\bo{d}_0. \]

Writing \(\bo{d}_0'\bo{Ad}_1=0\) as

(11.17)\[\begin{split} \label{eq:21} \begin{split} \bo{d}_0'\bo{Ad}_1 &= \bo{d}_0'\bo{Ar}_1 - \beta_{10}\bo{d}_0'\bo{Ad}_0\\ &= 0 \end{split} \end{split}\]

shows that

(11.18)\[ \label{eq:22} \beta_{10} = \frac{\bo{d}_0'\bo{Ar}_1}{\bo{d}_0'\bo{Ad}_0}. \]

In general,

(11.19)\[ \label{eq:23} \bo{r}_{i+1} = \bo{r}_i + \alpha_i\bo{Ad}_i, \]

and by induction, \(\bo{r}_{i+1}'\bo{d}_j=0\) provided \(\bo{d}_i'\bo{Ad}_j=0\) for \(j<i\). This can be achieved by using Grahm-Schmidt as

(11.20)\[ \label{eq:24} \bo{d}_i = \bo{r}_i - \sum_{j=0}^{i-1} \beta_{ij}\bo{d}_j, \]

where

(11.21)\[ \label{eq:25} \beta_{ij} = \frac{\bo{d}_j'\bo{Ar}_{i}}{\bo{d}_j'\bo{Ad}_j}. \]

Note that using the result \(\bo{r}_i\bo{d}_j=0\) for \(j<i\) in \eqref{eq:24} implies:

(11.22)\[ \label{eq:26} \bo{r}_i'\bo{r}_j=0\, \text{for}\, j<i. \]

Using (11.22) in (11.19), shows that \(\bo{r}_i'\bo{Ad}_j=0\) for \(j<i-1\). Thus, in (11.21), \(\beta_{ij}=0\) for \(j<i-1\), and the general expression for \(\bo{d}_i\) simplifies to

(11.23)\[ \label{eq:27} \bo{d}_i = \bo{r}_i - \beta_{i-1}\bo{d}_{i-1}, \]

where

(11.24)\[ \label{eq:28} \beta_{i-1} = \frac{\bo{d}_{i-1}'\bo{Ar}_i}{\bo{d}_{i-1}'\bo{Ad}_{i-1}}. \]

Writing \(\bo{r}_i\) as

(11.25)\[ \label{eq:29} \bo{r}_i = \bo{r}_{i-1} + \alpha_{i-1}\bo{Ad}_{i-1} \]

and pre-multiplying by \(\bo{r}_i'\) gives

(11.26)\[ \label{eq:30} \bo{r}_i'\bo{r}_i = \alpha_{i-1}\bo{r}_i'\bo{Ad}_{i-1}. \]

Pre-multiplying (11.25) by \(\bo{d}_{i-1}\) gives

(11.27)\[ \label{eq:31} \bo{d}_{i-1}'\bo{r}_{i-1} = -\alpha_{i-1}\bo{d}_{i-1}'\bo{Ad}_{i-1}. \]

Further, from (11.23), we can see that

(11.28)\[ \label{eq:32} \bo{d}_i'\bo{r}_i = \bo{r}_i'\bo{r}_i. \]

Using this is (11.27) gives

(11.29)\[ \label{eq:33} \bo{r}_{i-1}'\bo{r}_{i-1} = -\alpha_{i-1}\bo{d}_{i-1}'\bo{Ad}_{i-1}. \]

Using (11.26) and (11.29) in (11.24) gives

(11.30)\[ \label{eq:34} \beta_{i-1} = \frac{-\bo{r}_i'\bo{r}_i}{\bo{r}_{i-1}'\bo{r}_{i-1}}. \]

Finally, using (11.28) in (11.12) gives

(11.31)\[ \label{eq:35} \alpha_i = \frac{-\bo{r}_i'\bo{r}_i}{\bo{d}_i'\bo{Ad}_i}. \]

12. Preconditioned conjugate gradient method

In the PCCG method, the conjugate gradient method is applied to a transformed system of equations. The transformation of the system is based on a matrix \(\mathbf{M}\) that is approximately equal to \(\mathbf{A}\) and is easy to invert. A detailed explanation of PCCG is given in An Introduction to the Conjugate Gradient Method Without the Agonizing Pain by Jonathan Richard Shewchuk.

In PCCG, the solution at iteration \(n+1\) is:

(12.1)\[ \label{eq:36} \bo{x}_{n+1} = \bo{x}_n + \alpha_n\bo{d}_n, \]

where

(12.2)\[ \label{eq:37} \alpha_n = \frac{-\bo{r}_n'\bo{M}^{-1}\bo{r}_n} {\bo{d}_n'\bo{Ad}_n}, \]

(12.3)\[ \label{eq:38} \bo{r}_n = \bo{Ax}_n - \bo{b}, \]

(12.4)\[ \label{eq:39} \bo{d}_n = \bo{M}^{-1}\bo{r}_n - \beta_{n-1}\bo{d}_{n-1}, \]

(12.5)\[ \label{eq:40} \beta_{n-1} = \frac{-\bo{r}_n'\bo{M}^{-1}\bo{r}_n} {\bo{r}'_{n-1}\bo{M}^{-1}\bo{r}_{n-1}}. \]

As in the conjugate gradient method, the residual can be computed more efficiently as

(12.6)\[ \label{eq:41} \bo{r}_n = \bo{r}_{n-1} + \alpha_{n-1}\bo{Ad}_{n-1}. \]

However, it is recommended that (12.3) is used to every 50 iterations to avoid the accumulation of errors.

12.1. Exercise

Write a PCG function for the Conjugate Gradient Method algorithm described above.