The Conjugate Gradient Method
Solving Ax = b at massive scale, beating gradient descent's zig-zag.
On this page
Beginner: section 1.17's gradient descent diagram showed a real weakness: on a stretched, poorly-conditioned bowl, it zig-zags back and forth instead of heading straight for the minimum, needing many small steps. The conjugate gradient (CG) method is a smarter iterative algorithm for exactly this kind of problem — solving Ax = b (equivalently, minimizing a quadratic form, section 1.16) — that avoids this zig-zagging almost entirely.
Intermediate: the key idea is choosing each step direction to be conjugate (a specific kind of "A-orthogonal") to every previous direction, rather than just following the current negative gradient like plain gradient descent does. This guarantees CG never "undoes" progress it already made along an earlier direction.
Advanced: for a quadratic problem in n dimensions, this guarantee is exact and remarkable: conjugate gradient reaches the true minimum in at most n steps, in exact arithmetic — regardless of how poorly conditioned the problem is. In practice, on large, sparse, symmetric positive-definite systems (section 1.15), CG is run for far fewer than n iterations and still gets an excellent approximate answer, which is exactly why it's the standard method for enormous linear systems where forming or inverting the full matrix is completely infeasible.
Each new search direction d_k is built to be A-conjugate to every direction used so far — the algorithm only ever needs one matrix-vector product with A per iteration, never a full matrix inversion or decomposition.
Solving Ax = b for symmetric positive-definite A is equivalent to minimizing the quadratic form φ(x) = ½xᵀAx − bᵀx (section 1.16), since its gradient ∇φ(x) = Ax − b is exactly the residual, zero precisely at the solution. Given a current point xₖ and search direction dₖ, an exact line search picks the step size αₖ that minimizes φ(xₖ + αdₖ) along that one direction. Expand and differentiate with respect to α, then set to zero:
where rₖ = Axₖ − b is the current residual. Solving for α:
(the last simplification uses that dₖ is built to equal −rₖ plus a component along previous, A-conjugate directions, which vanishes against rₖ by the conjugacy property itself). This is exactly the alpha line in the code below — not a heuristic, but the exact 1D minimizer along the current search direction, which is precisely why each CG step never needs to backtrack or retry: it's already optimal for that direction by construction.
Where this is used: this exact-line-search argument is what distinguishes CG from generic gradient descent with a guessed learning rate — every step size in CG is derived, not tuned, which is a large part of why it needs no hyperparameter search to work well.
Red dashed = gradient descent, needing dozens of small zig-zagging steps. Green solid = conjugate gradient, reaching the exact minimum in at most 2 steps for this 2D problem.
This is a genuine, working large-scale linear solver in about a dozen lines — the same algorithm scales to systems with millions of variables when A is sparse (section 1.22).
import numpy as np
def conjugate_gradient(A, b, x0, tol=1e-8, max_iter=None):
x = x0.copy()
r = b - A @ x
d = r.copy()
max_iter = max_iter or len(b)
for _ in range(max_iter):
if np.linalg.norm(r) < tol:
break
Ad = A @ d
alpha = (r @ r) / (d @ Ad)
x = x + alpha * d
r_new = r - alpha * Ad
beta = (r_new @ r_new) / (r @ r)
d = r_new + beta * d
r = r_new
return x
n = 200
Q = np.random.randn(n, n)
A = Q.T @ Q + n * np.eye(n) # a large, well-behaved symmetric PD system (section 1.15)
b = np.random.randn(n)
x = conjugate_gradient(A, b, np.zeros(n))
print(np.allclose(A @ x, b, atol=1e-4)) # True, without ever forming A^-1- Natural gradient descent and other second-order optimization methods use CG to approximately solve the linear system involving the Fisher information matrix or Hessian (section 1.16), combined with Hessian-vector products (section 1.32) to avoid ever forming that matrix explicitly.
- Approximate Gaussian process inference uses CG to solve the large linear systems involving the kernel matrix (section 1.25) that exact Cholesky-based inference (section 1.15) would find too expensive at scale.
- Physics simulation and finite-element analysis rely on CG (and its preconditioned variants) as the standard workhorse for solving enormous sparse linear systems.
- Applying plain conjugate gradient to a non-symmetric or non-positive-definite matrix — the method (in its basic form) specifically requires A to be symmetric positive definite; related variants (e.g. GMRES, BiCGSTAB) exist for more general matrices.
- Running CG on a poorly conditioned system without preconditioning — while CG handles ill-conditioning far better than plain gradient descent, extremely large condition numbers still slow convergence in practice, and a good preconditioner is standard practice for genuinely large problems.
Going deeper
The "at most n steps" guarantee is a statement about exact arithmetic — in real floating-point computation, rounding error can slowly degrade the conjugacy property over many iterations, which is one reason CG is usually run as an iterative approximate method (stopping early once the residual is small enough) rather than insisting on running exactly n full steps.
At the master level: preconditioned conjugate gradient (PCG) transforms the system into an equivalent one with a much better condition number before running CG, typically using an approximate, cheap-to-invert version of A (like its diagonal, or an incomplete Cholesky factorization, section 1.15) — this is what makes CG practical on the genuinely huge, ill-conditioned systems that appear in real scientific computing and large-scale ML.
Why can conjugate gradient solve an n-dimensional quadratic problem exactly in at most n steps, while gradient descent generally cannot?
Each CG step direction is constructed to be A-conjugate to every previous direction, so progress made along one direction is never undone by a later step — after n conjugate directions, every dimension has been fully accounted for exactly. Gradient descent instead always follows the current gradient, which can repeatedly re-cross the same directions on ill-conditioned problems.