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

The fifth author

Four authors of a wrong answer have been named on this site 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.
36 min read 29 figures Whose fault is itWhere the algebra stops being the arithmeticMeasured, not assertedTwo routes to a number

Worth reading first: The exact answer to a nearby problem · An index that is a pair · Orthogonal is a number.

Three errors and one number is this field’s inventory. It sets out four authors of a wrong answer, the first two from the identity this site is built on and the second two from the fields that came after:

author what it is who chose it visible in a residual
the problem’s conditioning a property of the question nobody no, but estimable
the algorithm’s backward error a property of the code the library author yes, directly
the reformulation a property of a substitution the library author, once no
the approximation a number somebody typed the caller, every time no

Every one of those four is a statement about one computation: this answer, and how far it is from the truth, and whose fault the distance is.

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 = 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
Fig. 1 The identity as this field first drew it: two factors and one visible product, with two later terms in front of it.

There is a fifth, it does not fit the table’s shape, and the way it does not fit is the point of this essay.

The measurement

One vector of 4,096 binary64 numbers. One summation algorithm. One precision. Twenty-six ways of dividing the work across workers, which is what a parallel reduction is.

Twenty-one distinct answers.

Now apply the four columns. The problem’s conditioning is one number and it is the same for all twenty-six. Each run’s backward error is at the rounding level — every one of them is the exact sum of a vector within a rounding of the given one, which is the definition the exact answer to a nearby problem sets out. Nothing was reformulated: a sum is a sum. Nothing was approximated: no function was replaced, no degree chosen, no target set named.

Four authors, all silent, and twenty-one answers.

What the fifth author is

The implementation. Not the algorithm — the algorithm is fixed and is the same in every run — but the choices a machine makes about how to execute it: how many pieces a reduction is divided into, where the cuts fall, whether a multiply-add was fused, which kernel a length selected.

Set beside the other four it is different in every column:

conditioning backward error reformulation approximation the implementation
what it is a property of the question a property of the code a property of a substitution a number somebody typed a property of the machine
who chose it nobody the library author the library author the caller nobody, and no record is kept
when fixed fixed once every call every run
in a residual no, estimable yes no no no
about one run one run one run one run a set of runs

The last row is the one that makes it a fifth term rather than a special case of one of the others. Conditioning, backward error, reformulation and approximation are all quantities attached to a single computation and its single answer. The implementation’s contribution has no meaning for one run: a single answer is exactly as correct as the identity says it is, and the fifth author appears only when a second run exists to disagree with it.

Where it goes in the identity

The identity, as it stood after the two newest fields:

forward error  ⪅  ‖g − r‖  +  κ(linearised) × backward error.

The fifth term does not go in it, and saying why is more useful than forcing it.

Every term above bounds the distance between one computed answer and the truth. The fifth quantity bounds the distance between two computed answers, and the two statements have different subjects. What can be said, and is measured in a bound every answer satisfies, is a companion inequality:

|answer(machine A) − answer(machine B)|  ≈  0.25 · κ_sum · u,

where κ_sum is the summation condition number of the reductions involved — a quantity computable in one pass over the caller’s own data, stable to 2% over ten decades, and absent from every library’s output.

Two statements, then, rather than one. The identity says how wrong an answer may be. The companion says how much two answers may differ. The second is not implied by the first, because the first is satisfied by both answers with room to spare — 7,932 times the room, measured — and a bound covering both cannot separate them.

It is a backward error, and that is the reframing that helps

The most useful thing the field found is that the fifth author is not a new kind of quantity after all, once it is measured in the right units.

What a regression test can ask for measures the run-to-run variation on a solve at 3.2·10⁻¹² and the problem’s amplification at 4.3, so the machine’s disagreement is equivalent to a perturbation of the data of about 7.5·10⁻¹³ — a few thousand unit roundoffs, which is what several hundred amplified roundings in a 674-step iteration comes to.

So the fifth author perturbs the problem, exactly as a backward error does, and differs from the second author in one respect: the second is a property of the algorithm and is the same on every machine, while the fifth is drawn afresh from the implementation on every run.

That reframing pays for itself twice.

It says the magnitude is predictable. A backward error has a size, and this one’s is κ_sum·u.

And it says what a test can ask for. No procedure can distinguish two problems that differ by less than the backward error committed on them, so no regression test can catch a defect smaller than the implementation’s own perturbation. That is not a limitation of testing; it is the subject’s central bound, appearing in a place nobody expected it.

Three implementations, three mechanisms

