Field

Two errors, and whose fault they are

A good algorithm returns the exact answer to a nearby problem. So a wrong answer has two possible authors, and they can be measured separately: the backward error says how well the algorithm did, the condition number says how much the problem amplifies it, and only their product is what anyone sees.
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¹⁴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.

0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ∞(DA)cond(DA)Hilbert κ∞Hilbert condone system, two numbersκ∞ at no spread9.8κ∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it

The units the matrix is measured in

One linear system, written twice. The rows of the second are the rows of the first in different units, the solution is identical to the last bit, and the condition number has moved by eight orders of magnitude. One of those two numbers is a fact about the problem and the other is a fact about the notation.

10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ∞ · εcond(A,x) · εmeasuredboth bounds holdκ∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement

A condition number scaling cannot move

Skeel's componentwise condition number is invariant under any row scaling — exactly, before any norm is taken, because two diagonal factors cancel entry by entry. It is never larger than the normwise one and can be arbitrarily smaller, and the ratio between them is a diagnostic for which kind of ill-conditioning a matrix has.

the estimator maximises this quantity over the columns it visitscolumn 1 ‹visited›12column 2 ‹the answer›114column 311.4column 411.4column 511.4column 611.4column 711.4column 811.4column 911.4column 1011.4column 1111.4column 1211.4estimate 12.0a walk that stopped earlythe estimate returned12the true 1-norm114columns visited1products with the matrix5the walk's own stopping test firedand every column it could see was smaller

An estimate that can be fooled

Nobody computes a condition number, because forming an inverse costs more than the solve did. Every library estimates it instead, from four or five products with a factorisation already in hand. The estimate is exactly right on four random matrices out of five — and there is a matrix, three distinct entries wide, on which it returns a twentieth of the truth.

does this matrix look nearly singular?green: the test agrees with the truth · red: it does not · the bar under each number is its magnitude, over sixty-two decades|det A||det A|^(1/n)σₘᵢₙ1/κ = σₘᵢₙ/σₘₐₓ0.1·I at n = 40perfectly conditioned10⁻⁴⁰0.10.11κ = 10¹⁰, |det| = 1nearly singular1110·10⁻⁶10·10⁻¹¹Hilbert at n = 8nearly singular2.7·10⁻³³8.5·10⁻⁵1.1·10⁻¹⁰6.6·10⁻¹¹the two counterexamplesκ of the scaled identity1its determinant10⁻⁴⁰κ of the normalised matrix10¹⁰its determinant1det(cA) = cⁿ det(A)so a determinant carries the units n times over

The number that decides nothing

The determinant is the first scalar anybody attaches to a matrix and the last one worth consulting. A tenth of the identity has a determinant of 10⁻⁶⁰ and a condition number of exactly one. The Hilbert matrix's determinant stops being right at n = 13 and stops being a number at n = 29, and nothing in between reports either.

015304560759010512010⁻²10⁻¹110¹steprelative sizeleast error: 20discrepancy stop: 7errorresidualthe knob is an integerleast error, at step20error there0.14error at step 1206the residual falls at every stepthe error turns and keeps rising

The zero you are allowed to write

A deflation criterion sets a subdiagonal entry to zero because it is small. A drop tolerance discards an entry of a factor because it is small. A truncation discards a singular value because it is small. Three fields, three vocabularies, no shared arithmetic — and plotted as work saved against error accepted, one curve.

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

Deciding that a zero has arrived

The previous tolerances were offers — accept this much error, save this much work. A detection threshold is not an offer, because both directions are failures. One matrix here has three genuinely near-invariant subspaces, and the constant somebody typed decides which of them the recurrence stops at; at eight significand bits the same kind of constant produces a proof of something false.

10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³relative compression error, 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

An accuracy that is a backward error

Every backward error measured so far is something an algorithm produced and somebody then measured. This one is a line in the program. Solving with a compressed matrix gives a residual that is the compression's own error, at a slope of 1.000 over ten decades, so the knob that sets the storage sets the backward error directly.

