Conjugate gradient method

right|thumb|A comparison of the convergence of [[gradient descent with optimal step size (in green) and conjugate vector (in red) for minimizing a quadratic function associated with a given linear system. Conjugate gradient, assuming exact arithmetic, converges in at most n steps, where n is the size of the matrix of the system (here n = 2).]]

In mathematics, the conjugate gradient method is an algorithm for the numerical solution of particular systems of linear equations, namely those whose matrix is positive-semidefinite. The conjugate gradient method is often implemented as an iterative algorithm, applicable to sparse systems that are too large to be handled by a direct implementation or other direct methods such as the Cholesky decomposition. Large sparse systems often arise when numerically solving partial differential equations or optimization problems.

The conjugate gradient method can also be used to solve unconstrained optimization problems such as energy minimization. It is commonly attributed to Magnus Hestenes and Eduard Stiefel, who programmed it on the Z4, and extensively researched it.

The biconjugate gradient method provides a generalization to non-symmetric matrices. Various nonlinear conjugate gradient methods seek minima of nonlinear optimization problems.

Description of the problem addressed by conjugate gradients

Suppose we want to solve the system of linear equations

:<math>\mathbf{A}\mathbf{x} = \mathbf{b}</math>

for the vector <math>\mathbf{x}</math>, where the known <math>n \times n</math> matrix <math>\mathbf{A}</math> is symmetric (i.e., <math>\mathbf{A}^\mathsf{T} = \mathbf{A} </math>), positive-definite (i.e. <math>\mathbf{x}^\mathsf{T}\mathbf{Ax} > 0 </math> for all non-zero vectors <math>\mathbf{x}</math> in <math>\mathbb{R}^n</math>), and real, and <math>\mathbf{b}</math> is known as well. We denote the unique solution of this system by <math>\mathbf{x}_*</math>.

Derivation as a direct method

The conjugate gradient method can be derived from several different perspectives, including specialization of the conjugate direction method for optimization, and variation of the Arnoldi/Lanczos iteration for eigenvalue problems. Despite differences in their approaches, these derivations share a common topic—proving the orthogonality of the residuals and conjugacy of the search directions. These two properties are crucial to developing the well-known succinct formulation of the method.

We say that two non-zero vectors <math>\mathbf{u}</math> and <math>\mathbf{v}</math> are conjugate (with respect to <math>\mathbf{A}</math>) if

:<math> \mathbf{u}^\mathsf{T} \mathbf{A} \mathbf{v} = 0. </math>

Since <math>\mathbf{A}</math> is symmetric and positive-definite, the left-hand side defines an inner product

:<math>

\mathbf{u}^\mathsf{T} \mathbf{A} \mathbf{v} =

\langle \mathbf{u}, \mathbf{v} \rangle_\mathbf{A} :=

\langle \mathbf{A} \mathbf{u}, \mathbf{v}\rangle =

\langle \mathbf{u}, \mathbf{A}^\mathsf{T} \mathbf{v}\rangle =

\langle \mathbf{u}, \mathbf{A}\mathbf{v}\rangle.

</math>

Two vectors are conjugate if and only if they are orthogonal with respect to this inner product. Being conjugate is a symmetric relation: if <math>\mathbf{u}</math> is conjugate to <math>\mathbf{v}</math>, then <math>\mathbf{v}</math> is conjugate to <math>\mathbf{u}</math>. Suppose that

:<math>P = \{ \mathbf{p}_1, \dots, \mathbf{p}_n \}</math>

is a set of <math>n</math> mutually conjugate vectors with respect to <math>\mathbf{A}</math>, i.e. <math>\mathbf{p}_i^\mathsf{T} \mathbf{A} \mathbf{p}_j = 0</math> for all <math>i \neq j</math>.

Then <math>P</math> forms a basis for <math>\mathbb{R}^n</math>, and we may express the solution <math>\mathbf{x}_*</math> of <math>\mathbf{Ax} = \mathbf{b}</math> in this basis:

:<math>\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{p}_i \Rightarrow \mathbf{A} \mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{A} \mathbf{p}_i.</math>

Left-multiplying the problem <math>\mathbf{Ax} = \mathbf{b}</math> with the vector <math>\mathbf{p}_k^\mathsf{T}</math> yields

:<math>

\mathbf{p}_k^\mathsf{T} \mathbf{b}

= \mathbf{p}_k^\mathsf{T} \mathbf{A} \mathbf{x}_*

= \sum^{n}_{i=1} \alpha_i \mathbf{p}_k^\mathsf{T} \mathbf{A} \mathbf{p}_i

