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

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.
16 min read 6 figures Whose fault is itMeasured, not assertedExact ground truth

Suppose a solve of Ax = b and has produced an answer, and the question is whether the answer is right. The obvious check is to substitute it back: compute b − Ax̂ and see how small it is. It costs one matrix–vector product, it needs nothing that is not already to hand, and if the result is 10⁻¹⁷ then surely the answer is good.

It is not. The residual can be 10⁻¹⁷ and the answer can be wrong in its second significant figure, and this is not a contrived circumstance — it happens on a matrix anyone might write down, solved by the best available algorithm, in double precision.

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. 1 The case, without interpretation. A thirteen-by-thirteen Hilbert system, with a right-hand side built in exact rational arithmetic so that the true answer is the integers one to thirteen. The relative residual of the computed answer is 8.6·10⁻¹⁷. The eighth component is 0.049 where 8 belongs.

What the residual actually measures

Write the residual as r = b − Ax̂. Then A(x − x̂) = r, so the error in the answer is A⁻¹r.

That is the entire explanation. The residual is the error pushed through A, and getting the error back out requires pushing it through A⁻¹ — which multiplies it by as much as ‖A⁻¹‖. When A is nearly singular that is an enormous number, and a residual small enough to be invisible becomes an error large enough to be catastrophic.

The relationship, in the relative terms that are worth thinking in:

‖x − x̂‖/‖x‖ ⪅ κ(A) · ‖r‖/‖b‖

So the residual bounds the error only after multiplication by the condition number, and the condition number is precisely the quantity that was unknown when the decision was made that the residual was reassuring.

For the Hilbert system above, κ is 1.7·10¹⁸. A relative residual of 8.6·10⁻¹⁷ permits a relative error of about 150 — which is to say, no constraint at all. And indeed the answer has components that are wrong by 200%.

The residual measures the algorithm, and that is worth something

None of this makes the residual useless. It makes it a measurement of a different thing than most people take it for, and that thing is genuinely valuable.

A small residual says the computed answer is the exact answer to a nearby problem. It is, up to a scaling, the backward error, and it certifies the algorithm. A computed residual of 10⁻¹⁷ establishes that: whatever went wrong, the elimination did not do it. That is a real conclusion and it eliminates the entire class of explanations involving bugs, instability, and bad pivoting.

What it cannot do is certify the answer, and the reason it cannot is that certifying the answer requires knowing something about A that the residual does not contain.

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 two quantities, measured across twelve orders of magnitude of conditioning. The backward error — essentially the scaled residual — does not move. The forward error rises in step with κ. At the right-hand end they are eleven orders of magnitude apart, and everything about the computation that produced them is identical.

Three residuals that are all small

It helps to see the failure from the other side: cases where the residual is small and the answer is fine, and cases where it is small and the answer is not, with nothing distinguishing them except a number that has to be computed separately.

A well-conditioned system. κ = 100, residual 10⁻¹⁶, error 10⁻¹⁴. The residual was reassuring and it was right to be.

A moderately ill-conditioned system. κ = 10⁸, residual 10⁻¹⁶, error 10⁻⁸. Eight digits gone, and the residual is the same as in the first case.

The Hilbert system. κ = 10¹⁸, residual 10⁻¹⁷, error 2. The residual is the smallest of the three.

The residual is not correlated with the error at all. Across the whole sweep in the figure above, the residual varies by less than a factor of ten while the error varies by eleven orders of magnitude. Whatever the residual reports, it is not reporting the error.

Where the intuition comes from, and why it is wrong

The intuition behind trusting a residual is geometric and it is nearly right, which is what makes it durable. It goes: the true solution is the point where Ax − b is zero; my computed point makes Ax̂ − b very nearly zero; therefore my computed point is very nearly the true solution.

The missing step is that “very nearly zero” in the output of A does not mean “very nearly there” in the input of A. A maps a large region of input space onto a small region of output space whenever it has a small singular value, and the residual lives in output space.

The picture that fixes this is a badly scaled one. Imagine A is diagonal, with entries 1 and 10⁻¹². Then Ax − b is a two-component vector, and the second component is the second component of the error divided by 10¹². A point can be a trillion units away from the true answer in that direction and produce a residual of size one. The residual is not lying; it is reporting the size of the discrepancy in the space where the discrepancy was measured, and that space has been squashed.

