The same arithmetic at a different price

Every cost claim on this site so far has been a count of arithmetic. n³/3 multiply-adds for an elimination, 2mn² for a Householder QR, one matrix–vector product per conjugate gradient step. Those counts are correct, they are the ones every textbook gives, and they have not decided which of two implementations is faster for about thirty years.

The reason is not that the counts are wrong. It is that the arithmetic is not what a large computation spends its time on. A modern processor performs an arithmetic operation in a fraction of a nanosecond and fetches a word from main memory in something like a hundred of them, so a computation that fetches each of its numbers once per operation is idle for 99% of its life whatever the operation count says.

This essay makes that measurable, on the algorithm this site has drawn since its foundation phase.

The two orderings #

Unblocked is the textbook loop, and it is what matrix.js has done since the site was built. Find the pivot in column k, swap, scale the column below it, then update the whole trailing submatrix with a rank-one product. Repeat n times.

Blocked defers. Factor a panel of b columns — pivoting inside it, updating only its own columns — and then apply the whole panel’s worth of updates to the trailing submatrix at once, as a rank-b product rather than b rank-one products.

The arithmetic is identical. The same numbers are multiplied by the same numbers and added to the same accumulators; what changes is the order the accumulations happen in and, crucially, how many times each entry of the trailing block is visited. The unblocked version touches every trailing entry once per column of elimination. The blocked version touches it once per panel.

That is the entire idea, and stated that way it sounds like it must be a straightforward win.

Three counters, and two of them do not move #

The figure carries three and the argument is that only one of them differs.

Operations, compared as integers. At n = 48 both orderings perform exactly 72,568 of them. Not “the same to within a percent” — the same number, and the assertion compares them with === rather than to a tolerance, because they are counts and a tolerance on a count is a way of not noticing that two algorithms are doing different amounts of work.

The pivot sequence. Both choose the same one, which is asserted and is not obvious: a blocked factorisation pivots within the panel, with less of the matrix in view than the unblocked version has. On a general matrix that can pick a different sequence, which is why LAPACK’s dgetrf applies the panel’s row interchanges to the rest of the matrix afterwards. Here it does not, and the assertion says so — because if it did, the residual comparison below would be between two different factorisations and would mean nothing.

And the residual. ‖PA − LU‖/‖A‖ comes out at 4.487946226420872·10⁻¹⁶ from both. Every digit. That is a stronger statement than the assertion needs — it requires agreement to 10⁻¹⁵ — and it is worth printing at full length because it settles the only question a reader could reasonably raise about the comparison. Reassociating floating-point arithmetic is not in general innocent; this site has an essay about exactly that. Here the reassociation happens to fall out bit-identically, and the badge says so rather than the caption.

The words moved: 41,332 against 19,476. A factor of 2.12, on the identical computation.

What “a word moved” means here, and what it does not #

There is no timing anywhere in this field and there will not be. A measured runtime on the machine that built this page is a fact about that machine: it depends on the cache sizes, the associativity, the compiler, what else was running, and it cannot be checked against anything. This site’s whole proposition is that every number is measured against something that was not measured, and a duration has nothing to hold it against.

What is counted instead is a model. A fast memory of M words, a slow one holding the matrix, and a count of the times an entry is touched that is not already resident — with M named on every figure, because nothing in the model is optimal until M is named. That is the same discipline algorithms-data-structures.com uses for its external-memory field, and the boundary between the two sites is already drawn: that site owns cost counted in block transfers as a subject, and this one is licensed for blocked and tiled numerical algorithms, where the block size is chosen for a cache and the question is what the reordering does to the answer.

Which is why there is a third counter here that the other site does not have.

The cache is modelled as fully associative with least-recently-used eviction, and that is a choice worth stating rather than burying. Real caches are set-associative and their conflict misses are a large part of why tuned libraries pick the block sizes they do. Modelling that would make every number here depend on an associativity nobody stated. The effect this essay is about — that the order of the same arithmetic changes the traffic by a factor — is present in the idealised model, is the part that generalises, and makes every count here a lower bound on a real machine’s.

Deferring the update is not the saving #

