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

A bound every answer satisfies

The classical bound on a summation error is correct, it covers all twenty-six answers one vector produced, and it is 7,932 times larger than the difference between them. A statement true of every ordering cannot say which ordering you got.
32 min read 29 figures Where the algebra stops being the arithmeticMeasured, not assertedWhose fault is it

Worth reading first: The same program, twice · The order they are added in.

The same program, twice ends on a set of twenty-one distinct answers to one sum, and the natural next move is to reach for the bound. There is one, it is in every textbook on the subject, and for a sum of n floating-point numbers evaluated in any order it says

|ŝ − s| ≤ γₙ₋₁ · Σ|xᵢ|,     γₙ = nu / (1 − nu).

The bound is correct. Every one of the twenty-six runs satisfies it — that is asserted on every build of this site rather than assumed, because a bound nobody has fed a counterexample is a bound nobody has tested. And it is nearly useless for the question the previous essay ends on, for a reason that is worth stating carefully, because it is the reason this field is not a subsection of the arithmetic one.

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. 1 The measured disagreement between seven partitioned reductions, the quantity κu, and the classical bound, over ten decades of cancellation.

Three curves. Two of them are the same curve times a constant, and the third — the only one that is published — is nearly four orders above both.

What the bound is a statement about

Read the right-hand side. It contains n, the number of terms, and Σ|xᵢ|, the sum of their magnitudes. It does not contain the order.

That is not an oversight. The bound is derived by taking the worst case over every rounding the computation might commit, and the worst case is the same whichever order the terms arrive in — so the bound is a single number covering the whole family of orderings at once. It is a statement about the family.

The question a reader has after the previous essay is about one member of it. This run returned 1.0000000019790605; will the next one? A statement true of every member cannot separate two members, and no amount of tightening changes that: a tighter bound over all orderings is still one number for all of them.

This is the same distinction three walks and one bound draws in the arithmetic field, from the other side. There, three unrelated computations each have a standard bound linear in what they accumulate against, and all three attain a half of it — the bound is loose by a factor that is itself a measurement. Here the looseness is worse and the mismatch is different in kind: the quantity a reader wants is not smaller than the bound, it is a different quantity.

How much looser, measured

On the hero figure’s default vector, with 2,048 terms:

quantity value at κ = 10⁸
the measured spread between seven reductions 5.6·10⁻⁹ relative
κ · u 2.2·10⁻⁸
the classical bound γₙ₋₁Σ|xᵢ| 4.5·10⁻⁵ relative

The bound is 7,932 times the disagreement. That factor is not constant either: the bound carries a factor of n and the disagreement does not, so at 256 terms it is about three orders and at 8,192 it is about four. A reader who used the bound as a proxy for the spread would be wrong by a factor that depends on the length of a vector they may not have chosen.

And the factor between the spread and κu is 0.253 to 0.258 across ten decades, a spread of 1.02. One of these two is a description of the disagreement and the other is an envelope containing it.

The number that does describe it

So the useful quantity is κu, times a constant near a quarter, and the reason is the mechanism where the disagreement comes from sets out: the error of a reduction is a walk whose steps are the spacings of its partial sums, the partial sums reach about Σ|xᵢ|/√n in magnitude on a vector with mixed signs, and the difference between two such walks is a random variable of that scale.

Three things make it worth reporting rather than the bound.

It is computable before anything runs. One pass over the data gives Σ|xᵢ|, and the answer’s own magnitude is what the computation is for. A library that returned κ beside a residual norm would cost one pass and would tell a caller how much of the number to believe.

It is an estimate rather than an envelope. The bound answers a question about the worst case that nobody is in. κu answers the question a caller has, which is how much two runs will differ.

It is honest about being an estimate. The constant in front moves between about 0.2 and 0.5 with the vector, and a short vector widens it further. That is a factor of two, against a factor of ten thousand for the alternative.

Why the bound cannot be repaired

A reasonable objection: derive a sharper bound, one that knows about the order. The literature has those — a bound for pairwise summation carries log₂n rather than n, and a bound for a p-way reduction carries n/p + log₂p. Do they answer the question?

They do not, and the reason is instructive. Those bounds are still worst cases over the roundings within a fixed ordering, so each is a single number covering everything that ordering might do. Distinguishing two orderings would need, for each, a bound and a matching lower bound tight enough that the two intervals do not overlap — and the intervals are [0, γ·Σ|xᵢ|], which all contain zero, so they overlap by construction.

Every summation error bound in the literature has zero as its lower end, because a sum can be exact by luck. That is exactly what makes a family of them unable to separate two members. The refusal in this essay’s refutation is that statement, mechanised: an assertion is handed the claim that the bound excludes one of two measured answers, and it must reject.

