Inmathematics, theconjugate gradient method is analgorithm for thenumerical solution of particularsystems of linear equations, namely those whose matrix ispositive-semidefinite. The conjugate gradient method is often implemented as aniterative algorithm, applicable tosparse systems that are too large to be handled by a direct implementation or other direct methods such as theCholesky decomposition. Large sparse systems often arise when numerically solvingpartial differential equations or optimization problems.

The conjugate gradient method can also be used to solve unconstrainedoptimization problems such asenergy minimization. It is commonly attributed toMagnus Hestenes andEduard Stiefel,^[1][2] who programmed it on theZ4,^[3] and extensively researched it.^[4][5]

Thebiconjugate gradient method provides a generalization to non-symmetric matrices. Variousnonlinear conjugate gradient methods seek minima of nonlinear optimization problems.

Suppose we want to solve thesystem of linear equations

${\displaystyle \mathbf{A}\mathbf{x}=\mathbf{b}}$

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

The conjugate gradient method can be derived from several different perspectives, including specialization of the conjugate direction method for optimization, and variation of theArnoldi/Lanczos iteration foreigenvalue 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${\displaystyle \mathbf{u}}$ and${\displaystyle \mathbf{v}}$ are conjugate (with respect to${\displaystyle \mathbf{A}}$) if

${\displaystyle {\mathbf{u}}^{\mathsf{T}}\mathbf{A}\mathbf{v}=0.}$

Since${\displaystyle \mathbf{A}}$ is symmetric and positive-definite, the left-hand side defines aninner product

${\displaystyle {\mathbf{u}}^{\mathsf{T}}\mathbf{A}\mathbf{v}=⟨\mathbf{u},\mathbf{v}{⟩}_{\mathbf{A}}:=⟨\mathbf{A}\mathbf{u},\mathbf{v}⟩=⟨\mathbf{u},{\mathbf{A}}^{\mathsf{T}}\mathbf{v}⟩=⟨\mathbf{u},\mathbf{A}\mathbf{v}⟩.}$

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

${\displaystyle P=\{{\mathbf{p}}_1,\dots ,{\mathbf{p}}_n\}}$

is a set of${\displaystyle n}$ mutually conjugate vectors with respect to${\displaystyle \mathbf{A}}$, i.e.${\displaystyle {\mathbf{p}}_i^{\mathsf{T}}\mathbf{A}{\mathbf{p}}_j=0}$ for all${\displaystyle i\ne j}$. Then${\displaystyle P}$ forms abasis for${\displaystyle {\mathbb{R}}^n}$, and we may express the solution${\displaystyle {\mathbf{x}}_{\ast }}$ of${\displaystyle \mathbf{A}\mathbf{x}=\mathbf{b}}$ in this basis:

${\displaystyle {\mathbf{x}}_{\ast }=\sum _{i=1}^n{\alpha }_i{\mathbf{p}}_i\Rightarrow \mathbf{A}{\mathbf{x}}_{\ast }=\sum _{i=1}^n{\alpha }_i\mathbf{A}{\mathbf{p}}_i.}$

Left-multiplying the problem${\displaystyle \mathbf{A}\mathbf{x}=\mathbf{b}}$ with the vector${\displaystyle {\mathbf{p}}_k^{\mathsf{T}}}$ yields

${\displaystyle {\mathbf{p}}_k^{\mathsf{T}}\mathbf{b}={\mathbf{p}}_k^{\mathsf{T}}\mathbf{A}{\mathbf{x}}_{\ast }=\sum _{i=1}^n{\alpha }_i{\mathbf{p}}_k^{\mathsf{T}}\mathbf{A}{\mathbf{p}}_i=\sum _{i=1}^n{\alpha }_i{⟨{\mathbf{p}}_k,{\mathbf{p}}_i⟩}_{\mathbf{A}}={\alpha }_k{⟨{\mathbf{p}}_k,{\mathbf{p}}_k⟩}_{\mathbf{A}}}$

and so

