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

What determinism costs

Six ways to add up a vector, priced in operations per element and in accuracy. Nothing sits in the bottom left of the figure — an answer that is the same on every machine costs between three and twelve operations where an answer that is not costs one.
33 min read 28 figures Measured, not assertedWhere the algebra stops being the arithmeticTwo routes to a number

Worth reading first: The sum that cannot be wrong · The same program, twice · The order they are added in.

The sum that cannot be wrong establishes that order-independent summation exists. The question that follows immediately is what it costs, and the answer has two halves that point in different directions.

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. 1 Six summation policies on one vector: binary64 operations per element against relative error. Filled markers are the policies whose answer does not depend on the order.

The first half of the answer is on the figure. The second half is not, and is the more important of the two.

Reading the figure

Six policies, one vector of 2,048 numbers with a summation condition number of about 10⁶.

Runtime order — one add per element, one accumulator, whatever ordering the runtime picks. This is the baseline and it is the fastest thing on the figure.

A fixed partition — the same loop with the thread count pinned. Identical cost, identical accuracy, and reproducible only across runs that agree about the partition, which is why its marker is open rather than filled.

Pre-rounded, one pass — three operations an element, order-independent, and the least accurate policy on the figure at 2.3·10⁻⁷.

Pre-rounded, two passes — seven operations, order-independent, and on this vector exactly the correctly rounded answer.

Compensated — four operations, accurate, and not reproducible, which is the finding of the previous essay.

Exact accumulation — about twelve operations and a list of partials, correctly rounded, and order-independent by definition.

Nothing sits in the bottom left. There is no policy that is both as cheap as the baseline and independent of the order, and there cannot be one: the baseline’s cheapness is exactly its freedom to round wherever the arithmetic lands, and removing that freedom is the whole of what the other policies do.

The operation count is the half that misleads

Now the half of the answer that is not on the figure, and it is the one this site has a field about.

The same arithmetic at a different price measures two eliminations performing 72,568 operations each, choosing the same pivots, returning a factorisation identical to the last bit — and moving 41,332 and 19,476 words between fast and slow memory. The operation count is identical and the cost is not, by a factor of two.

A reduction is the extreme case of that. It reads n numbers and performs n adds, so it is memory-bound on every machine built in the last thirty years: the loop spends its time waiting for data and the arithmetic unit is idle most of the way. Adding two more operations per element to a loop that is waiting for memory costs approximately nothing.

So the honest reading of the figure’s horizontal axis is not three times the cost. It is:

  • the one-pass pre-rounded policy is nearly free on a memory-bound reduction, because the arithmetic it adds fits in time the loop was already spending;
  • the two-pass policy costs a second pass over the data, which on a vector that does not fit in cache is a genuine doubling — not because of the seven operations, but because of the second traversal;
  • exact accumulation costs what its state costs. Forty words of live state per accumulator and a data-dependent inner loop is the description of something that does not vectorise, and a reduction that does not vectorise gives up a factor that has nothing to do with its operation count.

Which reorders the figure. On flops, the ordering is 1, 3, 4, 7, 12. On what a reduction actually spends, it is closer to: the baseline and one-pass pre-rounding are the same speed, two-pass pre-rounding is twice that, and exact accumulation is a different kind of loop.

The state column, which is the one to read

Every policy on the figure carries a second number that the horizontal axis does not show: how much live state it needs per accumulator.

policy operations state order-independent
runtime order 1 1 no
fixed partition 1 8 at one partition count
pre-rounded, one pass 3 1 yes
pre-rounded, two passes 7 2 yes
compensated 4 2 no
exact accumulation 12 ~40 yes

One accumulator vectorises trivially — the compiler keeps four or eight of them in registers, which is the same trick the length that changes the kernel is about. Two is still fine. Forty is not, and the pattern of writes into them depends on the values, which defeats the vectoriser for a second reason.

That is the whole argument for the pre-rounded policies over the exact one, and it is an argument about registers rather than about arithmetic.

What each one is actually for

Prices are only useful next to purposes, so:

Runtime order is right whenever no comparison is downstream. A residual printed in a convergence log, a norm reported at the end, a diagnostic. Most reductions are in this category and this field is not an argument for changing them.

