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

The length that changes the kernel

A dot product's accuracy steps by a factor of 1.57 between 63 and 64 terms, on vectors drawn identically at both lengths. Nothing about the problem changes there. A library switches from one accumulator to four, at a constant in somebody else's source file.
33 min read 27 figures Where the algebra stops being the arithmeticMeasured, not assertedIdentical algebra, different arithmetic

Worth reading first: The same program, twice · The same arithmetic at a different price · The order they are added in.

Everything else in this field needs either a second worker or a compiler flag. This needs neither. It happens in one thread, in one process, from one binary, and it is caused by a constant somebody chose while tuning a library.

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. 1 Mean relative error of a dot product at each length from 58 to 70 terms, over forty draws each, with the vectors drawn identically at every length.

The error steps down by a factor of 1.57 between 63 terms and 64. The data is drawn from the same distribution on both sides. The arithmetic is the same precision. What changes at 64 is that the loop stops accumulating into one register and starts using four, because that is how a tuned kernel uses a vector unit.

Why four accumulators is the right thing to do

The switch is not a compromise. It is the correct implementation for two independent reasons and both are worth stating, because this essay is not an argument against it.

It is faster. A floating-point add has a latency of several cycles and a throughput of one per cycle, so a loop with a single accumulator is serialised on the dependency chain — each add waits for the previous one — and runs at a fraction of the machine’s rate. Four independent accumulators break the chain and let the pipeline fill. This is the same reasoning a block size is a property of the machine applies to a blocked elimination, one level down.

And it is more accurate, which is the direction of the step in the figure. Four accumulators is a four-way tree, each of whose walks is a quarter as long, which is exactly the √p improvement where the disagreement comes from derives — √4 = 2 predicted against 1.57 measured.

So the library is doing the right thing twice over. What it is also doing is putting a discontinuity in the accuracy of a routine at a length that has nothing to do with the caller’s problem.

The same shape one level up

Before leaving the dot product it is worth noting that this is the smallest instance of a pattern the site has already met at the level of a whole factorisation, and that the earlier essay reads differently in the light of it.

The same arithmetic at a different price measures a blocked and an unblocked elimination that perform 72,568 operations each, choose the same pivots, and return a factorisation identical to the last bit — while moving 41,332 and 19,476 words between fast and slow memory. Identical arithmetic, identical answer, half the traffic.

That identity is a property of that blocking: it re-orders which updates happen when, and each entry still receives the same sequence of individual multiply-subtracts. A blocking that accumulates a rank-b update before subtracting it — which is what a matrix-multiply-based factorisation does, and what every high-performance library actually runs — does not have that property. It sums b terms first and rounds once, which is a four-accumulator dot product wearing a larger hat.

So the fleet of blocked routines has a block size, the block size is tuned, and the answer depends on it in exactly the way this essay’s dot product depends on its cutoff. The earlier essay’s identical factorisations are the reassuring case, and it is worth knowing which case a library is in.

Why there is a cutoff at all

Below some length the four-accumulator loop is worse, for a reason that is entirely about bookkeeping: it needs a prologue to start four chains, an epilogue to combine them, and a scalar tail for the remainder. At sixteen elements that overhead is most of the work. So every tuned kernel has a threshold, and the threshold is chosen by measurement on the machine the library was tuned for.

Which means the constant is:

  • not a property of the problem — the caller’s data does not appear in the decision;
  • not a property of the mathematics — both loops compute the same sum in exact arithmetic;
  • a property of the library version and the machine it was tuned on, and therefore something that changes when either does.

What moves when the constant moves

The drag makes it concrete: moving the cutoff moves the step and changes nothing else. The discontinuity sits wherever the library was tuned to put it.

The consequence for a caller is specific and rarely noticed. A version bump that retunes this constant from 64 to 32 changes the answer of every dot product whose length lies between them — about a factor of two in the error, in either direction, on every such call. No release note has ever described that as a change of results, because from the library’s point of view it is a performance change: both kernels are correct, both are backward stable, and the accuracy of neither was specified.

And the caller’s symptom is the one this whole field produces: the answers moved and nothing in the program changed. Which is why algorithm selection belongs in this field rather than in the cost field, where its cause lives.

The refusal, which is the sharper statement

The refusal published with this essay is worth reading as a claim rather than as machinery, because it inverts the natural way to describe the figure.

