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

The vector that hides it

Every quick demonstration of a parallel sum uses positive numbers, and positive numbers are the one family where the effect is absent. Measured on six inner products this site already computes, the summation condition number runs from exactly 1 to 10¹⁷ — and the safe end is where nobody makes a decision.
33 min read 28 figures Measured, not assertedWhose fault is itWhere the algebra stops being the arithmetic

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

Three essays have argued from a vector that was built to argue from. A reader is owed the obvious question: does anything real look like that?

The answer is that the site has been computing sums like that for nineteen fields, and the ones that look nothing like it are the ones nobody makes a decision from.

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. 1 Six inner products taken out of computations elsewhere on this site, with their summation condition numbers measured rather than assumed.

Sixteen orders of magnitude between the top and the bottom of one figure, and the ordering is the argument.

The two at the bottom

rᵀr, a residual norm: κ = 1 exactly. Every term is a square, so no term can cancel another, and Σ|xᵢ| is Σxᵢ. There is nothing for a reduction order to disagree about beyond the last bit, and two machines summing this will agree to fourteen digits whatever they do.

pᵀAp, a curvature: κ = 4.1·10². A quadratic form in a positive definite matrix is positive, but its entrywise terms are not — the cross terms pᵢAᵢⱼpⱼ have both signs and cancel — so κ is a few hundred rather than one. Still small.

Those two are exactly the quantities the conjugate gradient method computes at every step. It is worth pausing on that, because it explains why the effect is so rarely noticed by people who work with iterative solvers every day: the method’s scalars are among the best-conditioned sums in the subject, and the disagreement they produce is a few hundred ulps rather than a few million. It is still enough to move an iteration count, which is what a stopping test is a race measures — but it is not enough to make anybody suspect the arithmetic.

The three at the top

A component of b − Ax: κ = 1.0·10¹⁷. The definition of a residual is a difference of two nearly equal things, so its terms cancel by construction, and the more accurate the answer the worse the cancellation. A residual of a converged solve is the worst-conditioned sum in the subject, and it gets worse as the solve gets better.

An orthogonality check qᵢᵀqⱼ: κ = 7.4·10¹⁵. Two orthonormal columns have an inner product of zero, so the terms sum to nothing at all while each is of size 1/n. The number the check returns is made entirely of rounding, and how much of it is made of which rounding is what a partitioning decides. Orthogonal is a number is the essay that turns the adjective into ‖QᵀQ − I‖; this is the observation that the number is one no two machines will agree about, and that the site’s own figures of it would move between builds if anything downstream of them were a comparison.

A nearly degenerate 2 × 2 determinant: κ = 3.6·10¹⁶. Two products of size 10¹⁶ whose difference is 1. One multiply the compiler removed follows that one all the way to its conclusion, which is that the answer is zero.

The reading

Sort the six by κ and the sorted list is very nearly a list of how the number is used.

The two well-conditioned sums are quantities that get reported: a residual norm goes into a convergence history, a curvature goes into a step length. Nobody compares either with a threshold in a way that changes the shape of an answer — a step length that differs in its eleventh digit gives a slightly different iterate, and the method carries on.

The three ill-conditioned sums are quantities that get compared. Is this residual below the tolerance? Is this inner product small enough to call the columns orthogonal? Is this determinant positive? Every one of them is a real number on its way to becoming a verdict, and every one of them is a sum whose value two machines will disagree about in its leading digits.

That is not a coincidence and it is not bad luck. A quantity is worth comparing with a threshold precisely when it is nearly zero, and a sum is ill conditioned precisely when it is nearly zero relative to its terms. The two conditions are the same condition. The reductions a program makes decisions from are, as a class, the reductions whose value is least well determined.

What this does to a benchmark

The immediate practical consequence is about measurement rather than about arithmetic, and it is the reason this effect has a reputation as a curiosity.

Run the hero experiment of the same program, twice on 4,096 positive numbers: the seven partitionings agree to a relative 1.6·10⁻¹⁵. Five of them still return distinct values — nothing about positive data makes addition associative — but the disagreement is at the last bit and no reasonable person would call it a problem.

