Conjugate Gradient Methods

Long Chen

Problem Setting

We introduce the conjugate gradient (CG) and and its preconditioned version PCG methods for solving

where $A$ is a symmetric and positive definite (SPD) operator defined on a Hilbert space $\mathbb V$ with $\dim \mathbb V = N$. When $\mathbb V = \mathbb R^N$, $A$ is an SPD matrix.

Equation $\eqref{maineq}$ can be derived from the following optimization problem

As $f$ is strongly convex, the global minimizer exists and unique and satisfies Euler equation $\nabla f(u) = 0$ which is exactly equation $\eqref{maineq}$ .

We shall derive CG from the $A$-orthogonal projection. We use $(\cdot,\cdot)$ for the standard inner product and $(\cdot,\cdot)_A$ for the inner product introduced by the SPD operator $A$. That is

When talking about orthogonality, we refer to the default $(\cdot,\cdot)$ inner product and emphasize the orthogonality in $(\cdot,\cdot)_A$ by $A$-orthogonal or conjugate.

Given a vector $x\in \mathbb V$, the $A$-orthogonal projection of $x$ to a subspace $S\subseteq \mathbb V$ is a vector, denoted by ${\rm Proj}_S^Ax\in S$, by the relation

CG Algorithm

Start from an initial guess $u_0$. Let $p_0 = r_0 = -\nabla f (u_0)$ . For $k=0,1,2,\ldots, n$, let $\mathbb V_k = {\rm span}\{p_0, p_1, \cdots , p_k\}$ be a subspace spanned by $A$-orthogonal basis, i.e. $(p_i, p_j)_A = 0$ for $i\neq j, i,j=0,\ldots, k$.

CG consists of three steps:

1. compute $u_{k+1}$ by the $A$-orthogonal projection of $u-u_0$ to $\mathbb V_k$.
2. add residual vector $r_{k+1}$ to $\mathbb V_k$ to get $\mathbb V_{k+1}$.
3. apply Gram-Schmit process to get $A-$orthogonal vector $p_{k+1}$ .

We now briefly explain each step and present recrusive formulae.

Remark 1. We treat $u, v$ as points and their difference as a vector. Spaces $\mathbb V, \mathbb V_k$ are vectors spaces and $u\in \mathbb V$ means $u-0\in \mathbb V$. If we mix the points with vectors, the space we are looking for $u_{k+1}$ is indeed the affine space $u_0+\mathbb V_k$. $\square$

Given a subspace $\mathbb V_k = {\rm span}\{p_0, p_1, \cdots , p_k\}$ spanned by $A$-orthogonal basis, we compute $u_{k+1}$ by

The projection can be computed without knowing $u$ as $Au=b$ is known. Indeed, by definition, $(u_{k+1} - u_0,p_i)_A = ({\rm Proj}_{V_k}^A (u - u_0), p_i)_A = (u - u_0, p_i)_A = (A(u - u_0), p_i) = (r_0, p_i)$ for $i=0,1,\ldots,k$, which leads to the formulae

With $u_{k+1}$, we can compute a new vector $r_{k+1} = b - Au_{k+1}$. If $r_{k+1} =0$, which means $u_{k+1}=u$ is the solution, then we stop. Otherwise expand the subspace to a larger one $\mathbb V_{k+1} = {\rm span}\{p_0, p_1, \cdots , p_k, r_{k+1}\}$.

Then use Gram-Schmit process to make new added vector $r_{k+1}$ $A$-orthogonal to others. The new conjugate direction is

The formulae $\eqref{updatep}$ is justified by the following orthogonality which will be proved later on.

Lemma 1. The residual $r_{k+1}$ is $(\cdot,\cdot)$-orthogonal to $\mathbb V_{k}$ and $A$-orthogonal to $\mathbb V_{k-1}$.

Last, we present recursive formula. Starting from an initial guess $u_0$ and $p_0= r_0$, for $k=0, 1, 2, \cdots$, we use three recursion formulae to compute

The method is called conjugate gradient method since the search direction is obtained by a correction of the negative gradient (the residual) direction and conjugate to all previous directions.