The natural description is a dot product of 60 terms has this accuracy. That is false. At 60 terms the mean relative error is one number if the library’s cutoff is 64 and a different number if it is 128, on identical data — so the accuracy at a given length is a property of which kernel that length selected, and the length is only a proxy for that.

Stated as a rule: a routine whose implementation is selected by the size of its input has an accuracy that is a step function of the size, and the steps are at sizes the caller cannot see.

How wide the family of such constants is

A dot product’s accumulator count is the smallest instance. Once the shape is recognised it is everywhere in numerical software:

A blocked factorisation’s block size, which decides how many updates are accumulated before being written back. The same arithmetic at a different price measures a blocked and an unblocked elimination performing identical operations and returning identical factors — that measurement is for a blocking that re-associates nothing, and a matrix-multiply-based update does re-associate, at a size threshold.

A crossover to a different algorithm entirely. A small eigenproblem solved by a direct formula and a large one by an iteration; a small sort by insertion and a large one by a merge; a small matrix multiplied conventionally and a large one by a recursive method. Each crossover is a size at which the answer changes because the method did.

A threshold that switches to a different precision. Mixed-precision libraries choose a working precision by size or by conditioning, which is the accuracy worth paying for’s subject — and there the change in the answer is large enough that nobody could miss it.

A dispatch on the machine. The same library, on two processors with different vector widths, runs different kernels with different accumulator counts. The binary is identical; the code path is not.

That last one is the version of this that is hardest to reason about, because there is no constant to look up. The answer depends on which instructions the processor supports.

What the step is worth, in the direction that matters

A reader who has followed the field this far will notice that the step in the figure is down — the four-accumulator kernel is more accurate — and may reasonably ask why a change that improves the answer is a problem.

It is a problem for exactly one reason and it is worth isolating: a discontinuity is not an improvement, whichever direction it points.

A caller who grows their problem from 63 unknowns to 64 sees the accuracy of one operation improve by 57%. If they were comparing against a stored result, the comparison fails. If they were extrapolating a convergence study across sizes, the extrapolation acquires a kink. If they were tuning a tolerance against measured behaviour, the tolerance was tuned on the wrong kernel.

And the same step, met from the other side, is a regression: a caller whose library retunes the cutoff upward finds the accuracy of their 100-term dot products getting worse by the same factor, from a version bump that mentioned only performance.

The site has a standing way of describing this and it applies here without modification. Accurate is not a property of a method makes the case that accuracy belongs to a method-and-matrix pair rather than to a method; this essay adds that on a tuned library it belongs to a method-and-matrix-and-size triple, and that the size dependence is not continuous.

What a caller can do about it

Less than in the rest of the field, and it is worth being straightforward about that.

Nothing about the arithmetic. The disagreement is between two kernels, both correct, chosen inside a library. No summation policy applied at the call site changes which one runs.

Pin the version, and record it. The constant belongs to the version, so a computation whose answer must be reproducible in a year needs the version recorded alongside it. This is ordinary software practice and this field’s contribution is a reason: not that the library might have bugs, but that a correct change to a tuning constant changes results.

Or take the reduction back. A caller who needs a bitwise-stable inner product can write one — three operations an element, pre-rounded, and no dispatch, which is the policy the sum that cannot be wrong sets out — and accept the loss of the vendor’s tuning. What determinism costs prices that trade, and on a memory-bound reduction it is cheaper than it sounds.

And, at minimum, do not test near a cutoff. A test suite whose problem size sits at 63 or 64 is sampling a discontinuity, and a version bump will move it across. Sizes chosen away from any power of two are less likely to sit on one.

The measurement is harder than it looks

Worth describing, because the figure took three attempts to get right and the two failures are instructive about how this class of defect hides.

The first attempt varied the length and the data together, drawing fresh vectors at each length from a fresh stream. That measures the step and the ordinary growth of error with n, and on a short sweep the two are the same size — so the figure showed a jagged line with a step somewhere in it, which is a picture of nothing.

The second attempt measured a single draw at each length, which is a sample of a random variable whose spread is larger than the effect. The step was invisible under the scatter. Forty draws per length is what makes the mean stable enough for a factor of 1.57 to be visible.

The third — the figure above — draws from one stream at every length and averages forty. What is held fixed is the distribution; what varies is the length and, at one length, the kernel.

The general lesson is the one this field keeps arriving at from different directions: an effect of this size is not visible without holding something fixed that nobody thinks to hold fixed. The partitioning effect needs one vector across many divisions; the contraction effect needs one matrix across two compilations; and this one needs one distribution across many lengths. Every one of the three is invisible to a measurement that varies the obvious thing.

