Chapter 91. Conjugate Gradient Method
Chapter 91. Conjugate Gradient Method
The conjugate gradient method is an iterative method for solving linear systems
$$ Ax=b $$
where (A) is symmetric positive definite. It is one of the most important Krylov subspace methods in numerical linear algebra.
The method combines three ideas:
| Idea | Role |
|---|---|
| Quadratic minimization | interprets (Ax=b) as an optimization problem |
| Conjugate directions | avoids undoing previous progress |
| Krylov subspaces | builds approximations using matrix-vector products |
For large sparse symmetric positive definite systems, conjugate gradient is often the standard first method to try. It avoids matrix factorization and requires mainly matrix-vector products, vector updates, and inner products. In exact arithmetic, it reaches the exact solution in at most (n) steps for an (n \times n) system, but in floating point arithmetic it is used as an iterative method and is usually stopped when the residual is small enough.
91.1 Symmetric Positive Definite Systems
The conjugate gradient method applies to systems
$$ Ax=b $$
where (A) is symmetric positive definite.
Symmetric means
$$ A^T=A. $$
Positive definite means
$$ x^TAx>0 $$
for every nonzero vector (x).
These assumptions are essential. They imply that (A) has positive eigenvalues, defines an inner product, and determines a strictly convex quadratic function.
Many applied problems produce symmetric positive definite matrices, including finite difference methods, finite element methods, least squares normal equations, graph Laplacian reductions, and quadratic optimization problems.
91.2 Linear Systems as Quadratic Minimization
If (A) is symmetric positive definite, solving
$$ Ax=b $$
is equivalent to minimizing the quadratic function
$$ \phi(x) = \frac{1}{2}x^TAx-b^Tx. $$
The gradient is
$$ \nabla \phi(x)=Ax-b. $$
Thus the stationary condition
$$ \nabla \phi(x)=0 $$
is exactly
$$ Ax=b. $$
Since (A) is positive definite, the quadratic function is strictly convex. Therefore it has a unique minimizer, and that minimizer is the unique solution of the linear system.
91.3 Residual and Gradient
At an approximate solution (x^{(k)}), define the residual
$$ r^{(k)}=b-Ax^{(k)}. $$
The gradient of (\phi) is
$$ \nabla\phi(x^{(k)})=Ax^{(k)}-b. $$
Therefore,
$$ r^{(k)}=-\nabla\phi(x^{(k)}). $$
The residual is the negative gradient of the quadratic objective. It points in the steepest descent direction.
Steepest descent uses the residual direction at each step. Conjugate gradient improves on steepest descent by choosing directions that are mutually conjugate with respect to (A).
91.4 The Energy Inner Product
For symmetric positive definite (A), define the (A)-inner product by
$$ \langle u,v\rangle_A=u^TAv. $$
The associated norm is
$$ |u|_A=\sqrt{u^TAu}. $$
This norm is called the energy norm.
Two vectors (p) and (q) are (A)-conjugate if
$$ p^TAq=0. $$
This is orthogonality measured in the geometry induced by (A), rather than ordinary Euclidean geometry.
91.5 Conjugate Directions
Suppose we have nonzero directions
$$ p^{(0)},p^{(1)},\ldots,p^{(m-1)} $$
that are mutually (A)-conjugate:
$$ (p^{(i)})^TAp^{(j)}=0 \qquad i\ne j. $$
If we minimize (\phi) successively along these directions, each step does not interfere with the progress made in previous directions.
This is the key advantage over steepest descent.
Steepest descent may zigzag across elongated level sets. Conjugate directions account for the geometry of the quadratic surface and avoid repeated correction in the same effective direction.
91.6 One-Dimensional Minimization
Given a current approximation (x^{(k)}) and a search direction (p^{(k)}), choose
$$ x^{(k+1)}=x^{(k)}+\alpha_k p^{(k)}. $$
We choose (\alpha_k) to minimize
$$ \phi(x^{(k)}+\alpha p^{(k)}). $$
Differentiate with respect to (\alpha):
$$ \frac{d}{d\alpha} \phi(x^{(k)}+\alpha p^{(k)}) = (p^{(k)})^T(Ax^{(k)}-b) + \alpha (p^{(k)})^TAp^{(k)}. $$
Since
$$ r^{(k)}=b-Ax^{(k)}, $$
setting the derivative to zero gives
$$ \alpha_k = \frac{(p^{(k)})^Tr^{(k)}}{(p^{(k)})^TAp^{(k)}}. $$
In the standard conjugate gradient recurrence, this becomes
$$ \alpha_k = \frac{(r^{(k)})^Tr^{(k)}}{(p^{(k)})^TAp^{(k)}}. $$
91.7 Residual Update
After choosing (\alpha_k), the new approximation is
$$ x^{(k+1)} = x^{(k)}+\alpha_kp^{(k)}. $$
The residual updates as
$$ r^{(k+1)} = b-Ax^{(k+1)}. $$
Substitute the update for (x^{(k+1)}):
$$ r^{(k+1)} = b-A(x^{(k)}+\alpha_kp^{(k)}) = b-Ax^{(k)}-\alpha_kAp^{(k)}. $$
Therefore,
$$ r^{(k+1)} = r^{(k)}-\alpha_kAp^{(k)}. $$
This update avoids recomputing (Ax^{(k+1)}) from scratch.
91.8 Constructing the Next Direction
The next search direction is formed from the new residual and the previous direction:
$$ p^{(k+1)} = r^{(k+1)}+\beta_kp^{(k)}. $$
The coefficient (\beta_k) is chosen so that
$$ (p^{(k+1)})^TAp^{(k)}=0. $$
For the standard conjugate gradient method,
$$ \beta_k = \frac{(r^{(k+1)})^Tr^{(k+1)}}{(r^{(k)})^Tr^{(k)}}. $$
Thus the method stores only the current approximation, residual, search direction, and one matrix-vector product.
91.9 The Basic Algorithm
Choose an initial vector (x^{(0)}).
Set
$$ r^{(0)}=b-Ax^{(0)}, \qquad p^{(0)}=r^{(0)}. $$
For
$$ k=0,1,2,\ldots $$
compute
$$ \alpha_k = \frac{(r^{(k)})^Tr^{(k)}}{(p^{(k)})^TAp^{(k)}}. $$
Then update
$$ x^{(k+1)} = x^{(k)}+\alpha_kp^{(k)}. $$
Update the residual:
$$ r^{(k+1)} = r^{(k)}-\alpha_kAp^{(k)}. $$
If the residual is small enough, stop.
Otherwise compute
$$ \beta_k = \frac{(r^{(k+1)})^Tr^{(k+1)}}{(r^{(k)})^Tr^{(k)}}. $$
Then update
$$ p^{(k+1)} = r^{(k+1)}+\beta_kp^{(k)}. $$
This is the standard three-term recurrence form of conjugate gradient.
91.10 Pseudocode
conjugate_gradient(A, b, x, tol, max_iter):
r = b - A * x
p = r
rsold = dot(r, r)
for k = 0 to max_iter - 1:
Ap = A * p
alpha = rsold / dot(p, Ap)
x = x + alpha * p
r = r - alpha * Ap
rsnew = dot(r, r)
if sqrt(rsnew) / norm(b) <= tol:
return x
beta = rsnew / rsold
p = r + beta * p
rsold = rsnew
return x
Each iteration requires:
| Operation | Count |
|---|---|
| Matrix-vector product | 1 |
| Dot products | 2 |
| Vector updates | several |
| Stored vectors | few |
The dominant cost for sparse matrices is usually the matrix-vector product.
91.11 Orthogonality of Residuals
In exact arithmetic, the residuals produced by conjugate gradient are mutually orthogonal:
$$ (r^{(i)})^Tr^{(j)}=0 \qquad i\ne j. $$
This means each iteration removes an error component in a new independent direction.
Since there can be at most (n) mutually orthogonal nonzero vectors in (\mathbb{R}^n), the method reaches zero residual in at most (n) steps in exact arithmetic. This is the finite termination property.
91.12 Conjugacy of Search Directions
The search directions are mutually (A)-conjugate:
$$ (p^{(i)})^TAp^{(j)}=0 \qquad i\ne j. $$
This is stronger than ordinary independence.
It means the directions are orthogonal under the energy inner product.
Because the objective function has Hessian (A), (A)-conjugacy is exactly the orthogonality needed to minimize the quadratic efficiently.
91.13 Krylov Subspace Interpretation
The (k)-th Krylov subspace is
$$ \mathcal{K}_k(A,r^{(0)}) = \operatorname{span}{r^{(0)},Ar^{(0)},A^2r^{(0)},\ldots,A^{k-1}r^{(0)}}. $$
The conjugate gradient iterate satisfies
$$ x^{(k)} \in x^{(0)}+\mathcal{K}_k(A,r^{(0)}). $$
Among all vectors in this affine Krylov space, (x^{(k)}) minimizes the error in the energy norm:
$$ |x-x^{(k)}|_A. $$
This Krylov interpretation explains why the method only needs matrix-vector products.
91.14 Polynomial Approximation View
Since
$$ e^{(k)}=x-x^{(k)}, $$
the error after (k) steps can be written as
$$ e^{(k)}=q_k(A)e^{(0)} $$
for a polynomial (q_k) of degree (k) satisfying
$$ q_k(0)=1. $$
The method chooses a polynomial that makes the error small over the eigenvalues of (A).
Thus convergence depends not only on the largest and smallest eigenvalues, but also on the distribution of all eigenvalues.
Clustered eigenvalues often lead to rapid convergence.
91.15 Convergence Bound
For symmetric positive definite (A), a standard bound is
$$ \frac{|x-x^{(k)}|_A}{|x-x^{(0)}|_A} \le 2 \left( \frac{\sqrt{\kappa_2(A)}-1} {\sqrt{\kappa_2(A)}+1} \right)^k. $$
Here
$$ \kappa_2(A) = \frac{\lambda_{\max}(A)}{\lambda_{\min}(A)} $$
is the spectral condition number.
This bound shows that convergence becomes slower as the condition number grows. It also explains why conjugate gradient is often much faster than Jacobi or Gauss-Seidel for well-preconditioned symmetric positive definite systems.
91.16 Effect of Eigenvalue Distribution
The condition number gives a useful worst-case estimate, but it does not tell the full story.
Conjugate gradient may converge rapidly when eigenvalues are clustered, even if the condition number is large.
Intuitively, each iteration builds a polynomial that is small at the eigenvalues of (A). If many eigenvalues lie in clusters, a low-degree polynomial can be small on most of them.
Thus two matrices with the same condition number may produce different convergence histories.
91.17 Preconditioned Conjugate Gradient
Preconditioning replaces the original system by an equivalent one that has better spectral properties.
Choose a symmetric positive definite matrix (M) that approximates (A) and is easy to solve with.
The preconditioned method uses
$$ Mz^{(k)}=r^{(k)}. $$
The vector (z^{(k)}) is the preconditioned residual.
The goal is to make the effective system better conditioned. In practice, preconditioning is often necessary for fast convergence of conjugate gradient.
91.18 Preconditioned Algorithm
Choose (x^{(0)}).
Set
$$ r^{(0)}=b-Ax^{(0)}. $$
Solve
$$ Mz^{(0)}=r^{(0)}. $$
Set
$$ p^{(0)}=z^{(0)}. $$
For (k=0,1,2,\ldots):
$$ \alpha_k = \frac{(r^{(k)})^Tz^{(k)}}{(p^{(k)})^TAp^{(k)}}. $$
$$ x^{(k+1)} = x^{(k)}+\alpha_kp^{(k)}. $$
$$ r^{(k+1)} = r^{(k)}-\alpha_kAp^{(k)}. $$
Solve
$$ Mz^{(k+1)}=r^{(k+1)}. $$
$$ \beta_k = \frac{(r^{(k+1)})^Tz^{(k+1)}}{(r^{(k)})^Tz^{(k)}}. $$
$$ p^{(k+1)} = z^{(k+1)}+\beta_kp^{(k)}. $$
91.19 Common Preconditioners
Common preconditioners include:
| Preconditioner | Main idea |
|---|---|
| Jacobi | use the diagonal of (A) |
| Incomplete Cholesky | approximate Cholesky factorization |
| SSOR | symmetric relaxation-based preconditioner |
| Multigrid | solve error components across scales |
| Domain decomposition | split problem into subdomains |
| Algebraic multigrid | build hierarchy from matrix structure |
A good preconditioner must balance two costs:
| Requirement | Meaning |
|---|---|
| Strong approximation | reduces iteration count |
| Cheap application | keeps each iteration inexpensive |
A preconditioner that is too expensive may erase the benefit of faster convergence.
91.20 Floating Point Behavior
In exact arithmetic, residuals remain mutually orthogonal and directions remain (A)-conjugate.
In floating point arithmetic, these properties gradually degrade.
Consequences include:
| Exact arithmetic property | Floating point behavior |
|---|---|
| finite termination in at most (n) steps | rarely exact termination |
| mutually orthogonal residuals | orthogonality is gradually lost |
| mutually (A)-conjugate directions | conjugacy is gradually lost |
| exact residual recurrence | recursive residual may drift |
For this reason, implementations sometimes recompute the true residual
$$ b-Ax^{(k)} $$
occasionally, rather than relying only on the recurrence.
91.21 Stopping Criteria
Common stopping rules use the relative residual:
$$ \frac{|r^{(k)}|_2}{|b|_2}\le \tau. $$
Other choices include:
| Criterion | Form |
|---|---|
| relative residual | (|r^{(k)}|/|b|) |
| preconditioned residual | ((r^{(k)})^Tz^{(k)}) |
| estimated error | bound using condition estimate |
| maximum iterations | fixed safety limit |
| stagnation | residual no longer decreases |
A small residual must still be interpreted with conditioning. For ill-conditioned systems, a small residual may correspond to a larger solution error.
91.22 Computational Cost
For a sparse matrix with (\operatorname{nnz}(A)) nonzero entries, one conjugate gradient iteration costs approximately
$$ O(\operatorname{nnz}(A)). $$
The method stores only a small number of vectors:
| Vector | Meaning |
|---|---|
| (x) | approximate solution |
| (r) | residual |
| (p) | search direction |
| (Ap) | matrix-vector product |
| (z) | preconditioned residual, if used |
This low storage cost makes the method suitable for very large sparse systems.
91.23 Comparison with Steepest Descent
Steepest descent uses
$$ p^{(k)}=r^{(k)} $$
at every step.
Conjugate gradient uses
$$ p^{(k+1)}=r^{(k+1)}+\beta_kp^{(k)}. $$
| Feature | Steepest descent | Conjugate gradient |
|---|---|---|
| Direction | residual only | residual plus previous direction |
| Geometry | Euclidean steepest direction | (A)-conjugate directions |
| Typical convergence | slow on elongated quadratics | much faster |
| Memory | very low | very low |
| Matrix requirement | SPD for standard form | SPD |
Conjugate gradient can be viewed as steepest descent with memory adjusted to the quadratic geometry.
91.24 Comparison with Direct Methods
Direct methods factor the matrix.
Conjugate gradient does not.
| Feature | Direct factorization | Conjugate gradient |
|---|---|---|
| Matrix type | many types | SPD only |
| Work | predictable | depends on convergence |
| Memory | may be large due to fill-in | low |
| Sparse systems | fill-in may dominate | often efficient |
| Accuracy | high if stable | depends on tolerance and conditioning |
| Reuse for many right-hand sides | strong | weaker unless preconditioner reused |
For a single large sparse SPD system, conjugate gradient is often preferable. For many right-hand sides with the same matrix, a direct factorization may become competitive.
91.25 Failure Modes
Conjugate gradient may fail or perform poorly when its assumptions are violated or when convergence is too slow.
| Failure mode | Cause |
|---|---|
| denominator nonpositive | (A) not positive definite |
| slow convergence | large condition number |
| stagnation | poor preconditioner or roundoff |
| loss of conjugacy | floating point error |
| residual drift | recursive residual differs from true residual |
| breakdown in PCG | preconditioner not SPD |
If (A) is symmetric indefinite, MINRES is often more appropriate. If (A) is nonsymmetric, GMRES or BiCGSTAB is usually used.
91.26 Practical Checklist
Before using conjugate gradient, check:
| Question | Reason |
|---|---|
| Is (A) symmetric? | required by standard CG |
| Is (A) positive definite? | required for positive denominators and convexity |
| Is (A) sparse or matrix-free? | favors CG |
| Is a good preconditioner available? | controls iteration count |
| Is the stopping tolerance meaningful? | residual does not equal error |
| Is the matrix scaled reasonably? | improves numerical behavior |
For many applications, the most important implementation decision is the preconditioner.
91.27 Summary
The conjugate gradient method solves symmetric positive definite systems
$$ Ax=b $$
by minimizing the quadratic
$$ \phi(x)=\frac{1}{2}x^TAx-b^Tx $$
over expanding Krylov subspaces.
Its basic recurrence is:
$$ r^{(0)}=b-Ax^{(0)}, \qquad p^{(0)}=r^{(0)}. $$
$$ \alpha_k = \frac{(r^{(k)})^Tr^{(k)}}{(p^{(k)})^TAp^{(k)}}. $$
$$ x^{(k+1)}=x^{(k)}+\alpha_kp^{(k)}. $$
$$ r^{(k+1)}=r^{(k)}-\alpha_kAp^{(k)}. $$
$$ \beta_k = \frac{(r^{(k+1)})^Tr^{(k+1)}}{(r^{(k)})^Tr^{(k)}}. $$
$$ p^{(k+1)}=r^{(k+1)}+\beta_kp^{(k)}. $$
In exact arithmetic, it terminates in at most (n) steps. In practice, it is used as a fast iterative method. Its effectiveness depends strongly on conditioning, eigenvalue distribution, and preconditioning.