More efficient formulae for $\alpha _k$ and $\beta _k$ can be derived using the orthogonality presented in Lemma 1.

We present the algorithm of the conjugate gradient method as follows.

​x
function u = CG(A,b,u,tol)
​
tol = tol*norm(b);
k = 1;
r = b - A*u; 
p = r;
r2 = r'*r;
while sqrt(r2) >= tol && k<length(b)
    Ap = A*p;
    alpha = r2/(p'*Ap);
    u = u + alpha*p;
    r = r - alpha*Ap;
    r2old = r2;
    r2 = r'*r;
    beta = r2/r2old;
    p = r + beta*p;
    k = k + 1;
end

Several remarks on the realization are listed below:

• The most time consuming part is the matrix-vector multiplication A*p. There is no need to form the matrix $A$ explicitly. A subroutine to compute the matrix-vector multiplication is enough. This is an attractive feature of Krylov subspace methods.
• The error is measured by the relative error of the residual $\|r\|\leq \text{tol} \|b\|$.
• A maximum iteration step is given to avoid a large loop.

Orthogonality and Optimality

We explore the orthogonality and optimality due to the $A$-projection to the subspace. By definition of $A$-projection, cf. $\eqref{Aproj}$ we have the orthogonality

Consequently

That is the shortest distance measured in the $A$-norm is given by the $A$-orthogonal projection.

We will use the $A$-orthogonality $\eqref{Aorth}$ to present a proof of Lemma 1.

Lemma 1. The residual $r_{k+1}$ is $(\cdot,\cdot)$-orthogonal to $\mathbb V_{k}$ and $A$-orthogonal to $\mathbb V_{k-1}$.

Proof. To simplify notation, we denote by $P_k = {\rm Proj}_{V_k}^A$. We write $u-u_{k+1} = u-u_0 - (u_{k+1}-u_0) = (I-P_{k}) (u-u_0)$. Therefore $u-u_{k+1} \bot _A \mathbb V_{k}$, which is equivalent to $r_{k+1} = A(u-u_{k+1})\bot \mathbb V_{k}$. The first orthogonality is thus proved.

Since at every step, the residual vector is added to expand the subspace, it can be easily proved by induction that $\mathbb V_{k} = {\rm span}\{p_0, p_1, \cdots, p_{k}\}= {\rm span}\{r_0, r_1, \cdots, r_{k}\}.$

By the recursive formula for the residual $r_{i+1} = r_{i} - \alpha _{i}Ap_{i}$, we get $Ap_{i} \in {\rm span}\{r_{i}, r_{i+1}\} \subset \mathbb V_{k}$ for $0\leq i\leq k-1$. Since we have proved $r_{k+1} \bot \mathbb V_{k}$, we get $(r_{k+1}, p_i)_A =(r_{k+1}, Ap_i) = 0$ for $0\leq i\leq k-1$, i.e. $r_{k+1}$ is $A$-orthogonal to $\mathbb V_{k-1}$. $\square$

In view of $\eqref{shortdist}$ , we have the optimality

As $f(v) - f(u) = \frac{1}{2}\|v - u\|_A^2$, we can obtain another optimality

We can also show

Exercise. Let $u_{k+1}$ be computed by $\eqref{uk1}$. Prove that

That is $u_{k+1}$ can be found by the steepest descent method starting from $u_k$ and moving along the searching direction $p_k$. From $\eqref{steep}$ we get another formulae $\alpha_k = (r_k, p_k)_A/(p_k, p_k)_A$.

PCG Algorithm

Let $B$ be an SPD matrix. We apply the Gram-Schmidt process to the subspace

to get another $A$-orthogonal basis. The three steps are still the same except in the second step, we add $Br_{k+1}$ instead of $r_{k+1}$.

1. compute $u_{k+1}$ by the $A$-orthogonal projection of $u-u_0$ to $\mathbb V_k$.
2. add residual vector $Br_{k+1}$ to $\mathbb V_k$ to get $\mathbb V_{k+1}$.
3. apply Gram-Schmit process to get $A-$orthogonal vector $p_{k+1}$ .

Starting from an initial guess $u_0$ and $p_0 = Br_0$. For $k=0, 1, 2, \cdots$, the three-term recursion formulae are

Proof of the following lemma is almost identical to that of Lemma 1.

Lemma 2. $Br_{k+1}$ is $A$-orthogonal to $\mathbb V_{k-1}$ and $p_{k+1}$ is $A$-orthogonal to $\mathbb V_k$.

We present the algorithm of preconditioned conjugate gradient method.

xxxxxxxxxx
function u = pcg(A,b,u,B,tol)
​
tol = tol*norm(b);
k = 1;
r = b - A*u;
rho =  1;
while sqrt(rho) >= tol && k<length(b)
    Br = B*r;
    rho = r'*Br;
    if k = 1 
       p = Br;
    else
       beta = rho/rho_old;
       p = Br + beta*p;
    end
    Ap = A*p;
    alpha = rho/(p'*Ap);
    u = u + alpha*p;
    r = r - alpha*Ap;
    rho_old = rho;
    k = k + 1;
end

Similarly we collect several remarks for the implementation

• If we treat the operator $A$ as a mapping from $\mathbb V$ to its dual space, i.e. $A: \mathbb V\to \mathbb V'$, then the residual is in the dual space $\mathbb V'$. Operator $B$ is the Reisz representation from $\mathbb V'$ to $\mathbb V$ in the inner product $(\cdot,\cdot)_B$. In the original CG, $B$ is the identity matrix.
• A general guiden principle for the design of preconditioners is to find the correct Hilbert space with the correct inner product such that the operator $A$ is continuous and stable. Then the corresponding Riesz representation of that inner product is a good preconditioner.
• For an effective PCG, the computation B*r is crucial. The matrix $B$ do not have to be formed explicitly. All we need is the matrix-vector multiplication which can be replaced by a subroutine.

Convergence Analysis

Recall that we have two bases for the subspace

Another basis is that

Lemma 3.

Proof. We prove it by induction. For $k = 0$, it is trivial. Suppose the statement holds for $k =i$. We prove it holds for $k=i+1$ by noting that

$\square$

The space ${\rm span}\{r_0, Ar_0, \cdots , A^{k}r_0\}$ is called Krylov subspace. The CG method belongs to a large class of Krylov subspace methods. In the sequel, we use $\mathcal P_k$ to denote the space of polynomials with degree $\leq k$.

Theorem 1. Let $u_k$ be the $k$th iteration in the CG method with an initial guess $u_0$. Then

Proof. For $v\in \mathbb V_{k-1}$, it can be expanded as

Let $p_k(x) = 1{\color{blue}+} \sum _{i=1}^k c_{i-1}x^i$. Then

The identity $\eqref{CGidentity}$ then follows from $\eqref{cg:inf}$.

Since $A^i$ is symmetric in the $A$-inner product, we have

which leads to the estimate $\eqref{CGratepoly}$. $\square$

The polynomial $p_k\in \mathcal P_k$ with constraint $p_k(0)=1$ will be called the residual polynomial. Various convergence results of CG method can be obtained by choosing specific residual polynomials.

Corollary 1. Let $\{\phi\}_{i=1}^N$ be the eigenvector basis of $A$ . Let $b \in {\rm span}\{\phi _{i_1}, \phi _{i_2}, \cdots , \phi _{i_k}\}$, $k\leq N$. Then CG method with $u_0=0$ will find the solution $Au=b$ in at most $k$ iterations. In particular, for a given $b\in \mathbb R^N$, the CG algorithm will find the solution $Au=b$ within $N$ iterations. Namely CG can be viewed as a direct method.

Proof. Denote by $\mathbb U_k = {\rm span}\{\phi _{i_1}, \phi _{i_2}, \cdots , \phi _{i_k}\}$ . Let $b = \sum_{l=1}^k b_l \phi_{i_l}$. Then $u = \sum _{l=1}^k u_l \phi _{i_l}$ with $u_l = b_l/\lambda_{i_l}$. Namely if $b\in \mathbb U_k$, so is $u$. Note that here the eigenvectors are introduced for analysis and unknown in the algorithm.

The subspace $\mathbb U_k$ is invariant under the action of $A$, i.e. $A\mathbb U_k = \mathbb U_k$. Therefore starting from $r_0=b\in \mathbb U_k$, the Krylov space $\mathbb V_i\subset \mathbb U_k$ for all $i\leq k-1$. If $r_{i+1}\neq 0$, the dimension $\mathbb V_{i+1}$ increases by 1. So at most $k$ iterations, $\mathbb V_k = \mathbb U_k$. The projection of $u\in \mathbb U_k$ to $\mathbb V_k = \mathbb U_k$ is itself.

Choose the polynomial $p(x) = \prod _{l=1}^k(1-x/\lambda _{i_l})$. Then $p$ is a residual polynomial of order $k$ and $p(0)=1, p(\lambda_{i_l})=0$. and

The assertion is then from the identity $\eqref{CGidentity}$.

Remark 2. The CG method can be also applied to symmetric and positive semi-definite matrix $A$. Let $\{\phi _i\}_{i=1}^k$ be the eigenvectors associated to $\lambda _{\min}(A)=0$. Then from Corollary \ref{ex:keigenvalue}, if $b \in {\rm range} (A) = {\rm span}\{\phi _{k+1}, \phi _{k+2}, \cdots, \phi _N\}$, then the CG method with $u_0\in {\rm range}(A)$ will find a solution $Au=b$ within $N-k-1$ iterations. $\square$

The CG is invented as a direct method but it is more effective to use as an iterative method. The rate of convergence depends crucially on the distribution of eigenvalues of $A$ and could converge to the solution within certain tolerance in steps $k\ll N$.

Theorem 2. Let $u_k$ be the $k$-th iteration of the CG method with $u_0$. Then

Proof. Let $b = \lambda _{\max}(A)$ and $a = \lambda _{\min}(A)$. Choose

where $T_k(x)$ is the Chebyshev polynomial of degree $k$ given by

It is easy to see that $p_k(0)=1$, and if $x\in [a, b]$, then

Hence

We set

Solving this equation for $\mathrm{e}^{\sigma}$, we have

with $\kappa (A)=b/a$. We then obtain

which complete the proof. $\square$

The estimate $\eqref{CGrate}$ shows that CG is in general better than the gradient method. Furthermore if the condition number of $A$ is close to one, the CG iteration will converge very fast. Even if $\kappa(A)$ is large, the iteration will perform well if the eigenvalues are clustered in a few small intervals; see the following estimate.

Corollary 2. Assume that $\sigma (A) = \sigma _0(A)\cup \sigma _1(A)$ and $l$ is the number of elements in $\sigma _0(A)$. Then

where

To perform the convergence analysis of PCG, we let $\tilde r_0 = Br_0$. Then one can easily verify that $\mathbb V_k = {\rm span}\{\tilde r_0, BA \tilde r_0, \cdots, (BA)^k \tilde r_0\}.$ Since the optimality $\eqref{cg:inf}$ still holds, and $BA$ is symmetric in the $A$-inner product, we obtain the following convergence rate of PCG.

Theorem 3. Let $A,B$ be SPD and let $u_k$ be the $k$th iteration in the PCG method with the preconditioner $B$ and an initial guess $u_0$. Then

The SPD matrix $B$ is called preconditioner. A good preconditioner should have the properties that the action of $B$ is easy to compute and that $\kappa (BA)$ is significantly smaller than $\kappa (A)$. The design of a good preconditioner requires the knowledge of the problem and finding a good preconditioner is a central topic of scientific computing.

An interesting fact is that an iterative method $u_{k+1} = u_k + B(f-Au_k)$ may not be convergent whereas $B$ can be a preconditioner. For example, the Jacobi method is not convergent for all SPD systems, but $B=D^{-1}$ can always be used as a preconditioner, which is often known as the diagonal preconditioner. Another popular preconditioner is obtained by an incomplete Cholesky factorization of $A$.