Where this leaves the field’s account

Three causes have now been measured, and they are genuinely different mechanisms:

cause what varies fixed by
the partitioning of a reduction how many pieces, and where the cuts are an order-independent policy
expression contraction whether a multiply-add was fused a compiler flag, or an explicit fma
algorithm selection which kernel a size selected pinning the version

The first is the field’s main subject and the one with a clean repair. The second is one instruction wide and has a repair that costs nothing. The third has no repair inside the caller’s program at all, and it is the one that will still be there after everything else has been made deterministic.

That ordering is also the ordering of how visible they are. A reduction’s disagreement can be measured from the call site by varying the thread count. A contraction’s can be measured by compiling twice. A kernel cutoff can only be measured by sweeping the problem size and looking for a step, which is not something anybody does unless they already suspect it.

What a library could publish

The field’s recommendations have so far been addressed to callers, because callers are the ones who can act. This one is addressed to the other side of the interface, and it is short.

Publish the cutoffs. A library that documents this routine uses a four-way accumulator above 64 elements has converted an invisible discontinuity into a fact a caller can design around: they can avoid sizes near it, or force one kernel, or record it beside a stored result. The constant is not a secret and it is already in the source; what is missing is a statement that it affects results.

And treat a change to one as a change to results. A release note that says improved dot product performance describes a change that moves numbers. Nothing about the change is wrong; what is wrong is the category it is announced in, and a note saying the accumulator threshold moved from 64 to 32; dot products of length 32 to 63 now return slightly different values would cost one line and would answer a class of bug report entirely.

Neither of these is a numerical proposal. Both are the same observation this field keeps making from different angles: the property that two runs return the same answer is real, it is currently nobody’s responsibility, and most of what would fix it is a matter of saying out loud what is already true.

