The answer that depends on the machine

Where the disagreement comes from

The error of a reduction is a walk whose step length is the spacing of the running total, not of the answer. That one sentence predicts the size of the disagreement to a factor of two, explains why dividing the work makes it smaller, and explains why the value cannot be predicted at all.

Worth reading first: The same program, twice · The order they are added in · Cancellation takes the answer, not a digit.

Two essays have now measured a disagreement without explaining it. This one explains it, in a sentence that turns out to carry everything the field needs:

A rounding is committed at every partial sum, and its size is the spacing of that partial sum, not of the answer.

Everything follows. The size of the disagreement, the fact that dividing the work reduces it, the fact that the value cannot be predicted, and the reason a vector of positive numbers shows nothing — all four are consequences of where the roundings are committed rather than of how many there are.

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. 1 The running error of two reductions of the same 2,048 numbers, against how many terms have been consumed. Both are exact statements about the arithmetic: the truth is carried in a compensated pair so that the instrument does not round.

Reading the picture

Two walks. They separate within the first handful of terms and never meet again, ending at 2.4·10⁻¹⁴ and 7.6·10⁻¹⁴ — a factor of three apart, on the same numbers, with the same algorithm.

Neither is a curve with noise on it. Each is a path: at every step the arithmetic commits one rounding of size at most half the spacing at the current running total, and adds it to what has accumulated. Sometimes the new rounding cancels part of what is there and the curve dips. Nothing about the path is deterministic in any useful sense, because it depends on the low bits of every partial sum.

The horizontal axis is where the mechanism is visible. This vector’s running total climbs to 52 on its way to an answer of 1.6·10⁻⁵. That is what a summation condition number of 10⁸ means in concrete terms: the arithmetic spends the entire computation handling numbers seven orders of magnitude larger than the one it returns, and rounds each of them at its scale rather than at the answer’s.

A rounding committed at a partial sum of 52 is worth about 5.8·10⁻¹⁵ in absolute terms. Relative to the answer of 1.6·10⁻⁵, it is worth 3.6·10⁻¹⁰. One rounding, in the middle of the loop, is already ten orders larger than the last bit of the result. There are 2,048 of them.

Why the answer’s own size is irrelevant

This is the point at which the field parts company with the ordinary reading of floating-point error, so it is worth stating twice.

The usual mental model is that each operation contributes a relative error of about u to the answer, so n operations contribute about nu. That model is right when the running total stays near the answer — which is what a sum of positive numbers does, and why the order they are added in can talk about a million reciprocals losing accuracy in a way that is proportional to how many of them there are.

When the terms cancel, the running total is nowhere near the answer, and the errors are committed at a scale the answer never sees. The right model is that each operation contributes about u·|running total| in absolute terms, and the answer’s smallness is what turns those absolutes into enormous relatives. That is the same mechanism cancellation takes the answer describes for one subtraction, applied n times to a running total instead of once to a pair.

The consequence is that the number of terms almost drops out. Fix κ and vary the length: the relative error of a single-accumulator reduction is 4.2·10⁻⁹ at 128 terms and 3.1·10⁻⁹ at 8,192, a fitted exponent in n of −0.083. No power of n at all, where the classical bound carries a full one — which a bound every answer satisfies works out in detail, and which comes out of exactly the two square roots this essay is about.

The two square roots

Written out once, because the cancellation of the two is what makes the whole field’s numbers small and stable.

The first √n is the roundings. They do not conspire. n roundings of similar size and independent sign accumulate like a random walk, reaching about √n times one of them rather than n times one of them. This is what the probabilistic error analysis of the last decade formalises, and it is why the deterministic bound’s factor of n is pessimistic by √n on any real data.

The second √n is the data. A vector of n mixed-sign terms with a fixed Σ|xᵢ| has entries of typical size Σ|xᵢ|/n, so its partial sums are themselves a walk and reach about √n · Σ|xᵢ|/n = Σ|xᵢ|/√n. The roundings are committed at that scale, which shrinks as the vector gets longer.

Multiply them: √n roundings, each of size u·Σ|xᵢ|/√n, gives u·Σ|xᵢ|. The length has gone. What is left is a quantity proportional to the sum of magnitudes, which divided by the answer is exactly κu — the constant this field reports, measured at 0.25 and derived here at 1.

The factor of four between the derivation and the measurement is honest and worth keeping: the derivation treats the roundings as independent and of full size, and they are neither. What the derivation buys is not the constant, it is the shape — no n, one κ, one u — and the shape is what survived ten decades of measurement at a spread of 1.02.

Why more workers is more accurate

Now the consequence that surprises everybody, including the author of the refusal this essay publishes.

