extended discussion of the conjugate gradient method

Suppose we wish to solve the system

$$A𝐱=𝐛$$(1)

where $A$ is a symmetric positive definite matrix. If we define the function

$$\mathrm{\Phi }(𝐱)=\frac{1}{2}{𝐱}^{T}A𝐱-{𝐱}^T𝐛$$(2)

we realise that solving (1) is equivalent to minimizing $\mathrm{\Phi }$. This is because if $𝐱$ is a minimum of $\mathrm{\Phi }$ then we must have

$$\nabla \mathrm{\Phi }(𝐱)=A𝐱-𝐛=\mathrm{𝟎}$$(3)

These considerations give rise to the steepest descent algorithm. Given an approximation $\widetilde{𝐱}$ to the solution $𝐱$, the idea is to improve the approximation by moving in the direction in which $\mathrm{\Phi }$ decreases most rapidly. This direction is given by the gradient of $\mathrm{\Phi }$ in $\widetilde{𝐱}$. Therefore we formulate our algorithm as follows

Given an initial approximation ${𝐱}^{(0)}$, for $k\in 1\mathrm{\dots }N$,

$${𝐱}^{(k)}={𝐱}^{(k-1)}+{\alpha }_k{𝐫}^{(k)}$$

where ${𝐫}^{(k)}=𝐛-A{𝐱}^{(k-1)}=-\nabla \mathrm{\Phi }({𝐱}^{(k-1)})$ and ${\alpha }_k$ is a scalar to be determined.

${𝐫}^{(k)}$ is traditionally called the residual vector. We wish to choose ${\alpha }_k$ in such a way that we reduce $\mathrm{\Phi }$ as much as possibile in each iteration, in other words, we wish to minimize the function $\varphi ({\alpha }_k)=\mathrm{\Phi }({𝐱}^{(k-1)}+{\alpha }_k{𝐫}^{(k)})$ with respect to ${\alpha }_k$

It’s possible to demonstrate that the steepest descent algorithm described above converges to the solution $𝐱$, in an infinite time. The conjugate gradient method improves on this by finding the exact solution after only $n$ iterations. Let’s see how we can achieve this.

We say that $\widetilde{𝐱}$ is optimal with respect to a direction $𝐝$ if $\lambda =0$ is a local minimum for the function $\mathrm{\Phi }(\widetilde{𝐱}+\lambda 𝐝)$.

In the steepest descent algorithm, ${𝐱}^{(k)}$ is optimal with respect to ${𝐫}^{(k)}$, but in general it is not optimal with respect to ${𝐫}^{(0)},\mathrm{\dots },{𝐫}^{(k-1)}$. If we could modify the algorithm such that the optimality with respect to the search directions is preserved we might hope that that ${𝐱}^{(n)}$ is optimal with respect to $n$ linearly independent directions, at which point we would have found the exact solution.

Let’s make the following modification

${𝐩}^{(k)}$$=$(4)

${𝐱}^{(k)}$$=$${𝐱}^{(k-1)}+{\alpha }_k{𝐩}^{(k)}$(5)

where ${\alpha }_k$ and ${\beta }_k$ are scalar multipliers to be determined

We choose ${\alpha }_k$ in the same way as before, i.e. such that ${𝐱}^{(k)}$ is optimal with respect to ${𝐩}^{(k)}$

$${\alpha }_k=\frac{{𝐫}^{(k)}{}^T𝐩^{(k)}}{{𝐩}^{(k)}{}^{T}A{𝐩}^{(k)}}$$(6)

Now we wish to choose ${\beta }_k$ such that ${𝐱}^{(k)}$, is also optimal with respect to ${𝐩}^{(k-1)}$. We require that

$${\frac{\partial \mathrm{\Phi }({𝐱}^{(k)}+\lambda {𝐩}^{(k-1)})}{\partial \lambda }|}_{\lambda =0}={[A{𝐱}^{(k)}-𝐛]}^T{𝐩}^{(k-1)}=0$$(7)

Since ${𝐱}^{(k)}={𝐱}^{(k-1)}+{\alpha }_k{𝐩}^{(k)}$, and assuming that ${𝐱}^{(k-1)}$ is optimal with respect to ${𝐩}^{(k-1)}$ (i.e. that ${[A{𝐱}^{(k-1)}-𝐛]}^T{𝐩}^{(k-1)}=0$) we can rewrite this condition as

$${[A({𝐱}^{(k-1)}+{\alpha }_k{𝐩}^{(k)})-𝐛]}^T{𝐩}^{(k-1)}={\alpha }_k{[{𝐫}^{(k)}-{\beta }_k{𝐩}^{(k-1)}]}^{T}A{𝐩}^{(k-1)}=0$$(8)

and therefore we obtain the required value for ${\beta }_k$,

$${\beta }_k=\frac{{𝐫}^{(k)}{}^{T}A{𝐩}^{(k-1)}}{{𝐩}^{(k-1)}{}^{T}A{𝐩}^{(k-1)}}$$(9)

Now we want to show that ${𝐱}^{(k)}$ is also optimal with respect to ${𝐩}^{(1)},\mathrm{\dots },{𝐩}^{(k-2)}$, i.e. that

$${[A{𝐱}^{(k)}-𝐛]}^T{𝐩}^{(j)}={𝐫}^{(k+1)}{}^T𝐩^{(j)}=0\mathit{\hspace{1em}\hspace{1em}}\forall j\in 1\mathrm{\dots }k-2$$(10)

