The bound that is never attained

The backward stability of Gaussian elimination rests on one quantity. Not the condition number, which belongs to the problem; not the multipliers, which pivoting bounds by one; but the growth factor — the ratio of the largest entry appearing anywhere in U to the largest entry of A.

The error analysis says: the backward error is bounded by roughly n²·ρ·u, where ρ is the growth factor. Everything else in that expression is known. n is the size of the matrix and u is the unit roundoff, and neither is a surprise. ρ is the one that has to be argued about.

Partial pivoting bounds it by 2ⁿ⁻¹, because with multipliers at most one, no single step can more than double the largest entry present. At n = 40 that bound is 5.5·10¹¹, which multiplied by u gives 1.2·10⁻⁴ — a backward error of one part in ten thousand, from a method described as backward stable.

The measured growth factor of a random forty-by-forty matrix is about three.

The one matrix that attains it #

The bound is not loose because it is badly derived. It is attained, exactly, by a matrix anyone can write down from memory: ones on the diagonal, −1 everywhere below it, and a column of ones down the right-hand edge.

For that matrix, at every size tested from four to twenty, the growth factor is exactly 2ⁿ⁻¹ — not approximately, not asymptotically. The check on every build is an equality to twelve digits rather than an inequality, because the value is known in closed form.

The detail that makes it interesting: partial pivoting performs no row interchanges at all. Every candidate pivot is already the largest entry in its column, so the strategy has nothing to choose, does exactly what it is supposed to, and the entries double anyway. This is not a case of the algorithm being defeated by a pathological input it fails to notice. It is a case of the algorithm working perfectly and the bound being tight.

All the growth is in the final column. Every other entry of U is unchanged from A, which is checked separately — it is the sort of claim that is easy to make from a picture and easy to be wrong about.

Why nobody worries #

Two reasons, and they are different in kind.

The empirical one. Growth factors above ten are essentially never observed on matrices arising from applications, and matrices constructed to produce large growth have to be constructed. Trefethen and Schreiber measured the average growth for random matrices as roughly n^(2/3), which at n = 1000 is a hundred — large, and eight orders of magnitude below 2⁹⁹⁹. The measurement in the figure above agrees: median growth rises slowly with n, from 1.13 at n = 5 to 3.23 at n = 40.

The structural one. For several important classes the growth factor is bounded by a small constant, provably. A symmetric positive definite matrix needs no pivoting at all and has ρ ≤ 1. A diagonally dominant matrix likewise. A totally nonnegative matrix has ρ = 1. Many real problems fall into one of those classes, and for them the worst case is not merely unlikely, it is impossible.

Neither reason is a proof about the general case, and that is the honest position: Gaussian elimination with partial pivoting has a worst-case bound nobody relies on and an observed behaviour everybody relies on, and the gap between them has been an open question for seventy years.

What to do with a useless bound #

This is the part worth thinking about, because it recurs everywhere in computing and is usually handled badly.

The bound is not useless in the sense of being wrong. It is tight — Wilkinson’s matrix proves it — and it cannot be improved without additional hypotheses. What is useless is using it to decide anything, because it is nowhere near what happens.

The three available responses:

Measure and report the growth. It costs one comparison per step to track max|u| while factorising, and it converts an untestable worry into a number. LAPACK does not return it by default, which is a pity; when it is available it is the first thing to look at if a solve behaves oddly.

Use a strategy with a better bound where it is required. Complete pivoting has growth bounded by roughly n^(¼ log n), which is polynomial. It costs O(n³) comparisons — the same order as the arithmetic — and essentially nobody pays that for a guarantee they will never invoke.

Rely on the empirical behaviour, and say so. This is what everyone does, and the only thing wrong with it is doing it silently. “Gaussian elimination with partial pivoting is backward stable in practice” is a true and useful sentence; dropping the last three words makes it false.

The measurement, and what it required #

The figure is a measurement rather than a citation, which meant deciding several things that a citation would have hidden.

Thirty matrices per size, not one. The growth factor of a single random matrix is a random variable, and its spread is not negligible — at n = 20 the thirty values run from about 1.4 to 2.8. Plotting one would have produced a jagged line that a reader would reasonably interpret as structure. The figure draws the median solid and the worst of the thirty dashed above it, so the spread is visible rather than averaged away.

Gaussian entries, and saying so. Growth statistics depend on the distribution the matrices are drawn from. These are independent standard normals, which is the usual choice; matrices from other distributions behave differently, and a figure that said only “random matrices” would be hiding a decision. One matrix is an anecdote is the seeded-generator habit this site carries from its sibling, and it is exactly what makes a claim like “about three” supportable.

A refusal. The measurement must be able to see large growth, or the finding that growth is always small is a property of the measurement rather than of the matrices. So on every build the growth routine is shown Wilkinson’s matrix at n = 20 and required to report a value above a hundred. Without that, a bug that clamped or mis-scaled the ratio would produce exactly the reassuring figure above and nothing would notice.

Growth and the backward error, connected #

It is worth closing the loop, because the growth factor is not interesting in itself — it is interesting because it appears in a bound on something that is.

Across the sweep of conditioning in a small residual is not a small error, the backward error of every solve sits below 10⁻¹⁴ regardless of the condition number. That flat line is the growth factor being small, expressed in the quantity that matters. If growth were behaving like 2ⁿ⁻¹, the line would rise with n and the whole framework would collapse: an algorithm whose backward error depends on the problem is not backward stable, and the separation between the algorithm’s fault and the problem’s would not exist.

