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

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.
15 min read 5 figures Exact ground truthTwo routes to a numberWhose fault is itMeasured, not asserted

There is a circularity at the centre of most writing about numerical error, and it is polite enough that it usually passes unnoticed.

The demonstration goes: here is a computation, here is the same computation done more carefully — in higher precision, or with a better algorithm — and the difference between them is the error. It is not the error. It is the difference between two approximations, and it is a lower bound on the error of the worse one only if the better one is genuinely better, which is exactly the thing being assumed.

For one family of problems the circularity can be cut, and this site leans on it heavily.

The Hilbert matrix has an integer inverse

The Hilbert matrix Hₙ has entries 1/(i+j−1). Every entry is a unit fraction; nothing about it looks remarkable. Its inverse has a closed form:

(H⁻¹)ᵢⱼ = (−1)^(i+j) (i+j−1) C(n+i−1, n−j) C(n+j−1, n−i) C(i+j−2, i−1)²

and every entry of it is an integer. That is the part that surprises: the inverse of a matrix made entirely of fractions has no fractions in it at all. For n = 8 the largest entry is 4,249,941,696. For n = 13 it is 106,518,477,825,760,000 — a number with eighteen digits, which a double cannot hold exactly and BigInt can hold without effort.

So the true solution of any Hilbert system with a rational right-hand side can be computed with no rounding whatsoever, and compared against a floating-point solve. The error is not estimated. It is known.

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 comparison. A thirteen-by-thirteen system whose right-hand side was constructed in BigInt rationals so that the exact answer is the integers one to thirteen. Beside it, what double-precision Gaussian elimination with partial pivoting returns. The number of correct digits is listed for each component; four of them have none.

What the comparison shows

Reading down that figure: the first three components are right to three or more digits. The eighth is 0.049 where 8 belongs. The ninth is 21.3 where 9 belongs. The tenth is −2.58 where 10 belongs — the wrong sign, on a problem whose exact answer is the positive integers.

And the backward error of that solve is 2.2·10⁻¹⁷.

Every part of that is worth holding at once. The elimination performed about 730 multiply-and- subtract operations and accumulated less rounding than a single operation’s worth. The answer it produced is the exact solution of a system within 2·10⁻¹⁷ of the one it was given. And it is wrong by 200% in four places.

κ(H₁₃) is 1.7·10¹⁸, which is above 1/u. When the condition number exceeds the reciprocal of the unit roundoff, the perturbation introduced by merely storing the matrix is enough to change the answer completely, and the computation is over before it starts.

Two exact routes, required to agree exactly

Having exact arithmetic available creates a temptation to trust it, and exact arithmetic can be implemented wrongly like anything else. So the two available exact routes are both computed and required to agree — not to a tolerance, exactly.

Route one: Gaussian elimination on the rational matrix, with every entry a BigInt fraction in lowest terms. Partial pivoting is retained, and it is worth noticing why: with no rounding, pivoting buys nothing in accuracy, and is needed only to avoid dividing by an exactly zero pivot. Separating those two jobs is something the floating-point version cannot do, and seeing the same code with the second job removed makes the point better than a paragraph.

Route two: multiply the right-hand side by the closed-form inverse above.

For an eight-by-eight, the two agree component by component, as exact fractions with identical numerators and denominators. There is no tolerance in that assertion because there is no rounding to tolerate. And the closed form itself is checked: H·H⁻¹ is computed in rationals and every entry is required to be exactly 1 or exactly 0 — not 1 ± 10⁻¹⁵, exactly the integer.

That is two routes to a number in its strongest available form, and it is one of the reasons this subject was worth a site: in most fields the second route is another approximation.

The conversion that nearly broke it

One implementation detail is worth recording because it is the kind of failure that produces an empty figure rather than an error message.

A rational with an eighteen-digit numerator converts to a double without difficulty. A rational arising from the intermediate stages of exact elimination on H₁₃ does not: the numerators and denominators grow, and by the last back-substitution step they run to several hundred digits. Number(n) / Number(d) on those gives Infinity / Infinity, which is NaN, which a plotting routine draws as nothing at all.

The fix is to strip the exponents before dividing — shift both BigInts down to about 200 bits, divide those, and multiply the scale back afterwards. It is four lines. What makes it worth mentioning is that the failure is silent in exactly the way this site keeps finding: no exception, no warning, a figure that renders and is empty. The assertion that now guards it converts a rational with a 400-digit numerator and requires the result to be finite and to equal what the arithmetic says it should be.

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 Why exact ground truth matters for the site’s central claim. This figure’s forward errors are measured against a known answer, constructed the same way. Without that, the rising line would be “the difference between this solve and a better solve”, and the flat line would be the only thing actually measured.