The first version of this file was wrong in an instructive way, and the wrongness is worth keeping in view because it is the mistake the idea invites.

Written the obvious way — defer the updates, then apply them in one nested loop over the trailing block — the blocked ordering moved 13,163 words against the unblocked version’s 13,018 at n = 32 with M = 64. Worse. The operation count was identical, the residual was identical, and the thing the whole exercise was for went the wrong way by 1%.

Deferring the update is not the saving. The saving is arranging the deferred update so that its working set fits. Three things have to be resident at once while the rank-b product runs: the tile of the trailing block being accumulated into, the strip of L beside it and the strip of U above it — about 3b² words. Written as one loop over the whole trailing block, the strips are far longer than that and every inner iteration evicts what the next one is about to read.

Tiling it — iterate over b×b tiles of the trailing block, and accumulate each tile completely before moving on — is a four-line change with no arithmetic in it and takes the same measurement to 19,476 against 41,332.

The relation 3b² ≈ M is where the next essay’s block size comes from, and it arrives here as a consequence of getting the loop right rather than as a tuning parameter. A blocked algorithm whose block does not fit is an unblocked algorithm with extra bookkeeping.

Why the flop count ever worked #

It is worth being fair to the count, because it was not always wrong and the reason it stopped being right is a fact about hardware rather than about algorithms.

On a machine where a memory access and an arithmetic operation cost about the same — which is what the machines these algorithms were analysed on were like — a count of operations is a count of everything, because the accesses are proportional to the operations and the constant is one. The n³/3 is then a genuine prediction.

What changed is that arithmetic got faster by about five orders of magnitude and memory latency by about two. The ratio between them is now roughly a hundred, and a computation that does one operation per fetched word runs at one per cent of the machine’s arithmetic rate. Everything about the design of numerical libraries since the 1980s follows from that number: the BLAS levels are a classification by operations per word, and level 3 exists because a matrix–matrix product does n³ operations on n² words and level 1 does n on n.

The blocked elimination above is that classification applied to elimination. The rank-one update is a level-2 operation, n² operations on n² words; the rank-b update is level 3, n²b on n². The whole of the difference in the figure is that ratio.

What this does not claim #

Not that blocking is a different complexity class. It is not, and the next essay measures the exponents: both orderings’ traffic grows like n³ — fitted at 3.182 and 3.185 — and blocking buys a constant of about 2.1 at every size drawn. The word “communication-avoiding” suggests something stronger than that, and for this algorithm at these sizes it is not stronger than that.

And not that the model is a machine. The 2.12 is the factor in an idealised memory with M = 144 words. On a real processor with three levels of cache, prefetching, and a compiler that unrolls the inner loop, the factor between a naive and a tuned LU is nearer ten, and almost none of the extra is in this model. What the model gets right is the sign and the mechanism, which is what an essay can be about.

The three counters are three different kinds of claim #

Worth separating, because the figure puts them on one axis and they are not one kind of thing.

The operation count is a combinatorial fact. It is decided by the algorithm before any machine is chosen, it is an integer, and it is the same integer on every machine that has ever existed. That is why it is compared with ===: a tolerance on it would be a tolerance on a fact about a loop nest.

The residual is a numerical fact. It depends on the arithmetic, on the order of the accumulations, and on the matrix. It has no reason to be identical between the two orderings and it is — which is a measurement rather than a theorem, and the badge prints all seventeen digits of it for that reason. Reassociating floating-point additions is not innocent in general; this site has an essay whose entire subject is a sum whose value depends on the order it is taken in. Here the reassociation happens to produce the identical rounding at every step, because the blocked update accumulates the same products into the same accumulator in the same sequence — it merely defers when that sequence runs.

The word count is a model. It depends on M, on the eviction policy, on the layout, and on the loop order — and it is the only one of the three that a reader could reasonably argue with. Which is why M appears in the title of the figure rather than in a footnote.

A field that had only the first two would be the elimination field. A field that had only the third would be algorithms-data-structures.com’s. Having all three on one figure is what makes the sentence the same arithmetic at a different price a measurement rather than a slogan.

What a reader should take from the factor of 2.12 #