${\displaystyle {\alpha }_k=\frac{⟨{\mathbf{p}}_k,\mathbf{b}⟩}{⟨{\mathbf{p}}_k,{\mathbf{p}}_k{⟩}_{\mathbf{A}}}.}$

This gives the following method^[4] for solving the equation${\displaystyle \mathbf{A}\mathbf{x}=\mathbf{b}}$: find a sequence of${\displaystyle n}$ conjugate directions, and then compute the coefficients${\displaystyle {\alpha }_k}$.

If we choose the conjugate vectors${\displaystyle {\mathbf{p}}_k}$ carefully, then we may not need all of them to obtain a good approximation to the solution${\displaystyle {\mathbf{x}}_{\ast }}$. So, we want to regard the conjugate gradient method as an iterative method. This also allows us to approximately solve systems where${\displaystyle n}$ is so large that the direct method would take too much time.

We denote the initial guess for${\displaystyle {\mathbf{x}}_{\ast }}$ by${\displaystyle {\mathbf{x}}_0}$ (we can assume without loss of generality that${\displaystyle {\mathbf{x}}_0=\mathbf{0}}$, otherwise consider the system${\displaystyle \mathbf{A}\mathbf{z}=\mathbf{b}-{\mathbf{A}\mathbf{x}}_0}$ instead). Starting with${\displaystyle {\mathbf{x}}_0}$ we search for the solution and in each iteration we need a metric to tell us whether we are closer to the solution${\displaystyle {\mathbf{x}}_{\ast }}$ (that is unknown to us). This metric comes from the fact that the solution${\displaystyle {\mathbf{x}}_{\ast }}$ is also the unique minimizer of the followingquadratic function

${\displaystyle f(\mathbf{x})=\frac{1}{2}{\mathbf{x}}^{\mathsf{T}}\mathbf{A}\mathbf{x}-{\mathbf{x}}^{\mathsf{T}}\mathbf{b},\mathbf{x}\in {\mathbf{R}}^n.}$

The existence of a unique minimizer is apparent as itsHessian matrix of second derivatives is symmetric positive-definite

${\displaystyle \mathbf{H}(f(\mathbf{x}))=\mathbf{A},}$

and that the minimizer (use${\displaystyle Df(\mathbf{x})=0}$) solves the initial problem follows from its first derivative

${\displaystyle \mathrm{\nabla }f(\mathbf{x})=\mathbf{A}\mathbf{x}-\mathbf{b}.}$

This suggests taking the first basis vector${\displaystyle {\mathbf{p}}_0}$ to be the negative of the gradient of${\displaystyle f}$ at${\displaystyle \mathbf{x}={\mathbf{x}}_0}$. The gradient of${\displaystyle f}$ equals${\displaystyle \mathbf{A}\mathbf{x}-\mathbf{b}}$. Starting with an initial guess${\displaystyle {\mathbf{x}}_0}$, this means we take${\displaystyle {\mathbf{p}}_0=\mathbf{b}-{\mathbf{A}\mathbf{x}}_0}$. The other vectors in the basis will be conjugate to the gradient, hence the nameconjugate gradient method. Note that${\displaystyle {\mathbf{p}}_0}$ is also theresidual provided by this initial step of the algorithm.

Let${\displaystyle {\mathbf{r}}_k}$ be theresidual at the${\displaystyle k}$th step:

${\displaystyle {\mathbf{r}}_k=\mathbf{b}-{\mathbf{A}\mathbf{x}}_k.}$

As observed above,${\displaystyle {\mathbf{r}}_k}$ is the negative gradient of${\displaystyle f}$ at${\displaystyle {\mathbf{x}}_k}$, so thegradient descent method would require to move in the directionr_k. Here, however, we insist that the directions${\displaystyle {\mathbf{p}}_k}$ must be conjugate to each other. A practical way to enforce this is by requiring that the next search direction be built out of the current residual and all previous search directions. The conjugation constraint is an orthonormal-type constraint and hence the algorithm can be viewed as an example ofGram-Schmidt orthonormalization. This gives the following expression:

(see the picture at the top of the article for the effect of the conjugacy constraint on convergence). Following this direction, the next optimal location is given by

