Numerical Linear Algebra
where the algebra stops being the arithmetic
The answer that depends on the machine

The tolerance that buys no agreement

Ask for four more orders of accuracy and you get them — the answers improve by a factor of 1.5 million. The ratio between the best and the worst run is 1.34, 1.48, 1.71 and 1.17 across the same sweep. The band falls and it does not close.
34 min read 27 figures Measured, not assertedWhose fault is itPrecision as the knob

Worth reading first: A stopping test is a race · The same program, twice · The zero you are allowed to write.

Every reader of a stopping test is a race has the same response, and it is the right response to have: if the runs disagree, converge further. The disagreement looks like the tail end of an answer that was not finished, and finishing it should remove it.

It does not, and the measurement is unambiguous.

Asking for more accuracy buys accuracy, and buys no agreement at allConjugate gradients on a 200 × 200 matrix with κ = 10⁴, solved at four tolerances, seven partition counts each. The bar at each tolerance spans the smallest and largest forward error the seven runs produced. The bars fall by 1.5·10⁶ across the sweep, which is the tolerance doing exactly what it is for. The ratio between the top and bottom of each bar is 1.34, 1.48, 1.71, 1.17, which does not fall with the tolerance. So the disagreement between machines is not a residue of an insufficiently converged answer that a tighter tolerance would remove; it is proportional to whatever accuracy was reached, and the runs stay a fixed factor apart all the way down.-12-10-8-610⁻¹³10⁻¹¹10⁻⁹10⁻⁷log₁₀ of the residual tolerance asked forforward error of the answer1.34×1.48×1.71×1.17×the band does not closetolerances swept4runs at each7accuracy gained1.5·10⁶ratio at 10⁻⁶1.3ratio at 10⁻¹²1.2the bars falland they keep their height
Fig. 1 The same 200 × 200 system solved at four tolerances, seven partition counts each. Each bar spans the best and worst forward error the seven runs produced.

The bars fall by 1.5·10⁶ across the sweep, which is the tolerance doing exactly what a tolerance is for: the worst run at 10⁻¹² is six orders more accurate than the best run at 10⁻⁶. The ratio between the top and the bottom of each bar is 1.34, 1.48, 1.71, 1.17 — no trend, and certainly no convergence towards one.

Why the intuition is wrong

The intuition treats the disagreement as a residue: two runs are near the answer, one is slightly nearer, and running longer takes both closer to the truth and therefore closer to each other. That picture is right for a method converging to a fixed point in exact arithmetic and wrong here for a reason that is one line.

The disagreement is not a distance from the answer. It is a distance between two trajectories that are both converging, and it is proportional to wherever they are. Two runs at 10⁻⁶ differ by about a third of 10⁻⁶; the same two runs at 10⁻¹² differ by about a fifth of 10⁻¹². The ratio is what stays fixed, because the mechanism producing it — a reduction disagreeing by κ_sum·u relative — is relative rather than absolute.

So converging further does buy something real: the absolute disagreement falls by six orders, exactly as the accuracy does. If a caller’s requirement is stated in absolute terms — the two machines must agree to 10⁻⁸ — tightening the tolerance meets it. If it is stated in relative terms, or as a bitwise comparison, no tolerance ever meets it.

That distinction is the practical content of this essay and it is worth stating as a rule: tightening a tolerance moves a disagreement, it does not remove one.

The same statement about precision

The identical argument applies to the other knob people reach for, and it is worth making explicit because the knob is this site’s founding one.

Where the disagreement comes from derives the size of the run-to-run variation as κu times a constant. It carries exactly one factor of u. So doubling the precision divides the disagreement by 2⁵³ — an enormous improvement in absolute terms — and leaves the ratio between two runs exactly where it was, and leaves the count of distinct answers exactly where it was.

A bitwise regression test fails as reliably at quadruple precision as at double. The answers are closer together; there are still many of them; and the test that asks whether they are equal gets the same answer.