Each accumulator in a p-way reduction sees n/p terms, so its own walk reaches about Σ|xᵢ|/(p√(n/p)) · √(n/p) — shorter by √p than the single accumulator’s. The roundings committed inside it are smaller in proportion. The final combination of p partial sums adds log₂p more roundings, at the scale of the answer, which is negligible when the answer is small.

So a divided reduction should be more accurate by about √p. Measured over eight vectors at 2,048 terms:

pieces mean absolute error mean peak partial sum
1 4.78·10⁻¹⁴ 52.1
4 3.74·10⁻¹⁴ 44.6
16 2.00·10⁻¹⁴ 37.4
64 8.15·10⁻¹⁵ 33.8

A factor of 5.9 from one accumulator to sixty-four, where √64 = 8. The prediction is right in direction and right within 25% in size, and the peak partial sum falls alongside it, which is the mechanism showing its work.

This is the same fact a reduction that changes the order measures in the communication field, where a tall-skinny QR computed as a tree of independent block factorisations returns a better backward error than the sweep it replaces — 1.65·10⁻¹⁵ against 9.95·10⁻¹⁵. That essay reads the result as a statement about communication being free of a stability cost. Read from here it is the same √p: the tree’s local factorisations each see less of the matrix.

So the folk claim that parallel arithmetic is less accurate is not merely unproven, it is backwards. What parallel arithmetic is, is different — and different is the whole problem, because the consumer of the number is a comparison rather than a reader.

Why the value is unpredictable

Everything above is a statement about the scale of the walk. None of it says where the walk ends.

That is not a gap in the analysis; it is the nature of the object. A random walk’s end point is not a function of its step sizes. Two walks with identical step-size distributions end in different places, and the difference between their end points has the same scale as either of them.

The practical consequence is the asymmetry this field lives with:

  • the size of the disagreement is computable in advance, from data the caller has, at a cost of one pass;
  • the answer any particular machine will return is not computable at all, short of running that machine.

Numerical analysis is usually in the opposite position — a worst case that is rarely attained, and an achieved error that is stable and reproducible. This is the one place on the site where the statistics are tight and the individual outcome is a coin flip.

Why sorting changes everything, and why nobody can

The derivation above assumed the partial sums wander. Break that assumption and the mechanism gives a different answer, which is the strongest evidence that the mechanism is the right one.

Sort a vector by magnitude and the partial sums stop wandering: they climb monotonically, so the walk is a line and the errors accumulate in a line too. That is why sorting matters at all — the order they are added in measures a factor of 470 from sorting a million reciprocals — and it is why the factor of n comes back for sorted data.

For a positive vector, sorting ascending is the classical repair. For a cancelling vector it is worse than useless: sorting by magnitude puts all the positives before all the negatives, drives the running total to its maximum possible value, and maximises every rounding on the way. The best ordering for a cancelling sum is one that keeps the running total small, which means interleaving the signs — and that is a decision that has to be made about the data, by somebody who has seen it.

Which is precisely what a parallel runtime cannot do. It divides the vector by index, because index is all it has. The ordering that would help is one no scheduler is in a position to choose, and the orderings it does choose are drawn from the family whose end points are unpredictable.

The order of the data, which the runtime cannot see either

There is a second ordering in play and it is easy to miss, because the field has spent three essays talking about the order the arithmetic combines things in. The terms also arrive in an order, and that one is the caller’s.

Take one vector and present it three ways to the same seven partitionings:

the terms, arranged peak partial sum error at one accumulator spread across seven reductions
shuffled 46.0 6.3·10⁻¹⁴ 5.7·10⁻⁹
sorted by magnitude 67.4 5.4·10⁻¹⁴ 5.2·10⁻⁹
all positives, then all negatives 838.9 4.2·10⁻¹³ 4.8·10⁻⁸

The third row is the one to read. Grouping by sign drives the running total to 839 — eighteen times the shuffled arrangement’s peak — because the sum climbs the whole way up and then comes the whole way down. Every rounding on that journey is committed eighteen times too large, and both the error and the run-to-run spread move by an order of magnitude.

Nothing about the arithmetic changed. The terms are the same multiset and the algorithm is the same algorithm; what changed is a property of the data’s layout, which in real code is decided by whatever assembled the vector — a sparse format, a sort somewhere upstream, the order the mesh was numbered in.

Two consequences worth carrying:

A benchmark on shuffled data understates the effect, which is the same warning the vector that hides it makes about positive data, one level down. Grouped signs are not exotic; they are what a physical problem with a source and a sink produces.

And this is why the ordering that would help cannot be chosen by the machine. Keeping the running total small means interleaving the signs, which requires having looked at the values. A runtime divides by index because index is all it has, so the one intervention that would shrink every rounding at once is available only to whoever built the vector — and available cheaply only to them, since a pass that inspects signs costs what the reduction costs.

What this predicts, checked