${\displaystyle {\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\alpha }_k{\mathbf{p}}_k}$

with

${\displaystyle {\alpha }_k=\frac{{\mathbf{p}}_k^{\mathsf{T}}(\mathbf{b}-{\mathbf{A}\mathbf{x}}_k)}{{\mathbf{p}}_k^{\mathsf{T}}\mathbf{A}{\mathbf{p}}_k}=\frac{{\mathbf{p}}_k^{\mathsf{T}}{\mathbf{r}}_k}{{\mathbf{p}}_k^{\mathsf{T}}\mathbf{A}{\mathbf{p}}_k},}$

where the last equality follows from the definition of${\displaystyle {\mathbf{r}}_k}$ .The expression for${\displaystyle {\alpha }_k}$ can be derived if one substitutes the expression forx_k+1 intof and minimizing it with respect to${\displaystyle {\alpha }_k}$

The above algorithm gives the most straightforward explanation of the conjugate gradient method. Seemingly, the algorithm as stated requires storage of all previous searching directions and residue vectors, as well as many matrix–vector multiplications, and thus can be computationally expensive. However, a closer analysis of the algorithm shows that${\displaystyle {\mathbf{r}}_i}$ is orthogonal to${\displaystyle {\mathbf{r}}_j}$, i.e.${\displaystyle {\mathbf{r}}_i^{\mathsf{T}}{\mathbf{r}}_j=0}$, for${\displaystyle i\ne j}$. And${\displaystyle {\mathbf{p}}_i}$ is${\displaystyle \mathbf{A}}$-orthogonal to${\displaystyle {\mathbf{p}}_j}$, i.e.${\displaystyle {\mathbf{p}}_i^{\mathsf{T}}\mathbf{A}{\mathbf{p}}_j=0}$, for${\displaystyle i\ne j}$. This can be regarded that as the algorithm progresses,${\displaystyle {\mathbf{p}}_i}$ and${\displaystyle {\mathbf{r}}_i}$ span the sameKrylov subspace, where${\displaystyle {\mathbf{r}}_i}$ form the orthogonal basis with respect to the standard inner product, and${\displaystyle {\mathbf{p}}_i}$ form the orthogonal basis with respect to the inner product induced by${\displaystyle \mathbf{A}}$. Therefore,${\displaystyle {\mathbf{x}}_k}$ can be regarded as the projection of${\displaystyle \mathbf{x}}$ on the Krylov subspace.

That is, if the CG method starts with${\displaystyle {\mathbf{x}}_0=0}$, then^[6]$${\displaystyle x_k={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{y\in {\mathbb{R}}^n}\left\{(x-y{)}^{\mathrm{\top }}A(x-y):y\in \mathrm{span}\left\{b,Ab,\dots ,A^{k-1}b\right\}\right\}}$$The algorithm is detailed below for solving${\displaystyle \mathbf{A}\mathbf{x}=\mathbf{b}}$ where${\displaystyle \mathbf{A}}$ is a real, symmetric, positive-definite matrix. The input vector${\displaystyle {\mathbf{x}}_0}$ can be an approximate initial solution or${\displaystyle \mathbf{0}}$. It is a different formulation of the exact procedure described above.

This is the most commonly used algorithm. The same formula for${\displaystyle {\beta }_k}$ is also used in the Fletcher–Reevesnonlinear conjugate gradient method.

We note that${\displaystyle {\mathbf{x}}_1}$ is computed by thegradient descent method applied to${\displaystyle {\mathbf{x}}_0}$. Setting${\displaystyle {\beta }_k=0}$ would similarly make${\displaystyle {\mathbf{x}}_{k+1}}$ computed by thegradient descent method from${\displaystyle {\mathbf{x}}_k}$, i.e., can be used as a simple implementation of a restart of the conjugate gradient iterations.^[4] Restarts could slow down convergence, but may improve stability if the conjugate gradient method misbehaves, e.g., due toround-off error.

The formulas${\displaystyle {\mathbf{x}}_{k+1}:={\mathbf{x}}_k+{\alpha }_k{\mathbf{p}}_k}$ and${\displaystyle {\mathbf{r}}_k:=\mathbf{b}-{\mathbf{A}\mathbf{x}}_k}$, which both hold in exact arithmetic, make the formulas${\displaystyle {\mathbf{r}}_{k+1}:={\mathbf{r}}_k-{\alpha }_k{\mathbf{A}\mathbf{p}}_k}$ and${\displaystyle {\mathbf{r}}_{k+1}:=\mathbf{b}-{\mathbf{A}\mathbf{x}}_{k+1}}$ mathematically equivalent. The former is used in the algorithm to avoid an extra multiplication by${\displaystyle \mathbf{A}}$ since the vector${\displaystyle {\mathbf{A}\mathbf{p}}_k}$ is already computed to evaluate${\displaystyle {\alpha }_k}$. The latter may be more accurate, substituting the explicit calculation${\displaystyle {\mathbf{r}}_{k+1}:=\mathbf{b}-{\mathbf{A}\mathbf{x}}_{k+1}}$ for the implicit one by the recursion subject toround-off error accumulation, and is thus recommended for an occasional evaluation.^[7]

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

In the algorithm,${\displaystyle {\alpha }_k}$ is chosen such that${\displaystyle {\mathbf{r}}_{k+1}}$ is orthogonal to${\displaystyle {\mathbf{r}}_k}$. The denominator is simplified from

${\displaystyle {\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}\mathbf{p}}_k}}$

