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

The roots are not the coefficients

A polynomial whose roots are the integers one to twenty, expanded exactly, handed to the routine every library uses. The computed roots are wrong in the third digit, the computation is backward stable for the matrix it factorised, and above degree eighteen the coefficients are not double-precision numbers at all.
34 min read 29 figures Whose fault is itMeasured, not assertedWhere the algebra stops being the arithmetic

Worth reading first: The condition number is an amplifier · The exact answer to a nearby problem · A matrix that depends on its own eigenvalue.

A scalar polynomial is a matrix polynomial of size one. Its companion matrix is the linearisation, the eigenvalues of that matrix are the roots, and everything the polynomial field says about linearising applies with n = 1.

Which means the field’s central claim can be tested at a size anybody can check by hand — and the test is the oldest example in this subject, sixty years old and still the sharpest.

The roots of ∏(x − i) up to degree 20, computed from the coefficients that define themA polynomial whose roots are the integers 1 to n, expanded exactly in BigInt and then handed to the route every library takes: the eigenvalues of the companion matrix. The measured relative error runs 1.711·10⁻¹³ at degree 6 to 0.007614 at degree 20 — the third digit of an answer whose truth is an integer. The other curve is the root condition number times the unit roundoff, which is a prediction made from the coefficients alone, and it tracks the measurement to within two orders at every degree. Above degree 19 the coefficients are integers larger than 2⁵³, so the polynomial handed to the routine is not the polynomial that was written down — the problem is damaged before any algorithm runs.57911131517192110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²degreerelative error of the worst rootmeasuredκ × ua polynomial is its coefficientsworst root, measured0.0076predicted, κ × u0.0023root condition number10¹³largest coefficient1.4·10¹⁹the roots are integersand the coefficients are not the roots
Fig. 1 The measured error in the roots of ∏(x − i) against degree, with the root condition number times the unit roundoff beside it. The prediction is made from the coefficients alone.

The polynomial

Take ∏(x − i) for i = 1 to n and expand it. The roots are the integers 1 to n, which is as known as an answer gets — no closed form to evaluate, no reference computation, just twenty integers.

The coefficients are integers too, and enormous: at n = 20 the constant term is 20! = 2.43·10¹⁸ and the coefficient of x¹⁹ is −210. Expanded exactly in BigInt, as they are here, so the only approximation anywhere in the demonstration is writing them down as doubles.

Then hand them to the route every numerical library takes: build the companion matrix, factorise it, read off the eigenvalues.

What comes back

degree worst relative error κ × u
6 1.7·10⁻¹³ 2.2·10⁻¹³
10 1.1·10⁻¹⁰ 1.4·10⁻¹⁰
14 7.2·10⁻⁷ 9.5·10⁻⁸
18 7.2·10⁻⁴ 7.9·10⁻⁵
20 7.6·10⁻³ 2.3·10⁻³

At degree twenty a root whose true value is an integer comes back wrong in the third significant figure. Not the fifteenth, not the eighth — the third.

The second column is the prediction, and it is made from the coefficients without reference to the computation. For a simple root r, perturbing coefficient a_k by a relative ε moves r by about ε|a_k||r|^k/|p′®|, so the relative sensitivity of r is the largest of those over k, divided by |r|. It reaches 1.0·10¹³ at degree twenty, at the root fourteen — not at the largest root and not at the smallest, but in the middle, where |r|^k is large and the derivative of the product is not.

Measured against predicted agrees to within a factor of nine at every degree and a factor of three at the last one. The error is the condition number times the unit roundoff, which is the condition number is an amplifier doing exactly what it says, on a problem where every other quantity is exactly known.

And the computation did nothing wrong

The temptation is to blame the eigenvalue routine, and the accounting refuses it.

The companion matrix at degree twenty has a norm of 2.4·10¹⁸, dominated by its first row. A real Schur factorisation of it is backward stable in the ordinary sense: the computed eigenvalues are exact for a companion matrix a relative rounding away. That theorem holds here, and it is the same one the form a real matrix can reach rests on.

But a perturbation of the matrix is not a perturbation of the polynomial. Change a single entry of the companion matrix by one unit in the last place and what results is the companion matrix of a polynomial whose coefficient has moved by one ulp — and one ulp of 2.43·10¹⁸ is about 500. The roots of a polynomial whose constant term has moved by 500 are not the integers.

