The exact answer to a nearby problem

Here is a computation that fails, and the interesting part is not that it fails.

Take the ten-by-ten Hilbert matrix — the matrix whose entry in row i and column j is 1/(i+j−1), which is as unremarkable a matrix as anyone could write down. Build a right-hand side so that the answer is exactly the integers one to ten. Solve it in double precision by Gaussian elimination with partial pivoting, which is the algorithm every library uses and the best available for a dense system.

The answer comes back wrong in its fourth decimal place. The relative error is 2.7·10⁻⁴, which is to say that of the sixteen significant digits double precision offers, twelve have gone.

Now the interesting part. Take the computed answer, substitute it back, and ask what problem it solves exactly. The answer is: a Hilbert system whose matrix and right-hand side differ from the ones given by a relative amount of 10⁻¹⁷.

The algorithm did not make a mistake. It solved a neighbouring problem perfectly, and the neighbouring problem’s answer is genuinely different from the one asked for.

Two definitions, and why the second one is the useful one #

Forward error is the obvious quantity: how far the computed answer is from the true answer, ‖x̂ − x‖/‖x‖. It is what the user cares about, and it is almost never computable — knowing it requires knowing the true answer, which is the thing being sought.

Backward error is the other direction: the smallest perturbation of the problem for which the computed answer would be exactly right. For a linear system, the smallest relative change to A and b such that (A + δA)x̂ = b + δb. It is computable, cheaply, from quantities already to hand — Rigal and Gaches gave the formula, and it amounts to the residual scaled by the sizes of the matrix, the answer and the right-hand side.

The reason backward error is the more useful of the two is that it isolates the algorithm. The forward error mixes together how well the algorithm did with how sensitive the problem is, and those are different things with different remedies. The backward error is purely the algorithm’s contribution, and it is measurable without knowing anything that is not already to hand.

An algorithm whose backward error is always at the level of the arithmetic’s rounding is called backward stable. That is the right standard to hold a numerical method to, and it is the strongest thing that can honestly be promised: no algorithm can give a small forward error on a problem where a small perturbation of the input changes the answer, because it cannot distinguish the given input from the perturbed one.

The identity that connects them #

The two are related by the thing that makes the problem hard or easy:

forward error ⪅ condition number × backward error

That is the whole framework. The condition number κ belongs to the problem — it is how much a relative perturbation of the input can be amplified into a relative perturbation of the output, and it exists whether or not any computer is involved. Multiplied by the backward error it gives a bound on what will be seen.

For the Hilbert solve above: κ = 1.6·10¹³, the backward error is 8·10⁻¹⁷, and the product is 1.3·10⁻³. The measured forward error is 2.7·10⁻⁴, comfortably inside the bound, as it should be — this is an inequality and the worst case is not always attained.

That flat line is the most important thing on this site. It says the elimination is not degrading. It never degrades. What degrades is the relationship between the problem as posed and its answer, and that had nothing to do with the elimination.

What this changes about debugging #

The practical consequence is a decision procedure, and it is worth stating baldly because it saves an enormous amount of time.

A computation has produced an answer nobody believes. Compute the backward error, which costs one matrix–vector product.

If it is large — much bigger than the unit roundoff times the size of the problem — the algorithm is at fault. Something is unstable, or there is a bug. This is fixable, and it is fixable by changing the algorithm.

If it is small — at the level of rounding — the algorithm did everything it could. The answer is wrong because the problem is sensitive, and the only remedies are more precision, a better conditioned formulation of the same question, or accepting that the quantity is not determined by the data to the accuracy that was wanted.

Those are entirely different situations and they are confused constantly. A great deal of effort has been spent making stable algorithms more stable, and a great deal of frustration comes from expecting an ill-conditioned problem to yield to a better solver.

How the backward error is measured here #

The formula used throughout this site is the normwise relative backward error,

η = ‖b − Ax̂‖ / (‖A‖·‖x̂‖ + ‖b‖),

which Rigal and Gaches showed is exactly the smallest relative perturbation of both A and b that makes x̂ exact. Three points about it are worth making, because each has caught someone.

The denominator matters. The residual alone, ‖b − Ax̂‖, is not a backward error: it has the units of b and scales with the size of the problem. Divide by the natural scale and it becomes a pure number that can be compared against the unit roundoff.

It is always computed in double. Where a figure shows a factorisation performed at sixteen or twenty-four bits, the measurement of its error is still done at fifty-three. Measuring a single-precision residual in single precision would report the error as zero, because the error is exactly the size of what single precision cannot see. This is the numerical equivalent of measuring a ruler with itself, and it is a single line in lib/matrix.js that keeps every figure on this site honest.

It has to be able to fail. A backward-error routine that returned something small no matter what would make every claim on this page vacuous, so it is fed a deliberately wrong answer on every build — a vector of ones, in place of a real solution — and required to report a large number. The thread assertions that reject is that discipline applied throughout.

