Elimination is a sequence of choices

Gaussian elimination is the first algorithm in the subject and it is normally presented as a recipe: subtract multiples of the first row from the others to clear the first column, repeat on what remains, back-substitute. Followed exactly, it produces the answer, and there is nothing in the description that looks like a choice.

There is one at every step, and it is the only thing standing between the algorithm and failure.

What is actually being computed #

The recipe hides what it produces. Elimination does not merely solve a system; it computes a factorisation, and the factorisation is more useful than the solution.

Every step subtracts a multiple of the pivot row from the rows below it. Record those multiples and they assemble into a unit lower triangular matrix L. What is left at the end is upper triangular, call it U. Then A = LU, or with row interchanges, PA = LU where P is the permutation that describes which rows were swapped.

This matters for three reasons. It separates the O(n³) work from the O(n²) work, so a second right-hand side costs a hundredth of the first. It makes the determinant free — the product of the pivots, with a sign for the swaps. And it makes the algorithm checkable in a way a bare solution is not: ‖PA − LU‖ is a number, and it must be at the level of rounding, and if it is not then something is wrong with the factorisation rather than with the problem.

That last point is the site’s rule and the reason it appears on every elimination figure here. A solution can be wrong for reasons that have nothing to do with the algorithm. A factorisation cannot.

The choice #

At step k the algorithm needs a nonzero entry in position (k,k) to divide by. If the entry sitting there is zero, it must swap in a row that has something nonzero in that column, or it cannot proceed.

That much is forced. What is not forced — and is the whole content of the choice — is what to do when the entry is nonzero but small.

Partial pivoting is the answer everyone uses: scan the column from row k downwards, find the entry of largest absolute value, and swap that row into position. It costs one pass over a column per step, which is nothing next to the elimination itself.

What it buys is a bound. The multiplier used to eliminate row i is U[i][k] divided by the pivot, and since the pivot was chosen to be the largest of them, every multiplier is at most 1 in absolute value. That single fact is the entire stability argument for Gaussian elimination, and it is worth seeing on the figure rather than taking on trust — the badge reports the largest multiplier, and it is required to be ≤ 1 on every build.

Five matrices, and five different sequences of choices #

The figure takes the matrix as its argument, so the sequence of choices can be watched rather than described. Five of them, at three sizes.

The distinction between those two figures is the one the section above draws in words. In the first, pivoting is optional and buys a smaller multiplier; in the second it is the difference between an algorithm and a division by zero. Every matrix in this collection sits somewhere between them, and nothing about the matrix announces which.

Three rows moved out of four, and every multiplier at or below 0.5, is the comfortable case: partial pivoting had a real choice at every step and took it. The next matrix is the one where the arithmetic finally leaves a mark.

Across the five the rows moved are 0, 2, 1, 3 and 3 and the largest multiplier is 0.267, 0.5, 10⁻⁸, 0.5 and 0.875 — every one of them at most one, which is the theorem, and none of them close to one except the largest problem. The residual is exactly zero on four of the five and 6.8·10⁻¹⁷ on the fifth: at these sizes an elimination in double precision reproduces its matrix bit for bit most of the time, and the one case where it does not is the one with the most arithmetic in it.

Without the swap, a multiplier can be anything at all. On the two-by-two in the swap that is not optional it is 10¹⁷, and the subsequent subtraction annihilates the useful part of the row it touches.

Why bounded multipliers are enough — almost #

The reasoning is short. Each entry of U is the original entry minus a sum of products of multipliers with earlier entries. If every multiplier is at most 1, no single step can more than double the largest entry present, so after n−1 steps the largest entry of U is at most 2ⁿ⁻¹ times the largest entry of A.

That ratio — max|u| over max|a| — is the growth factor, and the error analysis of Gaussian elimination says the backward error is bounded by roughly n² times the growth factor times the unit roundoff. So bounded multipliers give a bounded growth factor, which gives a bounded backward error, which is backward stability.

The “almost” is that 2ⁿ⁻¹ is a very large bound. At n = 40 it is 5.5·10¹¹, which would swamp double precision entirely. In practice the growth factor of a random matrix at n = 40 is about three. The bound that is never attained is about that gap, about the one matrix that does attain it, and about why an algorithm with a useless worst-case bound is nonetheless the one everybody uses.

What the picture shows that a description does not #

Three things, and they are the reason the figure is a matrix rather than a flowchart.

The zeros are made, not found. In the “after step 1” panel the first column below the pivot is grey, and it is grey because three subtractions put zeros there. The algorithm does not discover structure in the matrix; it destroys structure and records what it destroyed.

