Conjugate Gradient Method

The conjugate gradient method (CG) is an iterative Krylov subspace algorithm in computation and algorithms and numerical linear algebra for solving large, sparse systems of linear equations Ax = b whose coefficient matrix A ∈ R^n×n is real, symmetric (A^T = A), and positive-definite (x^TAx > 0 for all x ≠ 0). Developed independently and published jointly in 1952 by American mathematician Magnus R. Hestenes and Swiss mathematician Eduard Stiefel at the National Bureau of Standards, the algorithm bridges direct factorization methods (such as Cholesky Decomposition) and first-order iterative methods (such as Gradient Descent). In infinite-precision arithmetic, conjugate gradient reaches the analytic solution in at most n steps; in finite-precision scientific computing—following J. K. Reid's 1971 demonstration of its utility as an iterative solver—it computes high-accuracy approximations in k ≪ n iterations using only one sparse matrix-vector multiplication Ap_k and O(n) vector storage per step.

The Energy Norm and Steepest Descent Zig-Zagging

When A ∈ R^n×n is symmetric positive-definite (SPD), solving Ax = b is equivalent to finding the unique global minimum x^* = A⁻¹b of the strongly convex quadratic potential function:

ϕ(x) = (1/2) x^T A x − b^T x

Taking the gradient yields ∇ϕ(x) = Ax − b = −r(x), where r(x) = b − Ax is the linear residual vector. Minimizing ϕ(x) is equivalent to minimizing the A-norm (energy norm) of the error vector e_k = x_k − x^*:

‖e_k‖_A² = (x_k − x^*)^T A (x_k − x^*) = 2 [ ϕ(x_k) − ϕ(x^*) ]

In classical steepest descent, each iteration steps along the residual r_k with line-search step size α_k = (r_k^Tr_k)/(r_k^TAr_k). Because consecutive steepest-descent residuals are orthogonal in the Euclidean inner product (r_k+1^Tr_k = 0), steepest descent repeatedly steps back into subspaces it has already searched, zig-zagging across elongated ellipsoidal contours at a slow rate governed by (κ − 1)/(κ + 1), where κ = λ_max(A)/λ_min(A) is the spectral condition number.

A-Orthogonal (Conjugate) Search Directions

Hestenes and Stiefel eliminated this zig-zagging by requiring the search directions {p₀, p₁, ..., p_n−1} to be mutually A-orthogonal (or A-conjugate):

p_i^T A p_j = 0 for all i ≠ j

Under the coordinate transformation y = A^1/2x that maps the ellipsoidal level sets of ϕ(x) into spheres, A-conjugate vectors become orthogonal axes. Minimizing ϕ(x) sequentially along n mutually A-conjugate directions reduces the n-dimensional problem to n independent one-dimensional minimizations without spoiling optimality along previously searched directions.

Short Two-Term Recurrences

Because A is symmetric and each iterate x_k minimizes ‖x − x^*‖_A over the k-th Krylov subspace:

K_k(A, r₀) = span{ r₀, Ar₀, A²r₀, ..., A^k−1r₀ }

the new residual r_k+1 is automatically orthogonal to all preceding residuals (r_k+1^Tr_i = 0 for i ≤ k), and the new A-conjugate direction p_k+1 only requires orthogonalization against the preceding direction p_k.

Given an initial guess x₀ ∈ Rⁿ, initialize r₀ = b − Ax₀ and p₀ = r₀. For k = 0, 1, 2, ... until ‖r_k‖₂ / ‖b‖₂ ≤ ε:

1. Optimal Step Length: α_k = (r_k^T r_k) / (p_k^T A p_k)
2. State & Residual Update: x_k+1 = x_k + α_k p_k, and r_k+1 = r_k − α_k A p_k
3. Conjugacy Coefficient: β_k = (r_k+1^T r_k+1) / (r_k^T r_k)
4. Next Search Direction: p_k+1 = r_k+1 + β_k p_k

Each iteration requires one sparse matrix-vector product Ap_k, two inner products, and three vector updates, storing only four vectors of length n.

Square-Root Condition Number Acceleration