Every ill-conditioned matrix is that diagonal one in disguise. The singular value decomposition says so precisely: A = UΣVᵀ, so the squashing happens along the right singular vectors, with a factor for each singular value. Rank is a decision is about what happens when the smallest of those is small enough that the direction is effectively unresolvable.

The scaling trap, which is the same trap wearing a hat

There is a version of this that catches people who already know everything above.

The residual is usually reported relative to ‖b‖, which is sensible. But a system can be rescaled — multiply row i by any nonzero constant, and the solution is unchanged while the residual is not. A badly scaled system can therefore be made to show any residual at all without altering the answer at all.

This has two consequences worth carrying. First, a residual is only comparable across systems if the scaling is comparable, so “this solver achieves residuals of 10⁻¹⁵” is a statement about a benchmark suite’s scaling as much as about the solver. Second, the componentwise backward error — how much each entry of A and b would have to change, relative to its own size — is a stronger and often more honest quantity than the normwise one used here, and it is the right tool when the entries of A differ in magnitude by many orders.

This site uses the normwise measure throughout, because it is the one with the clean Rigal–Gaches characterisation and because the matrices here are deliberately not badly scaled. That is a choice, and stating it is part of the point: a residual is meaningless without saying which norm and which scaling produced it.

The one thing that makes it sufficient

There is a case where a small residual does bound the error usefully, and it is worth knowing because it is common: when an estimate is available of the condition number.

Estimating κ costs about as much as the solve itself if done honestly through the singular values, and much less through one of the standard estimators, which use a few solves with cleverly chosen right-hand sides to find a lower bound on ‖A⁻¹‖. LAPACK returns one on request; it is rcond, and almost nobody looks at it.

With κ in hand the residual becomes a genuine error bound. Without it, the residual is a statement about the algorithm and silence about the answer.

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. 3 What the missing factor looks like. A known perturbation is applied to the right-hand side in two hundred random directions and the amplification is measured; the worst reaches κ, and the median reaches nearly a third of it. This is the multiplier standing between a residual and an error, and it is a property of the matrix that no amount of looking at the residual will reveal.

Iterative refinement, and why it works at all

There is a classical trick that appears to contradict everything above, and understanding why it does not is worth the paragraph.

Iterative refinement: solve Ax̂ = b, compute the residual r = b − Ax̂, solve Ad = r using the same factorisation, and set x̂ ← x̂ + d. Repeat. It frequently recovers most of the lost digits.

If the residual carries no information about the error, how can correcting by it help?

The answer is that it does carry the information — A⁻¹r is the error, exactly. What it does not do is make the error *visible* without applying A⁻¹, and applying A⁻¹ is what the refinement step does. The reason refinement is not free is that the second solve has the same conditioning as the first, so the correction d is itself computed with a large relative error. What saves it is that d is small: a 10⁻⁴ relative error on a correction of size 10⁻⁴ is an absolute improvement of nearly the full amount.

The classical result is that refinement works if the residual is computed in higher precision than the solve. With the residual computed at the same precision, refinement improves the backward error but is limited in what it can do for the forward error. This is why mixed-precision iterative refinement — factorise in single, refine in double — has become interesting again on hardware where low precision is much faster, and it is one of the topics this site’s expansion is planned around.

What is measured here, and how

Two claims on this page are strong enough to deserve saying how they are checked.

“The backward error never leaves the rounding level.” Across the eight condition numbers in the figure, twenty seeded matrices each, the worst backward error observed is below 10⁻¹⁴. That tolerance is bracketed: the measured noise floor across all 160 solves is far below it, and a backward error of 10⁻⁶ would be a genuine failure of the algorithm, so the threshold sits strictly between something that cannot fire on correct code and something that would miss a real defect.

“The forward error is bounded by κ times the backward error.” Checked at every point, not at one, with the factor of safety written into the assertion rather than into a comment. The bound is an inequality, so the assertion is an inequality, and the figure would fail the build if any measured forward error exceeded it.

And the measurement of the residual is itself checked to be capable of failing: it is fed a deliberately wrong answer on every build and required to report a large number. Assertions that reject is the thread; the specific worry here is that a residual routine with a scaling bug would report small numbers for everything, and every figure on this page would look exactly the same.

What this means for testing numerical code

The advice that follows is unglamorous and it is the part that changes what people do.

A test that asserts a small residual is testing the factorisation, not the solver’s accuracy. That is worth having — it catches a great many real bugs, and it is fast — but it should be named for what it does. A test suite consisting entirely of residual checks will pass on a solver that is producing garbage on every ill-conditioned input.