So the answer is exact for a nearby polynomial, and a nearby polynomial has quite different roots. Which is the exact answer to a nearby problem with the second half of the sentence doing the damage, and it is the same structure as a backward-stable answer to a problem nobody asked one field over — there the mismatch is between a matrix and three coefficient matrices, and here between a matrix and a list of coefficients.

The coefficients are not binary64 numbers

There is a second failure underneath the first, and it arrives at a specific degree.

An integer is exactly representable in binary64 up to 2⁵³ ≈ 9.0·10¹⁵. The largest coefficient of the Wilkinson polynomial passes that at degree nineteen: 18! is 6.4·10¹⁵ and 19! is 1.2·10¹⁷.

So above degree eighteen the polynomial handed to the routine is not the polynomial that was written down. It has been rounded before any algorithm ran, by a relative 10⁻¹⁶ in a coefficient whose condition number is 10¹³, which accounts for three of the digits lost at degree twenty on its own.

That is worth separating from the conditioning argument, because it is a different kind of failure and has a different repair. Conditioning is a property of the question; representability is a property of the format, and it is the same boundary the units that overflow before the answer does measures for a matrix polynomial. Here it arrives at 2⁵³ rather than at 10³⁰⁸ because the quantity that has to be exact is an integer rather than merely finite.

What the site’s own machinery says about it

This site has an exact-arithmetic library, and it is worth saying plainly what it can and cannot do here.

exact.js expands the polynomial in BigInt, so the coefficients are exact at any degree — that is how the reference above is built. It can also evaluate the polynomial exactly at any rational point, so a root can be bracketed exactly: a sign change between two rationals is a certificate, and bisection on exact evaluations converges to any accuracy wanted with no rounding anywhere.

What it cannot do is make the double-precision route accurate, because the double-precision route’s problem is that the data it is given is already wrong. And what nobody can do is repair the conditioning: a polynomial whose roots move by 10⁻³ when a coefficient moves in its last bit has roots that are not determined by its coefficients to more than three digits, and no algorithm recovers information that the data does not contain.

That is the whole of what rank is a decision says about a different quantity, and it applies here word for word: the question “what are the roots of this polynomial, given these coefficients as doubles” has an answer good to three digits, and any routine claiming more is claiming something about a polynomial it was not given.

The perturbation that does the damage, written out

Wilkinson’s own account is worth reproducing because the arithmetic is small enough to check and because it makes the size of the effect concrete rather than asymptotic.

Perturb the coefficient of x¹⁹ — which is −210, an integer any format holds — by 2⁻²³, a relative change of about 10⁻⁹. The root at 20 moves by roughly

(2⁻²³ · 20¹⁹) / |p′(20)|,

and 20¹⁹ is 5.2·10²⁴ while p′(20) = 19! = 1.2·10¹⁷. The quotient is about 4·10⁻¹, so a change in the ninth digit of one coefficient moves a root by nearly half a unit. In Wilkinson’s original the pair at 16 and 17 leave the real axis altogether and become a complex conjugate pair — a polynomial with twenty distinct real roots, perturbed in one coefficient’s ninth digit, acquiring complex roots.

The measurement here is at the rounding level rather than at 2⁻²³, so the effect is smaller; what it is not is different in kind. And the imaginary parts are worth watching for the same reason: the computed spectrum of the companion matrix at degree twenty comes back with imaginary parts that are not zero, which is not merely inaccuracy but a change in the kind of answer — a real polynomial with real roots returning complex ones.

Where the site has met this shape before

Three times, and collecting them is the point of putting this essay in the error field rather than in the polynomial one.

The number that decides nothing is about a determinant: a scalar that is unrepresentable long before it is uninformative, and wrong sixteen sizes before it is absent. Here the coefficients are that scalar’s relatives — 20! is a determinant-sized number sitting in a coefficient — and they fail for the same arithmetic reason.

Two condition numbers of one matrix is about a quantity that is enormous in one norm and modest in another, so that the choice of measure decides whether a problem looks hard. The root condition number is the same kind of object: it is large because the monomial basis makes it large.

And the units the matrix is measured in is the repair that sometimes works — a scaling of the variable, x → αx, changes every coefficient and can improve the conditioning substantially. It is worth trying and it does not save this polynomial: the roots span a factor of twenty and no single α balances them.

Why anybody still does it this way

Given all that, it is fair to ask why the companion-matrix route is what every library ships.