This is the same shape as accuracy and agreement are different properties: every instrument that improves accuracy improves the absolute disagreement in proportion, and none of them touches the property that two runs return different bits.

Reading the four ratios

The four numbers are 1.34, 1.48, 1.71 and 1.17, and they are not monotone. That is worth an explanation rather than an average.

Seven runs is a small sample of a quantity that is itself a spread. The ratio between the largest and smallest of seven draws from a distribution is a noisy statistic — its own variability is comparable with the differences between the four values here — so the honest reading of the sequence is no trend, not rises then falls.

What the figure asserts whenever the figure is drawn is the weaker and safer claim: the last ratio is not smaller than the first by more than a factor of two, and every tolerance shows more than one distinct iteration count. Both are checked, and the second is what stops the figure from being a picture of four coincidences.

At n = 80 the same sweep gives 1.04, 4.83, 14.83, 1.36 — a much wider scatter, from a problem with fewer terms in its reductions and therefore a noisier disagreement. The claim survives that too, and it is the reason the assertion is written about the trend rather than about the flatness.

The iteration counts do the same thing

The figure plots forward errors, and the iteration counts underneath it tell the same story from a direction that is easier to act on.

At a tolerance of 10⁻⁶ the seven runs take 436 to 447 steps — six distinct counts, a spread of 2.5%. At 10⁻⁸: 558 to 572, six distinct, 2.5%. At 10⁻¹⁰: 676 to 690, six distinct, 2.1%. At 10⁻¹²: 800 to 803, four distinct, 0.4%.

Two readings, and the second is the more interesting.

The spread is roughly a fixed fraction of the work, not a fixed number of steps. That is the same relative-versus-absolute point the forward errors make: the disagreement scales with wherever the computation is rather than sitting at a fixed level.

And it narrows at the tightest tolerance, which is the one place the sweep shows something the simple story does not predict. At 10⁻¹² the method is close to the floor its arithmetic allows, the residual curve is flattening, and a flattening curve crosses a horizontal line over a narrower range of steps than a steeply falling one. So the step count becomes more reproducible while the answers stay a factor apart — the two quantities decouple, and a reader who was watching only the iteration count would conclude, wrongly, that convergence had fixed something.

That decoupling is worth carrying. An iteration count agreeing is not the answers agreeing, and of the two it is the cheaper one to observe and the less informative.

What a requirement should say

The practical residue, since somebody has to write the acceptance criterion.

A bitwise requirement is meetable and expensive. It needs order-independent reductions everywhere in the computation that reaches the compared quantity, and what determinism costs prices it at 1% to 122% of a solve depending on how cheap the matrix is. It is worth asking for when a verdict is downstream. It is not achieved by any tolerance.

A relative requirement is meetable and nearly free, and the number is κu. Two runs agree to within 10⁻¹⁰ relative is a statement the arithmetic can support if the reductions involved have κ_sum below about 10⁶, and cannot if they do not. Computing κ is one pass; the vector that hides it is how to know which reductions to compute it for.

An absolute requirement is meetable by tightening the tolerance, which is the one case where the instinct this essay is arguing against is right. If the requirement is that two runs agree to 10⁻⁸ absolute and the current disagreement is 10⁻⁷, converging two orders further meets it.

Three requirements, three different answers, and the mistake is to state one and act as though it were another.

Where the instinct does come from

It is worth taking the instinct seriously rather than only refuting it, because there is a class of computation in which it is exactly right.

A fixed-point iteration converging to a unique fixed point — a contraction, in the strict sense — forgets its history. Two runs that differ at step k are both within the contraction’s basin, and each subsequent step multiplies their difference by the contraction factor, so they converge to each other as well as to the answer. Run one long enough and two machines agree to the fixed point’s own accuracy, whatever they did on the way.

That is the picture the instinct comes from, and it describes a rate that is known in advance’s stationary methods exactly: a Jacobi iteration run to convergence returns the same answer from any starting vector and, near enough, from any arithmetic.