A fixed partition is the cheapest thing that helps and it is a partial answer. It removes the thread count as a variable and leaves the scheduler in place, so it fixes reproducibility within one machine’s configuration and not between two. It is worth doing because it is free, and worth not mistaking for a solution.

One-pass pre-rounding is for a reduction whose verdict matters and whose value does not. A rank threshold, a definiteness sign, a convergence test: the comparison needs the two machines to agree, and 10⁻⁷ of the largest element is far below any threshold worth setting. This is the case the policy was invented for and it is nearly free.

Two-pass pre-rounding is for a reduction whose value matters too, which is most residuals. Seven operations, two passes, correctly rounded on the vector measured here, and the same on every machine. If this field has one recommendation it is this one.

Compensation is for accuracy alone, on a machine that is not going to be compared with another one. It is the best accuracy-per-operation on the figure and it buys nothing this field is about.

Exact accumulation is for ground truth. It is what every measurement in this field is against, which is a use that justifies its cost completely: a reference implementation runs once and an inner loop runs 10⁹ times.

What it costs inside a solver, which is the number anybody wants

Operations per element is the right axis for a summation policy and the wrong axis for a decision. The decision is whether to spend it inside something, so here is the same price paid inside a conjugate gradient step.

One step is a matrix–vector product, three vector updates, and two inner products. Only the last of those is a reduction, so only the last is affected — and how much it matters depends entirely on how expensive the product is. Counting operations on a problem with n = 200:

the matrix one-pass two-pass compensated exact
dense, 200 per row 1.01× 1.03× 1.01× 1.05×
sparse, 27 per row 1.06× 1.19× 1.10× 1.35×
sparse, 5 per row 1.22× 1.67× 1.33× 2.22×

On a dense matrix the whole question is noise: making both reductions exactly accumulated costs five per cent of a step, and the answer stops depending on the machine. On a five-point stencil — the model problem this site’s iterative field is built on, at five nonzeros a row — exact accumulation doubles the step, and one-pass pre-rounding costs 22%.

That table is the practical content of the essay and its shape is worth stating: the cheaper the matrix, the dearer the reproducibility, because the reductions are a larger share of a step whose product is cheap. A solver on a dense matrix can have bitwise reproducibility for nothing. A solver on a sparse one is making a real choice.

And it points at the right intervention on the expensive side: one-pass pre-rounding at 22% on the sparse problem buys agreement in the stopping test, which is what a stopping test is a race shows costs 674 to 690 iterations — a spread of 2.4%. Spending 22% of a step to remove a 2.4% spread in the step count is a bad trade if the only goal is the count, and a good one if a verdict downstream needs two machines to agree.

What no summation policy buys

Three of this field’s results are outside the reach of everything on the hero figure, and a reader budgeting for reproducibility should know which they are before spending anything.

The kernel cutoff. A library that switches from one accumulator to four at a length threshold is changing algorithm, not ordering, and no policy applied to either kernel makes the two agree. The step in the length that changes the kernel is a property of a version number.

The contraction. Whether a multiply-add was fused is decided at compile time, and both forms are conforming. A reproducible summation over products that were themselves rounded differently is reproducibly wrong in two different ways. One multiply the compiler removed is that measurement.

The algorithm. Two libraries computing a QR by Householder and by a tall-skinny tree return different numbers, both backward stable, and the difference is not an ordering at all — it is a reduction that changes the order, where the tree is the better of the two.

So the recommendation this figure supports is narrower than make everything reproducible. It is: identify the reductions a verdict is read from, price them against the work around them, and spend there. The rest of the program’s variability has other causes and other repairs.

Where the second pass stops being free

One detail of the two-pass policy deserves its own paragraph, because it is the difference between a 3% cost and a 67% one in the table above, and it is not visible in either the operation count or the accuracy.

The second pass reads the vector again. If the vector fits in cache, that read is nearly free and the seven operations are what the policy costs. If it does not, the second pass is a second traversal of main memory, and on a memory-bound reduction that is a genuine doubling — the whole cost of the loop, paid twice, for a policy whose arithmetic is trivial.