The fifth author is one term in the inventory and three distinguishable mechanisms underneath it, and a diagnosis needs to know which.

The partitioning of a reduction. How many pieces, and where the cuts fall. Varies between machines, between runs on one machine with a work-stealing scheduler, and with an environment variable. It is the field’s main subject, its size is κ_sum·u, and it has a clean repair: the sum that cannot be wrong, at three operations an element.

Expression contraction. Whether a multiply-add was fused, decided at compile time from source that says nothing about it. One rounding wide, and enough to compute a determinant of one as zero. Repaired by writing the fused form explicitly, or by forbidding the contraction — and, measured over two hundred matrices, neither form is the correct one: a matrix that is definite on one machine finds the unfused build right seventeen times out of twenty-six.

Algorithm selection. Which kernel a size selected, which is a constant in a library’s source. Has no repair at the call site at all, and is the one that survives after everything else has been made deterministic.

The three are independent, they compose, and they are repaired by three different actions — a summation policy, a compiler flag, a pinned version. A diagnosis that names the fifth author without naming which of the three is a diagnosis that cannot be acted on.

The five, as a diagnostic

The practical use of an inventory is a sequence of questions, and the fifth adds one. Three errors and one number ends with four; here they are with the fifth, in the order worth asking them.

Was the problem sensitive? Its condition number. Rules out the answer most people reach for first.

Was the algorithm unstable? Its backward error against the object it factorised.

Was something reformulated? If a quadratic became a pencil or a least-squares problem became a normal-equations solve, the reformulated object’s conditioning is a third quantity.

Was anything approximated? If a function was replaced before the arithmetic began, the approximation error is a term in front of everything and is usually the largest.

And: does the answer change when nothing does? Run it again with a different thread count. If the answer moves, the fifth author is present, its size is κ_sum·u, and no amount of care in the other four will remove it.

The fifth question is the only one of the five that is answered by running the computation again rather than by computing a quantity, which is a fair description of how different it is.

Why it took nineteen fields to reach

The same question three errors and one number asks about the third and fourth authors, and the answer here is different from both.

The third and fourth were found by sweeping a parameter that was supposed to be free — a change of units that is exact, a target set that changes no eigenvalue — and watching something move. That technique needs a parameter, and the fifth author has none: there is no argument to a computation that says how many workers to use, and nothing in a program’s text mentions it.

What found this one instead was running the same thing twice, which is not a technique anybody applies to a deterministic program. A numerical analyst measures an error against a truth. A software engineer runs a test once and records whether it passed. Neither of those workflows has a step at which one computation is executed several ways and the answers compared with each other rather than with anything external — which is exactly the measurement every figure in this field is.

The instrument is trivial and it is not in either tradition’s toolkit. That is the whole reason a site with nineteen fields of error analysis had no term for a quantity that is present in every parallel program running today.

What it costs to have missed it

Two costs, and the second is the expensive one.

A misattributed diagnosis. A team that sees two machines disagree and has only four authors available will assign it to one of them: the problem must be ill conditioned, or the library must have a bug, or the data must have changed. All three are investigable, none is the cause, and the investigation ends without an explanation — which is where floating point is non-deterministic enters as folklore, correct and useless because it comes without a number.

And a test suite that is quietly disabled. A regression test that demands more agreement than the arithmetic can supply fails on correct builds, and a test that fails on correct builds is a test people stop believing and then stop running. The defect that arrives afterwards is caught by nothing. That sequence is the most damaging consequence of the missing term, and it is entirely preventable by measuring the band once.

What the four authors do to the fifth

The five are not a list of alternatives; they compose, and two of the compositions are worth naming because they are where the field’s numbers get large.

Conditioning amplifies it. The fifth author perturbs the problem by about κ_sum·u in the reductions and the problem’s condition number amplifies that perturbation into the answer, exactly as it amplifies a backward error. So a badly conditioned problem has a larger disagreement for the same reason it has a larger error, and a stopping test is a race measures both together: at κ = 100 every partitioning takes the same 21 steps, and at κ = 10⁴ they take between 674 and 690.

A chaotic algorithm amplifies it further. Conjugate gradients builds its search directions from its own residuals, so two runs that differ at step three build different subspaces from step four onward. The perturbation is 10⁻¹³ and the trajectories separate over hundreds of steps into a difference of sixteen iterations. That amplification belongs to the algorithm rather than to the problem, which means the second author and the fifth interact in a way the identity’s product form does not capture.

