The wall that elimination hits
In the earlier rung you learned the gold standard for solving A x = b: Gaussian elimination, which factors A = L U once and then solves cheaply by forward and back substitution. It is exact (up to rounding), it is robust with pivoting, and for one dense system of size n it costs about (2/3) n^3 floating-point operations. For n = 1000 that is roughly a billion flops — a fraction of a second. So far elimination looks unbeatable, and for small or moderate dense systems it genuinely is. This rung is not about replacing it; it is about the place where it quietly stops fitting.
Now picture where huge systems actually come from. Discretize the heat equation or a stress field on a 3D grid of 100 x 100 x 100 points and you get one unknown per grid point — a system with n = one million. The (2/3) n^3 cost of dense elimination is now about 7 x 10^17 flops, which even at a teraflop per second is roughly a week of computing for a single solve. Worse, storing the dense matrix A would need n^2 = 10^12 numbers, eight terabytes in double precision. Neither the time nor the memory is remotely affordable. The dense recipe did not get slightly slow; it fell off a cliff, because n^3 and n^2 grow brutally.
Sparsity: the structure elimination keeps destroying
There is a saving grace those grid matrices possess: they are sparse. A grid point only couples to its handful of immediate neighbours, so each row of A has maybe 5 or 7 nonzeros, not a million. The matrix that needed eight terabytes when stored densely needs only about 7 million numbers when you store just the nonzeros — a few dozen megabytes, easy. Sparsity is the gift that makes million-unknown problems imaginable at all. The catch is that you have to handle the matrix in a way that respects that gift, and naive elimination does not.
Here is the painful part. When elimination eliminates a variable, it adds multiples of one row to others, and positions that were zero can become nonzero. Those new nonzeros are called fill-in. A sparse matrix that started with 7 nonzeros per row can fill in until its L and U factors are dense, blowing the memory budget right back up. Clever reordering of the unknowns can hold fill-in down — this is the art behind a sparse direct solver, and for 2D problems it works beautifully. But for 3D grids the fill-in is unavoidable in the worst case, and the factors grow far past linear, dragging memory back toward the cliff you were trying to escape.
So the real tension is sharp. Direct elimination is exact and reliable, but it fights the sparsity — it keeps trying to fill in the zeros that made the problem tractable. We want a method that touches A only in ways that preserve sparsity. And there is exactly one operation on a sparse matrix that costs almost nothing and never creates fill-in: multiplying the matrix by a vector. That single observation is the doorway out.
The iterative idea: guess, measure, improve
An iterative method gives up on reaching the exact x in a fixed number of steps. Instead it starts from a guess x_0 — often just zeros — and produces a sequence x_1, x_2, x_3, ... that, if all goes well, marches toward the true solution. You stop when the current guess is good enough. This is the same spirit you already met for root-finding, where fixed-point iteration turned x = g(x) into x_{n+1} = g(x_n) and let a contraction pull you to the answer. Solving A x = b iteratively is that idea lifted from a single number to a whole vector.
How does the method know whether a guess is good without knowing the true x? It measures the residual r = b - A x_current. If x were exact, A x would equal b and the residual would be zero; the size of r tells you how far A x lands from b. Critically, computing the residual needs only one sparse matrix-vector product A x plus a subtraction — the cheap, sparsity-preserving operation we identified. Every iterative method in this rung is, at heart, a clever recipe for turning the residual into a correction that nudges x closer.
the skeleton every iterative solver shares:
x = x0 # initial guess (often all zeros)
repeat:
r = b - A*x # residual: one sparse mat-vec, no fill-in
if ||r|| / ||b|| < tol: # close enough? stop.
break
x = x + (correction from r) # the method's secret sauce lives here
cost per step is dominated by the single product A*x ,
which for a sparse A with ~7 nonzeros/row is O(n), not O(n^2).What you trade, and what you gain
The trade is real, so be honest about both sides. You give up exactness: an iterative solver returns an approximation, and you must decide how accurate is accurate enough by choosing a tolerance on the relative residual ||b - A x|| / ||b||. (We will treat residuals and stopping criteria carefully in a later guide, because a small residual is not automatically a small error — that gap is governed by the condition number, exactly as in the conditioning rung.) You also give up the lovely property of factor once, solve many: an LU factorization can solve A x = b for many right-hand sides b almost for free, whereas a basic iteration must restart its march for each new b.
What you gain is everything that lets the million-unknown problem fit. Memory: an iterative solver never builds L and U, so it needs nothing beyond A's nonzeros and a few work vectors of length n — sometimes it does not even store A explicitly, just a routine that returns A x, which is the idea behind a matrix-free method. Time: each step costs O(n) for a sparse A, and if the method converges in a number of steps that grows slowly with n, the whole solve can run in nearly O(n) — the near-linear growth rate we were chasing. That is the bargain: trade exactness and a fixed step count for memory you can actually afford and a cost that scales.
A first taste: the tiny stationary recipe, and where it goes
Let us make the abstraction touchable with the very first family of iterative methods, the stationary methods. They come from a trick called matrix splitting: chop A into an easy-to-invert part M and a remainder, A = M - (M - A). Then A x = b becomes M x = (M - A) x + b, which suggests the iteration x_{n+1} = M^{-1} ((M - A) x_n + b). You solve a cheap system with M each step instead of the hard one with A. The next guide builds this into the Jacobi and Gauss-Seidel methods, but you can already see the shape: pick an easy M, sweep, repeat.
But this is where honesty matters most, because stationary iterations do not always converge. Whether x_n marches toward the true x or wanders off forever is decided by the spectral radius of the iteration matrix M^{-1}(M - A) — the size of its largest eigenvalue in magnitude. If that number is below 1, the error shrinks by that factor every step and you win; if it is at or above 1, you lose, no matter how long you wait. So convergence is a genuine condition to check, not a guarantee — the direct analogue of the contraction condition you needed for fixed-point root-finding. And even when these simple methods converge, they can be slow, which is exactly what motivates the cleverer methods to come.
That road map is worth holding in your head as you climb the rest of this rung. The stationary methods are the warm-up. From them grows a deeper idea — building up the solution inside the Krylov subspace spanned by b, A b, A^2 b, ..., which leads to the conjugate gradient method for symmetric positive-definite systems and GMRES for general ones. And the practical key that turns 'converges, eventually' into 'converges fast' is the preconditioner: a cheap approximate inverse that reshapes the problem so the iteration races. Multigrid, glimpsed near the end, can even reach that near-O(n) ideal on the right problems. Every one of these still rests on the same humble foundation you just met — the sparse product A x, repeated.