The swap that is not optional

Here is a two-by-two system. The matrix has ε in its top-left corner and ones everywhere else; the right-hand side is (1, 2). Its exact solution, for any small ε, is x₁ = 1/(1−ε) and x₂ = (1−2ε)/(1−ε), which for ε = 10⁻¹⁷ is (1.000000, 1.000000) to every digit anyone would print.

Solve it by Gaussian elimination without a row interchange. The answer that comes back has 1 in its first component and 0 in its second.

Nothing failed. The pivot ε is nonzero, so there was no division by zero. No exception was raised, no condition was flagged, and the answer is a perfectly ordinary pair of numbers. It is just not the answer.

What went wrong, precisely #

The multiplier is 1/ε = 10¹⁷. Eliminating the second row subtracts 10¹⁷ times the first row from it, so the (2,2) entry becomes 1 − 10¹⁷ and the right-hand side becomes 2 − 10¹⁷.

Both of those are computed correctly. But 1 − 10¹⁷ is, in double precision, exactly −10¹⁷: the 1 is below half the gap between representable numbers at 10¹⁷, which is 16. So the information carried by that 1 — which is the entire content of the second equation — is gone. Likewise the 2.

Back-substitution then computes x₂ = (−10¹⁷)/(−10¹⁷) = 1, and x₁ = (1 − x₂)/ε = 0/ε = 0. The algorithm has solved, exactly, a system in which the two ones were absent.

Two things about this are worth separating. The subtraction 1 − 10¹⁷ is correctly rounded and introduces the smallest error it is permitted to. What destroyed the answer is that the operands were sixteen orders of magnitude apart, and they were sixteen orders apart because a multiplier of 10¹⁷ was used. The whole purpose of the row interchange is to make that impossible.

The growth factor is exactly 1/ε − 1 #

The unpivoted run does not merely produce a large entry — it produces a specific one, and knowing the value rather than a bound makes the claim checkable at every position of the slider.

max|a| is 1 and max|u| is |1 − 1/ε|, so the growth factor is 1/ε − 1 exactly. At ε = 10⁻³ it is 999; at ε = 10⁻¹⁷ it is 10¹⁷. That is asserted as an equality on every drag frame, not as an inequality, which means a change to the code that altered the arithmetic would be caught at every ε rather than at the one that happened to be tested.

The pivoted run’s growth factor is 1. It cannot be otherwise: with the swap, both multipliers are ε, and nothing grows.

The failure is graded, which is worse than a cliff #

Drag the slider on the figure and the most useful observation is what does not happen. There is no threshold. At ε = 10⁻² the unpivoted answer is fine. At 10⁻⁶ it has lost a few digits. At 10⁻¹² it has lost most of them. At 10⁻¹⁷ it has lost all of them. The transition is smooth and unmarked.

This is why unpivoted elimination survives in codebases. It works on the test cases, it works on the first few real problems, and then one day the data has a small leading entry and the answer is wrong. Nothing in the output distinguishes that day from the others.

A cliff would be better. A method that failed loudly at a known threshold could be checked against it. What this does instead is degrade at a rate that depends on data that was never inspected, which is the failure mode that costs the most to find.

Backward error says which, if it is computed #

The one thing that does distinguish the good run from the bad one is available and cheap.

The pivoted solve has a backward error of 0. Exactly zero — the computed answer is the exact solution of the system as given. The unpivoted solve has a backward error of 0.25.

That number is worth staring at. It says the answer returned is the exact solution of a system 25% away from the one that was asked about. Not 25% in some obscure component: a quarter of the size of the matrix, in the norm. No amount of ill-conditioning excuses that, because backward error is purely the algorithm’s contribution — the exact answer to a nearby problem is the essay about why that separation is available at all.

And here is the diagnostic value. The matrix in question has κ ≈ 2.6, which is about as well conditioned as a matrix can be. A well-conditioned problem with a huge backward error is unambiguous: the algorithm is at fault, and no amount of precision or reformulation is the answer. Swapping two rows is.

Why “just check for a zero pivot” is not enough #

The naive defence is to test whether the pivot is zero and swap only then. It fixes the case where the algorithm cannot proceed and leaves the case where it proceeds badly, which is the one that loses answers.

There is no threshold that repairs it either. Any rule of the form “swap if the pivot is smaller than δ” has to pick δ, and the right value depends on the sizes of the other entries in the column — which is exactly the comparison partial pivoting performs. Choosing the largest available entry is not a heuristic; it is the only scale-invariant rule, and it is the reason partial pivoting is not merely a strategy but the strategy.

Two stronger strategies exist. Complete pivoting searches the whole remaining submatrix rather than one column, which gives a much better growth bound — polynomial rather than exponential — at a cost of O(n³) comparisons, which is the same order as the arithmetic. Nobody uses it, because the growth bound partial pivoting fails to achieve is one it does not need in practice. Rook pivoting sits between them. The fact that the weakest of the three is universal is a small lesson in what bounds are worth.

The refusal #

An essay claiming that an algorithm fails owes a demonstration that the check could tell. Three assertions run on every build, and each is capable of firing.

The pivoted solve must be stable. Backward error below 10⁻¹⁵. If a change made pivoting stop working, this fires.

The unpivoted solve must be worse. Its forward error must exceed the pivoted one. If a change made the unpivoted version accidentally work, this fires — and that matters, because an essay whose central claim has quietly become false is worse than no essay.

The growth must be exactly 1/ε − 1. An equality, checked to nine digits, at all seventeen positions of the slider.