Not that a blocked LU is twice as fast, which is a claim about a machine this field does not make.

The transferable statement is about which quantity to instrument. A reader who profiles an elimination and finds it slower than n³/3 suggests has two hypotheses available: the operation count is wrong, or something other than operations is being paid for. The first is checkable in an afternoon and is almost always wrong; the second requires a different instrument, and the instrument is a count of memory traffic rather than a count of arithmetic.

This site has been separating the algorithm’s fault from the problem’s fault since its foundation phase, and the identity underneath it — forward error ⪅ condition number × backward error — exists because a wrong answer has two possible authors and they are separately measurable. The cost field is the same move applied to time rather than to accuracy. A slow computation has two possible authors, and the flop count measures exactly one of them.

Where the traffic actually goes #

The counters say the two orderings differ by 2.12. Reading why off the algorithm rather than off the number takes one more count, and it is the count that makes the block size predictable at all.

The unblocked update at step k touches every entry of an (n−k)² trailing block, and it does that for every k. Summed, that is

Σₖ (n − k)²  ≈  n³/3

entries touched, and — this is the part that matters — none of them are reused between steps. Step k reads a trailing entry, updates it, moves on, and step k+1 reads it again after having read everything else in between. With M much smaller than the trailing block, everything is evicted before it is wanted, so the touches are misses and the traffic is n³/3.

The blocked update touches the same entries, in tiles, b columns of elimination at a time. Each trailing entry is now read once and written once per panel rather than per column, so the count falls by a factor of b — provided the tile’s working set stays resident, which is the whole of the next essay.

n³/(3b) against n³/3. At b = 6 that predicts a factor of six, and the measurement is 2.12.

The gap between six and 2.12 is not a discrepancy to be explained away; it is the panel factorisation, which is unblocked by construction and is not accelerated at all. A blocked algorithm blocks the trailing update and leaves the panel alone, so the saving is bounded by how much of the work is in the trailing update — which is most of it asymptotically and rather less of it at n = 48. That is also why the ratio in the figure creeps up slightly with size.

The level the operation lives at #

There is a vocabulary for this and it is worth having, because it turns the essay’s measurement into a rule of thumb that transfers.

The BLAS — the interface every numerical library sits on — classifies its operations into three levels, and the classification is not by what they compute. It is by operations per word:

Levels 1 and 2 have a constant ratio. However large the problem, they perform a bounded number of operations per word fetched, so on a machine where a fetch costs a hundred operations they run at a few per cent of the arithmetic rate and there is nothing any implementation can do about it.

Level 3 is the only one whose ratio grows with the problem, and that is the entire reason it is a level. A matrix–matrix product touches 3n² words and performs 2n³ operations, so the arithmetic per word is proportional to n, and a large enough product spends as little of its life waiting as one likes.

The unblocked elimination’s inner loop is a rank-one update: a level-2 operation, ratio 2. The blocked one’s is a rank-b update: level 3, ratio proportional to b. The 2.12 in the figure is that substitution and nothing else, and it is why “cast the algorithm in terms of level-3 operations” is the single design rule that all of LAPACK is organised around.

Which also says where the rule stops. An algorithm with no level-3 formulation available — a sparse triangular solve, a Krylov iteration whose only operation is a matrix–vector product — cannot be blocked into one, and its cost genuinely is what the naive count says. That is not a failure of effort; it is a property of the computation.

What is left #

The block size, which is the next essay: the curve above is drawn at one b and the trade behind it has an interior optimum that the model predicts from M alone.

The other direction of reassociation, which is the third essay in this field. Blocking reorders an elimination and gets the same factorisation bit for bit. A tall-skinny QR reorders a reduction and gets a different one — a different sequence of roundings, the same backward stability — and the reason it is safe is a property this site has been measuring since its foundation phase.

The permutation this essay avoided. The matrices here are diagonally dominant by construction, so both orderings pick the identity permutation and the comparison is between two factorisations of the same matrix. On a general matrix a blocked factorisation can and does choose a different pivot sequence, and then the two residuals are of two different objects and the bit-identical agreement above is not available. What that costs in stability is a real question with a literature and it is not measured here.