The honest statement, then, is that the fifth author’s contribution is κ_sum·u and its effect is that times whatever the problem and the method amplify it by — which is the same two-factor structure the site was built on, applied to a perturbation nobody chose.

What the fifth author is not

Three exclusions, because a term this broad attracts things that belong elsewhere.

It is not the conditioning. Two machines are handed identical data. Nothing is perturbed on the way in; the perturbation is committed by the execution.

It is not randomness in the algorithm. A randomised method draws from a seed the program owns and can record, so its variation is reproducible by construction — a bound that holds with probability is where this site keeps that vocabulary, and the difference is that a seed is written down.

And it is not a bug. Every run is correct. The whole difficulty is that correctness, as the subject defines it, is a property each of the twenty-one answers has.

What a library would have to change

The inventory is only worth extending if the extension implies an action, so: five things, none of them research.

Return κ_sum beside a reduction. One pass, no estimator, and it turns the answers might disagree into a number.

Return the partitioning. One integer, known to the runtime and to nobody else, and enough to reproduce a result rather than merely to observe that it cannot be.

Document the cutoffs. Which kernel a size selects is already in the source; saying it affects results costs a line.

Offer an order-independent mode. Three operations an element for the cheap version, and the literature for it is decades old.

And say all three sentences in the specification rather than only the first — that the answer is the exact result for nearby data, that two runs may differ by about κ_sum·u, and that the reproducible mode makes the answer a function of the inputs alone.

Every one of those is available today. What was missing was a reason to think of them as one thing, which is what a name for the fifth author is for.

The one line

Four authors describe how far one answer is from the truth. The fifth describes how far two answers are from each other, it is a backward error the implementation commits rather than the algorithm, it is 0.25κ_sum·u, and every instrument this site has built reports zero for it.