Because the alternative is worse, and because the case above is rarer than it looks. A polynomial that arrives as a list of coefficients is unusual outside of textbooks — most polynomials in practice arrive as a characteristic polynomial nobody forms, a product of factors that can be kept factored, or an interpolant in a basis that is not the monomial one. The monomial basis is the bad one: its conditioning degrades exponentially with degree for roots away from the unit circle, and the same polynomial expressed in a Chebyshev basis, with a colleague matrix rather than a companion one, is far better behaved.

Which is the general form of the lesson, and it is the same one the polynomial field spends four essays on: the trouble is in the representation rather than in the algorithm. Change the units and a quadratic eigenvalue problem loses ten orders; change the basis and a scalar polynomial loses ten digits. In both cases the object is unchanged, every routine involved is backward stable, and what moved was the description.

What to do instead, priced

Three alternatives, and each is a real answer for a different situation.

Keep the factors. If the polynomial arrived as a product — a characteristic polynomial of a matrix, a product of linear factors from a model — never expand it. The roots are already there, and expanding is the step that destroys them. This is the same advice the inverse that is never formed gives about a matrix: the object that is easy to write down is not the object to compute with.

Change the basis. A polynomial expressed in Chebyshev coefficients on an interval containing its roots, with the corresponding colleague matrix, has a conditioning that does not degrade exponentially with degree. The routine is the same — build a matrix, factorise it — and the representation is the thing that changed.

Or work exactly. For integer or rational coefficients, exact evaluation and bisection give any accuracy wanted, at the cost of BigInt arithmetic whose numbers grow with the degree. On this polynomial that is entirely practical up to degree fifty, and it is how the reference for the measurement above is produced.

What none of the three does is make the double-precision monomial route accurate. That route answers a question about a polynomial that is not quite the one asked about, and the gap between the two is the whole of the error — which is, once more, the site’s oldest sentence with the emphasis on the last three words.

Why this belongs in the error field

The polynomial field would have taken this essay, and it is filed here instead for a reason worth making explicit.

Every quantity in it is one of the error field’s own. There is a condition number, computed from the data and not from the run. There is a backward error, at the rounding level, belonging to an object that is not the one posed. There is a forward error that is the product of the two to within a factor of nine at every degree. And there is an exactly known answer, which is what makes all three measurable rather than estimable.

The identity that runs through this site —

forward error ⪅ condition number × backward error

— is usually demonstrated on a linear system where all three terms are about the same matrix. Here the three terms are about different objects: the condition number is the polynomial’s, the backward error is the companion matrix’s, and the forward error is the roots’. Getting the attribution right requires noticing that, and the whole failure of intuition this essay records is the failure to notice it.

That is the error field’s subject rather than the polynomial field’s. What the polynomial field adds is that the same three-way mismatch occurs at every degree and every size, and that the repair — change the representation — is the same one.

What the measurement needed in order to exist

A closing note on method, since this essay’s numbers rest on three separate exactnesses and losing any of them would have produced a plausible wrong answer.

The roots are known. They are the integers, by construction, so the forward error is a difference from a known value rather than from a better computation.

The coefficients are exact. They are expanded in BigInt, so the polynomial being discussed is the polynomial that was meant, and the rounding that occurs when they are cast to doubles is identified as a separate, measurable effect rather than absorbed into the total.

And the prediction is independent. The root condition number is computed from the coefficients and the known roots, without reference to the computed answer, so “measured equals predicted” is a statement about two quantities rather than a restatement of one.

Take away the first and the error is unmeasurable. Take away the second and the polynomial drifts before the algorithm runs, with no way to say by how much. Take away the third and the agreement between measurement and prediction becomes circular. That is a lot of scaffolding for five numbers, and it is the reason those five numbers can be believed.

At other settings