What it takes to build the right-hand side

The construction deserves a paragraph because getting it wrong is easy and would invalidate everything.

The obvious approach: pick x = (1, 2, …, n), compute b = Hx in double precision, then solve Hx̂ = b in double and compare x̂ against x. That does not work, and the reason is subtle. Computing b in double rounds it, so the b that the solver is given is not Hx for the chosen x — it is H times some slightly different vector, and on a matrix with κ = 10¹⁸ that slightly different vector is enormously different. The comparison would then be measuring the rounding in the construction rather than the rounding in the solve.

So b is built in exact rational arithmetic: H is the exact rational Hilbert matrix, x is a vector of exact integers, and the product is a vector of exact rationals. Only then is it converted to double, and that single conversion is the one unavoidable rounding — the one that says “this is the problem a computer can be given”.

The distinction is not academic. It is the difference between a figure that shows what elimination does to a hard problem and a figure that shows what rounding does to a right-hand side, and the two would look nearly identical.

Where else exact truth is available

Three more places on this site, and it is worth cataloguing them because the technique is more widely applicable than it looks.

Determinants of integer matrices. The pivots of an LU factorisation multiply to the determinant; so does cofactor expansion, which on an integer matrix is exact. Every elimination figure here checks one against the other, and cofactor expansion on a four-by-four costs nothing.

Eigenvalues of a Jordan block with a corner entry. The characteristic polynomial of the n×n nilpotent shift with ε in the bottom-left is exactly λⁿ − ε, so the eigenvalues are the n-th roots of ε and their modulus is known in closed form. Symmetry is worth more than precision uses that as its ground truth, and checks it by evaluating the characteristic polynomial at the claimed root.

Eckart and Young’s theorem, which says the best rank-k approximation has error exactly σₖ₊₁ — not bounded by it, equal to it. That is an unusually sharp statement, and it makes the theorem itself usable as a second route: the measured approximation error and the next singular value are computed independently and required to agree to nine digits. The best approximation there is does that.

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. 3 The determinant case, in passing. This figure’s badge carries the factorisation residual, and behind it the build checks that the product of the pivots equals the cofactor expansion of the same integer matrix — two computations sharing no arithmetic, required to agree to twelve digits.

Why the Hilbert matrix, specifically

It has three properties that together are rare.

It is exactly representable as rationals, so exact arithmetic is available at all. A matrix of measured data is not.

It is symmetric positive definite, so it is not a pathological object dressed up as a matrix — it is as well behaved as a matrix can be in every respect except conditioning. Nobody can dismiss the result by pointing at a structural defect.

Its condition number grows explosively with size, roughly like e^(3.5n), so a small matrix can be made arbitrarily hard. κ(H₄) = 1.5·10⁴, κ(H₈) = 1.5·10¹⁰, κ(H₁₂) = 1.8·10¹⁶. That means a single family provides the whole range of the site’s horizontal axis, with everything else held fixed.

It also arises naturally, which is worth knowing: the Hilbert matrix is exactly the Gram matrix of the monomials 1, x, x², … on the interval [0,1]. So fitting a polynomial in the monomial basis by the normal equations is solving a Hilbert system, and the road that squares the problem is where that comes back. The most famously ill-conditioned matrix in the subject is not a curiosity; it is what comes out by doing the most obvious thing.

What “known” is worth

The value of having an exactly known answer is not that it makes the error smaller. It is that it removes the last place an argument can hide.

Consider the alternative version of the figure above: solve H₁₃ in double, solve it again in some higher precision, and plot the difference. Every number would be similar. But the claim would be “double disagrees with quadruple by this much”, and the reader would be entitled to ask whether quadruple is right — and to note that both computations used the same algorithm, the same pivoting strategy and the same code path, so a systematic error would cancel and be invisible.

With exact rationals, none of that is available. The truth was computed by a different algorithm, in a different number system, with no rounding, and checked against a closed form derived from binomial coefficients. If the floating-point answer is wrong, it is wrong against something that cannot be wrong in the same direction.

The cases where this is not available

Most of them, which is why the rest of the site does not rely on it.

For a general matrix there is no closed-form inverse, and exact elimination on a matrix of arbitrary doubles produces rationals whose numerators double in length at every step — the cost grows exponentially and it becomes impractical by about n = 20 even in BigInt.

So the sitewide strategy is a hierarchy. Where a closed form exists, use it and require exact agreement. Where it does not, construct the problem backwards — pick the answer, build the right-hand side — which gives a known solution at the cost of a rounding in the construction, and is what the seeded test matrices do. And where neither is available, fall back on measuring the backward error, which certifies the algorithm and says nothing about the answer, and say so.