At other settings

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 26 to 38 terms, over 40 draws each. Below 32 the kernel accumulates into one register; at 32 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.46× 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.2628303234363810⁻¹⁷10⁻¹⁶terms in the dot productmean relative errorthe kernel changes herea constant in a libraryone accumulator3.4·10⁻¹⁷four accumulators2.3·10⁻¹⁷step at the cutoff1.5cutoff32the problem did not changethe loop did
Fig. 2 The same measurement with the cutoff at 32, where the step has moved and nothing else has.
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 122 to 134 terms, over 40 draws each. Below 128 the kernel accumulates into one register; at 128 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.47× 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.12212412612813013213410⁻¹⁷10⁻¹⁶terms in the dot productmean relative errorthe kernel changes herea constant in a libraryone accumulator2.9·10⁻¹⁷four accumulators2·10⁻¹⁷step at the cutoff1.5cutoff128the problem did not changethe loop did
Fig. 3 And at 128 — which is where a dot product of 100 terms gets a different answer from the same source.
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 field’s main cause, which this one is independent of.
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. 5 The mechanism behind the step: four accumulators are four shorter 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. 6 The size of the disagreement a reduction produces, for comparison with this one.
A 6.25-bit block format against 8-bit E4M3, read two waysFour curves of relative error against the spread of the data. Two, for the block format, rise steeply when read as a median and gently when read as a norm; the other two are nearly flat.048121610⁻²10⁻¹1octaves of spread within a blockrelative errorblock, medianblock, 2-normE4M3, medianE4M3, 2-normthe same data, two readingsblock format, bits a value6.3norm across 16 octaves1.8median across 16 octaves62entries deleted at the wide end383640 values, blocks of 32a norm is dominated by what a block format keeps
Fig. 7 Another crossing between two implementations of one computation.
Words moved against the block size, n = 48, M = 100A U-shaped curve on logarithmic axes. Too small a block refactors the panel too often; too large a block does not fit in the fast memory and the tiled update evicts what it is about to read. The minimum is at b = 8, which is √M − 2. The three-line count says √(M/3) = 5, which has the right scaling and the wrong constant.110¹10⁴10⁵block size bwords movedthe count: √(M/3) = 5measured best: b = 8words movedat the best block1.6·10⁴at b = 13.9·10⁴at b = 243.9·10⁴derived from M with no measurement, and scannedthe two agree
Fig. 8 The cost model that puts a constant like this one in a library.
What each way of eliminating a constraint inherits, over five decades of κ(A)The same system solved twice, at 8 unknowns and 3 constraints with κ(H) = 100. The range-space method forms S = AH⁻¹Aᵀ and inherits κ(S), which rises from 12.99 to 3.981·10¹⁰ — the square of κ(A), for the reason the normal equations square it. The null-space method solves with the reduced Hessian ZᵀHZ, whose condition number is 21.13 at the start of the sweep and 21.13 at the end: it does not contain κ(A) at all. The two forward errors, measured against a solution computed in BigInt rationals, follow their own condition numbers: 5.314·10⁻⁶ against 5.788·10⁻¹² at the far end. Both methods are correct and one of them is usable.01234510⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹log₁₀ κ(A)condition number, and relative errorκ(S)κ(ZᵀHZ)range-space errornull-space erroragainst a BigInt answerκ(S) at κ(A) = 10⁵4·10¹⁰κ(ZᵀHZ), all stops21range-space forward error5.3·10⁻⁶null-space forward error5.8·10⁻¹²both are the same algebraand only one squares
Fig. 9 Two eliminations whose operations are identical and whose costs are not.
Communication for a 512×12 factorisation on 16 processorsTwo counts for the same factorisation. Rounds on the critical path: 48 for the column sweep, 4 for the reduction tree, 4 for Cholesky QR — a factor of n between the first and the other two. Words sent: 1170, 1170 and 2160 — the sweep and the tree send the same number, and the method with the fewest rounds sends the most.rounds on the critical pathHouseholder sweep48reduction tree4Cholesky QR4words sentHouseholder sweep1170reduction tree1170Cholesky QR2160two counts, two rankingsrounds, sweep ÷ tree12words, Cholesky ÷ tree1.8arithmetic, tree ÷ sweep1.5the rounds separate the threeand the words do not
Fig. 10 Why an implementation divides work 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. 11 The accuracy axis of the hero figure, in the essay that measured it as a choice.
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. 12 The policy that removes the field’s other cause and not this one.
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. 13 And what it would cost to take the reduction back from the library.
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. 14 The third cause, one instruction wide.
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. 15 What a change of kernel does when a comparison is downstream.
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. 16 And when it decides an integer.
The smallest defect a regression test can catch, on a machine that will not repeat itselfA 150 × 150 system with κ = 10⁴, solved by conjugate gradients to 10⁻¹⁰ on six partitionings. The shaded band is how far the answer moves across those six when nothing at all is wrong: 3.24·10⁻¹², relative. Each pair of dots is one defect — a single entry of the matrix changed by a relative δ, an assembly slip or a stale coefficient — solved on the same six machines, showing its nearest and furthest run. A defect is catchable only when its nearest run clears the band, because a test whose threshold sits inside the band fails a correct build. Everything at or below 10⁻¹³ overlaps and is invisible; 10⁻¹² is the smallest that separates, at 4.61·10⁻¹². So the tolerance exists, it is bracketed between 3.24·10⁻¹² and 4.61·10⁻¹² — a factor of 1.42 — and it is neither zero nor the 10⁻⁸ that usually gets typed.-15-13-11-910⁻¹³10⁻¹¹10⁻⁹10⁻⁷log₁₀ of the defect, relative to the entry it sits inhow far the answer movedthe machine, on its owna tolerance with two sidesmachine band3.2·10⁻¹²smallest caught10⁻¹²its nearest run4.6·10⁻¹²defects hidden3window, factor1.4below the band nothing is visibleand above it everything is
Fig. 17 The test that has to survive a version bump.
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. 18 Which inner products this matters for at all.
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. 19 An accuracy that has an optimum in a size, which is the shape a caller expects and does not get.
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. 20 The knob that is a caller’s, unlike the constant in this essay.
Multiplications in a hierarchical solve and in a dense factorisation, and where they crossBoth counted rather than estimated: the recursion carries a counter and reports what it actually did. The dense factorisation is n³/3 and rises at exactly three per doubling. The hierarchical solve rises at 2.13, 1.93, 1.74 — falling, because the cost is n log²n and its exponent is on its way to one. The two cross between n = 64 and n = 128: below it the format is the more expensive way to solve the system, at 1.48 times the dense count, and at n = 512 it is 6.2 times cheaper. Every point returns an answer at a backward error of about 1.4·10⁻¹⁰, so the comparison is between two ways of getting the same thing.567891010⁴10⁵10⁶10⁷10⁸10⁹log₂ nmultiplicationsdense factorisation, n³⁄3the recursion, countedwhere the format starts payingratio at n = 641.5ratio at n = 5120.16exponent, first doubling2.1exponent, last doubling1.7backward error1.4·10⁻¹⁰cheaper is a sizenot a property
Fig. 21 A cost model whose currency is not arithmetic, which is why the constant exists.
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 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. 23 An estimator’s spread, which is a different uncertainty about one number.
The residual conjugate gradients reports and the residual of the vector it holds, along one runA 50×50 diagonal matrix with one eigenvalue at 10⁻¹⁰ and the rest between 1 and 50, so every product with it is exact to a rounding and nothing below can be blamed on the operator. The lower curve is r ← r − αAp, updated by the recurrence, which is what a stopping test reads. The upper curve is ‖b − Ax‖/‖b‖, recomputed from the iterate. They start as the same vector and end 8.6·10⁸ apart, with the reported one at 6.89·10⁻²¹ — below the unit roundoff of 1.11·10⁻¹⁶, which is the shortest proof available that it is not the residual of anything.0122436486072849610⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best6.9·10⁻²¹the answer's, at its best5.9·10⁻¹²unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10⁹iterations drawn96the recurrence remembers every roundingand the stopping test is written in it
Fig. 24 Two residuals for one iteration, of which a stopping test sees one.
What the solver reported, and what it achievedFinal relative residual and final relative error against the working precision, on a logarithmic vertical axis. The residual sits at the stopping tolerance at every precision — the method converged, by every test available to it. The error runs from 1.4·10⁻¹³ at 53 bits to 0.00597 at 8, tracking the unit roundoff, which is drawn beside it.614223038465410⁻¹³10⁻¹⁰10⁻⁷10⁻⁴working significand bitsrelative sizeerrorresidualuthe gap the solver cannot seeerror ÷ residual at 24 bits1.5·10⁶error ÷ residual at 16 bits2.3·10⁷error ÷ residual at 8 bits2.4·10¹⁰every convergence test passesand the answer is wrong in proportion to u
Fig. 25 A residual that cannot see what it is asked about.
Numbers stored per unknown, against the size of the matrix, at 10⁻⁸The dense matrix stores n numbers per unknown, which is the straight line through the origin and doubles whenever n does. The hierarchical representations store 54, 79, 106, 133 and 50, 70, 94, 120, adding about 26 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 26 and 23 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.55, 1.42, 1.33, falling towards one and never arriving, which is what a logarithm looks like from inside.56789100108216324432log₂ nnumbers stored per unknowndenseweak: every off-diagonal blockstrong: only the admissible onesa constant per doublingstrong, n = 5126.8·10⁴weak, n = 5126.1·10⁴dense, n = 5122.6·10⁵per doubling26‖A − A_H‖ ⁄ ‖A‖3.8·10⁻¹⁰the dense line doublesand the other two add a constant
Fig. 26 A storage cost that follows from a decision about a threshold.
Conjugate gradients at κ = 104, against the bound κ permitsA semi-logarithmic plot of the relative A-norm error against iteration count. The measured curve falls below a smooth dashed curve showing the classical condition-number bound.0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417
Fig. 27 The convergence whose crossing point a tolerance sets.