A worked case, end to end #

It is worth doing one all the way through with the numbers visible, because the framework is abstract until it has been applied once.

The problem: solve H₁₀x = b where H₁₀ is the Hilbert matrix and b was constructed, in exact rational arithmetic, so that x is the integers one to ten.

What the arithmetic contributed. Partial pivoting performs about 330 multiply-and-subtract operations, each rounding at the level of 10⁻¹⁶. The accumulated effect, measured as ‖PA − LU‖/‖A‖, is 2·10⁻¹⁷ — smaller than a single rounding, because the roundings do not all point the same way.

What the algorithm’s answer is exact for. The backward error is 8·10⁻¹⁷. There exists a matrix within 8·10⁻¹⁷ of H₁₀, relatively, for which the computed vector is the exact answer.

What the problem does to that. κ(H₁₀) = 1.6·10¹³. A relative perturbation of 8·10⁻¹⁷ in the input can produce a relative perturbation of up to 1.3·10⁻³ in the output.

What was seen. A forward error of 2.7·10⁻⁴, which is a fifth of the bound.

Every one of those five numbers was computed rather than quoted, and three of them were computed twice by routes that share no arithmetic. The conclusion is not that double precision was insufficient — it is that the map from this b to this x magnifies relative error by thirteen orders of magnitude, and that fact is a property of the Hilbert matrix and would be true of a computation done by hand.

“Nearly a third” is measured at one κ, and a fraction stated once is a fraction that might be about that κ.

Seventy-seven of a hundred, and thirty-five of a hundred, are the two numbers to watch as the conditioning worsens — the question is whether they stay proportional or drift:

Two orders of magnitude of conditioning, and both numbers have moved by exactly two orders of magnitude — which is the first evidence that what the figure is really reporting is a pair of fractions rather than a pair of amplifications.

Both fractions are constants, and they settle by the fourth row. The worst of two hundred random directions reaches 0.761 κ and the median reaches 0.288 κ, to three digits, at every κ from 10⁵ to 10⁸ — and the first three rows are converging on the same two numbers from above. So the caption’s “nearly a third” is not a fact about κ = 10⁶; it is 0.288 everywhere, and the essay can say so.

The worst case is 76% of κ and not κ. Two hundred directions is enough to find three quarters of the available amplification and not enough to find all of it, which is what makes κ a bound rather than a prediction — and the shortfall is a stable 24% rather than something that closes as the problem gets worse.

And the spread between typical and worst is 2.64, flat. Not 100, not κ: a caller who samples one random perturbation and one adversarial one sees a factor of two and a half between them, whatever the conditioning. That is the number worth carrying, because it is the answer to “how unlucky do I have to be” — and the answer is barely.

When the backward error is genuinely large #

Backward stability is not automatic. Three algorithms on this site are not backward stable, and each was chosen for that reason.

Elimination without a row swap has a backward error of 0.25 on a two-by-two, which is to say the answer it returns is exact for a problem 25% away from the one it was given. That is not a small perturbation of anything.

The normal equations are not backward stable for least squares. The computed solution is the exact answer to a perturbed problem only if the perturbation is allowed to be of size κ·u rather than u, which for an ill-conditioned fit is a completely different problem.

Classical Gram–Schmidt produces a Q that is not orthogonal to anything like rounding, so any algorithm downstream that treats QᵀQ as the identity is solving a different problem from the one it believes it is solving.

In every one of those cases the fix is a different algorithm, it is known, and it costs nothing. That is what makes the backward-error question worth asking first.

Two routes to the same statement #

The claim that an algorithm is backward stable is a claim, and this site tries not to make claims it has not measured twice.

The first route is the residual computation above, applied to the computed answer. The second is to reconstruct the factorisation and measure it directly: for a solve by LU, ‖PA − LU‖/‖A‖ says how far the factors are from factorising the matrix they were given. For a QR, ‖A − QR‖/‖A‖ does the same job. Those two numbers arise from different arithmetic — one from a residual, one from a matrix product — and they must both be at rounding level or the claim is not supported.

The word “nearby” is doing work #

One subtlety, because it is where the framework is most often misapplied.

“The exact answer to a nearby problem” is only reassuring if a nearby problem is as good as the one asked about. For a system whose entries come from measurement, it is: A is already uncertain in its fourth digit, so an algorithm that perturbs it in its sixteenth has introduced nothing that was not already there. The backward error can be compared directly against the uncertainty in the data, and if it is smaller — which for a stable algorithm in double precision it overwhelmingly is — the computation has added nothing.

For a system whose entries are exact — integers, or a matrix defined by a formula — the argument is weaker, because there is no data uncertainty to hide behind. The Hilbert matrix is exactly this case, which is why it is used here so often: its entries are exact rationals, so a perturbation of 10⁻¹⁷ is a real perturbation and not merely a small one relative to what was known.