Conjugate gradients is not that. It is a direct method degraded by rounding rather than a contraction — its Krylov space is built from its own residuals, so two runs that differ build different subspaces and minimise over them, and there is no mechanism pulling the two trajectories together. An orthogonalisation nobody calls one is the essay about how far apart those subspaces get.

So the honest general rule is: the instinct is right for contractions and wrong for Krylov methods, which is unfortunate, because Krylov methods are what anybody solving a large system is running.

What is actually being paid for

A closing observation about where the money goes, because the sweep is also a cost measurement.

Going from a tolerance of 10⁻⁶ to 10⁻¹² costs the method about 80% more iterations — 441 to 801 on this problem, since conjugate gradients converges at a rate set by √κ and the extra digits are linear in the step count. That is a real price and it buys six orders of accuracy.

Going from an ordinary reduction to a pre-rounded one costs 1% to 22% of a step depending on the matrix and buys bitwise agreement.

So the two interventions are not competing for the same budget and they do not buy the same thing. The expensive one buys accuracy. The cheap one buys agreement. A team that has spent 80% more iterations trying to make two machines agree has bought the wrong thing at ten times the price, which is the specific waste this essay exists to prevent.

What to do when the requirement is somebody else’s

The awkward case, and the common one: the acceptance criterion was written before anybody knew about this, it says results must match to 10⁻¹², and it fails intermittently.

There are four honest responses and they are worth ranking.

Find out what the number is being compared with. A comparison against a stored answer is harder than a comparison against a second run — the stored answer was produced by one particular machine on one particular day, and matching it exactly requires reproducing that machine’s partitioning. A criterion phrased as any two runs agree is achievable with an order-independent reduction; one phrased as matches the reference requires the reference to have been produced the same way.

Measure the actual spread before renegotiating. Run the computation at several partition counts and report the observed ratio. On this problem it is 1.2 to 1.7 relative; on a problem whose reductions have κ = 10¹² it will be far larger. A number beats an argument, and the number is one afternoon’s work.

Then either widen the criterion to a relative band above the measured spread, or buy the bitwise property. Those are the only two, and which is right depends on whether anything downstream is a verdict. What a regression test can ask for measures how much room the first option leaves: on one problem the window between the machine’s own variation and the smallest defect worth catching is a factor of 1.42, which is not much.

And do not tighten the tolerance in the hope that it helps. It is the response that costs the most and delivers the least, and this essay is the measurement that says so. It is also the response that is hardest to argue against without one, because it appears to work: the failures become rarer, since a smaller absolute disagreement crosses a fixed absolute threshold less often, and a test that fails once a month instead of once a day reads as progress rather than as the same defect sampled less.

A note on the shape of the bars

One detail of the figure repays attention. The bars do not merely fall — their tops fall by the same factor as their bottoms, which is what makes the ratio constant. If the disagreement had an absolute floor, the bars would stop shrinking at some level and the ratio would blow up as the accuracy improved past it.

There is no such floor within the range measured, and the reason is that everything in sight scales with u: the reductions’ disagreement, the residual at which the method stops, and the forward error it corresponds to. A floor would appear if the tolerance were pushed to where the residual is comparable with the arithmetic’s own noise, which for this problem is around 10⁻¹⁵ — past the range of the sweep, and past the range at which a stopping test means anything, which is what a tolerance actually obtains’s subject.

The three knobs, and what each one moves

Since this field has now tested all three, a table is the honest summary. Each row is a knob a caller can turn, and the columns are what happens to the two properties that keep getting confused.

knob absolute disagreement ratio between runs distinct answers
tighter tolerance falls with the accuracy unchanged unchanged
higher precision falls by 2⁻ᵖ unchanged unchanged
pinned thread count unchanged within a machine unchanged between machines falls, not to one
order-independent reduction unchanged zero one

Only the last row changes the property in the last column, and it is the only row that does not improve the accuracy at all. That is the field’s whole shape in four lines: the instruments that make an answer better and the instrument that makes two answers the same are disjoint sets.