The remaining block is a smaller instance of the same problem. After step k the untouched region is an (n−k)×(n−k) matrix, and the next step is elimination applied to it. That recursive structure is what makes the operation count n³/3 and what makes blocked implementations possible.

The entries change size. This is the part that only a picture with numbers on it conveys. The entries of U are not the entries of A rearranged; they are new numbers, and how large they get is the quantity the stability of the whole thing depends on.

The row order, recorded #

The permutation is not a detail. The factorisation is of PA, not of A, and forgetting the P produces a solver that returns confident nonsense.

Concretely: the figure’s matrix ends with its rows in the order 4, 1, 3, 2. That order is the record of which row had the largest entry at each step, and it has to be applied to the right-hand side before the forward substitution. In an implementation the permutation is usually stored as a vector of indices rather than a matrix, because a matrix multiplication to permute n numbers is an extravagance.

The check that catches a permutation bug is the site’s rule again: ‖PA − LU‖. Get the permutation wrong and that residual is of order 1 rather than of order 10⁻¹⁶ — not slightly wrong, completely wrong — because L and U are then factors of a different matrix. It is the single most effective one-line test for a hand-written LU, and it is why the badge on these figures carries it rather than carrying the residual of the solve.

Two routes to the determinant #

The other quantity elimination produces for free is the determinant, and it is a good example of the site’s habit of computing everything twice.

Route one: the product of the diagonal entries of U, times −1 for each row interchange. This costs nothing beyond the factorisation.

Route two: cofactor expansion, which is the definition. For a four-by-four it is 24 terms of four factors each, it is exact on an integer matrix, and it shares no arithmetic whatever with the elimination.

They must agree, and the assertion is run on every build. What it catches is a sign error in the swap counting, which is otherwise invisible: a determinant of the wrong sign looks entirely plausible, and every other property of the factorisation is unaffected by it. Two routes to a number is the thread, and this is the cheapest instance of it on the site.

What elimination costs #

The operation count is n³/3 multiply-and-subtract pairs for the factorisation, and n² for each solve that follows. Three practical consequences follow, and they shape how the algorithm is used.

Reuse the factorisation. Solving with a second right-hand side costs 1% of the first at n = 100 and 0.01% at n = 10,000. Code that calls a black-box solve(A, b) in a loop over right-hand sides is doing a hundred times more work than it needs to.

Never form the inverse. Computing A⁻¹ costs three times as much as factorising, and then every solve is a matrix–vector product — the same n² as a triangular solve pair, but with worse accuracy, because the explicit inverse has been through more arithmetic. x = inv(A) * b is slower and less accurate than x = A \ b, and it is written constantly.

The n³ is why structure matters. A banded matrix, a sparse one, a symmetric positive definite one — each admits a cheaper factorisation, and for large problems that is the difference between possible and not. This site’s foundation is dense; sparsity and fill-in are in the expansion, and they are where the elimination order becomes a combinatorial question rather than a numerical one.

Why L is unit lower triangular #

A detail that looks like a convention and is a consequence.

The diagonal of L is all ones, always, and it is not stored. That is not a normalisation chosen for tidiness — it falls out of what L records. Row i of L holds the multipliers used to eliminate row i against each earlier pivot row, and row i’s relationship to itself is the identity, so the diagonal entry is 1 by construction.

The practical effect is that L and U together fit in exactly the storage the original matrix occupied. U goes in the upper triangle including the diagonal, L’s strictly lower part goes in the lower triangle, and the ones are implied. Every serious implementation factorises in place, and the permutation vector is the only extra storage required.

That in-place structure is also why the picture in the filmstrip is honest about what is happening: the greyed entries below the diagonal are not gone, they are where the multipliers now live. A figure showing them as zeros would be showing the mathematics and hiding the computation.

What symmetry buys #

For a symmetric positive definite matrix, the whole apparatus simplifies and one worry disappears entirely.

The factorisation becomes A = LLᵀ — Cholesky — with half the arithmetic, half the storage, and no pivoting at all. That is not a heuristic. Positive definiteness guarantees that the diagonal entry at each step is the largest in its column, so the pivot search would find it anyway, and the growth factor is bounded by 1.

The consequence for the bound that is never attained is worth noting: for this whole class of matrices there is no gap between the worst case and the behaviour, because the worst case is 1. Cholesky is the one dense factorisation in the subject whose stability requires no empirical argument whatsoever.

It also gives a test. Cholesky succeeds if and only if the matrix is positive definite — a negative number under the square root is a definitive answer, not a numerical accident. Attempting the factorisation is the standard way to check definiteness, and it is cheaper than computing eigenvalues.

The order the entries are touched in #