= \sum^{n}_{i=1} \alpha_i \left \langle \mathbf{p}_k, \mathbf{p}_i \right \rangle_{\mathbf{A

= \alpha_k \left \langle \mathbf{p}_k, \mathbf{p}_k \right \rangle_{\mathbf{A </math>

and so

:<math>\alpha_k = \frac{\langle \mathbf{p}_k, \mathbf{b} \rangle}{\langle \mathbf{p}_k, \mathbf{p}_k \rangle_\mathbf{A.</math>

This gives the following method of the algorithm shows that <math>\mathbf{r}_i</math> is orthogonal to <math>\mathbf{r}_j</math>, i.e. <math>\mathbf{r}_i^\mathsf{T} \mathbf{r}_j=0 </math>, for <math>i \neq j</math>. And <math>\mathbf{p}_i</math> is <math>\mathbf{A}</math>-orthogonal to <math>\mathbf{p}_j</math>, i.e. <math>\mathbf{p}_i^\mathsf{T} \mathbf{A} \mathbf{p}_j=0 </math>, for <math>i \neq j</math>. This can be regarded that as the algorithm progresses, <math>\mathbf{p}_i</math> and <math>\mathbf{r}_i</math> span the same Krylov subspace, where <math>\mathbf{r}_i</math> form the orthogonal basis with respect to the standard inner product, and <math>\mathbf{p}_i</math> form the orthogonal basis with respect to the inner product induced by <math>\mathbf{A}</math>. Therefore, <math>\mathbf{x}_k</math> can be regarded as the projection of <math>\mathbf{x}</math> on the Krylov subspace.

That is, if the CG method starts with <math>\mathbf{x}_0 = 0</math>, then<math display="block">x_k =

\mathrm{argmin}_{y \in \mathbb{R}^n}

{\left\{(x_*-y)^{\top} A(x_*-y): y \in \operatorname{span}\left\{b, A b, \ldots, A^{k-1} b\right\}\right\</math> where <math>x_*</math> is the solution to <math>\mathbf{A} \mathbf{x}= \mathbf{b}</math>.

The algorithm is detailed below for solving <math>\mathbf{A} \mathbf{x}= \mathbf{b}</math> where <math>\mathbf{A}</math> is a real, symmetric, positive-definite matrix. The input vector <math>\mathbf{x}_0</math> can be an approximate initial solution or <math>\mathbf{0}</math>. It is a different formulation of the exact procedure described above.

:<math>\begin{align}

& \mathbf{r}_0 := \mathbf{b} - \mathbf{A x}_0 \\

& \hbox{if } \mathbf{r}_{0} \text{ is sufficiently small, then return } \mathbf{x}_{0} \text{ as the result}\\

& \mathbf{p}_0 := \mathbf{r}_0 \\

& k := 0 \\

& \text{repeat} \\

& \qquad \alpha_k := \frac{\mathbf{r}_k^\mathsf{T} \mathbf{r}_k}{\mathbf{p}_k^\mathsf{T} \mathbf{A p}_k} \\

& \qquad \mathbf{x}_{k+1} := \mathbf{x}_k + \alpha_k \mathbf{p}_k \\

& \qquad \mathbf{r}_{k+1} := \mathbf{r}_k - \alpha_k \mathbf{A p}_k \\

& \qquad \hbox{if } \mathbf{r}_{k+1} \text{ is sufficiently small, then exit loop} \\

& \qquad \beta_k := \frac{\mathbf{r}_{k+1}^\mathsf{T} \mathbf{r}_{k+1{\mathbf{r}_k^\mathsf{T} \mathbf{r}_k} \\

& \qquad \mathbf{p}_{k+1} := \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k \\

& \qquad k := k + 1 \\

& \text{end repeat} \\

& \text{return } \mathbf{x}_{k+1} \text{ as the result}

\end{align}</math>

This is the most commonly used algorithm. The same formula for <math>\beta_k</math> is also used in the Fletcher–Reeves nonlinear conjugate gradient method.

Restarts

We note that <math>\mathbf{x}_{1}</math> is computed by the gradient descent method applied to <math>\mathbf{x}_{0}</math>. Setting <math>\beta_{k}=0</math> would similarly make <math>\mathbf{x}_{k+1}</math> computed by the gradient descent method from <math>\mathbf{x}_{k}</math>, i.e., can be used as a simple implementation of a restart of the conjugate gradient iterations.

A norm of the residual is typically used for stopping criteria. The norm of the explicit residual <math>\mathbf{r}_{k+1} := \mathbf{b} - \mathbf{A x}_{k+1}</math> provides a guaranteed level of accuracy both in exact arithmetic and in the presence of the rounding errors, where convergence naturally stagnates. In contrast, the implicit residual <math>\mathbf{r}_{k+1} := \mathbf{r}_k - \alpha_k \mathbf{A p}_k</math> is known to keep getting smaller in amplitude well below the level of rounding errors and thus cannot be used to determine the stagnation of convergence.

Computation of alpha and beta

In the algorithm, <math>\alpha_k</math> is chosen such that <math>\mathbf{r}_{k+1}</math> is orthogonal to <math>\mathbf{r}_{k}</math>. The denominator is simplified from

:<math>\alpha_k = \frac{\mathbf{r}_{k}^\mathsf{T} \mathbf{r}_{k{\mathbf{r}_{k}^\mathsf{T} \mathbf{A} \mathbf{p}_k} = \frac{\mathbf{r}_k^\mathsf{T} \mathbf{r}_k}{\mathbf{p}_k^\mathsf{T} \mathbf{A p}_k} </math>

since <math>\mathbf{r}_{k+1} = \mathbf{p}_{k+1}-\mathbf{\beta}_{k}\mathbf{p}_{k}</math>. The <math>\beta_k</math> is chosen such that <math>\mathbf{p}_{k+1}</math> is conjugate to <math>\mathbf{p}_{k}</math>. Initially, <math>\beta_k</math> is

:<math>\beta_k = - \frac{\mathbf{r}_{k+1}^\mathsf{T} \mathbf{A} \mathbf{p}_k}{\mathbf{p}_k^\mathsf{T} \mathbf{A} \mathbf{p}_k}</math>

using

:<math>\mathbf{r}_{k+1} = \mathbf{r}_{k} - \alpha_{k} \mathbf{A} \mathbf{p}_{k}</math>

and equivalently

<math> \mathbf{A} \mathbf{p}_{k} = \frac{1}{\alpha_{k (\mathbf{r}_{k} - \mathbf{r}_{k+1}), </math>

the numerator of <math>\beta_k</math> is rewritten as

:<math> \mathbf{r}_{k+1}^\mathsf{T} \mathbf{A} \mathbf{p}_k = \frac{1}{\alpha_k} \mathbf{r}_{k+1}^\mathsf{T} (\mathbf{r}_k - \mathbf{r}_{k+1}) = - \frac{1}{\alpha_k} \mathbf{r}_{k+1}^\mathsf{T} \mathbf{r}_{k+1} </math>

because <math>\mathbf{r}_{k+1}</math> and <math>\mathbf{r}_{k}</math> are orthogonal by design. The denominator is rewritten as

:<math> \mathbf{p}_k^\mathsf{T} \mathbf{A} \mathbf{p}_k = (\mathbf{r}_k + \beta_{k-1} \mathbf{p}_{k-1})^\mathsf{T} \mathbf{A} \mathbf{p}_k = \frac{1}{\alpha_k} \mathbf{r}_k^\mathsf{T} (\mathbf{r}_k - \mathbf{r}_{k+1}) = \frac{1}{\alpha_k} \mathbf{r}_k^\mathsf{T} \mathbf{r}_k </math>

using that the search directions <math>\mathbf{p}_k</math> are conjugated and again that the residuals are orthogonal. This gives the <math>\beta</math> in the algorithm after cancelling <math>\alpha_k</math>.

Example code in Julia (programming language)

using LinearAlgebra

"""

x = conjugate_gradient(A, b, x0 = zero(b); atol=length(b)*eps(norm(b))

Return the solution to `A * x = b` using the conjugate gradient method.

`A` must be a positive definite matrix or other linear operator.

`x0` is the initial guess for the solution (default is the zero vector).

`atol` is the absolute tolerance on the magnitude of the residual `b - A * x`

for convergence (default is machine epsilon).

Returns the approximate solution vector `x`.

"""

function conjugate_gradient(

A, b::AbstractVector, x0::AbstractVector = zero(b); atol=length(b)*eps(norm(b))

)

x = copy(x0) # initialize the solution

r = b - A * x0 # initial residual

p = copy(r) # initial search direction

r²old = r' * r # squared norm of residual

k = 0

while r²old > atol^2 # iterate until convergence

Ap = A * p # search direction

α = r²old / (p' * Ap) # step size

@. x += α * p # update solution

Update residual:

if (k + 1) % 16 == 0 # every 16 iterations, recompute residual from scratch

r .= b .- A * x # to avoid accumulation of numerical errors

else

@. r -= α * Ap # use the updating formula that saves one matrix-vector product

end

r²new = r' * r

@. p = r + (r²new / r²old) * p # update search direction

r²old = r²new # update squared residual norm

k += 1

end

return x

end

</syntaxhighlight>

Example code in MATLAB

function x = conjugate_gradient(A, b, x0, tol)

% Return the solution to `A * x = b` using the conjugate gradient method.

% Reminder: A should be symmetric and positive definite.

if nargin < 4

tol = eps;

end

r = b - A * x0;

p = r;

rsold = r' * r;

x = x0;

while sqrt(rsold) > tol

Ap = A * p;

alpha = rsold / (p' * Ap);

x = x + alpha * p;

r = r - alpha * Ap;

rsnew = r' * r;

p = r + (rsnew / rsold) * p;

rsold = rsnew;

end

</syntaxhighlight>

Numerical example

Consider the linear system Ax = b given by

:<math>\mathbf{A} \mathbf{x}= \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \end{bmatrix},</math>

we will perform two steps of the conjugate gradient method beginning with the initial guess

:<math>\mathbf{x}_0 = \begin{bmatrix} 2 \\ 1 \end{bmatrix}</math>

in order to find an approximate solution to the system.

Solution

For reference, the exact solution is

:<math> \mathbf{x} = \begin{bmatrix} \frac{1}{11} \\\\ \frac{7}{11} \end{bmatrix} \approx \begin{bmatrix} 0.0909 \\\\ 0.6364 \end{bmatrix}</math>

Our first step is to calculate the residual vector r0 associated with x0. This residual is computed from the formula r0 = b - Ax0, and in our case is equal to

:<math>\mathbf{r}_0 = \begin{bmatrix} 1 \\ 2 \end{bmatrix} -

\begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix}

\begin{bmatrix} 2 \\ 1 \end{bmatrix} =

\begin{bmatrix}-8 \\ -3 \end{bmatrix} = \mathbf{p}_0.</math>

Since this is the first iteration, we will use the residual vector r0 as our initial search direction p0; the method of selecting pk will change in further iterations.

We now compute the scalar using the relationship

:<math> \alpha_0 = \frac{\mathbf{r}_0^\mathsf{T} \mathbf{r}_0}{\mathbf{p}_0^\mathsf{T} \mathbf{A p}_0} = \frac{\begin{bmatrix} -8 & -3 \end{bmatrix} \begin{bmatrix} -8 \\ -3 \end{bmatrix{ \begin{bmatrix} -8 & -3 \end{bmatrix} \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} -8 \\ -3 \end{bmatrix} } =\frac{73}{331}\approx0.2205</math>

We can now compute x1 using the formula

:<math>\mathbf{x}_1 = \mathbf{x}_0 + \alpha_0\mathbf{p}_0 = \begin{bmatrix} 2 \\ 1 \end{bmatrix} + \frac{73}{331} \begin{bmatrix} -8 \\ -3 \end{bmatrix} \approx \begin{bmatrix} 0.2356 \\ 0.3384 \end{bmatrix}.</math>

This result completes the first iteration, the result being an "improved" approximate solution to the system, x1. We may now move on and compute the next residual vector r1 using the formula

:<math>\mathbf{r}_1 = \mathbf{r}_0 - \alpha_0 \mathbf{A} \mathbf{p}_0 = \begin{bmatrix} -8 \\ -3 \end{bmatrix} - \frac{73}{331} \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} -8 \\ -3 \end{bmatrix} \approx \begin{bmatrix} -0.2810 \\ 0.7492 \end{bmatrix}.</math>

Our next step in the process is to compute the scalar that will eventually be used to determine the next search direction p1.

:<math>\beta_0 = \frac{\mathbf{r}_1^\mathsf{T} \mathbf{r}_1}{\mathbf{r}_0^\mathsf{T} \mathbf{r}_0} \approx \frac{\begin{bmatrix} -0.2810 & 0.7492 \end{bmatrix} \begin{bmatrix} -0.2810 \\ 0.7492 \end{bmatrix{\begin{bmatrix} -8 & -3 \end{bmatrix} \begin{bmatrix} -8 \\ -3 \end{bmatrix = 0.0088.</math>

Now, using this scalar , we can compute the next search direction p1 using the relationship

:<math>\mathbf{p}_1 = \mathbf{r}_1 + \beta_0 \mathbf{p}_0 \approx \begin{bmatrix} -0.2810 \\ 0.7492 \end{bmatrix} + 0.0088 \begin{bmatrix} -8 \\ -3 \end{bmatrix} = \begin{bmatrix} -0.3511 \\ 0.7229 \end{bmatrix}.</math>

We now compute the scalar using our newly acquired p1 using the same method as that used for .

:<math> \alpha_1 = \frac{\mathbf{r}_1^\mathsf{T} \mathbf{r}_1}{\mathbf{p}_1^\mathsf{T} \mathbf{A p}_1} \approx \frac{\begin{bmatrix} -0.2810 & 0.7492 \end{bmatrix} \begin{bmatrix} -0.2810 \\ 0.7492 \end{bmatrix{ \begin{bmatrix} -0.3511 & 0.7229 \end{bmatrix} \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} -0.3511 \\ 0.7229 \end{bmatrix} } = 0.4122.</math>

Finally, we find x2 using the same method as that used to find x1.

:<math>\mathbf{x}_2 = \mathbf{x}_1 + \alpha_1 \mathbf{p}_1 \approx \begin{bmatrix} 0.2356 \\ 0.3384 \end{bmatrix} + 0.4122 \begin{bmatrix} -0.3511 \\ 0.7229 \end{bmatrix} = \begin{bmatrix} 0.0909 \\ 0.6364 \end{bmatrix}.</math>

The result, x2, is a "better" approximation to the system's solution than x1 and x0. If exact arithmetic were to be used in this example instead of limited-precision, then the exact solution would theoretically have been reached after n = 2 iterations (n being the order of the system).

Finite Termination Property

Under exact arithmetic, the number of iterations required is no more than the order of the matrix. This behavior is known as the finite termination property of the conjugate gradient method. It refers to the method's ability to reach the exact solution of a linear system in a finite number of steps—at most equal to the dimension of the system—when exact arithmetic is used. This property arises from the fact that, at each iteration, the method generates a residual vector that is orthogonal to all previous residuals. These residuals form a mutually orthogonal set.

In an n-dimensional space, it is impossible to construct more than n linearly independent and mutually orthogonal vectors unless one of them is the zero vector. Therefore, once a zero residual appears, the method has reached the solution and must terminate. This ensures that the conjugate gradient method converges in at most n steps.

To demonstrate this, consider the system:

<math>

A = \begin{bmatrix}

3 & -2 \\

-2 & 4

\end{bmatrix}, \quad

\mathbf{b} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}

</math>

We start from an initial guess <math>\mathbf{x}_0 = \begin{bmatrix} 1 \\ 2 \end{bmatrix}</math>. Since <math>A</math> is symmetric positive-definite and the system is 2-dimensional, the conjugate gradient method should find the exact solution in no more than 2 steps. The following MATLAB code demonstrates this behavior:

A = [3, -2; -2, 4];

x_true = [1; 1];

b = A * x_true;

x = [1; 2]; % initial guess

r = b - A * x;

p = r;

for k = 1:2

Ap = A * p;

alpha = (r' * r) / (p' * Ap);

x = x + alpha * p;

r_new = r - alpha * Ap;

beta = (r_new' * r_new) / (r' * r);

p = r_new + beta * p;

r = r_new;

end

disp('Exact solution:');

disp(x);

</syntaxhighlight>

The output confirms that the method reaches <math>\begin{bmatrix} 1 \\ 1 \end{bmatrix}</math> after two iterations, consistent with the theoretical prediction. This example illustrates how the conjugate gradient method behaves as a direct method under idealized conditions.

Application to Sparse Systems

The finite termination property also has practical implications in solving large sparse systems, which frequently arise in scientific and engineering applications. For instance, discretizing the two-dimensional Laplace equation <math>\nabla^2 u = 0</math> using finite differences on a uniform grid leads to a sparse linear system <math>A \mathbf{x} = \mathbf{b}</math>, where <math>A</math> is symmetric and positive definite.

Using a <math>5 \times 5</math> interior grid yields a <math>25 \times 25</math> system, and the coefficient matrix <math>A</math> has a five-point stencil pattern. Each row of <math>A</math> contains at most five nonzero entries corresponding to the central point and its immediate neighbors. For example, the matrix generated from such a grid may look like:

<math>

A = \begin{bmatrix}

4 & -1 & 0 & \cdots & -1 & 0 & \cdots \\

-1 & 4 & -1 & \cdots & 0 & 0 & \cdots \\

0 & -1 & 4 & -1 & 0 & 0 & \cdots \\

\vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\

-1 & 0 & \cdots & -1 & 4 & -1 & \cdots \\

0 & 0 & \cdots & 0 & -1 & 4 & \cdots \\

\vdots & \vdots & \cdots & \cdots & \cdots & \ddots

\end{bmatrix}

</math>

Although the system dimension is 25, the conjugate gradient method is theoretically guaranteed to terminate in at most 25 iterations under exact arithmetic. In practice, convergence often occurs in far fewer steps due to the matrix's spectral properties. This efficiency makes CGM particularly attractive for solving large-scale systems arising from partial differential equations, such as those found in heat conduction, fluid dynamics, and electrostatics.

Convergence properties

The conjugate gradient method can theoretically be viewed as a direct method, as in the absence of round-off error it produces the exact solution after a finite number of iterations, which is not larger than the size of the matrix. In practice, the exact solution is never obtained since the conjugate gradient method is unstable with respect to even small perturbations, e.g., most directions are not in practice conjugate, due to a degenerative nature of generating the Krylov subspaces.

As an iterative method, the conjugate gradient method monotonically (in the energy norm) improves approximations <math>\mathbf{x}_{k}</math> to the exact solution and may reach the required tolerance after a relatively small (compared to the problem size) number of iterations. The improvement is typically linear and its speed is determined by the condition number <math>\kappa(A)</math> of the system matrix <math>A</math>: the larger <math>\kappa(A)</math> is, the slower the improvement.

However, an interesting case appears when the eigenvalues are spaced logarithmically for a large symmetric matrix. For example, let <math>A = Q D Q^T</math> where <math>Q</math> is a random orthogonal matrix and <math>D</math> is a diagonal matrix with eigenvalues ranging from <math>\lambda_n = 1</math> to <math>\lambda_1 = 10^6</math>, spaced logarithmically. Despite the finite termination property of CGM, where the exact solution should theoretically be reached in at most <math>n</math> steps, the method may exhibit stagnation in convergence. In such a scenario, even after many more iterations—e.g., ten times the matrix size—the error may only decrease modestly (e.g., to <math>10^{-5}</math>). Moreover, the iterative error may oscillate significantly, making it unreliable as a stopping condition. This poor convergence is not explained by the condition number alone (e.g., <math>\kappa_2(A) = 10^6</math>), but rather by the eigenvalue distribution itself. When the eigenvalues are more evenly spaced or randomly distributed, such convergence issues are typically absent, highlighting that CGM performance depends not only on <math>\kappa(A)</math> but also on how the eigenvalues are distributed.

If <math>\kappa(A)</math> is large, preconditioning is commonly used to replace the original system <math>\mathbf{A x}-\mathbf{b} = 0</math> with <math>\mathbf{M}^{-1}(\mathbf{A x}-\mathbf{b}) = 0</math> such that <math>\kappa(\mathbf{M}^{-1}\mathbf{A})</math> is smaller than <math>\kappa(\mathbf{A})</math>, see below.

Convergence theorem

Define a subset of polynomials as

:<math>

\Pi_k^* := \left\lbrace \ p \in \Pi_k \ : \ p(0)=1 \ \right\rbrace \,,

</math>

where <math> \Pi_k </math> is the set of polynomials of maximal degree <math> k </math>.

Let <math> \left( \mathbf{x}_k \right)_k </math> be the iterative approximations of the exact solution <math> \mathbf{x}_* </math>, and define the errors as <math> \mathbf{e}_k := \mathbf{x}_k - \mathbf{x}_* </math>.

Now, the rate of convergence can be approximated as

:<math>

\begin{align}

\left\| \mathbf{e}_k \right\|_\mathbf{A}

&= \min_{p \in \Pi_k^*} \left\| p(\mathbf{A}) \mathbf{e}_0 \right\|_\mathbf{A}

&\leq \min_{p \in \Pi_k^*} \, \max_{ \lambda \in \sigma(\mathbf{A})} | p(\lambda) | \ \left\| \mathbf{e}_0 \right\|_\mathbf{A}

&\leq 2 \left( \frac{ \sqrt{\kappa(\mathbf{A})}-1 }{ \sqrt{\kappa(\mathbf{A})}+1 } \right)^k \ \left\| \mathbf{e}_0 \right\|_\mathbf{A}

&\leq 2 \exp\left(\frac{-2k}{\sqrt{\kappa(\mathbf{A})\right) \ \left\| \mathbf{e}_0 \right\|_\mathbf{A}

\,,

\end{align}

</math>

where <math> \sigma(\mathbf{A}) </math> denotes the spectrum, and <math> \kappa(\mathbf{A}) </math> denotes the condition number.

This shows <math>k = \tfrac{1}{2}\sqrt{\kappa(\mathbf{A})} \log\left(\left\| \mathbf{e}_0 \right\|_\mathbf{A} \varepsilon^{-1}\right)</math> iterations suffices to reduce the error to <math>2\varepsilon</math> for any <math>\varepsilon>0</math>.

Note, the important limit when <math> \kappa(\mathbf{A}) </math> tends to <math> \infty </math>

:<math>

\frac{ \sqrt{\kappa(\mathbf{A})}-1 }{ \sqrt{\kappa(\mathbf{A})}+1 }

\approx 1 - \frac{2}{\sqrt{\kappa(\mathbf{A})

\quad \text{for} \quad

\kappa(\mathbf{A}) \gg 1

\,.

</math>

This limit shows a faster convergence rate compared to the iterative methods of Jacobi or Gauss–Seidel which scale as <math> \approx 1 - \frac{2}{\kappa(\mathbf{A})} </math>.

No round-off error is assumed in the convergence theorem, but the convergence bound is commonly valid in practice as theoretically explained

:<math>\mathbf{r}_0 := \mathbf{b} - \mathbf{A x}_0</math>

:<math> \textrm{Solve:}\mathbf{M}\mathbf{z}_0 := \mathbf{r}_0</math>

:<math>\mathbf{p}_0 := \mathbf{z}_0</math>

:<math>k := 0 \, </math>

:repeat

::<math>\alpha_k := \frac{\mathbf{r}_k^\mathsf{T} \mathbf{z}_k}{\mathbf{p}_k^\mathsf{T} \mathbf{A p}_k}</math>

::<math>\mathbf{x}_{k+1} := \mathbf{x}_k + \alpha_k \mathbf{p}_k</math>

::<math>\mathbf{r}_{k+1} := \mathbf{r}_k - \alpha_k \mathbf{A p}_k</math>

::if rk+1 is sufficiently small then exit loop end if

::<math>\mathrm{Solve}\ \mathbf{M}\mathbf{z}_{k+1} := \mathbf{r}_{k+1}</math>

::<math>\beta_k := \frac{\mathbf{r}_{k+1}^\mathsf{T} \mathbf{z}_{k+1{\mathbf{r}_k^\mathsf{T} \mathbf{z}_k}</math>

::<math>\mathbf{p}_{k+1} := \mathbf{z}_{k+1} + \beta_k \mathbf{p}_k</math>

::<math>k := k + 1 \, </math>

:end repeat

:The result is xk+1

The above formulation is equivalent to applying the regular conjugate gradient method to the preconditioned system

:<math>\mathbf{E}^{-1}\mathbf{A}(\mathbf{E}^{-1})^\mathsf{T}\mathbf{\hat{x=\mathbf{E}^{-1}\mathbf{b}</math>

where

:<math>\mathbf{EE}^\mathsf{T}=\mathbf{M}, \qquad \mathbf{\hat{x=\mathbf{E}^\mathsf{T}\mathbf{x}.</math>

The Cholesky decomposition of the preconditioner must be used to keep the symmetry (and positive definiteness) of the system. However, this decomposition does not need to be computed, and it is sufficient to know <math>\mathbf{M}^{-1}</math>. It can be shown that <math>\mathbf{E}^{-1}\mathbf{A}(\mathbf{E}^{-1})^\mathsf{T}</math> has the same spectrum as <math>\mathbf{M}^{-1}\mathbf{A}</math>.

The preconditioner matrix <math>\mathbf{M}</math> has to be symmetric positive-definite and fixed, i.e., cannot change from iteration to iteration.

If any of these assumptions on the preconditioner is violated, the behavior of the preconditioned conjugate gradient method may become unpredictable.

An example of a commonly used preconditioner is the incomplete Cholesky factorization.

Using the preconditioner in practice

It is important to keep in mind that we don't want to invert the matrix <math>\mathbf{M}</math> explicitly in order to get <math>\mathbf{M}^{-1}</math> for use in the process, since inverting <math>\mathbf{M}</math> would take more time/computational resources than solving the conjugate gradient algorithm itself. As an example, let's say that we are using a preconditioner coming from incomplete Cholesky factorization. The resulting matrix is the lower triangular matrix <math>\mathbf{L}</math>, and the preconditioner matrix is:

<math>\mathbf{M}=\mathbf{LL}^\mathsf{T}</math>

Then we have to solve:

<math>\mathbf{Mz}=\mathbf{r}</math>

<math>\mathbf{z}=\mathbf{M}^{-1}\mathbf{r}</math>

But:

<math>\mathbf{M}^{-1}=(\mathbf{L}^{-1})^\mathsf{T}\mathbf{L}^{-1}</math>

Then:

<math>\mathbf{z}=(\mathbf{L}^{-1})^\mathsf{T}\mathbf{L}^{-1}\mathbf{r}</math>

Let's take an intermediary vector <math>\mathbf{a}</math>:

<math>\mathbf{a}=\mathbf{L}^{-1}\mathbf{r}</math>

<math>\mathbf{r}=\mathbf{L}\mathbf{a}</math>

Since <math>\mathbf{r}</math> and <math>\mathbf{L}</math> and known, and <math>\mathbf{L}</math> is lower triangular, solving for <math>\mathbf{a}</math> is easy and computationally cheap by using forward substitution. Then, we substitute <math>\mathbf{a}</math> in the original equation:

<math>\mathbf{z}=(\mathbf{L}^{-1})^\mathsf{T}\mathbf{a}</math>

<math>\mathbf{a}=\mathbf{L}^\mathsf{T}\mathbf{z}</math>

Since <math>\mathbf{a}</math> and <math>\mathbf{L}^\mathsf{T}</math> are known, and <math>\mathbf{L}^\mathsf{T}</math> is upper triangular, solving for <math>\mathbf{z}</math> is easy and computationally cheap by using backward substitution.

Using this method, there is no need to invert <math>\mathbf{M}</math> or <math>\mathbf{L}</math> explicitly at all, and we still obtain <math>\mathbf{z}</math>.

The flexible preconditioned conjugate gradient method

In numerically challenging applications, sophisticated preconditioners are used, which may lead to variable preconditioning, changing between iterations. Even if the preconditioner is symmetric positive-definite on every iteration, the fact that it may change makes the arguments above invalid, and in practical tests leads to a significant slow down of the convergence of the algorithm presented above. Using the Polak–Ribière formula

:<math>\beta_k := \frac{\mathbf{r}_{k+1}^\mathsf{T} \left(\mathbf{z}_{k+1}-\mathbf{z}_{k}\right)}{\mathbf{r}_k^\mathsf{T} \mathbf{z}_k}</math>

instead of the Fletcher–Reeves formula

:<math>\beta_k := \frac{\mathbf{r}_{k+1}^\mathsf{T} \mathbf{z}_{k+1{\mathbf{r}_k^\mathsf{T} \mathbf{z}_k}</math>

may dramatically improve the convergence in this case. This version of the preconditioned conjugate gradient method can be called flexible, as it allows for variable preconditioning.

The flexible version is also shown to be robust even if the preconditioner is not symmetric positive definite (SPD).

The implementation of the flexible version requires storing an extra vector. For a fixed SPD preconditioner, <math>\mathbf{r}_{k+1}^\mathsf{T} \mathbf{z}_{k}=0,</math> so both formulas for are equivalent in exact arithmetic, i.e., without the round-off error.

The mathematical explanation of the better convergence behavior of the method with the Polak–Ribière formula is that the method is locally optimal in this case, in particular, it does not converge slower than the locally optimal steepest descent method.

Vs. the locally optimal steepest descent method

In both the original and the preconditioned conjugate gradient methods one only needs to set <math>\beta_k := 0</math> in order to make them locally optimal, using the line search, steepest descent methods. With this substitution, vectors are always the same as vectors , so there is no need to store vectors . Thus, every iteration of these steepest descent methods is a bit cheaper compared to that for the conjugate gradient methods. However, the latter converge faster, unless a (highly) variable and/or non-SPD preconditioner is used, see above.

Conjugate gradient method as optimal feedback controller for double integrator

The conjugate gradient method can also be derived using optimal control theory. In this approach, the conjugate gradient method falls out as an optimal feedback controller,<math display="block">u = k(x, v):= -\gamma_a \nabla f(x) - \gamma_b v </math> for the double integrator system,<math display="block">\dot x = v, \quad \dot v = u </math> The quantities <math>\gamma_a</math> and <math>\gamma_b</math> are variable feedback gains.

A pathological example

This example is from

Let <math display="inline">t \in (0, 1)</math>, and define<math display="block">W= \begin{bmatrix}

t & \sqrt{t} & & & & \\

\sqrt{t} & 1+t & \sqrt{t} & & & \\

& \sqrt{t} & 1+t & \sqrt{t} & & \\

& & \sqrt{t} & \ddots & \ddots & \\

& & & \ddots & & \\

& & & & & \sqrt{t} \\

& & & & \sqrt{t} & 1+t

\end{bmatrix}, \quad b=\begin{bmatrix}

1 \\

0 \\

\vdots \\

\end{bmatrix}</math>Since <math>W</math> is invertible, there exists a unique solution to <math display="inline">W x = b </math>. Solving it by conjugate gradient descent gives us rather bad convergence:<math display="block">\|b- Wx_k\|^2 = (1/t)^{k}, \quad \|b- Wx_n\|^2 =0</math>In words, during the CG process, the error grows exponentially, until it suddenly becomes zero as the unique solution is found.

Conjugate gradient method

Description of the problem addressed by conjugate gradients

Derivation as a direct method

Restarts

Computation of alpha and beta

Example code in Julia (programming language)

Example code in MATLAB

Numerical example

Solution

Finite Termination Property

Application to Sparse Systems

Convergence properties

Convergence theorem

Using the preconditioner in practice

The flexible preconditioned conjugate gradient method

Vs. the locally optimal steepest descent method

Conjugate gradient method as optimal feedback controller for double integrator

A pathological example

See also

References

Further reading

External links