Run it on a vector with κ = 10¹⁰ and the disagreement is 4.9·10⁻⁷, which is three hundred million times larger.

A demonstration that draws its data from a uniform generator, or a benchmark that sums a vector of times, or a tutorial that adds up the first million integers, is a demonstration on the left-hand end of that range. It will show a discrepancy of one or two ulps, which reads as a rounding curiosity, and the reader will conclude that the phenomenon is real and unimportant. The counterexample this essay’s refusal is built on is exactly that vector: the assertion positive data disagrees with itself is fed 4,096 positive numbers and must reject.

Where else the site has been quietly living with it

Reading the census backwards through the collection turns up several places where a number this site publishes is one whose reproducibility was never checked. Naming them is more useful than not.

Every ‖QᵀQ − I‖ on the site. The orthogonality figures are the site’s signature measurement, and they are computed from inner products with κ ≈ 10¹⁶. The values drawn are stable because this site’s builds are serial and seeded — a seeded generator is what makes a figure byte-identical between builds — but a reader reproducing them in a threaded environment should expect the leading digits to move.

Every residual norm quoted at the rounding level. A figure that reports ‖PA − LU‖/‖A‖ = 4.49·10⁻¹⁶ is reporting a number whose own leading digits are a property of the reduction that computed it. The order of the residual is a fact about the factorisation; its digits are not.

Every trace estimate. The census puts zᵀAz at κ = 42, which is mild — but the choice of probe is a decision already shows that the variance across probes swamps everything, so the reduction order is not the dominant term there. It is worth knowing which term dominates rather than assuming.

None of that invalidates a figure. It qualifies a claim: a residual at 10⁻¹⁶ is a residual at 10⁻¹⁶, and its third digit belongs to the machine.

The census is not a ranking of difficulty

A tempting misreading is that the six are ordered by how hard they are to compute. They are not. Every one of them is a single loop, and the two at the top are the shortest — a 2 × 2 determinant is two products. What the census orders is how much of each answer is made of rounding, which is a question about the data and not about the work.

That distinction matters because it decides where an intervention goes. Making a hard computation easier is an algorithmic question and this site has nineteen fields of them. Making an ill-conditioned sum better conditioned is not available at all: κ is a property of the multiset, and no reordering, no compensation and no wider float changes it. What those buy is a smaller error against the same κ, which is a different and more modest thing.

So there are exactly three moves, and the field’s remaining essays are about the second and third:

  • change the data — compute a different quantity whose sum does not cancel, which is what forming rᵀr instead of Σrᵢ does and why the residual norm is at the bottom of the census while the residual component is at the top;
  • change the arithmetic — compensate, or accumulate exactly, which shrinks the error against a κ that has not moved;
  • change the summation — make it order-independent, which does not shrink the error at all and makes the answer the same everywhere.

The third is the one this field is here for, and it is the one with no entry in the standard vocabulary, because the vocabulary was built for the first two.

A sum this site could not have written before

One entry in the census deserves its own paragraph, because it is the reason the census is a measurement rather than a list.

The residual component sits at κ = 1.0·10¹⁷, which is above 1/u. A summation condition number larger than the reciprocal of the unit roundoff means the computed sum has no correct digits at all in the worst case — the amplification exceeds what the format can carry. That is not a pathological construction: it is one row of b − Ax for a system with a condition number of 10¹⁰, which is an ordinary badly-scaled problem of the kind the condition number is an amplifier is about.

The reason it does not destroy anything in practice is that nobody uses that number on its own. It goes into a norm, and a norm is a sum of squares, and the squaring puts κ back to one. The subject’s habit of working in norms is doing more numerical work than it is usually given credit for, and the census is where that becomes visible: the aggregate is well conditioned and every one of its components is not.