So the empirical smallness of the growth factor is load-bearing for this entire site. It is the reason the flat line is flat. It is also, uncomfortably, the one thing here that rests on observation rather than on proof.

The classes where the worst case is impossible #

The empirical argument is the weaker of the two reasons nobody worries, and it is worth setting out the stronger one because it covers a great deal of real computation.

Symmetric positive definite. Cholesky needs no pivoting and has growth factor at most 1. Every covariance matrix, every stiffness matrix from a finite element model, every Gram matrix falls here.

Diagonally dominant. If each diagonal entry exceeds the sum of the absolute values of the others in its row, the growth factor under no pivoting at all is at most 2. Discretised diffusion problems and many Markov chains are diagonally dominant by construction.

Totally nonnegative. Every minor nonnegative, growth factor exactly 1. Rarer, but it includes some interpolation matrices.

Upper Hessenberg. One subdiagonal, growth factor at most n. This matters because the QR algorithm reduces matrices to Hessenberg form first, so every step of the standard eigenvalue routine is operating in a regime with a linear bound.

For matrices in any of those classes the 2ⁿ⁻¹ worry is not merely unlikely, it cannot happen. The open question concerns only the general dense case, which in practice means matrices whose structure is unknown rather than absent.

How the growth factor enters the error bound #

The chain is worth writing out once, because the growth factor is otherwise a quantity with no visible consumer.

Wilkinson’s analysis says the computed factors satisfy L̂Û = A + E with ‖E‖ ≤ c·n·ρ·u·‖A‖, where c is a small constant. That E is exactly the perturbation the backward error measures: the factorisation is exact for a matrix E away from the one it was given.

So the backward error is proportional to n·ρ·u, and with ρ near 1 and n a few hundred, that is a few hundred roundings — 10⁻¹⁴ or so in double, which is what a small residual is not a small error measures across its whole sweep.

If ρ were 2ⁿ⁻¹ the same expression would give a backward error of order 1 for n above about 50, and Gaussian elimination would be unusable for any matrix larger than a small one. The entire practice of dense linear algebra rests on a quantity nobody can bound usefully and everybody has measured.

What to measure in practice #

Three numbers, each one line, and between them they cover most of what can go wrong with a dense solve.

The growth factor: track the largest absolute entry seen during elimination and divide by the largest entry of A. Above ten, look harder; above a thousand, something is unusual about the matrix.

The factorisation residual, ‖PA − LU‖/‖A‖, which is the site’s rule and which catches permutation bugs, indexing errors and the whole class of mistakes that produce a plausible-looking factorisation of the wrong matrix.

A condition estimate, which turns the residual into an error bound and is available from the factorisation for a fraction of its cost.

None of the three is returned by default by most libraries, and all three are cheap. The habit of computing them is most of what separates code that fails loudly from code that fails quietly, and this site’s rule — no decomposition without its residual — is that habit made compulsory for figures.

What this says about worst-case analysis #

A closing observation that is not really about matrices.

The 2ⁿ⁻¹ bound is a correct, tight, worst-case result, and it is a bad guide to behaviour. The matrix attaining it is measure-zero in any reasonable sense and does not arise from applications. An engineer who chose an algorithm by comparing worst-case bounds would reject the method the entire scientific computing world uses, in favour of one three times as expensive, on the strength of a scenario that has never occurred.

The lesson is not that worst cases do not matter — in a security setting they matter more than anything, because an adversary chooses the input. It is that a worst-case bound answers the question “what could an adversary do”, and most numerical computing has no adversary. When the inputs come from physics or from measurement rather than from someone trying to break the code, the distribution is the thing to measure, and the bound is a backstop.

Which is precisely why it is worth having both numbers on the same axes, which is what the figure at the top of this essay does. The bound and the behaviour are both true, they differ by eleven orders of magnitude, and knowing only one of them leads to a wrong decision in either direction.

What this field established #

Three essays, one algorithm, and a set of facts that the rest of the site leans on.

Elimination is a sequence of choices established what the factorisation records and why the multipliers are bounded by one. The swap that is not optional established what happens when they are not — a backward error of 0.25 on a well-conditioned problem, arrived at silently. And this essay establishes that the bound the whole argument rests on is never approached, which is the empirical fact everything else quietly assumes.

The thing that carries forward is the flat line. Every figure on this site that plots backward error against condition number shows a horizontal line at 10⁻¹⁶, and that line is the growth factor being near one, expressed in the quantity the reader cares about. It is what makes the separation between the algorithm’s fault and the problem’s available at all, and it is the reason a small residual is not a small error can make its argument without qualification.

The open question, stated plainly #

It is worth ending on, because it is one of the more surprising gaps in a mature field.

Nobody knows why the growth factor of Gaussian elimination with partial pivoting is small in practice. There is no theorem. There are results for structured classes, there is a large body of empirical measurement, and there is a probabilistic analysis for random matrices under specific distributions — and there is no general statement covering the matrices people actually factorise.

The most widely used numerical algorithm in existence rests, in the general case, on an observation. That is not a scandal; it is what an empirical science looks like when the theory has not caught up. But it is worth knowing, and it is worth knowing that the observation has been made carefully, many times, by people who would have been delighted to find a counterexample.

Which is also the argument for measuring it rather than assuming it. One line, one comparison per step, and the assumption becomes a number.