A small residual is not a small error is the essay about that last case, and the reason it comes before this one is that most real situations are that case. This essay is about the small number of situations where better is possible, and about the fact that the site’s central claim was established in one of them before being asserted anywhere else.

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 known answer, used as an axis. Eight precisions, one Hilbert system, and an error measured against the integers rather than against a better float. The prediction κ·u is marked; the measurement follows it over eleven orders of magnitude, which is a considerably stronger statement than “the error goes down as precision goes up”.

What is asserted

The exact inverse is required to be exact: H·H⁻¹ computed in rationals must be the integer identity in every entry, and every entry of the inverse must have denominator 1. The two exact routes must agree as fractions. The rational-to-double conversion must survive a 400-digit numerator and produce the right value. And the headline claim — that the solve is backward stable and that at least one component nonetheless has no correct digits — is checked as two separate assertions, either of which would fail the build alone.

That last pair is the shape of nearly every claim on this site: not “the answer is wrong”, which is easy, but “the answer is wrong and the algorithm is blameless”, which requires two measurements that could each independently refute the sentence.

The hierarchy, restated

Three levels of evidence about a numerical error, in decreasing order of strength, and it is worth knowing which one a given claim rests on.

The answer is known exactly. Available where a closed form exists — the Hilbert inverse, the Jordan block’s spectrum, Eckart–Young’s σₖ₊₁ — and conclusive when it is. Used wherever possible here.

The problem was constructed backwards. Pick the answer, build the right-hand side, and the truth is known up to one rounding in the construction. Available for any matrix, and what the seeded test matrices do; the caveat is that the construction’s rounding must be smaller than the effect being measured, which on an ill-conditioned matrix is not automatic.

Only the backward error is available. The general case. It certifies the algorithm completely and says nothing about the answer without a condition estimate.

Most published demonstrations sit at a fourth level below all three — comparison against a more careful computation — and the point of this essay is that the first level is available more often than people expect, and that it removes the last place an argument can hide.

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. 5 A fourth kind of known answer, in passing. Wilkinson’s matrix has growth factor exactly 2ⁿ⁻¹ by construction, so the measurement can be checked as an equality rather than a bound — and the same check confirms the growth routine can see a large value at all.

The habit generalises past the Hilbert matrix, and it is the one thing to take from this essay if nothing else survives: before measuring an error, ask what it is being measured against. If the answer is “a more careful version of the same computation”, the number is a difference rather than an error, and it is a lower bound at best. See a small residual is not a small error for the case where even a small honest number establishes nothing, and the condition number is an amplifier for the factor that turns one into the other.

Why the inverse is integral, briefly

The closed form is easy to state and easy to disbelieve, so it is worth a paragraph on why the inverse of a matrix of unit fractions contains no fractions.

The Hilbert matrix is the Gram matrix of the monomials on [0,1]: its (i,j) entry is the integral of xⁱ⁺ʲ, which is 1/(i+j+1). Inverting a Gram matrix amounts to expressing the monomials in terms of an orthogonal family for the same inner product — here the shifted Legendre polynomials, whose coefficients are integers. The determinant of Hₙ is the reciprocal of an integer, so Cramer’s rule puts an integer denominator underneath an integer cofactor and the fractions cancel exactly.

That the entries are integers is not a curiosity of the formula; it is what makes the check meaningful. Requiring every entry of H·H⁻¹ to be the integer 1 or the integer 0 is a much stronger demand than requiring it to be within 10⁻¹⁵, and it is available only because both factors are exact rationals. An implementation with a subtle error in a binomial coefficient would produce a nearly correct inverse and would fail the exact identity immediately.

It also explains the size of the numbers. The determinant of H₁₃ is about 10⁻¹⁶⁵, so its reciprocal — which is roughly the scale of the inverse’s entries — is enormous, and the largest entry is 1.07·10¹⁷. That number is the ill-conditioning made concrete: a matrix whose entries are all at most 1 has an inverse whose entries exceed 10¹⁷.

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. [ ε 1 ; 1 1 ] x = [ 1 ; 2 ], exact answer (1.000000, 1.000000)with partial pivoting1101U after elimination1.0000001.000000computed xbackward error 0forward error 0without10⁻¹⁷10-1·10¹⁷U after elimination0.0000001.000000computed xbackward error 0.25forward error 0.71no error is raisedgrowth 10¹⁷ The swap that is not optional Run elimination without a row interchange on a matrix that needs one and nothing announces a failure. There is no division by zero, no warning, and an answer of the right shape. It is simply wrong, and how wrong depends on a number you did not look at.