The repair is to fuse the two passes: snap, accumulate, and accumulate the residue, all in one traversal, at the price of a second accumulator and a second δ that has to be fixed in advance rather than computed from the residues. Fixing it in advance costs accuracy — the residue’s largest element is not known until the first pass has run — so the fused version is a policy between the two, and which of the three is right depends on whether the data fits in cache.

That is the same arithmetic at a different price in miniature, and it is the reason this essay’s table is a table of operation counts with a warning attached rather than a table of times: the operation count is the same on every machine and the answer to which policy is not.

The comparison nobody runs

A closing observation about why the figure is unusual rather than standard.

Every library that ships a reduction has measured its speed. Many have measured its accuracy against a reference. Almost none publish the third column, because the third column requires running the same input twice under different partitionings and comparing bits — which is not a benchmark, it is a test, and it lives somewhere else entirely from a timing.

The consequence is that the trade-off in the hero figure is not visible from any documentation. A caller choosing between two routines can find out which is faster and which is more accurate, and cannot find out which returns the same answer twice. That is the gap what a regression test can ask for measures from the consumer’s side, where the same missing column turns into a tolerance nobody can choose.

And it is the reason this field’s recommendation is stated in terms of the caller’s own data rather than in terms of a library’s. κ of the caller’s sum is computable in one pass; which reduction the library performed is not discoverable at all.

A note on the axis that is missing

The figure has cost and accuracy on it and not the third property, because the third property is not a number: a policy either returns the same bits or it does not. That is why the reproducible ones are drawn as filled markers rather than as a third axis, and it is worth saying plainly, because the temptation to make it continuous is strong and produces nonsense.

How reproducible is not a quantity. The count of distinct answers over four hundred permutations — 303, 72, 119, 1, 1 in the sum that cannot be wrong — looks like one, and it is not: it is a property of the sample rather than of the policy, and a longer sample moves every number above one and none of the ones equal to it. The only honest reading of that column is binary, and the only honest way to draw it is a mark.

The same applies to the spread in ulps that the earlier essays report. It is a useful size, it predicts what a tolerance has to clear, and it is not a measure of reproducibility — a policy with a spread of two ulps is exactly as irreproducible as one with a spread of two million, and it will fail a bitwise test just as reliably. What the small spread buys is that a tolerant test can be written, which is the whole subject of what a regression suite can ask for.