Which is also the warning. A method that looks at residual components rather than at a norm — a componentwise stopping test, an entrywise backward error, a check for which equation is worst satisfied — is a method that has stepped from the bottom of this census to the top, and accuracy that is not a property of a method is the essay that argues those methods are the sharper ones. They are. They are also the ones whose answers two machines will disagree about.

Two things this does not say

It does not say ill-conditioned sums are wrong. Every one of the runs in the census is backward stable, and the residual with κ = 10¹⁷ is the exact residual of a nearby system. Cancellation exposes error that was already there, which is the reading cancellation takes the answer insists on and which is as true here as it is for a single subtraction.

And it does not say the ill-conditioned sums should be avoided. They are the sums worth computing. A residual is the most informative quantity a solver produces; an orthogonality check is the measurement this whole site is built on. The recommendation is not to stop computing them, it is to know that their reproducibility is not free — which is what the sum that cannot be wrong prices, and what what determinism costs puts a number on.

How to use the census

The operational version fits on a line: before worrying about reproducibility, compute κ of the sums involved. One pass, no estimator, no matrix.

If κ is near one, the arithmetic will agree between machines to the last bit or two, and any disagreement observed has another cause — a different library version, a different algorithm, a different problem — which is a much more useful thing to know than a vague suspicion of the hardware.

If κ is 10⁸ and the number is on its way to a comparison, the comparison is the thing to look at rather than the arithmetic, and the rest of this field is about what to do.

And if κ is 10¹⁶, the quantity is made of rounding and its leading digits are not a property of the problem at all. That is worth knowing before publishing them.

One last note on the census’s own honesty. Every κ in it is measured on a specific instance — one Q, one system, one probe — and κ is a property of a vector rather than of a kind of vector, so a different instance moves the numbers. What does not move is the ordering, and the ordering is what the essay is about: sums of squares at one, cross terms in the hundreds, differences of nearly equal things at 10¹⁵ and above. Those three bands come from the signs of the terms, and the signs of the terms are a property of what the quantity means.

That is the reason the census is worth having at all rather than a table of formulas: it is a measurement of the site’s own arithmetic, taken with the same instrument as everything else, and it would move if the site’s libraries changed underneath it.

The version of this to remember

Three sentences, and the third is the one that is not obvious.

A sum’s reproducibility is decided by its summation condition number — the total size of its terms divided by the size of its answer — and by nothing else about it: not its length, not its precision, not the number of workers it was divided across.

That number is exactly one for a sum of squares and above 10¹⁵ for a difference of nearly equal things, and it costs one pass to compute.

And the sums a program compares with a threshold are, as a class, the sums at the bad end — because a quantity is worth comparing with a threshold when it is nearly zero, and a sum is ill conditioned when it is nearly zero. The two conditions are the same condition, which is why this is a field and not a footnote.