The factor of n is not merely pessimistic — it is absent

The bound’s looseness has a shape, and measuring it is more useful than noting it.

Hold the cancellation fixed at κ = 10⁸ and vary the length. The bound’s right-hand side carries γₙ₋₁ and therefore grows linearly with n: from 128 terms to 8,192 it rises by a factor of 64. The measured error does this:

terms mean relative error of a single-accumulator reduction
128 4.18·10⁻⁹
512 3.15·10⁻⁹
2,048 2.18·10⁻⁹
8,192 3.12·10⁻⁹

Twenty-four vectors at each length, and a fitted exponent in n of −0.083. Not a smaller power of n. Not √n, which is what the modern probabilistic bounds replace n with. No power of n at all.

The explanation is the walk again and it is worth writing out, because two square roots cancel and that is not obvious.

A vector of n mixed-sign terms with Σ|xᵢ| fixed has entries of typical size Σ|xᵢ|/n, so its partial sums are a random walk that reaches about √n · Σ|xᵢ|/n = Σ|xᵢ|/√n. Each of the n roundings is therefore about u·Σ|xᵢ|/√n, and n of them accumulate as a walk rather than in a line, giving √n · u · Σ|xᵢ|/√n = u · Σ|xᵢ|. The first √n is the improvement a probabilistic bound already knows about — roundings do not conspire — and the second is one nobody writes down, because it comes from the data being a walk rather than from the arithmetic.

So the ratio between the bound and the disagreement grows linearly with the length: 2.8·10³ at 256 terms, 8.1·10³ at 2,048, 1.5·10⁵ at 8,192. A reader using the bound as an estimate is not out by a constant. They are out by a factor proportional to a length they may not have chosen.

The honest caveat: this is a statement about vectors whose partial sums wander like a walk, which is what mixed signs and no structure produce. A vector sorted by magnitude, or one with a systematic sign pattern, has partial sums that go somewhere rather than wandering, and there the n comes back — which is exactly the case the order they are added in measures when it sorts a million reciprocals and buys a factor of 470.

The asymmetry this leaves

The result is a strange epistemic position, and naming it is most of the value of this essay.

The size of the disagreement is predictable. 0.25κu, computable in one pass, stable across ten decades.

Which answer arrives is not predictable at all. It is the end point of a walk, and knowing the step sizes gives the scale of the end point and nothing about its value. The hero figure of the previous essay is not a curve with noise on it — it is noise, and the shape it has is the shape that particular vector and that particular scheduler produced.

Numerical analysis is generally in the opposite position. A condition number is a worst case, the achieved error is usually far below it, and the useful skill is knowing when the worst case is reached. Here the statistics are tight and the individual outcome is a coin flip, which is the situation the randomised field is in and no other: a bound that holds with probability is where this site first met a guarantee that changes kind, and the parallel is close enough to be worth stating and not close enough to be the same thing. There the randomness is drawn deliberately, from a seed the program owns. Here it is drawn by the scheduler, and nothing records what it drew.

What a lower bound would need

There is one honest way to bracket a computed answer, and it is the interval field’s: a limit is not a computation and the essays around it compute in directed rounding, so every operation returns an interval containing the true value. Applied to a reduction, that gives an interval containing the exact sum whatever the order — and therefore containing every ordering’s answer.

That is the correct object, and its width is the honest report. It costs a factor of about two in arithmetic and it needs rounding-mode control, which is exactly the machinery the direction the error leans measures. What it does not do is name which answer arrived: it says the truth is in here and so is every run, which is the same shape of statement as the bound and a much tighter version of it.

For a reader who needs the interval, that is the tool. For a reader who needs two machines to agree, the interval is not the answer either — the answer is to make the computation order-independent, which is the sum that cannot be wrong.

What a library could print, and why it does not

The gap between the bound and the disagreement is not an argument for a new theorem. It is an argument for a different return value, and the change is three lines.

A routine that computes a sum or an inner product already touches every term. Accumulating Σ|xᵢ| alongside costs one add and one absolute value per element — under a factor of two on a memory-bound loop, and nothing at all on one that is compute-bound. Dividing by the answer gives κ. Multiplying by u gives an estimate of how much two machines will disagree, good to a factor of two.

Compare that with what the neighbouring field had to build. An estimate that can be fooled is about the condition number of a matrix, which cannot be computed cheaply at all: the honest route is a singular value decomposition costing more than the solve it is describing, so every library ships an estimator instead, and the estimator can be fooled by a matrix built to fool it. κ of a sum has none of those problems. It is exact, it is one pass, and no matrix can be constructed to defeat it, because there is nothing to defeat — it is a ratio of two numbers the loop already has.