since${\displaystyle {\mathbf{r}}_{k+1}={\mathbf{p}}_{k+1}-{\beta }_k{\mathbf{p}}_k}$. The${\displaystyle {\beta }_k}$ is chosen such that${\displaystyle {\mathbf{p}}_{k+1}}$ is conjugate to${\displaystyle {\mathbf{p}}_k}$. Initially,${\displaystyle {\beta }_k}$ is

${\displaystyle {\beta }_k=-\frac{{\mathbf{r}}_{k+1}^{\mathsf{T}}\mathbf{A}{\mathbf{p}}_k}{{\mathbf{p}}_k^{\mathsf{T}}\mathbf{A}{\mathbf{p}}_k}}$

using

${\displaystyle {\mathbf{r}}_{k+1}={\mathbf{r}}_k-{\alpha }_k\mathbf{A}{\mathbf{p}}_k}$

and equivalently

${\displaystyle \mathbf{A}{\mathbf{p}}_k=\frac{1}{{\alpha }_k}({\mathbf{r}}_k-{\mathbf{r}}_{k+1}),}$

the numerator of${\displaystyle {\beta }_k}$ is rewritten as

${\displaystyle {\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}}$

because${\displaystyle {\mathbf{r}}_{k+1}}$ and${\displaystyle {\mathbf{r}}_k}$ are orthogonal by design. The denominator is rewritten as