At other settings

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. 2 The identity at a smaller size, where the fifth term is the same size and the others are smaller.
The same 4096 numbers, added up 26 ways, 21 different answersEach dot is one reduction of one vector of 4096 binary64 numbers whose summation condition number is 1.01·10⁸: the vector is cut into p contiguous pieces, each piece summed left to right, and the pieces combined in a tree. The vertical position is the distance from the exactly rounded sum, relative. Filled dots are the static split a fixed thread count gives; open dots are unequal splits, which is what a work-stealing scheduler produces at the same p. 21 of the 26 runs returned distinct values, spanning 1.55·10⁻¹³ — 2.3·10⁷ ulps of the answer. The dashed line is the classical bound γ₍ₙ₋₁₎Σ|xᵢ|, which every one of them satisfies with 2.58·10⁴ to spare, and which is the same number for all of them: it contains n and Σ|xᵢ| and nothing about the order.110¹10⁻¹¹10⁻⁸10⁻⁵pieces the vector was divided intodistance from the exact sum, relativethe published boundκ · uone vector, one algorithmdistinct answers21runs26spread, in ulps2.3·10⁷κ of the sum10⁸bound ÷ worst error2.6·10⁴nobody chose pand no answer is the answer
Fig. 3 The measurement: one algorithm, one precision, twenty-one answers.
What the runs disagree by, and what the bound says they mightOver ten decades of summation condition number, on vectors of 2048 numbers: the measured spread between seven partitioned reductions of the same vector, the quantity κu, and the classical error bound γ₍ₙ₋₁₎Σ|xᵢ| relative to the sum. The measured spread and κu are one curve times a constant — the ratio runs from 0.253 to 0.258 across the whole sweep, a spread of 1.021 while both quantities move by ten orders. The bound sits 7932 to 8100 above the disagreement and is identical for every ordering, which is what makes it useless for the question a reader actually has: not how wrong is this, but will the next run say the same.10²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹κ of the sumrelative sizethe boundκ · umeasured spreadtwo curves and one constantspread ÷ κu, low0.25spread ÷ κu, high0.26bound ÷ spread, low7932decades swept10the spread is computablethe bound cannot see the order
Fig. 4 The companion inequality, and the bound that cannot substitute for it.
The summation condition number of six inner products this site already computesκ = Σ|xᵢ| ÷ |Σxᵢ| for six reductions taken from computations elsewhere on this site, at 64 terms: an orthogonality check between two columns of a Q, one component of a residual b − Ax on an ill-conditioned system, a nearly degenerate 2 × 2 determinant, a Hutchinson trace probe zᵀAz, a residual norm rᵀr, and a curvature pᵀAp. The range is 1 to 1.01·10¹⁷. The two at the safe end are sums of squares, where no term can cancel another; the three above 10¹⁰ are the orthogonality check, the residual and the determinant — which is to say every quantity a stopping test, an orthogonality test or an orientation test is written in. The reductions that are reproducible are the ones nobody makes a decision from.κ of the suma component of b − Ax1.01·10¹⁷ad − bc, near-degenerate3.6·10¹⁶qᵢᵀqⱼ, an orthogonality check7.39·10¹⁵pᵀAp, a curvature409zᵀAz, a trace probe41.8rᵀr, a residual norm1measured, not assumedhighest10¹⁷lowest1above 10¹⁰3terms64sums of squares are safeand nobody decides anything from one
Fig. 5 The quantity the companion is written in, measured on this site’s own reductions.
The error of a reduction is a walk, and the number of pieces changes which walkRunning error of two reductions of the same 2048 numbers, against how many terms have been consumed: one accumulator, and 8 pieces combined at the end. Each is the distance between what the arithmetic has accumulated so far and the exact sum of the same terms, carried in a compensated pair so the instrument does not round. The two separate within the first few terms and never meet again, ending at 7.56·10⁻¹⁴ and 2.41·10⁻¹⁴. The mechanism is in the horizontal axis: a rounding is committed at every partial sum and its size is that partial sum's spacing, and this vector's running total reaches 32.1 on the way to an answer of 1.63·10⁻⁵. Dividing the work shortens each walk, so the mean error over eight vectors falls from 4.78·10⁻¹⁴ at one piece to 8.14·10⁻¹⁵ at sixty-four: more workers is more accurate, which is the opposite of what the word *parallel* suggests here.051210241536204810⁻¹⁷10⁻¹⁴terms consumederror accumulated so farone accumulator8 piecesthe steps are the partial sumspeak partial sum32the answer1.6·10⁻⁵error, one piece7.6·10⁻¹⁴error, 8 pieces2.4·10⁻¹⁴mean at p = 648.1·10⁻¹⁵the walk sets the sizeand nothing sets the value
Fig. 6 Where it comes from.
The smallest defect a regression test can catch, on a machine that will not repeat itselfA 150 × 150 system with κ = 10⁴, solved by conjugate gradients to 10⁻¹⁰ on six partitionings. The shaded band is how far the answer moves across those six when nothing at all is wrong: 3.24·10⁻¹², relative. Each pair of dots is one defect — a single entry of the matrix changed by a relative δ, an assembly slip or a stale coefficient — solved on the same six machines, showing its nearest and furthest run. A defect is catchable only when its nearest run clears the band, because a test whose threshold sits inside the band fails a correct build. Everything at or below 10⁻¹³ overlaps and is invisible; 10⁻¹² is the smallest that separates, at 4.61·10⁻¹². So the tolerance exists, it is bracketed between 3.24·10⁻¹² and 4.61·10⁻¹² — a factor of 1.42 — and it is neither zero nor the 10⁻⁸ that usually gets typed.-15-13-11-910⁻¹³10⁻¹¹10⁻⁹10⁻⁷log₁₀ of the defect, relative to the entry it sits inhow far the answer movedthe machine, on its owna tolerance with two sidesmachine band3.2·10⁻¹²smallest caught10⁻¹²its nearest run4.6·10⁻¹²defects hidden3window, factor1.4below the band nothing is visibleand above it everything is
Fig. 7 And what it costs a test suite.
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. 8 The two-factor identity, complete for every field before the polynomial one.
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. 9 The first author, which is the same for every run.
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. 10 The standing all five are measured against.
One quadratic eigenvalue problem in 9 systems of units: what the solver reports and what the answer is worthAn overdamped chain of 6 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 7.1·10⁻¹⁶ at γ = 1 and 2.57·10⁻¹² at γ = 108 — it moves by a factor of 7022 while the other two move by 7.37·10¹¹. The backward error for the QUADRATIC, which is what the person who posed the problem is entitled to, grows by 7.03·10¹¹ across the same sweep, and the forward error follows it: 9.14·10⁻¹⁵ to 0.00674. 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), worst2.6·10⁻¹²η(quadratic), worst4.6·10⁻⁴forward error, worst0.0067coefficient spread4.2·10¹⁵the solver is right at every stopabout a problem nobody asked
Fig. 11 Two backward errors for one eigenvalue, which is where the third author was found.
Three residuals, and the one a solver returns is the one about nothingFor a rational approximant with m poles, on a target set reaching to 0.3: ‖T̃(λ)x‖ — the residual against the problem the eigensolver was handed — stays at 1.2·10⁻¹⁵ to 1.5·10⁻¹² — it does not fall with m, and past a point it slowly rises, because each added pole makes the fit's own basis worse conditioned; ‖T(λ)x‖, which costs one further evaluation of γ√(λ + c), falls with the approximation; and the forward error against the closed form falls with it, staying 20.5 to 24.8 times larger. The free number rises by 1306 across the sweep while the answer improves by 885. So the free number says nothing and the nearly-free number says almost everything. Evaluate the residual against the function you asked about, not against the one you handed over is the whole practical content of this field.23456710⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²poles in the approximantresiduals and errorforward error‖T(λ)x‖‖T̃(λ)x‖one extra evaluationagainst the approximant1.5·10⁻¹²against the problem asked1.9·10⁻⁶forward error4.8·10⁻⁵‖g − r‖ there8.5·10⁻⁵the free residual is flatand the answer is not
Fig. 12 Three residuals for one answer, which is where the fourth was.
The error that happens before any arithmetic: |g − r| on the target setγ√(λ + c) minus its approximant, plotted against the distance from the branch point at −0.4. The target set runs from -0.2 to 8.99, and its left end is 0.2 from the branch point — a choice made before anything is computed. The rational approximant with 8 poles clustered on the cut reaches 8.46·10⁻⁴; the polynomial of degree 9, which linearises to a matrix of the same size, reaches 0.00676. Both errors are committed at this stage. No residual computed later contains either of them, because a solver's residual is taken against whichever of these two functions it was handed.110¹10⁻⁷10⁻⁵10⁻³distance from the branch point, λ + c|g − r|the eigenvaluesdegree 98 poleschosen, not computedreach0.06left end, from the cut0.2rational, worst8.5·10⁻⁴polynomial, worst0.0068linearisation size, both54committed before the solveand invisible to it
Fig. 13 The term the free residual cannot contain.
Six ways of computing one spectrum, at a change of units of 106Three linearisations — the first companion form, the second, and the symmetric member of the DL(Q) family — each reduced to a standard eigenvalue problem in two ways: by inverting the leading coefficient, whose eigenvalues are λ, and by inverting the trailing one, whose eigenvalues are 1/λ. All six have exactly the eigenvalues of the quadratic in exact arithmetic. Measured against the closed form at γ = 106, the best is first/leading at 5.473·10⁻⁶ and the worst second/trailing at 2.234·10⁻⁴, a spread of 40.81. Under the leading reduction the symmetric form and the first companion form are THE SAME MATRIX — −A₁⁻¹A₀ is [[−M⁻¹C, −M⁻¹K], [I, 0]] for both — so the distinction between them is one the inversion discards, and it survives only under the trailing one.first · leading5.47·10⁻⁶first · trailing4.62·10⁻⁵second · leading4.27·10⁻⁵second · trailing2.23·10⁻⁴symmetric · leading5.47·10⁻⁶symmetric · trailing4.05·10⁻⁵all six are the same algebrabest route5.5·10⁻⁶worst route2.2·10⁻⁴spread across the six41condition of the linearisation8.3·10¹²the spectra agreeand the arithmetic does not
Fig. 14 Six routes to one spectrum, differing by a factor a residual cannot see.
One system, one tolerance, and 11 different amounts of workConjugate gradients on a 200 × 200 symmetric positive definite matrix with κ = 10⁴, stopped when the relative residual falls below 10⁻¹⁰, run at each of 13 partition counts. The only thing that differs between the runs is how many pieces the two global inner products were summed in — which on a real machine is the number of workers. Iteration counts run from 674 to 690, 11 of them distinct, a spread of 2.4% of the work. Every run converged and every answer is right to the accuracy asked for; the forward errors span 6.31·10⁻¹¹ to 1.08·10⁻¹⁰, a factor of 1.72. What differs is not a digit at the end of an answer — it is the bill, and it is an integer.110¹670673676679682685688691694pieces the inner products were summed initerations to the tolerance674, the cheapest run690, the dearestthe same solve, pricedpartitionings run13distinct counts11spread, per cent2.4best forward error6.3·10⁻¹¹worst1.1·10⁻¹⁰one matrix, one toleranceand the cost is the machine's
Fig. 15 The fifth author deciding an integer.
The numerical rank of one matrix is 10, 12, or somewhere betweenSingular values of one 60 × 14 matrix whose spectrum falls by a fixed factor per index, so there is no gap anywhere and the threshold decides the rank. Each curve is the same matrix, its Gram matrix formed with the inner products summed in a different number of pieces. The curves lie on top of each other for the first several values and separate below about 7.3·10⁻¹³, which is where forming AᵀA has put the rounding. The dashed line is the threshold, σ₁ · 10⁻¹⁴. Counting the values above it gives 12 at p = 1, 12 at p = 2, 12 at p = 4, 10 at p = 8, 10 at p = 16, 11 at p = 32, 11 at p = 60 — the rank of one matrix, as a function of how many workers were available. Not a digit of an answer: the number of columns a model built from this matrix would have.13579111310⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹index ksingular valueσ₁ · 10⁻¹⁴, the thresholdone matrix, three rankspartitionings7lowest rank10highest rank12threshold7.3·10⁻¹⁴σ₁7.3the curves separate in the noiseand the threshold is drawn through it
Fig. 16 And deciding the shape of a model.
Two conforming builds, one matrix, and 26 different verdicts in 200 triesForty 6 × 6 Gram matrices at each of five conditionings, each formed exactly from a tall matrix whose last column is nearly the sum of its first two. Every one of them is asked *are you positive definite* by a Cholesky, once with the inner products accumulated with a single rounding per term and once with each product rounded before it is subtracted. The bar is how many of the forty came back with different verdicts: 7 at 10⁻⁸, 2 at 10⁻⁹, 4 at 10⁻¹⁰, 10 at 10⁻¹¹, 3 at 10⁻¹². The two shaded parts split those by which build was right, settled by the signs of the leading principal minors computed in BigInt rationals — 9 to the fused form and 17 to the rounded one. So neither is the correct implementation. The verdict is a property of the build, and the answer to which build is right is a property of the matrix.-12-11-10-9-8024681012log₁₀ of how nearly dependent the columns areverdicts that disagreed, of 40724103which one is correctmatrices tested200verdicts disagreed26fused was right9unfused was right17lower part: thefused build was righta sign has no last digitso a verdict has nowhere to hide
Fig. 17 And deciding a verdict, with an exact adjudicator saying neither build is the right one.
A determinant of one, computed as zero by the expression that is written downThe matrix [[x+2, x+1], [x+3, x+2]] at x = 2ᵏ has determinant exactly 1 at every k. The two products are x² + 4x + 4 and x² + 4x + 3, so they differ by one part in x², and above k = 26 they round to the same binary64 number. Relative error of two evaluations: ad − bc as written, which is exact until k = 27 and then returns zero at every larger size; and Kahan's form, w = bc, e = bc − w recovered by the fusion, (ad − w) − e, which is exact throughout. The difference between them is one rounding of one product — a relative 10⁻¹⁶ — and it is the whole of the answer, because the subtraction that follows removes everything else. Both forms are IEEE-754 conforming and a compiler may emit either from the same source.2022242628303210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²k, where the entries are near 2ᵏrelative error of the determinantas writtenfused: exactproducts need 54 bitsone rounding, the whole answertrue determinant1naive, k = 300fused, k = 301first wrong at k27sizes returning 06both forms conformand the source does not say which
Fig. 18 The smallest instance of it, one instruction wide.
Five summation policies, 400 permutations each, and how many answers each returnedThe same 1,024 numbers, summation condition number 10⁸, presented in 400 different orders. The bar is the number of distinct binary64 values the policy returned; the figure beside it is the worst relative error it made. A single accumulator returns 303 values and eight pieces return 72; Kahan's compensated loop is the most accurate of the three at 4.96·10⁻¹⁰ and still returns 119. Pre-rounded summation returns one value, at an error of 4.91·10⁻⁶ — four orders worse than the compensated loop and the same on every machine. Exact accumulation returns one value and the right one. The comparison is bitwise because the claim is bitwise: a tolerance here would pass everything.distinct answersone accumulator303eight pieces72compensated119pre-rounded1exact1worst error 1.17·10⁻⁸worst error 4.61·10⁻⁹worst error 4.96·10⁻¹⁰worst error 4.91·10⁻⁶worst error 0bitwise, or not at allpermutations400one accumulator303pre-rounded1its error4.9·10⁻⁶compensated error5·10⁻¹⁰accuracy and agreement are different propertiesand the accurate one is not the agreed one
Fig. 19 The property that removes it.
Six ways to add up a vector, priced against what they returnBinary64 operations per element on the horizontal axis and relative error on the vertical, for one vector of 2,048 numbers with a summation condition number of 1.05·10⁶. Filled markers are the policies whose answer does not depend on the order the terms arrive in; open markers are the ones whose answer does. The cheapest order-independent policy costs 3 operations an element — a divide, a round and an add — and is 2.29·10⁻⁷ wrong; the most accurate policy here is pre-rounded, two passes at 7. Nothing sits in the bottom left. The cost of an answer that is the same on every machine is between three and twelve times the cost of an answer that is not, and the two-pass pre-rounded policy is the one worth knowing about: seven operations, order-independent, and accurate.02468101210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸binary64 operations per elementrelative errorruntime orderfixed partitionpre-rounded, one passpre-rounded, two passescompensated (Kahan)exact accumulationwhat determinism costsreproducible policies3cheapest, ops/element3its error2.3·10⁻⁷two passes10⁻¹⁷ordinary, 8 pieces2.7·10⁻¹¹nothing is cheap and exactthe middle of the figure is the answer
Fig. 20 And what that property costs.
Every column's 1-norm on a 12×12 matrix built to stop Hager's walk one column short12 bars, one per column, each the sum of the absolute values in that column — the quantity the estimator maximises, whose largest is the matrix's 1-norm. The walk visits 1 of them and returns 12.00; column 2 has 1-norm 114.00 and is never visited, because its entries alternate in sign and cancel against every sign vector the walk stands on.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
Fig. 21 An estimator for the first author, which exists.
What Hager's estimator returns, as a share of the truth, against the size of the matrix built to defeat itThe estimate over the true 1-norm against n, both axes logarithmic. It falls from 0.526 at n = 4 to 0.0658 at n = 32, along the line 1/t where t is the construction's own multiplier — which is bounded only by the size of the matrix, so the ratio has no floor. A line at one marks a correct estimate.10¹10⁻²10⁻¹1size of the matrixestimate ÷ true 1-norma correct estimatewhat it returns1 / tno floorratio at n = 40.53ratio at n = 320.066products, either size5the estimate is always a lower boundwhich is the direction that flatters the matrix
Fig. 22 And the case built to fool it, which is the nearest the other four come to being unknowable.
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. 23 The knob that moves four of the five and leaves the fifth’s ratio where it was.
Four pivot rules on [[ε, 1], [1, 1]] with ε = 10⁻¹², before and after its first row is scaled by 1/εEight bars of forward error. On the matrix as given, no pivoting loses the answer entirely and all three pivot rules return it to rounding. Multiply the first row by 1/ε — which does not change the solution — and partial pivoting makes no interchange at all and returns the same wrong answer that no pivoting does. Scaled partial pivoting and complete pivoting are unmoved.forward error, relative to a solution of exactly (1, 1)no pivoting · as given1.3·10⁻⁴ 0 interchangesno pivoting · rows scaled1.2·10⁻⁴ 0 interchangespartial · as given2.2·10⁻¹⁶ 1 interchangepartial · rows scaled1.2·10⁻⁴ 0 interchangesscaled partial · as given2.2·10⁻¹⁶ 1 interchangescaled partial · rows scaled2.2·10⁻¹⁶ 1 interchangecomplete · as given1.1·10⁻¹⁶ 1 interchangecomplete · rows scaled2.2·10⁻¹⁶ 1 interchangethe same problem twicepartial, as given2.2·10⁻¹⁶partial, rows scaled1.2·10⁻⁴its relative residual1.2·10⁻¹⁶complete, rows scaled2.2·10⁻¹⁶the two systems have the same solutionand one pivot rule cannot see it
Fig. 24 A decision that looks like an algorithm and is a scaling.
The accuracy worth paying for in a preconditioner, against the condition number of the problemFour complete cost sweeps, each reduced to the ε that minimised the total work, against the κ of the problem it was run on. The shift that moves κ over 2.7 decades leaves every distance between every pair of clusters exactly where it was, so the ranks the partition asks for stay at 10, 9, 7, 5 throughout and nothing geometric is moving. The optimum walks 0.5 → 0.5 → 10⁻⁶ → 10⁻⁸ — six decades of accuracy, bought because the problem got harder and for no other reason. On the easiest problem the best hierarchical preconditioner in the comparison has rank one; on the hardest it is the tightest one on the sweep. How accurate an approximate inverse should be is a question with an answer, and the answer is not in the matrix's structure.10¹10²10³10⁴10⁵-9-7-5-3-11condition number of the problemlog₁₀ ε at the cheapest totalrank 1rank 1rank 6rank 5the knob answers to κκ = 210.5κ = 1.1·10⁴10⁻⁸rank at the easy end1rank at the hard end5ranks at 10⁻⁸, spread5the geometry did not moveand the right accuracy did
Fig. 25 An accuracy with an optimum in a size, which is the shape a caller expects and does not always get.
Classical and modified Gram–Schmidt: the same subtraction, in a different orderTwo panels of pseudocode differing in one argument, with the resulting pairwise dot products of the computed columns listed beneath each.for each previous column i, subtract the projection of column j onto qᵢclassicalr[i][j] = qᵢ · a[j] ↑ the ORIGINAL columnv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 7 and 8: 1columns 6 and 8: 0.13columns 6 and 7: 0.13modifiedr[i][j] = qᵢ · v ↑ what is LEFT of itv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 1 and 8: 4.4·10⁻⁷columns 2 and 8: 2.7·10⁻⁷columns 3 and 8: 2.4·10⁻⁸The two R factors agree to 1.2·10⁻⁶ relative. The two Q factors do not.the 8×8 Hilbert matrixone word, eight orders
Fig. 26 Two orderings of one orthogonalisation, in the field that made the point first.
The proved bound against the error it bounds, at 24 significand bitsTwo curves against the Hilbert size on a logarithmic vertical axis: the error the computation actually made, which is known because the exact answer is the integers by construction, and the bound proved for it in interval arithmetic. The bound is above the error at every size, by 9.949 to 23.12 times. Past n = 5 the method refuses to return a bound at all.345678910111210⁻⁵10⁻³10⁻¹Hilbert size nerrorrefusedrefusedrefusedrefusedrefusedthe proved boundthe actual errorbound ÷ errorslack at n = 313slack at n = 423slack at n = 59.9the error is known, not estimatedthe exact answer is the integers 1 … n
Fig. 27 A bound measured against what it bounds.
Error of a difference quotient for J(x)v against ε, on the Bratu problem at n = 64The Jacobian of this problem is a formula, so the error of each quotient is measured against a derivative that is exact rather than against a better quotient. The forward difference falls with a slope of 1.01 — first order — reaches 1.28·10⁻¹⁰ at ε = 10^-6, and rises again with a slope of -1.00 as cancellation takes over. The central difference falls with a slope of 2.00 and bottoms at 1.11·10⁻¹² for two residual evaluations instead of one. Neither gets near the unit roundoff at any ε.10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice
Fig. 28 A floor under a computed difference, which is where a verdict stops meaning anything.
The accuracy asked for, against the accuracy obtained, for two kernels on one partitionε is applied to each block against that block's own largest singular value; the error is then reported against the whole matrix's norm. The two are not the same number and the dashed diagonal is where they would be. For 1/r the obtained error is 29 times smaller than the tolerance asked for; for log r on the identical partition it is 13 times smaller. Both curves are straight and parallel to the diagonal, so the knob does what a knob should — a decade in buys a decade out — and neither of them sits on it. The factor is how much of the matrix's mass lives off the diagonal, which is a property of the kernel; it is not the number of blocks, which is 112 and would put the curves on the other side of the diagonal.10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³ε asked for, per block‖A − A_H‖ ⁄ ‖A‖ obtainedasked = obtainedlog r1 ⁄ rtwo numbers, not one1/r, obtained at 10⁻⁸3.4·10⁻¹⁰log r, obtained at 10⁻⁸7.8·10⁻¹⁰1/r, obtained ⁄ asked0.034log r, obtained ⁄ asked0.078blocks in the partition112the tolerance is per blockand the error is per matrix
Fig. 29 What a tolerance actually obtains, measured rather than assumed.

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 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. 110¹10⁻¹¹10⁻⁸10⁻⁵pieces the vector was divided intodistance from the exact sum, relativethe published boundκ · uone vector, one algorithmdistinct answers21runs26spread, in ulps2.3·10⁷κ of the sum10⁸bound ÷ worst error2.6·10⁴nobody chose pand no answer is the answer The same program, twice One vector of 4,096 numbers, one summation algorithm, one precision, twenty-six runs — and twenty-one different answers. Nothing in the program chose between them, every one of them satisfies the textbook bound, and the exactly rounded answer is not among them. 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.

What links here

Computed from the collection, not written here: the essays that point at this one.

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.

Algorithm selectionApproximation before linearisationBackward errorBitwise reproducibilityCondition numberForward errorLinearisationReduction orderResidualRun-to-run variation