The Conjugate Gradient Method

Conjugate gradient solves symmetric positive-definite systems by minimizing a quadratic along mutually conjugate search directions, using only matrix-vector products.

Minimization view of a linear system

For a symmetric positive-definite matrix A, solving A x = b is equivalent to minimizing the quadratic (1/2) x^T A x - b^T x, whose gradient is the residual A x - b. The conjugate gradient method (CG) descends this quadratic, but instead of steepest descent it chooses A-conjugate directions that do not undo previous progress.

import numpy as np
def cg(A, b, x, tol=1e-10, nmax=1000):
    r = b - A @ x; p = r.copy(); rs = r @ r
    for _ in range(nmax):
        Ap = A @ p
        alpha = rs / (p @ Ap)
        x += alpha*p; r -= alpha*Ap
        rs_new = r @ r
        if rs_new**0.5 < tol: break
        p = r + (rs_new/rs)*p; rs = rs_new
    return x

Why it is efficient

Each iteration needs just one matrix-vector product and a few vector operations, so it never forms or factors A and exploits sparsity fully. In exact arithmetic CG converges in at most n steps, but in practice it reaches high accuracy far sooner, especially when A's eigenvalues are clustered.

Convergence and conditioning

The convergence rate depends on the condition number of A: the error after k steps shrinks by roughly ((sqrt(kappa)-1)/(sqrt(kappa)+1))^k. Ill-conditioned systems converge slowly, which is why CG is nearly always paired with a preconditioner that clusters the eigenvalues and cuts the iteration count sharply.

Scope

CG requires A to be symmetric positive definite. For symmetric indefinite systems, MINRES applies; for general nonsymmetric systems, GMRES or BiCGSTAB are used. Within its domain, preconditioned CG is the standard method for the large sparse SPD systems from discretized elliptic PDEs.

Preconditioned CG solves the symmetric positive-definite systems arising from potential and diffusion operators in breeder Hyperion field models.