110¹10²10⁻¹²10⁻⁸10⁻⁴110⁴ncondition number, and residual reachedmarks above: κ of the step · dashes: its closed formmiddle: a fit from a random startbelow: the same fit started at the answera decomposition that is ill-conditionedκ at n = 1283.3·10⁴its closed form3.3·10⁴cosine of the terms1from a random start0.012from the answer1.4·10⁻¹⁰the answer existsand cannot be found

A tensor that cannot be decomposed

Every member of a certain sequence is exactly a sum of two rank-one terms, and both terms are written down in closed form. A three-hundred-sweep fit from a random start does not find them, and stalls at the same one per cent however far the sequence goes — while a fit started at the answer loses digits exactly as 2n² says it should.

-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μcondition number, and relative errorκ₂, augmentedcomponentwise, condensedcomponentwise, augmentederror, augmentedone matrix, two numbersκ₂ at μ = 10⁻¹²3·10¹³componentwise, same matrix13their ratio2.3·10¹²measured error9.4·10⁻¹⁶every library prints the top lineand the error obeys the third

Two condition numbers of one matrix

κ₂ is a worst case over perturbations of a given norm, and a normwise perturbation may put its whole budget on the smallest entry. The componentwise number is a worst case over perturbations proportional to the entries, which is what a backward-stable factorisation actually makes. On one matrix they are 3·10¹³ and 13.3, and the error obeys the second.

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

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.

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

Three errors and one number

This site's identity has two factors and a division of blame between them. Two fields have now added a third party and a fourth, and only one of the four is a property of anything — the others are decisions, made before the arithmetic, reported by nothing.

the problem you posedA = H12b = A·(1, 2, …, 12)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 1.8·10⁻¹⁷the answer you wantedx = (1, 2, …, 12), exactlythe answer you gotx̂, wrong by 0.02 relativebackward error 1.8·10⁻¹⁷forward error 0.02κ = 1.8·10¹⁶κ · η = 0.33, and the measured forward error is 0.02.The algorithm is not at fault. The problem is.H12, LU with partial pivotingresidual and error differ

The fifth author

Four authors of a wrong answer have been named already and each is a statement about one computation. The fifth is not: it is what separates two computations that are both correct, it is a backward error of measurable size, and no residual, bound or condition number contains it.

10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²-1-0.500.51relative backward error acceptedfraction of the work not doneQR deflation criterionincomplete Cholesky droplow-rank truncationthree fields, no shared arithmeticslope, deflation criterion0.039slope, drop tolerance0.051slope, rank truncation0.23widest apart, as a ratio5.7a tolerance is an offerand the three offers are one offer

A tolerance is priced by the problem

Three tolerances from three fields sit on one pair of axes and agree to within a factor of 5.74. That factor is the ratio of the two curves that cannot move. Change the only problem in the comparison and the third curve's fitted slope swings from 0.188 to 0.040 while the printed spread does not shift by a digit.

each bar is a percentage — of the sample, or of the true condition numberexactly right86.0%inside 10%91.0%inside a factor of 291.0%worst in the sample, ×10035.7%the constructed matrix, ×1007.7%usually exactexact share0.86worst of the sample0.36the constructed matrix0.077a routine that is right most of the timeand never wrong in the safe direction

The tail a sample never reaches

Hager's estimator is exactly right on four random matrices in five, and that share is stable — between 80.5 and 87.5 per cent across nine sizes. The worst underestimate is not stable at all: it falls every time more matrices are drawn, from 0.746 at sixty to 0.377 at four hundred, and the matrix built to defeat the estimator sits five times below anything four hundred draws found.

LAPACK's walk, exact83.0% · 5.2 productsLAPACK's walk, worst × 10037.7% · 5.2 productsblock of one, exact84.0% · 4.4 productsblock of one, worst × 10037.7% · 4.4 productsblock of two, exact96.5% · 8.5 productsblock of two, worst × 10059.6% · 8.5 productsblock of four, exact100.0% · 16.8 productsblock of four, worst × 100100.0% · 16.8 productseach bar a percentage — of the sample, or of the truthtwo random vectors close most of the tail

