The school trick, grown up
You already know how to solve a small system A x = b by hand: pick an equation, use it to wipe out one variable from all the others, repeat until only one variable is left, then back your way out. That schoolbook ritual has a grown-up name, Gaussian elimination, and it is the single most-run algorithm in scientific computing. The whole point of this rung is to see it not as a one-off procedure for one b, but as a machine you can crank cleanly, count exactly, and reuse — which is where the LU factorization comes in.
Let us crank it on a tiny 3-by-3 system so the moving parts are visible. Picture the augmented array — the matrix A with the right-hand side b tacked on as a final column. The algorithm sweeps left to right, one pivot column at a time. In each column it picks the diagonal entry as the pivot, then for every row below it subtracts a multiple of the pivot row chosen exactly so the entry under the pivot becomes zero. That multiplier — the factor you scaled the pivot row by — is the one number you must not throw away. Most people, doing this by hand, compute it and forget it. Remembering it is the whole secret of LU.
[ 2 1 1 | 5 ] [ 2 1 1 | 5 ]
[ 4 3 3 | 13] --> [ 0 1 1 | 3 ] row2 -= 2*row1
[ 6 7 9 | 22] [ 0 4 6 | 7 ] row3 -= 3*row1
[ 2 1 1 | 5 ]
--> [ 0 1 1 | 3 ]
[ 0 0 2 | -5 ] row3 -= 4*row2
multipliers used: 2, 3, 4 (keep these!)Two triangles fall out for free
Look at what elimination left behind in the example. The transformed matrix is upper triangular — all zeros below the diagonal — call it U. And the multipliers we were told to keep (2, 3, 4) slot exactly into the below-diagonal positions of a lower triangular matrix L with ones on its diagonal. The astonishing fact is that these two pieces multiply back to the original: A = L U. Elimination did not just solve a system; it silently factored A into a lower-triangular times an upper-triangular matrix. That product of two triangular matrices is the LU factorization.
Why is a product of two triangles such a prize? Because triangular systems are trivially easy to solve. With L lower triangular, L y = b can be read off from the top down — the first equation has one unknown, the second has two but you already know the first, and so on; that is forward substitution. With U upper triangular, U x = y peels apart from the bottom up by back substitution. So once you have A = L U, solving A x = b splits into two cheap triangular sweeps. The next guide in this rung dwells on exactly how these substitutions work and why they cost so little.
- Start from A x = b and substitute A = L U, so the system reads (L U) x = b. Regroup the parentheses as L (U x) = b.
- Name the inner part y = U x. The outer system is now L y = b, a lower-triangular system — solve it top down by forward substitution.
- With y in hand, solve U x = y, an upper-triangular system — solve it bottom up by back substitution. The x you get satisfies the original A x = b.
Factor once, solve many
Now the payoff that makes LU more than a tidy restatement of Gaussian elimination. The expensive part — the elimination sweep that produces L and U — touches A alone and never looks at b. So if you must solve the SAME matrix A against many different right-hand sides b_1, b_2, b_3, ... you pay the elimination cost only once, then handle each new b with a pair of dirt-cheap triangular solves. This is the factor-once-solve-many idea, and it is everywhere: a structural engineer pushes a hundred different load cases through one stiffness matrix; an animation pins the same physics and re-solves every frame.
The numbers say why this matters. The elimination that builds L and U for a dense n-by-n matrix costs about (2/3) n^3 operations — an O(n^3) job, the dominant expense. Each triangular solve costs only about n^2 operations, an O(n^2) job. For n = 1000 that is a billion-ish operations to factor versus a million to solve: the second right-hand side is roughly a thousand times cheaper than the first. The flop count is not pedantry; it is the difference between an overnight job and an instant one, and it is why nobody re-runs elimination per b.
Where naive elimination breaks
Everything above assumed the pivot — the diagonal entry we divide by — was never zero. But the very first matrix you try might have a zero sitting on the diagonal. Then the multiplier is a division by zero and the algorithm dies on the spot. Swapping that row with a lower one that has a nonzero entry there fixes it instantly, and that row swap is the seed of pivoting. So pivoting is not an optional refinement bolted on later; sometimes plain elimination simply cannot proceed without it.
The subtler trouble is a pivot that is not zero but merely tiny. Dividing by a small number makes the multipliers huge, and huge multipliers magnify the round-off error already lurking in floating-point arithmetic — remember that every operation rounds, and that 0.1 has no exact binary form, so the machine is never working with your exact numbers. A tiny pivot can amplify those rounding errors until the computed answer is nonsense, even though the underlying problem was perfectly well-behaved. That gap between a sound problem and a method that wrecks it is the difference between a problem's condition number and an algorithm's stability.
The cure is partial pivoting: before eliminating in a column, scan it for the entry with the largest magnitude and swap that row up to be the pivot. This keeps every multiplier at most 1 in size, which keeps round-off under control. With partial pivoting in place, Gaussian elimination becomes a backward-stable algorithm in practice — it returns the exact answer to a problem within round-off of yours, which is the most an honest method can promise. The next guide, why we pivot, unpacks the stability story; for now just hold the rule: real solvers always pivot, and what they actually compute is P A = L U, where P records the row swaps.
Structure, and what comes next
That (2/3) n^3 price tag assumes a dense matrix, every entry possibly nonzero. But the matrices that come out of real models are often sparse — mostly zeros, with a few nonzeros in a predictable pattern. A heat equation on a grid, for instance, gives a banded matrix whose nonzeros hug the diagonal. A sparse direct solver never stores or touches the zeros, so an LU factorization that would cost O(n^3) dense can drop to nearly O(n) for a narrow band. The one-dimensional special case is so clean it has its own name, the Thomas algorithm, an O(n) sweep for tridiagonal systems.
There is an honest catch, though. Elimination can turn a zero into a nonzero — a position that was empty in A becomes occupied in L or U. That is fill-in, and on a badly ordered sparse matrix it can flood the factors with so many new nonzeros that the sparsity advantage evaporates. The remedy is to reorder the rows and columns first, choosing an elimination order that creates as little fill as possible. This fill-in and reordering tradeoff is the heart of sparse direct solving, and the final guide of this rung returns to cost and sparsity in earnest.