Even so the conclusion holds, and the reason can be measured rather than argued, because an answer that is known supplies both Hilberts — the exact rationals and the doubles a program actually holds. Solving Hx = 1 three ways, and measuring the second and third against the first:

The middle column is what a perfect solver returns given the matrix a program can hold — an exact rational solve of the rounded entries, with no elimination error at all.

An algorithm that made no rounding error whatsoever would be wrong by between a fifth and all of what the ordinary one is wrong by, and it is within a factor of 4.6 of LU at every size. So on this matrix the elimination is not where the error comes from. Converting 1/(i+j+1) to a double is, and that happens before any algorithm runs.

The n = 4 row deserves a sentence of its own. There the ordinary double solve is more accurate than an exact solve of the matrix it was given — 1.01·10⁻¹³ against 1.41·10⁻¹³. That is not a virtue of the algorithm. It is two errors of opposite sign partly cancelling, and a reader who saw only that row would draw exactly the wrong conclusion about which stage to improve.

So the framework’s reassurance survives on an exactly specified matrix, for a reason that has nothing to do with data uncertainty and is stronger than “higher precision would not help enough”. Removing the algorithm’s error entirely, at the precision in hand, changes the answer by less than a factor of five. There is no algorithm to blame, because the perturbation the framework is about had already happened when the matrix was written down.

That last point generalises further than the Hilbert matrix, and it is the reason to state it as a column of a table rather than as a remark. Any problem specified by a formula whose values are not binary fractions — a covariance from a kernel, a stiffness matrix from an integral, a Vandermonde on non-dyadic nodes — is perturbed at the moment it is assembled, by a relative amount the assembly cannot avoid. A backward-error analysis of the solver then accounts for the second perturbation and is silent about the first, and the two are the same size for a stable algorithm.

Which means the honest budget has three terms rather than two. The data’s own uncertainty, where there is any; the perturbation the representation introduces, which there always is and which no algorithm can undo; and the algorithm’s own backward error, which is the only one a solver’s documentation talks about. On the Hilbert matrix the middle term is most of the total, and the measurement above is the only way to know that, because every quantity a solver reports lumps the first two together into “the problem”.

Why this is not the usual framing #

Most introductions to the subject present the condition number first, as a property of a matrix, and then treat numerical error as a separate topic about rounding. That ordering is defensible and it produces a predictable confusion: readers come away believing that ill-conditioned problems are ones where floating point misbehaves, and that a sufficiently careful implementation would fix them.

Putting backward error first inverts it. Rounding becomes the small, well-understood, bounded thing it actually is; conditioning becomes a property of the question being asked; and the interesting work is in the multiplication between them. It also makes several later results land differently.

The bound that is never attained is about the growth factor, which is the one place where Gaussian elimination’s backward stability is not guaranteed — the bound permits an exponential blow-up that random matrices miss by eleven orders of magnitude. That is a statement about the algorithm, and it is interesting precisely because everything else in this field is a statement about the problem.

Symmetry is worth more than precision is the opposite extreme: an eigenvalue problem where a perturbation of 10⁻¹⁶ moves the answer by 10⁻², with a backward-stable algorithm and no arithmetic anywhere at fault. The perturbation there is not even rounding error introduced by a computation — it is the rounding that occurred when the matrix was stored.

What to take from this #

Three sentences, and the rest of this field is their consequences.

A wrong answer has two possible authors, the algorithm and the problem, and they are separately measurable. The algorithm’s contribution is the backward error, it is cheap to compute, and for every method used on this site except the three named above it sits at the level of rounding no matter how hard the problem is. The problem’s contribution is the condition number, it belongs to the mathematics and not to any computation, and no algorithm can remove it.

A small residual is not a small error is the same statement approached from the direction people usually arrive from, which is having computed a residual and found it reassuring. It is not reassuring on its own, and the essay is about what it takes to make it so.

One last distinction #

Backward stability is a property of an algorithm, not of a run. Saying “this solve was backward stable” is shorthand for “this algorithm has small backward error on all inputs, and this run is an instance”.

That matters when a measured backward error comes out small on a method that is not backward stable. Elimination without pivoting has a backward error of zero on most matrices — it is only the ones with a small leading pivot where it fails, and they are a minority. Measuring one run and concluding the method is fine is exactly the error the essay is about.

So the two uses of the number are different. On a method, it is evidence to be gathered across inputs. On a run, it is a diagnostic that says which half of the problem to look at, and that is what it is used for here.

Two routes, and only one of them earns this sentence #

The property this page defines is not automatic and it is not implied by accuracy. Two routes to the same solution, from the same factorisation, can have forward errors within a factor of fifty and backward errors twelve orders of magnitude apart.