Where this goes next

10²10³10⁴10⁵matrix size ncountoperations, bothwords, unblockedwords, blocked (b = 6)the answer does not move‖PA − LU‖/‖A‖, unblocked2.8·10⁻¹⁶‖PA − LU‖/‖A‖, blocked2.8·10⁻¹⁶difference between them0the dashed curve is both orderings' operation countthe solid pair is what they cost The same arithmetic at a different price A blocked and an unblocked elimination perform 72,568 operations each — the same operations, associated differently — choose the same pivots, and return a factorisation identical to the last bit: ‖PA − LU‖/‖A‖ = 4.487946226420872·10⁻¹⁶ in both. One of them moves 41,332 words between fast and slow memory and the other moves 19,476. 110¹10⁴10⁵block size bwords movedthe count: √(M/3) = 5measured best: b = 8words movedat the best block1.6·10⁴at b = 13.9·10⁴at b = 243.9·10⁴derived from M with no measurement, and scannedthe two agree A block size is a property of the machine Three lines of counting say the best block size is √(M/3). Scanned over every integer at five fast memories, the measured optimum is √M − 2 — exactly, at all five. The count has the right scaling and the wrong constant, low by a factor of 1.56, and the wrong form: the answer is affine in √M rather than proportional to it. 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.

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.

Algorithm selectionBitwise reproducibilityBlocked algorithmError accumulationLoop orderMemory hierarchyReduction orderRun-to-run variation