So the reason it is not printed is not cost and not difficulty. It is that the quantity has not been asked for, because the question it answers — will the next run agree — is not one the existing vocabulary can express. The whole of the previous essay is that observation and this essay is its consequence: the bound answers the question the vocabulary has, and the vocabulary is missing a question.

There is one more thing such a routine could return and it is free: the number of pieces the reduction was actually performed in. That single integer is what separates two runs, it is known to the runtime and to nobody else, and a result that carried it would let a caller reproduce a number rather than merely observe that they cannot.

The reading for the rest of the site

This essay is the reason the field is separate rather than a section of the arithmetic one, and the argument compresses to three lines.

The arithmetic field’s question is how large is the error, and the bound answers it.

This field’s question is how large is the disagreement between two correct runs, and the bound does not answer it, because it is the same number for both of them.

The quantity that does answer it — κ of the sum, times u, times a constant near a quarter — is not in any library’s output, is one pass over the data away, and is a property of the caller’s own data rather than of anybody’s code.

One further consequence, for anybody writing the documentation rather than the code. The sentence this routine is backward stable is true of every run in the hero figure and is heard by readers as a promise about repeatability, which it has never been. The two statements a caller needs are different and both are short: the answer is the exact sum of a vector a rounding away from this one, and two runs of it may differ by about κu. The first is the guarantee the literature gives. The second is the one that decides whether a test suite can be written, and it is missing from every manual page this site has read.

At other settings