There is a fourth in lib/matrix.js, which is the version of this system with a hard zero in the corner: with pivoting, the answer is exact; without, the factorisation breaks down and the routine must report it rather than returning infinities. A solver that returns NaN is a solver whose output a figure will draw as an empty picture, which is the failure this site has been caught by before. Assertions that reject is the thread.

The scaling that defeats the rule #

Partial pivoting chooses the largest entry in the column, and there is a case where that is the wrong entry — worth knowing because it is the one exception, and because the fix is not more pivoting.

Multiply the first row of the two-by-two by 10¹⁸. The system is unchanged mathematically: the same equation, scaled. But now the entry in the corner is 10¹⁸ε rather than ε, and if ε is 10⁻¹⁷ that is 10, which is larger than the 1 below it. Partial pivoting looks at the column, sees 10 against 1, and declines to swap — reproducing exactly the failure it exists to prevent.

The problem is that “largest entry in the column” is only a sensible criterion when the rows are comparably scaled. Scaled partial pivoting fixes it by comparing each candidate against the largest entry in its own row rather than in absolute terms, and it is what a careful implementation does when the rows may have wildly different magnitudes.

The deeper point is that row scaling is a property of how the problem was written down rather than of the problem, and a pivoting rule that depends on it is depending on an arbitrary choice. Equilibrating the matrix first — scaling rows and columns so the entries are comparable — makes the question go away and usually reduces the condition number as well, which is the recommendation in the condition number is an amplifier.

Two more pivoting strategies, and why nobody uses them #

Complete pivoting searches the entire remaining submatrix at each step and swaps both a row and a column, bringing the largest available entry to the pivot position. Its growth factor is bounded by roughly n^(¼ log n) — polynomial rather than exponential, and a genuine improvement on partial pivoting’s guarantee.

It costs O(n³) comparisons across the factorisation, which is the same order as the arithmetic, so it is roughly twice as slow. It also destroys any structure the matrix had — a banded matrix does not stay banded — and it requires tracking a column permutation as well as a row one.

Rook pivoting searches alternately along rows and columns until it finds an entry that is largest in both, which usually terminates quickly. Its growth bound is between the other two and its expected cost is close to partial pivoting’s.

Neither is used in general-purpose libraries, and the reason is the subject of the bound that is never attained: partial pivoting’s guarantee is useless and its behaviour is excellent, so paying for a better guarantee buys nothing observable. That is a decision the field made empirically and has never had to revisit.

What a pivot search costs #

Worth quantifying, because “pivoting is free” is asserted more often than it is measured.

The search at step k examines n − k entries, so the total across the factorisation is about n²/2 comparisons, against n³/3 multiply-and-subtract pairs for the arithmetic. At n = 100 that is 5,000 comparisons against 333,000 arithmetic operations: about 1.5%, and comparisons are cheaper than multiplications.

The row interchange itself is the more interesting cost. Swapping two rows of a large matrix moves n numbers and, on a machine with a cache, may move them a long way. In a blocked implementation the swaps are deferred and applied to a whole panel at once for exactly this reason, and on distributed memory the interchange can require communication, which is why alternative schemes such as tournament pivoting exist at scale.

For anything that fits in memory the answer is unambiguous: pivoting costs a couple of percent and buys the difference between an answer and a number.

What this generalises to #

The pattern here appears three more times on this site and it is worth naming.

An algorithm has two variants. They are algebraically identical — the same derivation produces both, and in exact arithmetic they return the same answer. One of them forms a quantity from operands of wildly different magnitude and the other does not. The second is stable and the first is not, the difference is invisible in the derivation, and it is visible in one measured number.

Two Gram–Schmidts is the same shape: one word changed in a line of pseudocode, and eight orders of magnitude in ‖QᵀQ − I‖. The road that squares the problem is the same shape again, with a whole method rather than a line.

In each case the lesson is not “be careful”. It is that the choice exists, that it has a right answer, and that the right answer costs nothing — one comparison per step here, one variable’s worth of bookkeeping there. The expensive part is knowing the choice was there at all, which is what running the failing version is for.

The same shape, three more times #

The pattern of this essay recurs, and recognising it is worth more than any of its instances.

An algorithm has two variants that a derivation cannot tell apart. One of them, at some point, forms a quantity from operands of wildly different magnitude. That variant is not backward stable, the other is, the cost of choosing correctly is nil, and the failure is silent.

Two Gram–Schmidts differs by which vector a projection coefficient is computed against — one word — and produces a Q that is not orthogonal on a matrix where the other produces one that is. The road that squares the problem differs by whether AᵀA is formed, and turns a condition number of 10⁶ into 10¹². Cancellation takes the answer is the same pattern in a single expression, where (1 − cos x)/x² and 2 sin²(x/2)/x² are the same function and differ by sixteen orders of magnitude.

In every case the derivation is silent, the failure is quiet, and the fix is free. That is why this site runs the algorithms rather than describing them: the difference is not visible in the algebra, and it is completely visible in one measured number.

What to do in practice #

Use a library. LAPACK’s dgesv pivots, has done for forty years, and is faster than anything hand-written.

Where a hand-written version is unavoidable — and there are reasons to, including the one this site is built on — then the checklist is short. Pivot, always, without a threshold. Record the permutation and apply it to the right-hand side. Print ‖PA − LU‖/‖A‖ once and confirm it is at rounding level. Track the largest multiplier and confirm it is at most one, which is a one-line check that catches a pivot search comparing the wrong quantity.

And when a solve surprises, compute the backward error before doing anything else. It costs a matrix–vector product and it says which half of the problem to look at — a small residual is not a small error is about what that number does and does not establish, which is the other half of the same skill.