At other settings

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 106. 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 1.91·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⁻¹¹binary64 operations per elementrelative errorruntime orderfixed partitionpre-rounded, one passpre-rounded, two passescompensated (Kahan)exact accumulationwhat determinism costsreproducible policies3cheapest, ops/element3its error1.9·10⁻¹¹two passes10⁻¹⁷ordinary, 8 pieces2.6·10⁻¹⁵nothing is cheap and exactthe middle of the figure is the answer
Fig. 2 At κ = 100 the ordering of the six barely changes and every error is small, which is the case in which none of this is worth doing.
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 0.00238 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⁻⁸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
Fig. 3 And at 10¹⁰, where the one-pass policy’s fixed error stops being the worst on the figure.
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. 4 The property being bought, 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. 5 And the property being avoided.
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. 6 How large the thing being bought is, as a function of the data.
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. 7 Which of this site’s own reductions would justify the spend.
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. 8 What the policies remove: a walk with a rounding at every step.
The accuracy of a dot product, either side of a length nobody in the program choseMean relative error of a dot product of two vectors whose entries are drawn identically at every length, from 58 to 70 terms, over 40 draws each. Below 64 the kernel accumulates into one register; at 64 and above it uses four independent accumulators and combines them, which is how a tuned library uses a vector unit. The error steps down by 1.57× at the cutoff and is flat either side of it. Nothing about the problem changes there: the vectors are drawn from the same distribution, the arithmetic is the same precision, and the cutoff is a constant in somebody else's source file. A user whose problem grows across it sees the answer move, and there is nothing in their program to look at.5860626466687010⁻¹⁷10⁻¹⁶terms in the dot productmean relative errorthe kernel changes herea constant in a libraryone accumulator3.2·10⁻¹⁷four accumulators2.1·10⁻¹⁷step at the cutoff1.6cutoff64the problem did not changethe loop did
Fig. 9 A cost decision inside a library that changes an answer, taken for exactly the reasons this essay’s horizontal axis understates.
A 6.25-bit block format against 8-bit E4M3, read two waysFour curves of relative error against the spread of the data. Two, for the block format, rise steeply when read as a median and gently when read as a norm; the other two are nearly flat.048121610⁻²10⁻¹1octaves of spread within a blockrelative errorblock, medianblock, 2-normE4M3, medianE4M3, 2-normthe same data, two readingsblock format, bits a value6.3norm across 16 octaves1.8median across 16 octaves62entries deleted at the wide end383640 values, blocks of 32a norm is dominated by what a block format keeps
Fig. 10 Another accumulator whose width is the whole design.
Communication for a 512×12 factorisation on 16 processorsTwo counts for the same factorisation. Rounds on the critical path: 48 for the column sweep, 4 for the reduction tree, 4 for Cholesky QR — a factor of n between the first and the other two. Words sent: 1170, 1170 and 2160 — the sweep and the tree send the same number, and the method with the fewest rounds sends the most.rounds on the critical pathHouseholder sweep48reduction tree4Cholesky QR4words sentHouseholder sweep1170reduction tree1170Cholesky QR2160two counts, two rankingsrounds, sweep ÷ tree12words, Cholesky ÷ tree1.8arithmetic, tree ÷ sweep1.5the rounds separate the threeand the words do not
Fig. 11 Why a reduction is divided at all, counted in messages rather than in operations.
Multiplications in a hierarchical solve and in a dense factorisation, and where they crossBoth counted rather than estimated: the recursion carries a counter and reports what it actually did. The dense factorisation is n³/3 and rises at exactly three per doubling. The hierarchical solve rises at 2.13, 1.93, 1.74 — falling, because the cost is n log²n and its exponent is on its way to one. The two cross between n = 64 and n = 128: below it the format is the more expensive way to solve the system, at 1.48 times the dense count, and at n = 512 it is 6.2 times cheaper. Every point returns an answer at a backward error of about 1.4·10⁻¹⁰, so the comparison is between two ways of getting the same thing.567891010⁴10⁵10⁶10⁷10⁸10⁹log₂ nmultiplicationsdense factorisation, n³⁄3the recursion, countedwhere the format starts payingratio at n = 641.5ratio at n = 5120.16exponent, first doubling2.1exponent, last doubling1.7backward error1.4·10⁻¹⁰cheaper is a sizenot a property
Fig. 12 A cost model whose currency is not arithmetic, in the field that made the argument first.
Relative error of three summation algorithms in binary32A log–log plot of relative error against the number of terms for naive, pairwise and compensated summation, each measured against the exactly rounded sum.10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸
Fig. 13 The accuracy axis of this figure, in the essay that measured it.
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. 14 What exact accumulation returns, and what it is for.
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. 15 Buying accuracy back at the end rather than during, which is the other shape of this trade.
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. 16 The knob that moves accuracy without moving agreement.
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. 17 The distinction that keeps a reproducible answer from being a correct one.
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. 18 The comparison that makes the spend worth making.
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. 19 And the verdict that would stop moving.
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. 20 What the test suite gets in exchange.
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. 21 A reproducibility failure that costs nothing to fix and is not a summation policy at all.
What one symmetric permutation does to the storage, on a matrix it does not changeThe same 256 × 256 matrix, its rows and columns renumbered by one permutation and its columns by the same one. The condition number is 24.3948 either way, to eight digits; the Frobenius norm is 6139.964 either way, to twelve. Nothing a norm can see has moved. Clustered, the partition stores 27,008 numbers. Shuffled, the admissibility test finds no admissible pair anywhere — every cluster of a shuffled numbering spans the whole interval, so every q is infinite — and the format degenerates to dense storage exactly. The rule with no test to fail does worse than that: it compresses every off-diagonal block regardless, gets ranks up to 119 out of 128, and stores 118,208 numbers — 1.80 times the matrix it was compressing. A rank-119 factorisation of a 128-column block is a more expensive way to write down the block than the block.numbers stored, 256 × 256clustered, strong27,008clustered, weak24,064the dense matrix65,536shuffled, strong65,536shuffled, weak118,208the same matrix, twiceκ, clustered24κ, shuffled24‖A‖_F, clustered6140‖A‖_F, shuffled6140shuffled weak ⁄ dense1.8the compressibility is in the numberingand the numbering is not in the matrix
Fig. 22 A choice with no numerical content and a large cost, which is the shape this field’s decisions have.
How close Hager's estimate is to the true κ₁, over 200 seeded 8×8 matricesFive bars. The estimate is exactly the true condition number on 81% of the sample and inside ten per cent on 87%; the worst underestimate in the whole sample returns 38% of the truth. The last bar is the matrix built to defeat it, at 7.7%, well below anything the random sample reached.each bar is a percentage — of the sample, or of the true condition numberexactly right80.5%inside 10%87.0%inside a factor of 287.0%worst in the sample, ×10037.7%the constructed matrix, ×1007.7%usually exactexact share0.81worst of the sample0.38the constructed matrix0.077a routine that is right most of the timeand never wrong in the safe direction
Fig. 23 An estimator’s spread, which is a different uncertainty about one number.
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. 24 Two residuals for one iteration, of which a stopping test sees one.
What the solver reported, and what it achievedFinal relative residual and final relative error against the working precision, on a logarithmic vertical axis. The residual sits at the stopping tolerance at every precision — the method converged, by every test available to it. The error runs from 1.4·10⁻¹³ at 53 bits to 0.00597 at 8, tracking the unit roundoff, which is drawn beside it.614223038465410⁻¹³10⁻¹⁰10⁻⁷10⁻⁴working significand bitsrelative sizeerrorresidualuthe gap the solver cannot seeerror ÷ residual at 24 bits1.5·10⁶error ÷ residual at 16 bits2.3·10⁷error ÷ residual at 8 bits2.4·10¹⁰every convergence test passesand the answer is wrong in proportion to u
Fig. 25 A residual that cannot see what it is asked about.
Numbers stored per unknown, against the size of the matrix, at 10⁻⁸The dense matrix stores n numbers per unknown, which is the straight line through the origin and doubles whenever n does. The hierarchical representations store 54, 79, 106, 133 and 50, 70, 94, 120, adding about 26 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 26 and 23 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.55, 1.42, 1.33, falling towards one and never arriving, which is what a logarithm looks like from inside.56789100108216324432log₂ nnumbers stored per unknowndenseweak: every off-diagonal blockstrong: only the admissible onesa constant per doublingstrong, n = 5126.8·10⁴weak, n = 5126.1·10⁴dense, n = 5122.6·10⁵per doubling26‖A − A_H‖ ⁄ ‖A‖3.8·10⁻¹⁰the dense line doublesand the other two add a constant
Fig. 26 A storage cost that follows from a decision about a threshold.
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. 27 The convergence whose crossing point a tolerance sets.
The residual basis of conjugate gradients, at κ = 106A semi-logarithmic plot against iteration count showing the loss of orthogonality among the residual vectors rising while the relative residual falls.05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here
Fig. 28 The loss of orthogonality that amplifies a small perturbation over hundreds of steps.

Where this goes next

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 The sum that cannot be wrong Snap every addend to a common multiple before adding, and every partial sum is exact — so the order stops mattering, by construction rather than by luck. Four hundred permutations return one value where an ordinary reduction returns three hundred and three. 10²10³10⁴10⁵matrix size ncountoperations, bothwords, unblockedwords, blocked (b = 6)the answer does not move‖PA − LU‖/‖A‖, unblocked2.8·10⁻¹⁶‖PA − LU‖/‖A‖, blocked2.8·10⁻¹⁶difference between them0the dashed curve is both orderings' operation countthe solid pair is what they cost The same arithmetic at a different price A blocked and an unblocked elimination perform 72,568 operations each — the same operations, associated differently — choose the same pivots, and return a factorisation identical to the last bit: ‖PA − LU‖/‖A‖ = 4.487946226420872·10⁻¹⁶ in both. One of them moves 41,332 words between fast and slow memory and the other moves 19,476. 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 reproducibilityData movementExact accumulationFlop countKahan summationMemory hierarchyPre-rounded summationReproducible summation