What the runs disagree by, and what the bound says they mightOver ten decades of summation condition number, on vectors of 256 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.090 to 0.127 across the whole sweep, a spread of 1.414 while both quantities move by ten orders. The bound sits 2008 to 2839 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.09spread ÷ κu, high0.13bound ÷ spread, low2008decades swept10the spread is computablethe bound cannot see the order
Fig. 2 At 256 terms the constant in front of κu is 1.4 times wider, because a short reduction commits few enough roundings that one of them is a share of the total.
What the runs disagree by, and what the bound says they mightOver ten decades of summation condition number, on vectors of 8192 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.049 to 0.056 across the whole sweep, a spread of 1.155 while both quantities move by ten orders. The bound sits 1.5·10⁵ to 1.7·10⁵ 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.049spread ÷ κu, high0.056bound ÷ spread, low1.5·10⁵decades swept10the spread is computablethe bound cannot see the order
Fig. 3 And at 8,192 the bound has moved up with n while the disagreement has not.
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 answers the bound covers, all twenty-six of them.
The same 4096 numbers, added up 26 ways, 15 different answersEach dot is one reduction of one vector of 4096 binary64 numbers whose summation condition number is 100: 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. 15 of the 26 runs returned distinct values, spanning 1.56·10⁻¹³ — 22 ulps of the answer. The dashed line is the classical bound γ₍ₙ₋₁₎Σ|xᵢ|, which every one of them satisfies with 2.47·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 answers15runs26spread, in ulps22κ of the sum100bound ÷ worst error2.5·10⁴nobody chose pand no answer is the answer
Fig. 5 And the same picture where κ is a hundred, which is where a bound and a spread agree about there being nothing to say.
The error of a reduction is a walk, and the number of pieces changes which walkRunning error of two reductions of the same 2048 numbers, against how many terms have been consumed: one accumulator, and 8 pieces combined at the end. Each is the distance between what the arithmetic has accumulated so far and the exact sum of the same terms, carried in a compensated pair so the instrument does not round. The two separate within the first few terms and never meet again, ending at 7.56·10⁻¹⁴ and 2.41·10⁻¹⁴. The mechanism is in the horizontal axis: a rounding is committed at every partial sum and its size is that partial sum's spacing, and this vector's running total reaches 32.1 on the way to an answer of 1.63·10⁻⁵. Dividing the work shortens each walk, so the mean error over eight vectors falls from 4.78·10⁻¹⁴ at one piece to 8.14·10⁻¹⁵ at sixty-four: more workers is more accurate, which is the opposite of what the word *parallel* suggests here.051210241536204810⁻¹⁷10⁻¹⁴terms consumederror accumulated so farone accumulator8 piecesthe steps are the partial sumspeak partial sum32the answer1.6·10⁻⁵error, one piece7.6·10⁻¹⁴error, 8 pieces2.4·10⁻¹⁴mean at p = 648.1·10⁻¹⁵the walk sets the sizeand nothing sets the value
Fig. 6 The walk the bound is a worst case over.
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 Where κ actually sits, for the reductions this site computes elsewhere.
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. 8 The count no bound predicts: how many distinct answers four hundred orderings return.
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. 9 And what it costs to make that count one.
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. 10 A bound and what it bounds, in the arithmetic field, where the same looseness was first measured.
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. 11 The error a bound on summation is about.
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. 12 The amplifier κ of a sum is the summation-shaped version of.
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. 13 The distinction every bound on this site is measured against.
Backward error, forward error and the condition number, with measured valuesTwo boxes at the top — the problem posed and the nearby problem the algorithm answered exactly — and two answers below them, with the distances between all four labelled by numbers from a Hilbert solve.the problem you posedA = H10b = A·(1, 2, …, 10)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 2.3·10⁻¹⁷the answer you wantedx = (1, 2, …, 10), exactlythe answer you gotx̂, wrong by 2.7·10⁻⁴ relativebackward error 2.3·10⁻¹⁷forward error 2.7·10⁻⁴κ = 1.6·10¹³κ · η = 3.6·10⁻⁴, and the measured forward error is 2.7·10⁻⁴.The algorithm is not at fault. The problem is.H10, LU with partial pivotingresidual and error differ
Fig. 14 Two authors of error and a bound over their product.
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. 15 The exactly known answer that makes a measured error a measurement.
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. 16 Where the sum of magnitudes and the magnitude of the sum part company.
The width of a guaranteed bracket, and the error it guarantees, at 24 bitsTwo rising curves on logarithmic axes. The upper one, the bracket width, rises faster than the lower one, the actual error.10¹10²10³10⁴10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³terms summedrelative widthbracket widthactual errora bound with no tolerance in itbracket contains the truth7width at 10,000 terms4.4·10⁻⁴error it bounds2.4·10⁻⁶every sum computed twice, up and downtrue at every size, useful at none of the large ones
Fig. 17 The honest bracket, computed rather than derived, in the field that does it properly.
The randomised SVD against the optimum it cannot beatA semi-logarithmic plot of approximation error against target rank. A shaded band shows the spread across seeds, a solid line the optimal error from the exact singular values, and a dashed line the published probabilistic bound well above both.04812162010⁻¹10⁻⁰.⁵1target rank k‖A − A_k‖₂published boundrandomisedσₖ₊₁, optimalhow far apart the three areworst seed spread1.6bound / median at k = 125.9median / optimum at k = 121.960×60, 6 seeds, oversampling p = 5band is best to worst
Fig. 18 A guarantee that changes kind, which is the nearest thing on this site to a spread rather than an error.
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. 19 What the disagreement costs when a comparison is downstream of 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. 20 And when the comparison decides an integer.
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. 21 Another bound this site measured against what it bounds, and found loose by a factor with a shape.
Two perturbation bounds and the error that was measured, on a 8×8 matrix spread over 8 decades of unitsThree curves against the size of an entrywise relative perturbation. The normwise bound κ_∞·ε is a valid bound and sits 4·10⁷ times above the componentwise one cond(A, x)·ε, which is also a bound and is nearly attained by the worst of forty random perturbations at each size.10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ_∞ · εcond(A,x) · εmeasuredboth bounds holdκ_∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement
Fig. 22 A bound that describes an error better because it asks a different question about it.
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. 23 A discontinuity no bound in the literature contains, because it is a property of a library version.
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. 24 Buying accuracy back after the fact, which buys none of this field’s property.
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. 25 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. 26 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. 27 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. 28 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. 29 The convergence whose crossing point a tolerance sets.

Where this goes next

110¹10⁻¹¹10⁻⁸10⁻⁵pieces the vector was divided intodistance from the exact sum, relativethe published boundκ · uone vector, one algorithmdistinct answers21runs26spread, in ulps2.3·10⁷κ of the sum10⁸bound ÷ worst error2.6·10⁴nobody chose pand no answer is the answer The same program, twice One vector of 4,096 numbers, one summation algorithm, one precision, twenty-six runs — and twenty-one different answers. Nothing in the program chose between them, every one of them satisfies the textbook bound, and the exactly rounded answer is not among them. 110¹10²10³110¹10²10³the accumulating quantity, relative to its first valuethe error, relative to its first valuethe bounds: slope 1what all three do: slope ½three mechanisms, one exponenta left-to-right sum0.49a chain of rotations0.55a residual recurrence0.51every bound's slope1spread of the three0.067a bound is a sum of the roundingsand the roundings have signs Three walks and one bound A left-to-right sum, a chain of three thousand rotations and a conjugate gradient residual recurrence share no arithmetic and no vocabulary. Each has a standard bound that is linear in whatever it accumulates against. All three come out at a half — 0.486, 0.554 and 0.507 — and nothing is rescaled. 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 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.

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.

A-priori boundError accumulationLower boundReduction orderRun-to-run variationSummation condition numberUnit roundoffWorst-case analysis