A mechanism earns its place by predicting something it was not built to explain. Three predictions, all checked every time this page is produced:

A positive vector should show almost nothing. Its partial sums climb to the answer and stay there, so every rounding is at the answer’s scale and the relative error is a few nu. Measured: 1.6·10⁻¹⁵ across seven partitionings, against 4.9·10⁻⁷ on a cancelling vector at κ = 10¹⁰.

Compensated summation should be accurate and still irreproducible. Kahan’s loop removes most of each rounding but not all, and what is left is still committed at the running total’s scale. So the error should fall a long way and the count of distinct answers should not fall to one. Measured over four hundred permutations: 119 distinct answers, against 303 for a plain loop — three times better and not remotely reproducible.

The disagreement should be immune to a wider float. Every quantity in the derivation carries one factor of u, so doubling the precision moves everything down together and changes nothing about the shape. That is the prediction most often assumed away — use doubles is the standard advice, and it buys accuracy while buying no agreement, exactly as the tolerance that buys no agreement measures for a tolerance rather than for a format.

What it does not explain

Two things, named because a mechanism that explained everything would be suspicious.

It does not explain the constant. The derivation gives u·Σ|xᵢ| and the measurement gives about a quarter of it. The gap is the independence assumption — consecutive roundings are correlated, because a partial sum that has just been rounded up is more likely to round down next — and recovering the constant would need a model of that correlation. The site’s habit is to report the measured number and to say which part of it was derived, which is what the direction the error leans does for a rounding mode: the size of one rounding is set by the precision, and how ten thousand of them combine is set by something else, with fitted exponents of 0.47 and 1.01 for two modes on identical data.

It does not explain the tail. The measurements here are means over seeds; the worst of twenty runs is about twice the mean and the analysis says nothing about how much worse a rare one can be. For a walk the answer is a Gaussian tail, and the honest version of that claim needs a distribution rather than a mean — which is the shape of statement the randomised field already makes, and a bound that holds with probability is where this site keeps its vocabulary for it.

The one line

The error of a reduction is a walk over its own partial sums. What a machine changes is which walk, the length of the vector nearly cancels out of the answer, dividing the work shortens each walk and improves the average, and none of that says where any particular walk ends.