The roots of ∏(x − i) up to degree 10, computed from the coefficients that define themA polynomial whose roots are the integers 1 to n, expanded exactly in BigInt and then handed to the route every library takes: the eigenvalues of the companion matrix. The measured relative error runs 1.711·10⁻¹³ at degree 6 to 1.127·10⁻¹⁰ at degree 10 — the third digit of an answer whose truth is an integer. The other curve is the root condition number times the unit roundoff, which is a prediction made from the coefficients alone, and it tracks the measurement to within two orders at every degree. Every coefficient here is still exactly representable, which is what makes the failure a conditioning one rather than a storage one.5791110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²degreerelative error of the worst rootmeasuredκ × ua polynomial is its coefficientsworst root, measured1.1·10⁻¹⁰predicted, κ × u1.4·10⁻¹⁰root condition number6.1·10⁵largest coefficient1.3·10⁷the roots are integersand the coefficients are not the roots
Fig. 2 Degree ten, where the coefficients are still small integers and the error is 10⁻¹⁰.
The roots of ∏(x − i) up to degree 14, computed from the coefficients that define themA polynomial whose roots are the integers 1 to n, expanded exactly in BigInt and then handed to the route every library takes: the eigenvalues of the companion matrix. The measured relative error runs 1.711·10⁻¹³ at degree 6 to 7.19·10⁻⁷ at degree 14 — the third digit of an answer whose truth is an integer. The other curve is the root condition number times the unit roundoff, which is a prediction made from the coefficients alone, and it tracks the measurement to within two orders at every degree. Every coefficient here is still exactly representable, which is what makes the failure a conditioning one rather than a storage one.57911131510⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²degreerelative error of the worst rootmeasuredκ × ua polynomial is its coefficientsworst root, measured7.2·10⁻⁷predicted, κ × u9.5·10⁻⁸root condition number4.3·10⁸largest coefficient3.9·10¹¹the roots are integersand the coefficients are not the roots
Fig. 3 Fourteen, where it passes 10⁻⁶.
The roots of ∏(x − i) up to degree 22, computed from the coefficients that define themA polynomial whose roots are the integers 1 to n, expanded exactly in BigInt and then handed to the route every library takes: the eigenvalues of the companion matrix. The measured relative error runs 1.711·10⁻¹³ at degree 6 to 0.3333 at degree 22 — the third digit of an answer whose truth is an integer. The other curve is the root condition number times the unit roundoff, which is a prediction made from the coefficients alone, and it tracks the measurement to within two orders at every degree. Above degree 19 the coefficients are integers larger than 2⁵³, so the polynomial handed to the routine is not the polynomial that was written down — the problem is damaged before any algorithm runs.5791113151719212310⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²degreerelative error of the worst rootmeasuredκ × ua polynomial is its coefficientsworst root, measured0.33predicted, κ × u0.066root condition number3·10¹⁴largest coefficient6.8·10²¹the roots are integersand the coefficients are not the roots
Fig. 4 And twenty-two, past the point where the coefficients can be stored.
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. 5 The amplifier itself, in the essay that introduced it.
The backward error of a hierarchical solve, against the error of the representation it usedSolve A_H x = b exactly and the residual against the matrix that was wanted is b − Ax = (A_H − A)x, so ‖b − Ax‖ ⁄ ‖A‖‖x‖ cannot be anything but the representation's own error. Measured across six accuracies spanning ten decades, the two track at a slope of 0.994 and sit a constant 6.7× apart, which is the difference between a norm of a matrix and a norm of that matrix applied to one vector. The reading is the one this field is built on: the accuracy knob is not an accuracy, it is a backward error chosen in advance. Everywhere else on this site a backward error is something an algorithm produced and somebody then measured; here it is a line in the program, and its size is known before the solve starts.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.99representation at 10⁻⁸1.2·10⁻⁹backward error there1.3·10⁻¹⁰representation at 10⁻¹²1.2·10⁻¹³backward error there2.5·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem
Fig. 6 The backward error, which is what the routine here is entitled to be judged on.
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. 7 A small residual and a large error, on a linear system.
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. 8 The per-eigenvalue sensitivity of a matrix, which the root condition number is the scalar case of.
One quadratic eigenvalue problem in 9 systems of units: what the solver reports and what the answer is worthAn overdamped chain of 8 masses, with λ replaced by γμ so that the coefficients become (γ²M, γC, K). That substitution is exact in both directions and divides the spectrum by γ exactly, so the closed form is still available and every error here is measured against it. The backward error of the eigenpair for the LINEARISED MATRIX — the residual a solver's own error analysis is about — is 6.76·10⁻¹⁶ at γ = 1 and 7.59·10⁻¹³ at γ = 108 — it moves by a factor of 1928 while the other two move by 1.68·10¹⁰. The backward error for the QUADRATIC, which is what the person who posed the problem is entitled to, grows by 8.76·10¹⁰ across the same sweep, and the forward error follows it: 7.49·10⁻¹⁴ to 0.00126. Nothing went wrong with the solver at any stop.0246810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ γ, the change of unitsrelative errorforward errorη, the quadraticη, the linearisationagainst a closed formη(linearisation), worst7.6·10⁻¹³η(quadratic), worst1.2·10⁻⁴forward error, worst0.0013coefficient spread4.2·10¹⁵the solver is right at every stopabout a problem nobody asked
Fig. 9 The same argument for a matrix polynomial, where the substitution rather than the basis is at fault.
Where the missing eigenvalues went: an 6 × 6 quadratic with 1 of its masses removedA chain of 6 masses with the first 1 of them set to zero, so M is singular of rank 5. det Q(λ) is a polynomial of degree at most 12; interpolated exactly in BigInt rationals at 13 nodes it has degree 11, so 1 of the 12 eigenvalues are at infinity — the same object a descriptor pencil has, arriving here because a degree of freedom with no inertia is an algebraic constraint. The float route counts the singular values of M judged to be zero and reaches 1, backed by a gap of ∞ between consecutive singular values. One integer, two routes, and only the second of them is a decision.finite eigenvalues (degree of det Q)11at infinity (2n − degree)1at infinity, by the rank of M12n, if M were nonsingular12a degree, not a decisiondegree of det Q11at infinity1by the rank of M1singular-value gap∞the count is a degreeand the other route is a judgement
Fig. 10 A characteristic polynomial computed exactly, which is how this essay’s reference is built.
log₁₀|det H_n| by three routes, to n = 30The exact rational determinant of the Hilbert matrix falls to 10^-518 at n = 30. The accumulated logarithm of the pivots follows it to n = 13 and then departs, reaching 10^-352 — wrong by 166 orders of magnitude and still an ordinary-looking number. The product of the pivots is a separate curve that ends at n = 29, where it underflows to exactly zero.26101418222630-550-450-350-250-150-50nlog₁₀|det H_n|exact rationalΣ log|u_ii|Π u_iithree routes, one theoremexact at the largest n-518accumulated logarithm-352decades of disagreement166smallest pivot at that n1.2·10⁻¹⁷every pivot is a normal numberat every size on this axis
Fig. 11 The determinant, which is the same object’s value rather than its roots.
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. 12 An answer known in closed form, and what it is for.
Wilkinson's 12×12 matrix and its upper triangular factorThe matrix on the left has ones on the diagonal, minus ones below it and a column of ones at the right. On the right, its U factor, whose last column doubles down the rows to sixty-four.1··········1-11·········1-1-11········1-1-1-11·······1-1-1-1-11······1-1-1-1-1-11·····1-1-1-1-1-1-11····1-1-1-1-1-1-1-11···1-1-1-1-1-1-1-1-11··1-1-1-1-1-1-1-1-1-11·1-1-1-1-1-1-1-1-1-1-111-1-1-1-1-1-1-1-1-1-1-11A1··········1·1·········2··1········4···1·······8····1······16·····1·····32······1····64·······1···128········1··256·········1·512··········11024···········2048U‖PA − LU‖/‖A‖0growth factor2048the 2ⁿ⁻¹ bound2048Row interchanges performed: 0.Partial pivoting had nothing to choose.every entry is 0, 1 or −1the bound is attained here
Fig. 13 Wilkinson’s other matrix, and the growth it produces under elimination.
Where 24 perturbations of size 10⁻⁸ put the eigenvalues of two 6×6 matrices with the same spectrumTwo clouds of eigenvalues in the complex plane. Both matrices have every eigenvalue at 0.8 exactly; both were perturbed by the same 24 random matrices of norm 10⁻⁸. The normal matrix's eigenvalues stay within 8.86·10⁻⁹ of where they were — the size of the perturbation — and the bidiagonal's spread out to 0.0701, a factor of 7.9·10⁶ further on the same data.100.1750.35real partimaginary part0.8bidiagonalnormalone spectrum, two matricesthe perturbation10⁻⁸normal, furthest moved8.9·10⁻⁹bidiagonal, furthest0.07ratio7.9·10⁶the two matrices have identical eigenvaluesand one of them says so under perturbation
Fig. 14 What a perturbation does to a spectrum when the matrix is not normal.
The singular values of a sum of two rank-4 blocks, and the 4 a truncation has to discardA rank-4 block times a vector is a rank-4 block times a vector. A rank-4 block times a rank-4 block is a rank-4 block. A rank-4 block PLUS a rank-4 block is a rank-8 block, exactly, and the 8 bars here are why: the sum of two 4-dimensional spaces is generally 8-dimensional, and none of the 8 singular values is small. Truncating back to 4 costs 61.9 per cent of the block. Below the 8th the values are the unit roundoff, which is the check that the doubling is exact rather than approximate. Every product, every factorisation and every Schur complement inside this format is a chain of these, and there is nothing else to do: without the truncation the ranks double at every level and the format is dense by the bottom.σ ⁄ σ₁ of the sum, 64 × 64σ1, kept1σ2, kept0.922σ3, kept0.812σ4, kept0.785σ5, discarded0.778σ6, discarded0.758σ7, discarded0.658σ8, discarded0.571σ92.63·10⁻¹⁶σ102.3·10⁻¹⁶the operation that is not closedrank of each term4rank of the sum8truncated back to4cost of the truncation0.62the best there is0.62two planesmake a four-space
Fig. 15 The integers a double can hold exactly, which is where the second failure comes from.
Where binary32 stops being able to hold the linearisation, and where the scaled problem does notForming the first companion linearisation needs γ²M to be representable. The upper row is whether every entry of (γ²M, γC, K) is finite in binary32, at twelve changes of units; the lower row is the same question after Fan–Lin–Van Dooren scaling. The boundary is √(largest finite value ÷ largest entry of M) = 1.8447·10¹⁹, computed from the format's parameters and nothing measured — and the sweep agrees with it at every stop: the last change of units that survives is 10¹⁹ and the first that does not is 10²⁰. The scaled row never fails, because δγ²M has norm about one whatever γ was. In a narrow format the scaling is not an accuracy device: it is what makes the problem exist.02468101201change of units, by exponent10⁰10¹10²10³10⁴10⁶10¹⁰10¹⁶10¹⁹10²⁰10⁴⁰10¹⁵⁰as writtenafter scalinga range questionlargest finite value3.4·10³⁸predicted boundary γ1.8·10¹⁹last γ that forms10¹⁹stops where scaling fails0not a poor answerno answer at all
Fig. 16 And the same boundary for a matrix polynomial’s coefficients.
The roots of ∏(x − i) up to degree 18, computed from the coefficients that define themA polynomial whose roots are the integers 1 to n, expanded exactly in BigInt and then handed to the route every library takes: the eigenvalues of the companion matrix. The measured relative error runs 1.711·10⁻¹³ at degree 6 to 7.174·10⁻⁴ at degree 18 — the third digit of an answer whose truth is an integer. The other curve is the root condition number times the unit roundoff, which is a prediction made from the coefficients alone, and it tracks the measurement to within two orders at every degree. Every coefficient here is still exactly representable, which is what makes the failure a conditioning one rather than a storage one.579111315171910⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²degreerelative error of the worst rootmeasuredκ × ua polynomial is its coefficientsworst root, measured7.2·10⁻⁴predicted, κ × u7.9·10⁻⁵root condition number3.5·10¹¹largest coefficient3.4·10¹⁶the roots are integersand the coefficients are not the roots
Fig. 17 Degree eighteen, the last one whose coefficients are exactly representable in binary64.
The roots of ∏(x − i) up to degree 20, computed from the coefficients that define themA polynomial whose roots are the integers 1 to n, expanded exactly in BigInt and then handed to the route every library takes: the eigenvalues of the companion matrix. The measured relative error runs 1.711·10⁻¹³ at degree 6 to 0.007614 at degree 20 — the third digit of an answer whose truth is an integer. The other curve is the root condition number times the unit roundoff, which is a prediction made from the coefficients alone, and it tracks the measurement to within two orders at every degree. Above degree 19 the coefficients are integers larger than 2⁵³, so the polynomial handed to the routine is not the polynomial that was written down — the problem is damaged before any algorithm runs.57911131517192110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²degreerelative error of the worst rootmeasuredκ × ua polynomial is its coefficientsworst root, measured0.0076predicted, κ × u0.0023root condition number10¹³largest coefficient1.4·10¹⁹the roots are integersand the coefficients are not the roots
Fig. 18 And twenty, where the largest is 2.4·10¹⁸ and is not.
Fitting the same degree-11 polynomial in two basesOn the left, the data and two fitted curves that lie on top of each other. On the right, the condition numbers of the two design matrices, ten orders of magnitude apart.00.10.20.30.40.50.60.70.80.910.20.40.60.81xyboth fits, drawn on top of one anothercondition number of the design matrixmonomial1.2·10⁸Chebyshev2.5largest fitted coefficientmonomial113Chebyshev0.51rms residual: 3.9·10⁻⁶ and 2.9·10⁻⁷ — the data is fitted either way.30 points, degree 11, single precisionthe basis is part of the problem
Fig. 19 The conditioning of a monomial basis against degree, which is what the whole failure is about.
Fitting the same degree-10 polynomial in two basesOn the left, the data and two fitted curves that lie on top of each other. On the right, the condition numbers of the two design matrices, ten orders of magnitude apart.00.10.20.30.40.50.60.70.80.910.20.40.60.81xyboth fits, drawn on top of one anothercondition number of the design matrixmonomial2.1·10⁷Chebyshev2.3largest fitted coefficientmonomial176Chebyshev0.51rms residual: 4.2·10⁻⁶ and 3.1·10⁻⁷ — the data is fitted either way.30 points, degree 10, single precisionthe basis is part of the problem
Fig. 20 And one degree lower.
The 16 eigenvalues of an 8 × 8 quadratic eigenvalue problem, computed and in closed formλ²M + λC + K for a chain of 8 masses with C = 0.3M + 0.1K. The crosses are the closed form — one scalar quadratic per eigenvalue of K, whose roots are known exactly — and the discs are the eigenvalues a real Schur factorisation returns from the 16 × 16 first companion linearisation. There are 16 of them for a matrix with 8 rows, of which 16 are complex and arrive in conjugate pairs, so the eigenvectors cannot be independent: 9 vectors in 8 dimensions never are. The worst disagreement between the two routes is 2.12·10⁻¹⁵, and the routes share nothing but the three coefficient matrices.00.0939569-2-1.24498-0.4899610.2650581.020081.7751real partimaginary parttwo routes, one spectrumeigenvalues16rows8complex16against the closed form2.1·10⁻¹⁵n rows and 2n eigenvaluesso the eigenvectors are not a basis
Fig. 21 The matrix case of the same object, where n is eight rather than one.
One quadratic eigenvalue problem in 7 systems of units: what the solver reports and what the answer is worthAn overdamped chain of 8 masses, with λ replaced by γμ so that the coefficients become (γ²M, γC, K). That substitution is exact in both directions and divides the spectrum by γ exactly, so the closed form is still available and every error here is measured against it. The backward error of the eigenpair for the LINEARISED MATRIX — the residual a solver's own error analysis is about — is 6.76·10⁻¹⁶ at γ = 1 and 1.26·10⁻¹³ at γ = 106 — it moves by a factor of 320 while the other two move by 7.31·10⁷. The backward error for the QUADRATIC, which is what the person who posed the problem is entitled to, grows by 1.86·10⁸ across the same sweep, and the forward error follows it: 7.49·10⁻¹⁴ to 5.47·10⁻⁶. Nothing went wrong with the solver at any stop.012345610⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ γ, the change of unitsrelative errorforward errorη, the quadraticη, the linearisationagainst a closed formη(linearisation), worst1.3·10⁻¹³η(quadratic), worst2.6·10⁻⁷forward error, worst5.5·10⁻⁶coefficient spread4.2·10¹¹the solver is right at every stopabout a problem nobody asked
Fig. 22 And the matrix case of the same argument about backward stability.
Relative error of two algebraically identical expressions for (1 − cos x)/x², in binary64A log–log plot of relative error against x. The expression written as it reads loses accuracy as x falls and is eventually wrong in every digit; the rearranged form stays at rounding level.10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at allbinary64 throughoutone function, two spellings · zero below 1.5·10⁻⁸
Fig. 23 Cancellation in a sum, which is what expanding a product of factors does to them.
How far |r_nn| sits above σ_min on Kahan's matrix, against the size and the parameter3 curves of |r_nn| ÷ σ_min against n, one per Kahan parameter. Every curve rises without turning over, reaching 7·10⁴ at n = 30, c = 0.5. Column pivoting makes no interchange at any point on any of them, so the failure is not a poor choice — there is nothing to choose.813182328110¹10²10³10⁴10⁵10⁶size of the matrix|r_nn| ÷ σ_minthe two agreec = 0.2c = 0.35c = 0.5no ceilingratio at n = 1021ratio at n = 307·10⁴interchanges, anywhere0the greedy rule never had a choiceand the gap grows with every row
Fig. 24 Compensated summation, and the digits it buys back.
Backward error of Cramer's rule and of elimination on 11974 2×2 systems at 24 bitsThe share of systems whose normwise backward error exceeds a given multiple of the unit roundoff. Elimination never exceeds 1.2u. Cramer's rule reaches 395u on the same systems, on a matrix whose condition number is 3.39·10⁶ and which elimination solved to 0.3u. The control curve, over uniformly random 2×2 systems, stops at 4.9u — so the failure belongs to the near-parallel family and not to a choice of scale.10⁻²10⁻¹110¹10²10³10⁴10⁻⁴10⁻³10⁻²10⁻¹1backward error, in units of ushare of systems above ituCramereliminationCramer, control24-bit arithmeticworst Cramer, in u395worst elimination, in u1.2control, worst Cramer4.9κ of the worst system3.4·10⁶one derivation, two computationsand only one of them is stable
Fig. 25 A formula that is correct and unusable, from the elimination field.
Where conjugate gradients certifies that a positive definite matrix is indefiniteEvery matrix in this grid is 30×30 and positive definite by construction — its spectrum is κ^(−i/(n−1)), so the smallest eigenvalue is 1/κ and none of them is negative. A filled cell is a run in which pᵀAp came out non-positive and the iteration produced a direction it would report as a proof of indefiniteness. 14 of the 72 runs did. The region is a staircase whose top edge is at twelve significand bits and whose left edge, at eight, is at κ = 10⁵. Nothing rounded incorrectly anywhere: every comparison was performed exactly as written, on a number that was computed as accurately as the format allows.significand bits10³10⁴10⁵10⁶10⁷10⁸10⁹10¹⁰κ(A)87892634456329625911109435761199127514162024every matrix positive definiteruns producing a false certificate14of runs in total72never above, in significand bits12first κ at eight bits10⁵the comparison was correctand what it proved was not true
Fig. 26 A certificate that is not one, which is what a small residual is here.
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. 27 And an answer known exactly, which is the only way any of this is measurable.
The roots of ∏(x − i) up to degree 14, computed from the coefficients that define themA polynomial whose roots are the integers 1 to n, expanded exactly in BigInt and then handed to the route every library takes: the eigenvalues of the companion matrix. The measured relative error runs 1.711·10⁻¹³ at degree 6 to 7.19·10⁻⁷ at degree 14 — the third digit of an answer whose truth is an integer. The other curve is the root condition number times the unit roundoff, which is a prediction made from the coefficients alone, and it tracks the measurement to within two orders at every degree. Every coefficient here is still exactly representable, which is what makes the failure a conditioning one rather than a storage one.57911131510⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²degreerelative error of the worst rootmeasuredκ × ua polynomial is its coefficientsworst root, measured7.2·10⁻⁷predicted, κ × u9.5·10⁻⁸root condition number4.3·10⁸largest coefficient3.9·10¹¹the roots are integersand the coefficients are not the roots
Fig. 28 Degree fourteen, where the error passes 10⁻⁶.
Fitting the same degree-10 polynomial in two basesOn the left, the data and two fitted curves that lie on top of each other. On the right, the condition numbers of the two design matrices, ten orders of magnitude apart.00.10.20.30.40.50.60.70.80.910.20.40.60.81xyboth fits, drawn on top of one anothercondition number of the design matrixmonomial2.1·10⁷Chebyshev2.3largest fitted coefficientmonomial176Chebyshev0.51rms residual: 4.2·10⁻⁶ and 3.1·10⁻⁷ — the data is fitted either way.30 points, degree 10, single precisionthe basis is part of the problem
Fig. 29 And the basis whose conditioning decides all of it.

Where this goes next

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. 00.0939569-2-1.24498-0.4899610.2650581.020081.7751real partimaginary parttwo routes, one spectrumeigenvalues16rows8complex16against the closed form2.1·10⁻¹⁵n rows and 2n eigenvaluesso the eigenvectors are not a basis A matrix that depends on its own eigenvalue A damped structure does not produce Ax = λx. It produces (λ²M + λC + K)x = 0, where the matrix whose null vector is wanted is a function of the number being solved for — so there is nothing to factorise, an n × n problem has 2n answers, and the eigenvectors cannot be a basis. 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.

Reads more easily once this is understood

Essays that name this one as worth reading first.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Backward errorCharacteristic polynomialCompanion formCondition numberExact ground truthMatrix polynomialRepresentable numbersWilkinson polynomial