At other settings

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 16 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 9.28·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 − Ax9.28·10¹⁶ad − bc, near-degenerate3.6·10¹⁶qᵢᵀqⱼ, an orthogonality check8.57·10¹⁵pᵀAp, a curvature86.6zᵀAz, a trace probe10.8rᵀr, a residual norm1measured, not assumedhighest9.3·10¹⁶lowest1above 10¹⁰3terms16sums of squares are safeand nobody decides anything from one
Fig. 2 The same six at sixteen terms, where the two sums of squares are still exactly one and the orthogonality check has less to cancel.
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 200 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 2.98·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 − Ax2.98·10¹⁷qᵢᵀqⱼ, an orthogonality check1.09·10¹⁷ad − bc, near-degenerate3.6·10¹⁶zᵀAz, a trace probe269pᵀAp, a curvature137rᵀr, a residual norm1measured, not assumedhighest3·10¹⁷lowest1above 10¹⁰3terms200sums of squares are safeand nobody decides anything from one
Fig. 3 And at two hundred, where it has more.
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. 4 What the field’s first figure looks like on data from the safe end of the census.
The same 4096 numbers, added up 26 ways, 20 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. 20 of the 26 runs returned distinct values, spanning 1.55·10⁻¹³ — 3.8·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 answers20runs26spread, in ulps3.8·10¹¹κ of the sum10¹²bound ÷ worst error2.6·10⁴nobody chose pand no answer is the answer
Fig. 5 And on data from the other end.
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 The line between the two, which is where every one of the six sits.
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. 7 Why a sum whose terms cancel commits its roundings at the wrong scale.
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. 8 The same mechanism in one subtraction, which is where this site first measured it.
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 And in a sum, where the order is a decision an author makes.
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. 10 The site’s signature measurement, computed from an inner product at κ ≈ 10¹⁶.
QᵀQ from classical Gram–Schmidt and from Householder on the 8×8 Hilbert matrixTwo eight-by-eight tables of QᵀQ. The upper one has ones on the diagonal and entries as large as one off it; the lower one is the identity to three decimal places everywhere.A = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle
Fig. 11 The same quantity printed entry by entry, which is where the cancellation is visible.
Classical and modified Gram–Schmidt: the same subtraction, in a different orderTwo panels of pseudocode differing in one argument, with the resulting pairwise dot products of the computed columns listed beneath each.for each previous column i, subtract the projection of column j onto qᵢclassicalr[i][j] = qᵢ · a[j] ↑ the ORIGINAL columnv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 7 and 8: 1columns 6 and 8: 0.13columns 6 and 7: 0.13modifiedr[i][j] = qᵢ · v ↑ what is LEFT of itv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 1 and 8: 4.4·10⁻⁷columns 2 and 8: 2.7·10⁻⁷columns 3 and 8: 2.4·10⁻⁸The two R factors agree to 1.2·10⁻⁶ relative. The two Q factors do not.the 8×8 Hilbert matrixone word, eight orders
Fig. 12 Two orderings of one orthogonalisation, in the field that first made the point.
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 A residual, which the census puts at the top of the range.
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. 14 A trace probe, where the variance across probes is the term that dominates.
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. 15 An estimator’s spread, which is a different kind of uncertainty about one number.
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. 16 The matrix version of the amplifier the census measures for a sum.
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. 17 The standing that makes each of these a measurement rather than a comparison of two guesses.
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. 18 What the reductions at the top of the census do under four hundred orderings.
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. 19 And what it would cost to make them agree.
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. 20 The conjugate gradient method, whose two reductions are the two at the bottom of the census.
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. 21 A Gram matrix, whose entries are the sums in the middle of it.
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. 22 And the determinant, which is the shortest sum in the census and the worst.
Four pivot rules on [[ε, 1], [1, 1]] with ε = 10⁻¹², before and after its first row is scaled by 1/εEight bars of forward error. On the matrix as given, no pivoting loses the answer entirely and all three pivot rules return it to rounding. Multiply the first row by 1/ε — which does not change the solution — and partial pivoting makes no interchange at all and returns the same wrong answer that no pivoting does. Scaled partial pivoting and complete pivoting are unmoved.forward error, relative to a solution of exactly (1, 1)no pivoting · as given1.3·10⁻⁴ 0 interchangesno pivoting · rows scaled1.2·10⁻⁴ 0 interchangespartial · as given2.2·10⁻¹⁶ 1 interchangepartial · rows scaled1.2·10⁻⁴ 0 interchangesscaled partial · as given2.2·10⁻¹⁶ 1 interchangescaled partial · rows scaled2.2·10⁻¹⁶ 1 interchangecomplete · as given1.1·10⁻¹⁶ 1 interchangecomplete · rows scaled2.2·10⁻¹⁶ 1 interchangethe same problem twicepartial, as given2.2·10⁻¹⁶partial, rows scaled1.2·10⁻⁴its relative residual1.2·10⁻¹⁶complete, rows scaled2.2·10⁻¹⁶the two systems have the same solutionand one pivot rule cannot see it
Fig. 23 A decision that looks like an algorithm and is a scaling.
The accuracy worth paying for in a preconditioner, against the condition number of the problemFour complete cost sweeps, each reduced to the ε that minimised the total work, against the κ of the problem it was run on. The shift that moves κ over 2.7 decades leaves every distance between every pair of clusters exactly where it was, so the ranks the partition asks for stay at 10, 9, 7, 5 throughout and nothing geometric is moving. The optimum walks 0.5 → 0.5 → 10⁻⁶ → 10⁻⁸ — six decades of accuracy, bought because the problem got harder and for no other reason. On the easiest problem the best hierarchical preconditioner in the comparison has rank one; on the hardest it is the tightest one on the sweep. How accurate an approximate inverse should be is a question with an answer, and the answer is not in the matrix's structure.10¹10²10³10⁴10⁵-9-7-5-3-11condition number of the problemlog₁₀ ε at the cheapest totalrank 1rank 1rank 6rank 5the knob answers to κκ = 210.5κ = 1.1·10⁴10⁻⁸rank at the easy end1rank at the hard end5ranks at 10⁻⁸, spread5the geometry did not moveand the right accuracy did
Fig. 24 An accuracy with an optimum in a size, which is the shape a caller expects and does not always get.
The proved bound against the error it bounds, at 24 significand bitsTwo curves against the Hilbert size on a logarithmic vertical axis: the error the computation actually made, which is known because the exact answer is the integers by construction, and the bound proved for it in interval arithmetic. The bound is above the error at every size, by 9.949 to 23.12 times. Past n = 5 the method refuses to return a bound at all.345678910111210⁻⁵10⁻³10⁻¹Hilbert size nerrorrefusedrefusedrefusedrefusedrefusedthe proved boundthe actual errorbound ÷ errorslack at n = 313slack at n = 423slack at n = 59.9the error is known, not estimatedthe exact answer is the integers 1 … n
Fig. 25 A bound measured against what it bounds.
Error of a difference quotient for J(x)v against ε, on the Bratu problem at n = 64The Jacobian of this problem is a formula, so the error of each quotient is measured against a derivative that is exact rather than against a better quotient. The forward difference falls with a slope of 1.01 — first order — reaches 1.28·10⁻¹⁰ at ε = 10^-6, and rises again with a slope of -1.00 as cancellation takes over. The central difference falls with a slope of 2.00 and bottoms at 1.11·10⁻¹² for two residual evaluations instead of one. Neither gets near the unit roundoff at any ε.10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice
Fig. 26 A floor under a computed difference, which is where a verdict stops meaning anything.
The accuracy asked for, against the accuracy obtained, for two kernels on one partitionε is applied to each block against that block's own largest singular value; the error is then reported against the whole matrix's norm. The two are not the same number and the dashed diagonal is where they would be. For 1/r the obtained error is 29 times smaller than the tolerance asked for; for log r on the identical partition it is 13 times smaller. Both curves are straight and parallel to the diagonal, so the knob does what a knob should — a decade in buys a decade out — and neither of them sits on it. The factor is how much of the matrix's mass lives off the diagonal, which is a property of the kernel; it is not the number of blocks, which is 112 and would put the curves on the other side of the diagonal.10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³ε asked for, per block‖A − A_H‖ ⁄ ‖A‖ obtainedasked = obtainedlog r1 ⁄ rtwo numbers, not one1/r, obtained at 10⁻⁸3.4·10⁻¹⁰log r, obtained at 10⁻⁸7.8·10⁻¹⁰1/r, obtained ⁄ asked0.034log r, obtained ⁄ asked0.078blocks in the partition112the tolerance is per blockand the error is per matrix
Fig. 27 What a tolerance actually obtains, measured rather than assumed.
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. 28 A bound whose worst case is never attained, in the field that measured it first.

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. A = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle Orthogonal is a number "Q is orthogonal" is a claim about a measurable quantity, ‖QᵀQ − I‖, and on the eight-by-eight Hilbert matrix two standard algorithms return 10⁻¹⁵ and 1 for it. The one that returns 1 still reconstructs the matrix perfectly, which is why nothing warns you. 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.

CancellationCatastrophic cancellationOrthogonality lossReduction orderResidualSummationSummation condition numberTrace estimation