Writing the error at step k as e_k = P_k(A)e₀ over polynomials of degree k with P_k(0) = 1 and bounding the discrete spectrum [λ₁, λ_n] by a scaled Chebyshev polynomial of the first kind T_k yields the error bound:

‖x_k − x^*‖_A ≤ 2 [ (√κ − 1) / (√κ + 1) ]^k ‖x₀ − x^*‖_A, where κ = λ_n / λ₁

Replacing κ in steepest descent with √κ in conjugate gradient reduces the iteration count needed for accuracy ε from O(κ log(1/ε)) to O(√κ log(1/ε)). Moreover, if A has only m < n distinct eigenvalues, CG converges in at most m steps.

Preconditioned Conjugate Gradient (PCG)

To accelerate convergence when κ(A) ≫ 1, Preconditioned Conjugate Gradient (PCG) introduces an easily invertible SPD preconditioner M ≈ A (such as Incomplete Cholesky M = L̃L̃^T or Algebraic Multigrid [AMG]), solving Mz_k = r_k at each step and replacing r_k^Tr_k with r_k^Tz_k so convergence scales with √κ(M⁻¹A) ≪ √κ(A).

Fletcher–Reeves and Polak–Ribière Formulas

In 1964, Roger Fletcher and Colin M. Reeves extended CG to smooth nonlinear objectives f(x): Rⁿ → R (such as potential energy minimization in protein folding and molecular dynamics):

• Fletcher–Reeves (1964): β_k^FR = ‖∇f(x_k+1)‖₂² / ‖∇f(x_k)‖₂²
• Polak–Ribière (1969): β_k^PR = [ ∇f(x_k+1)^T (∇f(x_k+1) − ∇f(x_k)) ] / ‖∇f(x_k)‖₂²

Whereas both formulas coincide on quadratic functions under analytic line searches, Polak–Ribière automatically resets β_k → 0 whenever consecutive gradients become nearly equal (x_k+1 ≈ x_k), restarting along the steepest-descent direction −∇f(x_k+1).

Why is Conjugate Gradient used as an iterative method if it finishes in n steps in theory?

For sparse 3D PDE grids with n = 10⁷ to 10⁹ unknowns, taking n steps or storing a dense Cholesky Decomposition is infeasible; Preconditioned Conjugate Gradient reaches machine accuracy in k ≈ 50–500 iterations (k ≪ n).

What happens if Conjugate Gradient is applied to a non-symmetric or indefinite matrix?

If A^T ≠ A, the short two-term recurrence loses orthogonality over earlier Krylov vectors; if A is symmetric indefinite, p_k^TAp_k can vanish or turn negative. Symmetric indefinite systems use MINRES, while non-symmetric systems use GMRES or BiCGSTAB.

How does Incomplete Cholesky preconditioning speed up Conjugate Gradient?

Incomplete Cholesky (IC(0)) factors A ≈ L̃L̃^T while dropping fill-in outside the sparsity pattern of A, preserving O(nnz(A)) memory while clustering the eigenvalues of (L̃L̃^T)⁻¹A near 1.

1. Hestenes, M. R., & Stiefel, E. (1952). Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 49(6), 409–436. https://doi.org/10.6028/jres.049.044
2. Fletcher, R., & Reeves, C. M. (1964). Function minimization by conjugate gradients. The Computer Journal, 7(2), 149–154. https://doi.org/10.1093/comjnl/7.2.149
3. Polak, E., & Ribière, G. (1969). Note sur la convergence de méthodes de directions conjuguées. Revue Française d'Informatique et de Recherche Opérationnelle, Série Rouge, 3(16), 35–43. https://doi.org/10.1051/m2an/196903R100351
4. Meijerink, J. A., & van der Vorst, H. A. (1977). An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix. Mathematics of Computation, 31(137), 148–162. https://doi.org/10.1090/S0025-5718-1977-0438681-4
5. Saad, Y. (2003). Iterative Methods for Sparse Linear Systems (2nd ed.). Society for Industrial and Applied Mathematics (SIAM). https://doi.org/10.1137/1.9780898718003