One more thing the filmstrip shows that the recipe does not, and it is the reason the algorithm as written above is not the algorithm in a library.

The version drawn here touches the matrix a row at a time, which is how it is derived and how it is taught. A production implementation reorganises the same arithmetic into blocks — factor a panel of columns, then update the remaining submatrix with a single large matrix multiplication — because the arithmetic is identical and the memory traffic is not.

The natural next sentence is that the rounding is therefore not identical — reordering the operations changes which sums are accumulated in which order, and the order they are added in is exactly about what that does — and that blocked LU is typically slightly more accurate, because the inner updates are matrix multiplications performed with a longer accumulation.

That is worth measuring rather than repeating, and it does not come out that way. Right-looking LU with partial pivoting at three block sizes, 64×64, backward error ‖PA − LU‖/‖A‖:

The first three columns are not merely equal — they are bit-identical. Every entry of U and every pivot, at every block size and every condition number, ten comparisons out of ten. Blocking on its own changes no rounding whatsoever.

The reason is that a longer accumulation buys nothing when the accumulator is the same width as the storage. Keeping a running sum in a register across eight terms and storing it back eight times are the same sequence of roundings, because storing a double is exact. What a genuinely wider accumulator buys is the fourth column: compensated summation in the trailing update improves the backward error by 4% to 11% — real, and about a tenth of a digit.

The fifth column is there because the fused multiply–add is the mechanism usually offered next, and it removes one rounding per product. It does not help: slightly worse at four of the five rows. The roundings that dominate LU’s backward error live in the sums, not in the products, which is the same lesson the order they are added in draws about a dot product.

So the conclusion below stands and the reason for it changes. A bit-for-bit comparison between two LU implementations is not a meaningful test — but not because blocking changes the answer, since it demonstrably does not. It is because implementations differ in accumulator width and summation order, and two that agree on those agree exactly, as these do. That is why the checks on this site compare residuals and invariants rather than entries.

assertBlockingChangesNothingAndWidthChangesALittle requires the bit-identity at every block size and requires the compensated version to be better by a few per cent rather than by a digit.

What is asserted here #

Every claim in this essay is checked when the figure is generated, which means on every build.

The filmstrip’s last frame must be exactly the U that the factorisation routine returns, entry by entry, to 10⁻¹². That sounds pedantic and it is the assertion that catches a figure showing one computation while its caption describes another — the elimination in the generator is written out separately, so that it can record the intermediate states, and nothing but this check keeps the two implementations in agreement.

‖PA − LU‖/‖A‖ must be below 10⁻¹³. The largest multiplier must be at most 1. The determinant from the pivots must equal the determinant from cofactors to 10⁻¹².

And the growth measurement is fed something it must not miss: Wilkinson’s matrix, whose growth factor is exactly 2ⁿ⁻¹. A measurement routine that reported small growth for everything would make the bound that is never attained a meaningless essay, so it is shown the one matrix that breaks it and required to notice.

Where elimination sits in the collection #

Elimination is the first algorithm in the subject and the last one to be fully understood, and the three essays of this field take it in that order.

This one is about what the algorithm computes and what the choice at each step is. The swap that is not optional runs the version that makes the choice badly, and measures a backward error of 0.25 on a matrix with κ ≈ 2.6 — a well conditioned problem answered wrongly, which is the one situation where the algorithm is unambiguously to blame. The bound that is never attained is about the quantity the whole stability argument rests on, and about the fact that its bound is eleven orders of magnitude away from its behaviour.

Read together, they are a case study in the distinction the exact answer to a nearby problem sets up. The same algorithm, on the same kind of matrix, is backward stable or not depending on one comparison per step; and when it is stable, everything that remains wrong with the answer belongs to the problem.

The one thing to take away #

Anyone writing an LU factorisation by hand should, print ‖PA − LU‖/‖A‖ once and look at it.

It costs two matrix multiplications on a matrix already in memory, it takes five minutes, and it distinguishes a working factorisation from every common way of getting one wrong: a permutation applied in the wrong direction, an off-by-one in the multiplier loop, a forgotten sign, an L whose diagonal was stored when it should have been implied. All of those produce a factorisation that looks entirely normal and reconstructs a different matrix.

That is the reason the rule this site runs on is about residuals rather than about answers. An answer can be wrong for reasons no code change would fix. A factorisation that does not reconstruct its matrix is a bug, every time, and it is one number away from being visible.

A rule with no choices in it #

Every step of an elimination compares numbers and picks one, and that is what keeps its intermediates bounded. A closed-form rule has no such step — its products are named in advance — and there is nowhere in it to be careful.