Two columns see what one walk cannot

The condition estimator every library ships walks from the all-ones vector, and a matrix whose largest column cancels against that vector hides from it: at n = 24 it reports five per cent of the truth. The block estimator behind MATLAB's condest walks with two vectors, the second random. On the same matrix at three sizes it is exact on every one of twenty seeds. On four hundred random 8 × 8 matrices it is exact on 96.5 per cent where the single walk is exact on 83.0, and its worst case, 0.596, is reached in the first fifty draws and not lowered by the next 1,550. The single walk's worst was still falling at 1,600. Four vectors are exact on all 400.

largest eigenvalue erroras given, order 960.043balanced, order 964.4·10⁻⁴symmetrised, order 962.2·10⁻¹²016324864809611210⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹order nlargest |computed − exact| eigenvaluean error of one: the integers are no longer told apartas givenbalancedsymmetrisedopen dots: complex pairs returned for a real spectrumevery entry is an integer, stored exactly

Balanced is not symmetric

The Sylvester–Kac matrix is made of small integers, so a double holds it exactly, and its eigenvalues are the integers from −(n − 1) to n − 1 in steps of two. Every digit an eigensolver loses on it is therefore the solver's own, and it loses them at exactly the rate first-order perturbation theory predicts: the median error is half the prediction across 1,568 eigenvalues. Balancing, the preprocessing libraries apply for this kind of matrix, divides every condition number by about sixty and leaves their growth untouched, and at order 112 the unbalanced solver returns eighteen complex eigenvalues for a spectrum of integers.

each steering column half the hidden oneone random vector, 34 columns: share fooled0.77block of two, 34 columns: share fooled0.64block of four, 34 columns: share fooled0.390510152025303500.20.40.60.81steering columnsshare of seeds fooledone random vectorblock of twoblock of fourdashed: the single vector's share squared and to the fourthmore columns, more seeds, same factor

Columns that steer a random start

A condition estimator whose first vector is random cannot be built against, and concentrating the hidden column did not hide it: its own entries steer the first sign vector onto it. Give the steering to other columns — each equal on both rows of every pair the hidden column alternates over — and the random start is fooled after all. The prediction was a factor of one over the square root of the number of steering columns, on a share of seeds that falls as they lose to the hidden column, and the block of two fooled as often. The factor is whatever the builder chooses; the share, at a fixed factor of one half, rises with every column added, from 9 per cent with four to 77 with thirty-four. And the block of two is not fooled as often: it is fooled on the square of the single start's share, and the block of four on its fourth power — independent guesses, which is the one thing that still works.

exact share, points apartsize 8: random start less all-ones start0.025size 16: random start less all-ones start0.00570%80%90%100%share of matrices on which the estimate is exactsize 8size 16LAPACK's walk83.0% · 5.2 productsone walk from all ones84.0% · 4.4 productsone walk from random signs86.5% · 4.3 productsblock of two96.5% · 8.5 productsLAPACK's walk83.5% · 5.3 productsone walk from all ones83.5% · 4.4 productsone walk from random signs83.0% · 4.3 productsblock of two96.0% · 8.5 productsbars start at seventy per centthe first vector's sign pattern costs nothing

A first vector nobody can build against

The matrix built to fool a condition estimator is built against one vector, the all-ones vector its walk starts from, and the block estimator escaped it by adding a second, random one. Starting the single walk from random signs instead escapes it on every one of forty seeds at every size from 8 to 48, for the same 4.3 products, and loses nothing on random matrices — 86.5 per cent exact at size 8 against 84.0 from all ones, with tails that cross between sizes. The obvious way to build against a random start, a hidden column on few rows whose signs a random vector cancels half the time, fails on every seed: the hidden column writes itself into the walk's first product and turns the walk towards it. What the block of two's second vector buys is ten points of exact share, not the escape.

