What a slow iteration was secretly building
In the last guide you watched the Jacobi and Gauss-Seidel iterations creep toward the answer, converging only when the spectral radius of their iteration matrix sits below 1 — and often creeping painfully slowly. They share one humble habit worth staring at: from the current guess they form the residual vector r = b - A x, and the only expensive thing they ever do with the matrix is multiply by it. They never factor A, never even look inside it. That single multiply, A times a vector, is the whole budget.
Now ask a sly question. Start from x_0 = 0, so the first residual is just r_0 = b. After one matrix multiply you can also see A b. After two, A^2 b. After k multiplies, the only vectors you could possibly have formed by combining what you have seen are weighted sums of b, A b, A^2 b, up to A^{k-1} b. That growing span has a name: the Krylov subspace of dimension k, written K_k. Every iteration that only multiplies by A is, whether it knows it or not, confined to this staircase of subspaces — each rung one dimension taller than the last.
Turning A x = b into rolling downhill
Conjugate gradient works on a special but enormously common case: A is symmetric positive definite (SPD) — symmetric, and x^T A x > 0 for every nonzero x. These are exactly the matrices you get from springs, from least-squares normal equations, from discretized Laplacians in heat and electrostatics. For an SPD matrix there is a magic trick: solving A x = b is identical to finding the bottom of a bowl. Define the energy phi(x) = (1/2) x^T A x - b^T x. Because A is SPD this surface is a perfect upward bowl with a single lowest point, and at that point its gradient is exactly A x - b. The gradient vanishes — the bowl bottoms out — precisely when A x = b.
So linear solving becomes minimization, and the residual r = b - A x is literally the downhill direction (the negative gradient). The most obvious idea, steepest descent, is to repeatedly step along r. It works, but for a long thin bowl it zig-zags maddeningly: each step undoes part of the progress of the last, because the new downhill direction is not aware of where the old steps already pointed. The number of zig-zags scales with the condition number of A — a stretched, ill-conditioned bowl is a long narrow valley, and naive descent bounces between its walls forever.
The conjugate trick: search directions that never interfere
The fix that gives conjugate gradient its name is to choose search directions that, in the geometry warped by A, are mutually perpendicular — called A-conjugate: two directions p and q are A-conjugate when p^T A q = 0. The payoff is gorgeous. If you minimize exactly along a set of A-conjugate directions, optimizing along a new one never spoils the optimality you already won along the old ones. No undoing, no zig-zag. Each step is a clean, permanent piece of progress, and in exact arithmetic you reach the exact answer of an n-by-n system in at most n such steps.
Here is the part that feels like a conjuring trick. You might fear that building a whole set of A-conjugate directions needs you to remember and re-orthogonalize against all previous ones — expensive and memory-hungry. For SPD matrices it does not. A short three-term recurrence does the whole job: the next direction is the new residual plus a single scalar times the previous direction. CG stores just a handful of vectors and does one matrix-vector product per iteration, yet each new direction comes out automatically A-conjugate to every earlier one. This near-miraculous economy is exactly what the symmetry of A buys you; lose symmetry and it vanishes, which is why the next guide needs the heavier GMRES.
Conjugate Gradient for SPD system A x = b (start x = 0)
r = b # residual = negative gradient of the bowl
p = r # first search direction
repeat:
Ap = A * p # the ONE matrix-vector product per step
alpha = (r . r) / (p . Ap) # exact 1-D minimum along p
x = x + alpha * p # take the step
r_new = r - alpha * Ap # update residual (no extra A multiply)
if ||r_new|| <= tol * ||b||: break # relative stopping test
beta = (r_new . r_new) / (r . r)
p = r_new + beta * p # new direction, auto A-conjugate to all past
r = r_newHow fast, and when to stop
The 'at most n steps' guarantee is theoretically lovely but practically beside the point — for a system with millions of unknowns you never want n steps. The real reason CG is treasured is that it usually reaches a good enough answer in far fewer. The honest convergence estimate says the error shrinks roughly like ((sqrt(kappa) - 1)/(sqrt(kappa) + 1))^k after k steps, where kappa is the condition number of A. Compare that to the condition-number dependence of steepest descent: CG depends on sqrt(kappa), not kappa, which can be the difference between thousands of iterations and dozens.
That formula carries a blunt warning. Convergence speed is governed by kappa, so an ill-conditioned A still makes CG crawl. This is the entire reason the fifth guide exists: a preconditioner M approximates A but is cheap to solve, and applying CG to the better-conditioned transformed system collapses the iteration count. Plain CG on a nasty Laplacian might need thousands of steps; with a good preconditioner, dozens. Never reach for CG on a hard problem without planning to precondition it.
Because you stop early, you need an honest stopping criterion, and CG hands you one for free: it updates the residual every step, so you can watch ||r|| fall and halt when ||b - A x|| <= tol times ||b|| — a relative test, since absolute size is meaningless without scale. Two honest cautions, though. First, a small residual does not guarantee a small error in x: the two are linked by exactly the condition number again, so on an ill-conditioned system you can have a tiny residual while x is still noticeably off. Second, in real floating-point arithmetic the beautiful A-conjugacy slowly erodes through round-off, so the textbook 'exact in n steps' never quite holds; CG is used as an honest iterative method, run until the residual is small, not as a finite direct one.
Why this matters and where it leads
Step back and weigh the bargain. A direct Cholesky factorization of a dense n-by-n SPD matrix costs O(n^3) work and O(n^2) storage, and on a sparse matrix it can suffer fill-in that bloats both. CG costs one matrix-vector product per iteration and stores only a few vectors. If A has just a handful of nonzeros per row — the rule for a discretized PDE — that product is O(n), and a few hundred iterations beats O(n^3) by an astronomical margin once n runs into the millions. This is precisely the 'why iterate instead of factor' promise from the first guide of this rung, now made sharp.
There is a deeper reward hiding in 'A is touched only through A times a vector'. You never need A as an explicit array of numbers at all — only a routine that, given a vector v, returns A v. This is the matrix-free idea, and it is liberating: the action of a giant discretized operator can be computed on the fly from a stencil, with the matrix never assembled or stored. Krylov methods are the natural partners of matrix-free computing, which is how solvers reach problem sizes that could never fit in memory as an explicit matrix.