At other settings

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 2 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 1.13·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 accumulator2 piecesthe steps are the partial sumspeak partial sum32the answer1.6·10⁻⁵error, one piece7.6·10⁻¹⁴error, 2 pieces1.1·10⁻¹³mean at p = 648.1·10⁻¹⁵the walk sets the sizeand nothing sets the value
Fig. 2 Two pieces: the shortest division there is, and already a different walk.
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 128 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 1.86·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 accumulator128 piecesthe steps are the partial sumspeak partial sum32the answer1.6·10⁻⁵error, one piece7.6·10⁻¹⁴error, 128 pieces1.9·10⁻¹⁴mean at p = 648.1·10⁻¹⁵the walk sets the sizeand nothing sets the value
Fig. 3 And a hundred and twenty-eight, where no accumulator sees more than sixteen terms.
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. 4 The end points of twenty-six of these walks.
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. 5 Their scale, which the mechanism predicts, against the bound, which it does not.
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. 6 How far the running totals climb in the reductions this site actually computes.
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 The prediction about compensated summation, checked: accurate, and 119 answers.
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. 8 And the policies whose walks are not walks at all.
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. 9 The accuracy this mechanism is usually described in terms of.
Relative error of two algebraically identical expressions for (1 − cos x)/x², in binary32A log–log plot of relative error against x. The expression written as it reads loses accuracy as x falls and is eventually wrong in every digit; the rearranged form stays at rounding level.10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻⁹10⁻⁶10⁻³1xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at allbinary32 throughoutone function, two spellings · zero below 3.5·10⁻⁴
Fig. 10 One subtraction, doing to a pair what a reduction does to a running total.
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. 11 The amplifier, in its matrix form.
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. 12 A residual and an error, which is the distinction a walk’s end point sits inside.
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. 13 The instrument every measurement here is against.
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. 14 A loss of orthogonality, which is a walk of exactly this kind in a matrix.
‖QᵀQ − I‖ at every step of a Lanczos run, n = 40The departure from orthogonality of the Lanczos basis, against the step, on a logarithmic vertical axis. It sits at the level of rounding for 13 steps and then climbs by a factor of about seventeen a step for ten steps running, twelve orders of magnitude, before saturating. The dashed line marks the step at which the first Ritz value converged, which is step 11.051015202530354010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹step‖QᵀQ − I‖first Ritz value converges√uno reorthogonalisationfullthe two events are one eventstep orthogonality crosses √u14step the first Ritz value converges11consecutive steps at ≥5× growth10a drift would grow like a square rootthis grows like a geometric series and then stops
Fig. 15 And one that stops being a walk when the recurrence loses its footing.
Growth factor under partial pivoting to n = 40: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8
Fig. 16 Another quantity whose worst case carries a power of n it never attains.
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. 17 The looseness measured rather than described.
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 What an unpredictable end point costs when a comparison reads it.
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 when it decides an integer.
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. 20 A change in the walk’s shape at a length chosen by a library.
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. 21 The knob that moves every quantity in the derivation together, and therefore changes none of its shape.
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. 22 Buying accuracy back after the fact, which is the other thing a wider float is usually asked for.
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. 23 The loss of orthogonality that amplifies a small perturbation over hundreds of steps.
The spectrum of C⁻¹T against T's own, n = 64, ρ = 0.8Two sorted spectra on a linear vertical axis. The unpreconditioned Toeplitz matrix spreads its eigenvalues across its whole range; the preconditioned operator puts 97% of them within a tenth of 1, with 0 below zero.08162432404856640123456789index, sortedeigenvalueT aloneC⁻¹Teigenvalues, not singular valuesfraction within 0.1 of 10.97eigenvalues below zero0λ_min of the preconditioner0.11the dashed line is 1, where the cluster formsthe solid line is zero
Fig. 24 A spectrum whose clustering decides a rate, computed from sums of exactly this kind.
Two condition numbers of one 8×8 system, as its rows are put into different unitsFour curves against the spread of the row units, in decades. κ_∞ of the scaled matrix rises from 9.83 to 1.9·10¹⁰ while the componentwise condition number stays at 6.98 throughout — the same system, the same solution, and one of the two numbers is a fact about the units. Hilbert's two numbers are drawn flat beside them at 3.4·10¹⁰ and 1.2·10¹⁰: a matrix whose sensitivity no scaling repairs.0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it
Fig. 25 Units, which are an author of error wearing its usual disguise.
Backward and forward error of two routes to x, on 30×30 systems across eight decades of κBoth routes start from the same LU factorisation. The LU solve's backward error is 2.2·10⁻¹⁷ at every conditioning — a flat line at the unit roundoff. Multiplying by the explicitly formed inverse gives 4.5·10⁻⁵ at κ = 10¹⁴, a slope of 0.94 against κ. The two forward errors, drawn above them, are 2.8·10⁻⁴ and 0.015 — within a factor of 54, which is why the difference between the two methods is invisible to anyone measuring how wrong the answer is.10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A)relative errorforward, A⁻¹bforward, A\bbackward, A⁻¹bbackward, A\bat κ = 10¹⁴η, LU solve2.2·10⁻¹⁷η, via the inverse4.5·10⁻⁵forward, LU solve2.8·10⁻⁴forward, via the inverse0.015one factorisation, two ways to use itand one of them forfeits the backward error
Fig. 26 Two routes to one answer, differing by orders in the backward error.
What 12 members of a sequence cost, by what changes between themEvery bar is 12 solves of an 120×120 system, costed in multiplications with a dense Cholesky at n³/3 and a solve at 4n². When nothing changes the answer is reused and the cost is one factorisation. When only the right-hand side changes, one factorisation serves every member — 16.5 per cent of the independent cost, with every solve still at a relative residual of 1.54·10⁻¹³. When the matrix drifts, the factorisation is reusable for a while and the while has to be measured: 4 factorisations for 12 members. When everything changes, nothing carries.nothing8.3%the right-hand side16.5%the matrix, slowly85.9%everything100.0%what can be reused: the answerwhat can be reused: the factorisationwhat can be reused: the factorisation, for a whilewhat can be reused: nothingshare of the cost of a sequence that shares nothingwhat 12 members cost in factorisationssame: factorisations1rhs: factorisations1drift: factorisations4independent: factorisations12what changes between the membersdecides what may be carried
Fig. 27 The cost of solving the same problem again, where a spread in step counts is paid repeatedly.
160 single-probe estimates of one 40×40 trace, from the two standard probe distributionsTwo clouds of points, each one estimate of the trace from one random vector. Both are centred on the true value of 171.14 — the estimator has no bias at all — and the dashed bands are the standard deviation each distribution is known in closed form to have: 59.14 for the ±1 probe and 72.01 for the normal one. Measured over these draws they come out 56.68 and 68.57.027548110813544119.306194.613269.919345.226420.532drawzᵀAz from one probenormal±1two routes to one spreadthe trace171±1 spread, predicted59±1 spread, measured57normal ÷ ±11.2no bias in either cloudand one of them is narrower for a reason
Fig. 28 A spread that comes from a deliberate random draw rather than from a scheduler.

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.

Error accumulationPairwise summationRandom walkReduction orderReduction treeRun-to-run variationSummation condition numberUnit roundoff