${\displaystyle {\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}$

using that the search directions${\displaystyle {\mathbf{p}}_k}$ are conjugated and again that the residuals are orthogonal. This gives the${\displaystyle \beta }$ in the algorithm after cancelling${\displaystyle {\alpha }_k}$.

"""    conjugate_gradient!(A, b, x)Return the solution to `A * x = b` using the conjugate gradient method."""functionconjugate_gradient!(A::AbstractMatrix,b::AbstractVector,x::AbstractVector;tol=eps(eltype(b)))# Initialize residual vectorresidual=b-A*x# Initialize search direction vectorsearch_direction=copy(residual)# Compute initial squared residual normnorm(x)=sqrt(sum(x.^2))old_resid_norm=norm(residual)# Iterate until convergencewhileold_resid_norm>tolA_search_direction=A*search_directionstep_size=old_resid_norm^2/(search_direction'*A_search_direction)# Update solution@.x=x+step_size*search_direction# Update residual@.residual=residual-step_size*A_search_directionnew_resid_norm=norm(residual)# Update search direction vector@.search_direction=residual+(new_resid_norm/old_resid_norm)^2*search_direction# Update squared residual norm for next iterationold_resid_norm=new_resid_normendreturnxend

Consider the linear systemAx =b given by

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

${\displaystyle {\mathbf{x}}_0=\left[\begin{array}{c}2\\ 1\end{array}\right]}$

in order to find an approximate solution to the system.

For reference, the exact solution is

${\displaystyle \mathbf{x}=\left[\begin{array}{c}\frac{1}{11}\\ \\ \frac{7}{11}\end{array}\right]\approx \left[\begin{array}{c}0.0909\\ \\ 0.6364\end{array}\right]}$

Our first step is to calculate the residual vectorr₀ associated withx₀. This residual is computed from the formular₀ =b -Ax₀, and in our case is equal to

Since this is the first iteration, we will use the residual vectorr₀ as our initial search directionp₀; the method of selectingp_k will change in further iterations.

We now compute the scalarα₀ using the relationship

We can now computex₁ using the formula

${\displaystyle {\mathbf{x}}_1={\mathbf{x}}_0+{\alpha }_0{\mathbf{p}}_0=\left[\begin{array}{c}2\\ 1\end{array}\right]+\frac{73}{331}\left[\begin{array}{c}-8\\ -3\end{array}\right]\approx \left[\begin{array}{c}0.2356\\ 0.3384\end{array}\right].}$

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

Our next step in the process is to compute the scalarβ₀ that will eventually be used to determine the next search directionp₁.

Now, using this scalarβ₀, we can compute the next search directionp₁ using the relationship

${\displaystyle {\mathbf{p}}_1={\mathbf{r}}_1+{\beta }_0{\mathbf{p}}_0\approx \left[\begin{array}{c}-0.2810\\ 0.7492\end{array}\right]+0.0088\left[\begin{array}{c}-8\\ -3\end{array}\right]=\left[\begin{array}{c}-0.3511\\ 0.7229\end{array}\right].}$

We now compute the scalarα₁ using our newly acquiredp₁ using the same method as that used forα₀.

Finally, we findx₂ using the same method as that used to findx₁.

${\displaystyle {\mathbf{x}}_2={\mathbf{x}}_1+{\alpha }_1{\mathbf{p}}_1\approx \left[\begin{array}{c}0.2356\\ 0.3384\end{array}\right]+0.4122\left[\begin{array}{c}-0.3511\\ 0.7229\end{array}\right]=\left[\begin{array}{c}0.0909\\ 0.6364\end{array}\right].}$

The result,x₂, is a "better" approximation to the system's solution thanx₁ andx₀. If exact arithmetic were to be used in this example instead of limited-precision, then the exact solution would theoretically have been reached aftern = 2 iterations (n being the order of the system).

The conjugate gradient method can theoretically be viewed as a direct method, as in the absence ofround-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 aniterative method, the conjugate gradient method monotonically (in the energy norm) improves approximations${\displaystyle {\mathbf{x}}_k}$ 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 thecondition number${\displaystyle \kappa (A)}$ of the system matrix${\displaystyle A}$: the larger${\displaystyle \kappa (A)}$ is, the slower the improvement.^[8]

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

Define a subset of polynomials as

${\displaystyle {\mathrm{\Pi }}_k^{\ast }:=\left\{p\in {\mathrm{\Pi }}_k:p(0)=1\right\},}$

where${\displaystyle {\mathrm{\Pi }}_k}$ is the set ofpolynomials of maximal degree${\displaystyle k}$.

Let${\displaystyle {\left({\mathbf{x}}_k\right)}_k}$ be the iterative approximations of the exact solution${\displaystyle {\mathbf{x}}_{\ast }}$, and define the errors as${\displaystyle {\mathbf{e}}_k:={\mathbf{x}}_k-{\mathbf{x}}_{\ast }}$.Now, the rate of convergence can be approximated as^[4][9]

where${\displaystyle \sigma (\mathbf{A})}$ denotes thespectrum, and${\displaystyle \kappa (\mathbf{A})}$ denotes thecondition number.

This shows${\displaystyle k=\frac{1}{2}\sqrt{\kappa (\mathbf{A})}\log \left({\Vert {\mathbf{e}}_0\Vert }_{\mathbf{A}}{\epsilon }^{-1}\right)}$ iterations suffices to reduce the error to${\displaystyle 2\epsilon }$ for any${\displaystyle \epsilon >0}$.

Note, the important limit when${\displaystyle \kappa (\mathbf{A})}$ tends to${\displaystyle \mathrm{\infty }}$

${\displaystyle \frac{\sqrt{\kappa (\mathbf{A})}-1}{\sqrt{\kappa (\mathbf{A})}+1}\approx 1-\frac{2}{\sqrt{\kappa (\mathbf{A})}}\text{for}\kappa (\mathbf{A})\gg 1.}$

This limit shows a faster convergence rate compared to the iterative methods ofJacobi orGauss–Seidel which scale as${\displaystyle \approx 1-\frac{2}{\kappa (\mathbf{A})}}$.

Noround-off error is assumed in the convergence theorem, but the convergence bound is commonly valid in practice as theoretically explained^[5] byAnne Greenbaum.

If initialized randomly, the first stage of iterations is often the fastest, as the error is eliminated within the Krylov subspace that initially reflects a smaller effective condition number. The second stage of convergence is typically well defined by the theoretical convergence bound with$\sqrt{\kappa (\mathbf{A})}$, but may be super-linear, depending on a distribution of the spectrum of the matrix${\displaystyle A}$ and the spectral distribution of the error.^[5] In the last stage, the smallest attainable accuracy is reached and the convergence stalls or the method may even start diverging. In typical scientific computing applications indouble-precision floating-point format for matrices of large sizes, the conjugate gradient method uses a stopping criterion with a tolerance that terminates the iterations during the first or second stage.

In most cases,preconditioning is necessary to ensure fast convergence of the conjugate gradient method. If${\displaystyle {\mathbf{M}}^{-1}}$ is symmetric positive-definite and${\displaystyle {\mathbf{M}}^{-1}\mathbf{A}}$ has a better condition number than${\displaystyle \mathbf{A},}$ a preconditioned conjugate gradient method can be used. It takes the following form:^[10]

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

${\displaystyle \text{Solve:}\mathbf{M}{\mathbf{z}}_0:={\mathbf{r}}_0}$

${\displaystyle {\mathbf{p}}_0:={\mathbf{z}}_0}$

${\displaystyle k:=0}$

repeat

${\displaystyle {\alpha }_k:=\frac{{\mathbf{r}}_k^{\mathsf{T}}{\mathbf{z}}_k}{{\mathbf{p}}_k^{\mathsf{T}}{\mathbf{A}\mathbf{p}}_k}}$

${\displaystyle {\mathbf{x}}_{k+1}:={\mathbf{x}}_k+{\alpha }_k{\mathbf{p}}_k}$

${\displaystyle {\mathbf{r}}_{k+1}:={\mathbf{r}}_k-{\alpha }_k{\mathbf{A}\mathbf{p}}_k}$

ifr_k+1 is sufficiently smallthen exit loopend if

${\displaystyle \mathrm{S}\mathrm{o}\mathrm{l}\mathrm{v}\mathrm{e}\mathbf{M}{\mathbf{z}}_{k+1}:={\mathbf{r}}_{k+1}}$

${\displaystyle {\beta }_k:=\frac{{\mathbf{r}}_{k+1}^{\mathsf{T}}{\mathbf{z}}_{k+1}}{{\mathbf{r}}_k^{\mathsf{T}}{\mathbf{z}}_k}}$

${\displaystyle {\mathbf{p}}_{k+1}:={\mathbf{z}}_{k+1}+{\beta }_k{\mathbf{p}}_k}$

${\displaystyle k:=k+1}$

end repeat

The result isx_k+1

The above formulation is equivalent to applying the regular conjugate gradient method to the preconditioned system^[11]

${\displaystyle {\mathbf{E}}^{-1}\mathbf{A}({\mathbf{E}}^{-1}{)}^{\mathsf{T}}\hat{\mathbf{x}}={\mathbf{E}}^{-1}\mathbf{b}}$

where

${\displaystyle {\mathbf{E}\mathbf{E}}^{\mathsf{T}}=\mathbf{M},\hat{\mathbf{x}}={\mathbf{E}}^{\mathsf{T}}\mathbf{x}.}$

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${\displaystyle {\mathbf{M}}^{-1}}$. It can be shown that${\displaystyle {\mathbf{E}}^{-1}\mathbf{A}({\mathbf{E}}^{-1}{)}^{\mathsf{T}}}$ has the same spectrum as${\displaystyle {\mathbf{M}}^{-1}\mathbf{A}}$.

The preconditioner matrixM 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 usedpreconditioner is theincomplete Cholesky factorization.^[12]

It is important to keep in mind that we don't want to invert the matrix${\displaystyle \mathbf{M}}$ explicitly in order to get${\displaystyle {\mathbf{M}}^{-1}}$ for use it in the process, since inverting${\displaystyle \mathbf{M}}$ 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${\displaystyle \mathbf{L}}$, and the preconditioner matrix is:

${\displaystyle \mathbf{M}={\mathbf{L}\mathbf{L}}^{\mathsf{T}}}$

Then we have to solve:

${\displaystyle \mathbf{M}\mathbf{z}=\mathbf{r}}$

${\displaystyle \mathbf{z}={\mathbf{M}}^{-1}\mathbf{r}}$

But:

${\displaystyle {\mathbf{M}}^{-1}=({\mathbf{L}}^{-1}{)}^{\mathsf{T}}{\mathbf{L}}^{-1}}$

Then:

${\displaystyle \mathbf{z}=({\mathbf{L}}^{-1}{)}^{\mathsf{T}}{\mathbf{L}}^{-1}\mathbf{r}}$

Let's take an intermediary vector${\displaystyle \mathbf{a}}$:

${\displaystyle \mathbf{a}={\mathbf{L}}^{-1}\mathbf{r}}$

${\displaystyle \mathbf{r}=\mathbf{L}\mathbf{a}}$

Since${\displaystyle \mathbf{r}}$ and${\displaystyle \mathbf{L}}$ and known, and${\displaystyle \mathbf{L}}$ is lower triangular, solving for${\displaystyle \mathbf{a}}$ is easy and computationally cheap by usingforward substitution. Then, we substitute${\displaystyle \mathbf{a}}$ in the original equation:

${\displaystyle \mathbf{z}=({\mathbf{L}}^{-1}{)}^{\mathsf{T}}\mathbf{a}}$

${\displaystyle \mathbf{a}={\mathbf{L}}^{\mathsf{T}}\mathbf{z}}$

Since${\displaystyle \mathbf{a}}$ and${\displaystyle {\mathbf{L}}^{\mathsf{T}}}$ are known, and${\displaystyle {\mathbf{L}}^{\mathsf{T}}}$ is upper triangular, solving for${\displaystyle \mathbf{z}}$ is easy and computationally cheap by usingbackward substitution.

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

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 thePolak–Ribière formula

${\displaystyle {\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}}$

instead of theFletcher–Reeves formula

${\displaystyle {\beta }_k:=\frac{{\mathbf{r}}_{k+1}^{\mathsf{T}}{\mathbf{z}}_{k+1}}{{\mathbf{r}}_k^{\mathsf{T}}{\mathbf{z}}_k}}$

may dramatically improve the convergence in this case.^[13] This version of the preconditioned conjugate gradient method can be called^[14]flexible, as it allows for variable preconditioning. The flexible version is also shown^[15] 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,${\displaystyle {\mathbf{r}}_{k+1}^{\mathsf{T}}{\mathbf{z}}_k=0,}$ so both formulas forβ_k are equivalent in exact arithmetic, i.e., without theround-off error.

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

In both the original and the preconditioned conjugate gradient methods one only needs to set${\displaystyle {\beta }_k:=0}$ in order to make them locally optimal, using theline search,steepest descent methods. With this substitution, vectorsp are always the same as vectorsz, so there is no need to store vectorsp. Thus, every iteration of thesesteepest 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-SPDpreconditioner is used, see above.

The conjugate gradient method can also be derived usingoptimal control theory.^[17] In this approach, the conjugate gradient method falls out as anoptimal feedback controller,$${\displaystyle u=k(x,v):=-{\gamma }_a\mathrm{\nabla }f(x)-{\gamma }_{b}v}$$ for thedouble integrator system,$${\displaystyle \dot{x}=v,\dot{v}=u}$$ The quantities${\displaystyle {\gamma }_a}$ and${\displaystyle {\gamma }_b}$ are variable feedback gains.^[17]

The conjugate gradient method can be applied to an arbitraryn-by-m matrix by applying it tonormal equationsA^TA and right-hand side vectorA^Tb, sinceA^TA is a symmetricpositive-semidefinite matrix for anyA. The result isconjugate gradient on the normal equations (CGN orCGNR).

A^TAx =A^Tb

As an iterative method, it is not necessary to formA^TA explicitly in memory but only to perform the matrix–vector and transpose matrix–vector multiplications. Therefore, CGNR is particularly useful whenA is asparse matrix since these operations are usually extremely efficient. However the downside of forming the normal equations is that thecondition number κ(A^TA) is equal to κ²(A) and so the rate of convergence of CGNR may be slow and the quality of the approximate solution may be sensitive to roundoff errors. Finding a goodpreconditioner is often an important part of using the CGNR method.

Several algorithms have been proposed (e.g., CGLS, LSQR). TheLSQR algorithm purportedly has the best numerical stability whenA is ill-conditioned, i.e.,A has a largecondition number.

The conjugate gradient method with a trivial modification is extendable to solving, given complex-valued matrix A and vector b, the system of linear equations${\displaystyle \mathbf{A}\mathbf{x}=\mathbf{b}}$ for the complex-valued vector x, where A isHermitian (i.e., A' = A) andpositive-definite matrix, and the symbol ' denotes theconjugate transpose. The trivial modification is simply substituting theconjugate transpose for the realtranspose everywhere.

The advantages and disadvantages of the conjugate gradient methods are summarized in the lecture notes by Nemirovsky and BenTal.^{[18]: Sec.7.3}

This example is from^[19]Let$t\in (0,1)$, and defineSince${\displaystyle W}$ is invertible, there exists a unique solution to$Wx=b$. Solving it by conjugate gradient descent gives us rather bad convergence:$${\displaystyle \Vert b-W{x}_k{\Vert }^2=(1/t{)}^k,\Vert b-W{x}_n{\Vert }^2=0}$$In words, during the CG process, the error grows exponentially, until it suddenly becomes zero as the unique solution is found.

• Atkinson, Kendell A. (1988). "Section 8.9".An introduction to numerical analysis (2nd ed.). John Wiley and Sons.ISBN 978-0-471-50023-0.
• Avriel, Mordecai (2003).Nonlinear Programming: Analysis and Methods. Dover Publishing.ISBN 978-0-486-43227-4.
• Golub, Gene H.; Van Loan, Charles F. (2013). "Chapter 11".Matrix Computations (4th ed.). Johns Hopkins University Press.ISBN 978-1-4214-0794-4.
• Saad, Yousef (2003-04-01)."Chapter 6".Iterative methods for sparse linear systems (2nd ed.). SIAM.ISBN 978-0-89871-534-7.
• Gérard Meurant: "Detection and correction of silent errors in the conjugate gradient algorithm", Numerical Algorithms, vol.92 (2023), pp.869-891. url=https://doi.org/10.1007/s11075-022-01380-1
• Meurant, Gerard; Tichy, Petr (2024).Error Norm Estimation in the Conjugate Gradient Algorithm. SIAM.ISBN 978-1-61197-785-1.

• "Conjugate gradients, method of",Encyclopedia of Mathematics,EMS Press, 2001 [1994]

• Numerical linear algebra
• Gradient methods

• CS1 German-language sources (de)
• CS1 maint: multiple names: authors list
• Articles with short description
• Short description is different from Wikidata
• Articles with example MATLAB/Octave code