We do this by strong induction on $k$, assuming that for all $\mathrm{\ell }\in 1\mathrm{\dots }k-1,j\in 1\mathrm{\dots }\mathrm{\ell }$, ${𝐱}^{(\mathrm{\ell })}$ is optimal with respect to ${𝐩}^{(j)}$ or equivalently that

$${[A{𝐱}^{(\mathrm{\ell })}-𝐛]}^T{𝐩}^{(j)}={𝐫}^{(\mathrm{\ell }+1)}{}^T𝐩^{(j)}=0$$(11)

Noticing that ${𝐫}^{(k+1)}={𝐫}^{(k)}+{\alpha }_{k}A{𝐩}^{(k)}$, we want to show

$${𝐫}^{(k)}{}^T𝐩^{(j)}+{\alpha }_k{𝐩}^{(k)}{}^{T}A{𝐩}^{(j)}=0$$(12)

and therefore since ${𝐫}^{(k)}{}^T𝐩^{(j)}=0$ by inductive hypothesis, it suffices to prove that

$${𝐩}^{(k)}{}^{T}A{𝐩}^{(j)}=0\mathit{\hspace{1em}\hspace{1em}}\forall j\in 1\mathrm{\dots }k-2$$(13)

But, again by the definition of ${𝐩}^{(k)}$, this is equivalent to proving

$${𝐫}^{(k)}{}^{T}A{𝐩}^{(j)}-{\beta }_k{𝐩}^{(k-1)}{}^{T}A{𝐩}^{(j)}=0$$(14)

and since ${𝐩}^{(k-1)}{}^{T}A{𝐩}^{(j)}=0$ by inductive hypothesis, all we need to prove is that

$${𝐫}^{(k)}{}^{T}A{𝐩}^{(j)}=0$$(15)

To proceed we require the following lemma.

Let $V_{\mathrm{\ell }}=\mathrm{span}\{A^{\mathrm{\ell }-1}{𝐫}^{(1)},\mathrm{\dots },A{𝐫}^{(1)},{𝐫}^{(1)}\}$.

• ${𝐫}^{(\mathrm{\ell })},{𝐩}^{(\mathrm{\ell })}\in V_{\mathrm{\ell }}$

• If $𝐲\in V_{\mathrm{\ell }}$ then ${𝐫}^{(k)}⟂𝐲$ for all $k>\mathrm{\ell }$.

Since ${𝐩}^{(j)}\in V_j$, it follows that $A{𝐩}^{(j)}\in V_{j+1}$, and since ,

$${𝐫}^{(k)}{}^T(A{𝐩}^{(j)})=0$$(16)

To finish let’s prove the lemma. The first point is by induction. The base case $\mathrm{\ell }=1$ holds. Assuming that ${𝐫}^{(\mathrm{\ell }-1)},{𝐩}^{(\mathrm{\ell }-1)}\in V_{\mathrm{\ell }-1}$, we have that

$${𝐫}^{(\mathrm{\ell })}={𝐫}^{(\mathrm{\ell }-1)}+{\alpha }_{\mathrm{\ell }-1}A{𝐩}^{(\mathrm{\ell }-1)}\in V_{\mathrm{\ell }}\mathit{\hspace{1em}\hspace{1em}}{𝐩}^{(\mathrm{\ell })}={𝐫}^{(\mathrm{\ell })}-{\beta }_{\mathrm{\ell }}{𝐩}^{(\mathrm{\ell }-1)}\in V_{\mathrm{\ell }}$$(17)

For the second point we need an alternative characterization of $V_{\mathrm{\ell }}$. Since ${𝐩}^{(1)},\mathrm{\dots },{𝐩}^{(\mathrm{\ell })}\in V_{\mathrm{\ell }}$,

$$\mathrm{span}\{{𝐩}^{(1)},\mathrm{\dots },{𝐩}^{(\mathrm{\ell })}\}\subseteq V_{\mathrm{\ell }}$$(18)

By (11), we have that if for some $s\in 1\mathrm{\dots }\mathrm{\ell }$

$${𝐩}^{(s)}={\lambda }_{s-1}{𝐩}^{(s-1)}+\mathrm{\cdots }+{\lambda }_1{𝐩}^{(1)}$$(19)

then ${𝐩}^{(s)}{}^T𝐫^{(s)}=0$, but we know that

$${𝐩}^{(s)}{}^T𝐫^{(s)}={[{𝐫}^{(s)}-{\beta }_s{𝐩}^{(s-1)}]}^T{𝐫}^{(s)}={𝐫}^{(s)}{}^T𝐫^{(s)}$$(20)

but were this zero, we’d have ${𝐫}^{(s)}=\mathrm{𝟎}$ and we would have solved the original problem. Thus we conclude that ${𝐩}^{(s)}\notin \mathrm{span}\{{𝐩}^{(1)},\mathrm{\dots },{𝐩}^{(s-1)}\}$, so the vectors ${𝐩}^{(1)},\mathrm{\dots },{𝐩}^{(s)}$ are linearly independent. It follows that $\dim \mathrm{span}\{{𝐩}^{(1)},\mathrm{\dots },{𝐩}^{(\mathrm{\ell })}\}=\mathrm{\ell }$, which on the other hand is the maximum possible dimension of $V_{\mathrm{\ell }}$ and thus we must have

$$V_{\mathrm{\ell }}=\mathrm{span}\{{𝐩}^{(1)},\mathrm{\dots },{𝐩}^{(\mathrm{\ell })}\}$$(21)

which is the alternative characterization we were looking for.Now we have, again by (11), that if $𝐲\in V_{\mathrm{\ell }}$, and , then

$${𝐫}^{(k)}{}^T𝐲=0$$(22)

thus the second point is proven.