The Conjugate Gradient Method
The optimal short-recurrence Krylov solver for symmetric positive-definite systems, minimizing energy error with constant memory per step.
The idea
Conjugate gradient (CG) solves Ax = b when A is symmetric positive definite. It can be viewed as minimizing the quadratic energy function (1/2) x^T A x minus b^T x, whose minimum is exactly the solution. Rather than steepest descent, which zigzags, CG searches along directions that are A-conjugate (orthogonal in the inner product defined by A), so it never undoes progress from earlier steps.
The short recurrence
Because A is symmetric, the Lanczos process gives a three-term recurrence, and CG inherits it. Each iteration needs one matrix-vector product, a couple of inner products, and vector updates. Memory is constant: only a few vectors are stored regardless of how many iterations run. This is why CG scales to enormous problems.
import numpy as np
def cg(A, b, x0, tol=1e-10, maxit=1000):
x = x0.copy()
r = b - A @ x
p = r.copy()
rs = r @ r
for k in range(maxit):
Ap = A @ p
alpha = rs / (p @ Ap)
x += alpha * p
r -= alpha * Ap
rs_new = r @ r
if np.sqrt(rs_new) < tol:
break
p = r + (rs_new / rs) * p
rs = rs_new
return x, k
Convergence
In exact arithmetic CG converges in at most n steps for an n-by-n system, but its practical value is that it converges much faster when the eigenvalues are clustered. A standard bound shows the error decreases at a rate governed by the square root of the condition number: reducing the condition number by preconditioning directly accelerates convergence.
Preconditioned CG
In practice one solves M^(-1) A x = M^(-1) b with a symmetric positive-definite preconditioner M that approximates A but is cheap to invert. Incomplete Cholesky and algebraic multigrid are common choices. Preconditioned CG is the default solver for the elliptic pressure and potential equations that appear throughout plasma and fluid codes.