The Conjugate Gradient Method

Solve a symmetric positive-definite system in at most n steps by searching along mutually conjugate directions.

Beyond gradient descent

Steepest descent on an ill-conditioned quadratic zig-zags and converges slowly. Conjugate gradient (CG) fixes this by choosing search directions that are A-orthogonal (conjugate), so progress in one direction is never undone by later steps. In exact arithmetic it reaches the solution of an n-by-n system in at most n steps.

import numpy as np
A=np.array([[4.0,1.0],[1.0,3.0]]); b=np.array([1.0,2.0])
x=np.zeros(2); r=b-A@x; p=r.copy(); rs=r@r
for i in range(10):
    Ap=A@p; alpha=rs/(p@Ap)
    x=x+alpha*p; r=r-alpha*Ap
    rs_new=r@r
    if rs_new<1e-20: break
    p=r+(rs_new/rs)*p; rs=rs_new
print(np.round(x,6),'in',i+1,'steps')

How it works

Each step moves the optimal distance alpha along the current direction p, then builds the next direction from the new residual plus a correction that keeps it conjugate to all previous directions. Remarkably this only needs the last direction and residual - no long history - thanks to the properties of the underlying Krylov subspace.

Convergence in practice

Rounding means CG rarely finishes in exactly n steps, but it converges fast when eigenvalues are clustered. The error after k steps is bounded by a factor involving sqrt(kappa), a big improvement over gradient descent's kappa. Preconditioning - solving M^-1 A x = M^-1 b with M approximating A - clusters the spectrum and is essential for hard problems.

Why it dominates

CG needs only matrix-vector products and a handful of vectors of storage, so it solves the enormous sparse SPD systems from discretized PDEs that direct methods cannot touch. It is the default inner solver in implicit fluid, structural, and plasma-equilibrium codes.