The third row is worth its own note because it is the one people actually do. Pinning the thread count removes one of the three sources in the same program, twice — the number of pieces — and leaves the other two, the split points a scheduler chooses and the order the partial sums arrive in. Measured, that takes the count of distinct answers down and not to one, which is exactly the outcome that sustains the belief that the problem is nearly fixed.

One line

A tolerance decides how accurate an answer is and has no effect at all on how many answers there are, which are two properties the same word has been carrying for as long as anybody has been writing acceptance criteria.

At other settings

Asking for more accuracy buys accuracy, and buys no agreement at allConjugate gradients on a 80 × 80 matrix with κ = 10⁴, solved at four tolerances, seven partition counts each. The bar at each tolerance spans the smallest and largest forward error the seven runs produced. The bars fall by 5.4·10⁶ across the sweep, which is the tolerance doing exactly what it is for. The ratio between the top and bottom of each bar is 1.04, 4.83, 14.83, 1.36, which does not fall with the tolerance. So the disagreement between machines is not a residue of an insufficiently converged answer that a tighter tolerance would remove; it is proportional to whatever accuracy was reached, and the runs stay a fixed factor apart all the way down.-12-10-8-610⁻¹³10⁻¹¹10⁻⁹10⁻⁷log₁₀ of the residual tolerance asked forforward error of the answer1.04×4.83×14.83×1.36×the band does not closetolerances swept4runs at each7accuracy gained5.4·10⁶ratio at 10⁻⁶1ratio at 10⁻¹²1.4the bars falland they keep their height
Fig. 2 At eighty unknowns the ratios scatter much more widely, from a smaller sample of a noisier quantity.
Asking for more accuracy buys accuracy, and buys no agreement at allConjugate gradients on a 300 × 300 matrix with κ = 10⁴, solved at four tolerances, seven partition counts each. The bar at each tolerance spans the smallest and largest forward error the seven runs produced. The bars fall by 8.7·10⁵ across the sweep, which is the tolerance doing exactly what it is for. The ratio between the top and bottom of each bar is 1.27, 1.36, 1.11, 1.38, which does not fall with the tolerance. So the disagreement between machines is not a residue of an insufficiently converged answer that a tighter tolerance would remove; it is proportional to whatever accuracy was reached, and the runs stay a fixed factor apart all the way down.-12-10-8-610⁻¹³10⁻¹¹10⁻⁹10⁻⁷log₁₀ of the residual tolerance asked forforward error of the answer1.27×1.36×1.11×1.38×the band does not closetolerances swept4runs at each7accuracy gained8.7·10⁵ratio at 10⁻⁶1.3ratio at 10⁻¹²1.4the bars falland they keep their height
Fig. 3 And at three hundred, where more iterations mean more reductions to disagree about.
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. 4 The spread this essay is testing a repair for.
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. 5 And what can be asked for instead.
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. 6 The intervention that does buy agreement, priced.
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. 7 And the property it buys, tested bitwise.
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. 8 The disagreement at the bottom of the chain.
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. 9 Its size, which carries one factor of u and therefore moves with every knob at once.
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. 10 What a stopping value obtains, in the essay that measured it before the machine was in the picture.
Conjugate gradients at κ = 104, against the bound κ permitsA semi-logarithmic plot of the relative A-norm error against iteration count. The measured curve falls below a smooth dashed curve showing the classical condition-number bound.0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417
Fig. 11 The convergence whose crossing point the tolerance sets.
The residual conjugate gradients reports and the residual of the vector it holds, along one runA 50×50 diagonal matrix with one eigenvalue at 10⁻¹⁰ and the rest between 1 and 50, so every product with it is exact to a rounding and nothing below can be blamed on the operator. The lower curve is r ← r − αAp, updated by the recurrence, which is what a stopping test reads. The upper curve is ‖b − Ax‖/‖b‖, recomputed from the iterate. They start as the same vector and end 8.6·10⁸ apart, with the reported one at 6.89·10⁻²¹ — below the unit roundoff of 1.11·10⁻¹⁶, which is the shortest proof available that it is not the residual of anything.0122436486072849610⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best6.9·10⁻²¹the answer's, at its best5.9·10⁻¹²unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10⁹iterations drawn96the recurrence remembers every roundingand the stopping test is written in it
Fig. 12 The two residuals a stopping test chooses between, before any machine is involved.
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. 13 The other knob, which does the same thing this one does.
Iterative refinement from a 24-bit factorisation, κ = 10⁴A semi-logarithmic plot of forward error against refinement step. One curve falls steeply to the level of a double-precision solve; the other is nearly flat.012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹refinement step‖x − x*‖ / ‖x*‖a full double-precision solveresidual in24-bitresidual indoubleone argument apartκ·u of the factorisation6·10⁻⁴double residual, final3.2·10⁻¹³same-precision, final1.3·10⁻⁴30×30, κ = 10⁴, same factors in both runsidentical cost
Fig. 14 Buying accuracy back at the end, which also buys no agreement.
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. 15 The number a relative requirement should be written against.
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. 16 The distinction between what the test reads and what the requirement is about.
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. 17 A verdict no tolerance can stabilise.
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. 18 And another.
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. 19 A disagreement with no tolerance in it at all.
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. 20 The standing that makes a forward error a measurement rather than an estimate.
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. 21 The floor the bars would meet if the sweep went far enough.
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. 22 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. 23 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. 24 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. 25 A bound measured against what it bounds.
Loss of orthogonality against condition number, in binary64A log–log plot of the norm of Q-transpose-Q minus the identity against condition number. Classical Gram–Schmidt rises steeply, modified Gram–Schmidt rises gently, and Householder is flat.110²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classicalmodifiedHouseholderκ²uκu8×8, eight seeds per κ, binary64all three reconstruct A
Fig. 26 Orthogonality as a number, computed from an inner product at the wrong end of the census.
QᵀQ from classical Gram–Schmidt and from Householder on the 8×8 Hilbert matrixTwo eight-by-eight tables of QᵀQ. The upper one has ones on the diagonal and entries as large as one off it; the lower one is the identity to three decimal places everywhere.A = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle
Fig. 27 The same quantity entry by entry, where the cancellation that makes it fragile is visible.