Build at least one test case whose answer is known exactly. Not a case where the answer was computed by a more careful method: a case where the answer is known by construction. Multiply a matrix by a chosen integer vector to make the right-hand side, and the answer is that vector. An answer that is known does this in exact rational arithmetic so that even the right-hand side is free of rounding, which closes the last hole in the argument.

Report the condition number in the test output. If a test on an ill-conditioned matrix starts failing at 10⁻⁶ where it used to pass at 10⁻⁸, the useful question is whether the matrix got harder or the code got worse, and κ answers it in one line.

Forward error of a 6×6 Hilbert solve at eight precisionsA bar for each significand width from 12 to 53 bits showing the relative error in the computed solution, with the condition number times the unit roundoff marked as a prediction.κ = 1.5·10⁷ · the exact answer is (1, 2, …, 6)12 bits3.316 bits0.8620 bits4.324 bits0.09330 bits7.6·10⁻⁴36 bits6.2·10⁻⁶43 bits2·10⁻⁷53 bits4.5·10⁻¹¹dashed: κ · unit roundoffone matrix, eight arithmeticsmeasured against a known answer
Fig. 4 The same system solved at eight precisions, with the error measured against the exact answer rather than against a residual. The bars follow κ·u across eleven orders of magnitude — the prediction, not a fit. A residual-only test would show essentially the same small number for every one of these eight runs.

Two other residuals on this site

The word residual is used in two more places here, and they are worth distinguishing because they answer different questions.

The factorisation residual, ‖PA − LU‖/‖A‖ or ‖A − QR‖/‖A‖, asks whether the factors multiply back to the matrix they came from. This is the quantity the site’s rule is about, and it is printed on every figure that draws a decomposition. It is a check on the factorisation and not on anything downstream — orthogonal is a number contains the sharpest example of that gap, where a factorisation reconstructs its matrix to 10⁻¹⁷ while the Q it produced is not orthogonal in any useful sense.

The least-squares residual, ‖b − Ax‖ at the minimising x, is not an error measure at all — it is the answer. Its size says how well the model fits the data, and a large one means the data is not close to the column space, which may be entirely correct. The valley with no bottom is about the fact that this residual is insensitive to the coefficients, which is exactly why the coefficients are hard to determine.

The three uses share a name and almost nothing else. Keeping them apart is worth the effort.

The practical rule

With only the residual in hand, what is known is the algorithm and nothing about the answer.

Knowing about the answer requires a second number, and there are three places to get one. Estimate the condition number, which is the general answer and costs a fraction of the solve. Perturb the input and re-solve, which measures the sensitivity directly and costs a solve per direction — the condition number is an amplifier is about doing exactly that. Or construct a problem whose answer is known, which is only available in special cases but is completely conclusive when it is — an answer that is known is that case.

What is not available is inferring the error from the residual alone, and the reason it is worth being firm about that is that the residual is so cheap, so natural, and so nearly always the only check performed.

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. 5 The relationship in one picture, with the numbers from a ten-by-ten solve. The residual lives on the top edge, and it is small. The error lives on the bottom edge, and it is not. Reading one off the other requires the number down the right-hand side, and it is not in the residual.

This is also the reason the rule this site is named for is about residuals rather than about errors. No decomposition is drawn without its residual printed is a claim about the factorisation — that it factorises the matrix it was given — and that is exactly what a residual can establish. It says nothing about whether the answer downstream is any good, and every figure that shows an answer prints the second number too.

And the corollary for figures

The rule this site runs on follows directly, and it is worth stating as a consequence rather than as a preference.

No decomposition is drawn without its residual printed — because the residual is exactly what a factorisation figure can honestly claim. It says the factors factorise the matrix they were given, which is a complete statement about the thing being drawn.

Every figure that shows an answer carries a second number as well, and it is a different number in each case: a forward error where the truth is known, a condition number where it is not, a backward error where the point is the algorithm. That asymmetry is deliberate, and it is the shortest summary of this essay.

Growth factor under partial pivoting: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8
Fig. 6 The quantity behind the flat backward-error line. Growth stays near three for random matrices where the bound permits 5.5·10¹¹, and that empirical smallness is what keeps the residual small on every matrix — including the ones where the answer is worthless.

The two neighbouring essays complete the argument. The condition number is an amplifier supplies the missing factor, measured by perturbing rather than defined; and the bound that is never attained is about the one quantity that could make the flat line stop being flat, and about the fact that nobody can prove it will not.

Where this goes next

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 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. 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.