612182430364210⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1degree nworst relative root errormonomial, integersmonomial, Chebyshev pointsChebyshev, integersChebyshev, Chebyshev pointsred: monomial; blue and green: Chebyshevthe basis has to match the roots

A basis that has to know where the roots are

Wilkinson's polynomial loses its roots in the third figure at degree twenty when they are computed from its monomial coefficients, and the repair offered was to change the basis: in Chebyshev coefficients, with the colleague matrix, the conditioning 'does not degrade exponentially with degree'. Expanded exactly in that basis and rounded once, the integers' roots come back ten orders better at degree twenty — 4.4·10⁻¹³ against 7.6·10⁻³ — and still degrade exponentially, at a third of a decade a degree, to 2.8·10⁻⁴ at forty. Place the same number of roots at the interval's Chebyshev points instead and they stay at 10⁻¹⁴ at every degree, while the monomial route loses them as fast as the integers. The basis does not repair the problem; it has to match where the roots are.

1015202530354010⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³degree n, roots at 1 … nworst relative root errortwo factorsunbalancedbalancedconditioning alonedashed: the two predicted linesbalancing removes the second factor

The second factor was the last row

Wilkinson's polynomial in Chebyshev coefficients loses its roots exponentially from the colleague matrix, and two factors accounted for the loss: the roots' own conditioning, and the size of the matrix's last row against the rest of it. A diagonal similarity by powers of two — balancing, which changes no eigenvalue and rounds nothing — removes the second factor entirely: the balanced error sits within a factor of two and a half of the conditioning at every degree from fourteen to forty, three and a half thousand times better than unbalanced at forty. What is left is the conditioning, still exponential. Balancing the monomial companion matrix takes its norm down by forty-five orders and gains nothing, because there the error was never the matrix's.

16243240485664728010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴order nerrormiddle eigenvalues, errormiddle eigenvectorsworst eigenvector, componentwiseworst eigenvector, normwisefloored at 10⁻¹⁷the middle is the best of the spectrum

The middle eigenvectors come back exact

The Sylvester–Kac matrix's eigenvectors are known exactly, and the middle ones were predicted to be lost before their eigenvalues, since the gap between neighbouring eigenvalues is only two. Started from the all-ones vector, inverse iteration seems to confirm it: half the middle eigenvectors are wholly wrong at order 80 while their eigenvalues keep five digits. But the all-ones vector has, in exact arithmetic, no component at all along every other eigenvector. From a random start the middle eigenvectors come back exact to fourteen digits at every order measured, while their eigenvalues lose up to eleven. What is lost is at the ends of the spectrum, in the smallest entries of eigenvectors whose entries span twenty decades — and never by more than three units of roundoff times that span.

16243240485664728010⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³orderworst over the middleeigenvalue erroreigenvector, Schur routeeigenvector, inverse iterationmiddle eigenvalues, |λ| ≤ 5the route is the answer

The vector moves with its eigenvalue

The Sylvester–Kac matrix's middle eigenvectors came back exact to fourteen digits by inverse iteration while their eigenvalues lost eleven, and the route LAPACK uses — back substitution in the computed Schur form, then the accumulated orthogonal transformation — was predicted to lose them, at about the eigenvalue's error divided by the gap of two: 10⁻⁵ at order 80. It does, at 2.3·10⁻⁵, so the route is the whole of the difference: nine orders of magnitude between two ways of computing the same vectors from the same eigenvalues. The constant is not one over the gap. On all 274 eigenvectors whose eigenvalues are off by more than 10⁻¹¹, at every order from 40 to 80, the Schur route's error is 0.76 to 0.84 times its eigenvalue's error, and 0.79 is a number the exact vectors carry: how fast a unit Krawtchouk vector turns as its eigenvalue moves. The Schur route returns the right eigenvector of a wrong eigenvalue.

All essays