Where this goes next

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 A stopping test is a race One matrix, one right-hand side, one tolerance, thirteen partition counts — and eleven different iteration counts between 674 and 690. Every run converged, every answer is right to the accuracy asked for, and what differs is the bill. 0123456789101110⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110³Newton steptolerance asked for, and iterations paiditerations paidtolerance asked forouter residualthe adaptive policy, step by stepNewton steps10inner iterations, total1009first step's cost1last step's cost271final outer residual3.4·10⁻¹¹the rule reads the last two residualsand asks for nothing it cannot use A tolerance that reads its own residual The cheapest constant forcing term costs 980 inner iterations and arrives with a hundred times the forward error of the dearest, which costs 9,358. A rule that sets each step's tolerance from the ratio of the last two residuals costs 1,009 and arrives with neither problem — and it is not a constant, so it does not appear on the curve the constants are compared on. 02468101210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²binary64 operations per elementrelative errorruntime orderfixed partitionpre-rounded, one passpre-rounded, two passescompensated (Kahan)exact accumulationwhat determinism costsreproducible policies3cheapest, ops/element3its error0.0024two passes10⁻¹⁷ordinary, 8 pieces3.6·10⁻⁷nothing is cheap and exactthe middle of the figure is the answer Accuracy and agreement are different properties The most accurate policy on this site's summation figure returns 119 different answers, and the one that returns a single answer is four orders less accurate. Neither property implies the other, and the vocabulary has one word for both.

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.

Bitwise reproducibilityConjugate gradientsConvergence rateReduction orderRun-to-run variationStopping criterionToleranceUnit roundoff