Numerical Linear Algebra
where the algebra stops being the arithmetic
Two errors, and whose fault they are

The exact answer to a nearby problem

A good algorithm does not give an approximate answer to your problem. It gives the exact answer to a problem very close to yours — and once that is the definition, a wrong result has two possible authors and they can be measured apart.
16 min read 6 figures Whose fault is itWhere the algebra stops being the arithmeticMeasured, not assertedTwo routes to a number

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.

Backward error, forward error and the condition number, with measured valuesTwo boxes at the top — the problem posed and the nearby problem the algorithm answered exactly — and two answers below them, with the distances between all four labelled by numbers from a Hilbert solve.the problem you posedA = H10b = A·(1, 2, …, 10)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 2.3·10⁻¹⁷the answer you wantedx = (1, 2, …, 10), exactlythe answer you gotx̂, wrong by 2.7·10⁻⁴ relativebackward error 2.3·10⁻¹⁷forward error 2.7·10⁻⁴κ = 1.6·10¹³κ · η = 3.6·10⁻⁴, and the measured forward error is 2.7·10⁻⁴.The algorithm is not at fault. The problem is.H10, LU with partial pivotingresidual and error differ
Fig. 1 The two errors and the number connecting them, with values from that solve. Across the top: how far the problem the algorithm answered is from the one it was given. Down the right: how much this particular problem amplifies a perturbation. Along the bottom: the product, which is what anyone actually sees.

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.

Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 2 The identity as a measurement rather than a formula. Matrices are built to a prescribed condition number, from ten to ten to the thirteenth, and solved. The backward error is a flat line at the level of rounding — the algorithm is exactly as good on the hardest problem as on the easiest. The forward error rises in step with κ. Twenty matrices at each point, median drawn solid, worst of the twenty dashed.

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.

The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 3 The second situation at its most extreme. A thirteen-by-thirteen Hilbert system whose exact answer is the integers one to thirteen; the computed one has 0.05 where 8 belongs, 21.3 where 9 belongs and −2.6 where 10 belongs. The backward error is 2.2·10⁻¹⁷. There is no better algorithm to reach for.

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.

How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 4 The amplification, measured by doing it. The right-hand side of a system is perturbed in two hundred random directions by a known relative amount, and the amplification of the answer is recorded. No direction exceeds κ — it is a worst case — but the median direction reaches nearly a third of it, so ill-conditioning is not a rare misfortune you have to be unlucky to encounter.

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.

Gaussian elimination on a 4×4, one step at a timeFour copies of the same matrix: as given, and after each of the three elimination steps. The pivot in use is outlined and the entries reduced to zero are greyed.21-13-3-121-212-443-12as givenrows in the order 1 2 3 443-1201.251.252.502.51.5-30-0.5-0.52after step 1pivot 443-1202.51.5-3000.5400-0.21.4after step 2pivot 2.543-1202.51.5-3000.540003after step 3pivot 0.5‖PA − LU‖/‖A‖0largest multiplier0.75row order 4 3 2 1the pivot is chosen
Fig. 5 The second route, printed. Elimination step by step on a four-by-four, with ‖PA − LU‖/‖A‖ on the figure. This is the rule the whole site runs on: no decomposition is drawn without its residual, so a factorisation that has gone wrong says so in the picture rather than in a caption nobody reads.

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, because the sensitivity is so extreme that nothing representable in double precision would help. An answer that is known is about exactly that case, and about the fact that the true answer for it can be computed in exact rational arithmetic, so the two contributions can be separated without any interpretation at all.

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.

How far a perturbation of size ε moves an eigenvalue, 8×8A log–log plot of eigenvalue movement against perturbation size. The symmetric case lies on a line of slope one; the non-symmetric case lies on a line of slope one eighth, and at a perturbation of ten to the minus sixteen it has already moved by a hundredth.10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹size of the perturbation ‖δA‖how far the eigenvalues moveJordan block, ε^(1/8)symmetric, ≤ ‖δA‖rounding error alone moves it to 10⁻²six seeds per symmetric point; Jordan is closed formsymmetry beats precision
Fig. 6 Where the framework’s second factor has no upper limit worth having. For a symmetric matrix, an eigenvalue moves by at most the size of the perturbation — a bound with no constant in it. For an eight-by-eight Jordan block it moves by the eighth root, so a perturbation at the level of storing the matrix has already moved the spectrum by a hundredth.

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.

Where this goes next

110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method A small residual is not a small error Substituting the answer back and finding that it fits is the most natural check there is, and it verifies the wrong thing. A residual of 10⁻¹⁷ is entirely compatible with an answer whose second digit is wrong. 110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ The condition number is an amplifier κ is usually introduced as a definition and then quoted. It is a measurement: perturb the input by a known amount, look at how much the output moves, and the largest ratio you can find is the number. H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known An answer that is known Almost every demonstration of numerical error estimates the error by computing the same thing more carefully. The Hilbert matrix does not need that: its inverse is a closed form in integers, so the true answer is available exactly and the error is measured rather than approximated.