A block size is a property of the machine

The previous essay draws one curve for the blocked ordering, at one block size, and takes the block size for granted. This essay is about where that number comes from, and about the fact that it is not a matter of taste.

A trade with a cliff on one side #

The scan runs over every integer block from 1 to 24. At n = 96 with M = 144 words of fast memory:

It is not a symmetric U. The left-hand side is a slope — 4.35× at b = 1, falling smoothly through 2.67, 1.43 and 1.04 — and the right-hand side is a cliff: 109,202 at b = 10 and 294,445 at b = 11, a factor of 2.7 for one more column. Past it the curve is flat at about 2.9× and stays there all the way to b = 24.

That asymmetry is the mechanism, drawn. Too small a block is a gradual waste, because the panel factorisation is unblocked and its share of the work shrinks smoothly as b grows. Too large a block is a discrete failure: the working set either fits in M or it does not, and one column past the boundary every inner iteration evicts what the next one is about to read. There is no partial credit on the right-hand side.

Which means a block size chosen conservatively — too small — costs a factor of 1.4 and one chosen optimistically costs a factor of 2.9, on the same machine, for the same reason a reader might round either way. The assertion requires the largest block in the scan to cost more than 1.5 times the best; it measures 2.91.

Every integer, because the grid was lying #

The first version of this scan sampled 1, 2, 3, 4, 6, 8, 12, 16, 24, 32 — a reasonable-looking geometric grid — and it gets the answer wrong. On this matrix it reports the optimum at 6 and the true optimum is 10.

The reason is visible in the table above and is not noise. The grid’s next sample after 6 is 8, and 8 is worse than 6: 122,413 against 113,249. So the coarse scan sees the curve turn up and stops, having found a local minimum at 6 on a curve whose real minimum is at 10 and 4% lower. The non-monotonicity between 6 and 10 — 113,249, 132,963, 122,413, 115,629, 109,202 — is an integer effect, the panel width dividing into 96 or not, and a grid that samples two of those five points can land on either.

That mattered more than it should have, because the number being extracted is compared against a prediction. A prediction confirmed by a scan that cannot resolve it is not confirmed by anything, and the coarse grid’s answer of 6 agreed with √(M/3) exactly.

What the count says, and what the scan says #

The tiled update holds three things at once, so about 3b² words, so

b  ≈  √(M/3)  =  0.577 √M

which is three lines with no measurement in them anywhere. At M = 144 it gives 6.

The scan, over every integer, at n = 96:

b = √M − 2, exactly, at all five.

Drawn rather than tabulated, the five memories are five scans, and the column of twos is the distance from each minimum to the dashed line the count puts on the same axes.

The dashed line is not lagging by a fixed fraction, which is what a wrong constant would look like. Over the five memories the prediction reads 4, 5, 6, 8 and 9 against a measured 6, 8, 10, 12 and 14, so its ratio wanders — 1.50, 1.60, 1.67, 1.50, 1.56 — while its difference is two every time. A constant ratio and a constant difference are different relations, and the figure shows which one is on the page.

That last number is the one to keep from the sweep, and it is a scaling rather than a constant. The Hong–Kung floor is n³/√M, so at a fixed matrix the traffic a perfect schedule moves falls as M⁻⁰·⁵ — quadruple the cache and halve the words. The scan’s own optima fall from 158,119 to 88,877 across exactly that quadrupling, which is M⁻⁰·⁴¹⁶. Fast memory is bought at a discount here: each doubling returns about 1.33× rather than 1.41×, and the shortfall is not a constant overhead that a bigger cache would amortise, because a constant overhead would make the exponent approach −0.5 from below as M grows and it does not.

The second half of the title is checkable at all five memories rather than at the one the earlier figure draws. Every scan above reports its optimum unchanged at a matrix half again as large — still 6 at M = 64, still 10 at 144, still 12 at 196, still 14 at 256 — so the invariance measured at M = 100 is not a coincidence of that memory.

So the count is wrong twice over. Its constant is low by a mean factor of 1.56, and its form is wrong: the measurement is affine in √M rather than proportional to it, and no choice of constant in front of √M reproduces a column of twos.

Working backwards from the measurement, the count that fits is b² + 2b = M — a tile of b² words plus a row strip and a column strip of b each, rather than three tiles — which gives b = √(M+1) − 1 and is within one of the measurement everywhere. The original overcounted the working set by a factor of three by charging tile-sized rents for strip-sized things.

modelledBlock is left as it was rather than corrected, and the dashed line on the figure still marks √(M/3). Replacing it would remove the only interesting thing about it. A derivation that gets the scaling right and the constant wrong is a normal and instructive outcome, and this collection’s habit is to print both numbers rather than to publish whichever one survived.

And it does not move with the problem #

The other half of the title. Fix M and grow the matrix: at M = 100 the optimum is 8 at n = 48, 56, 64, 72, 80 and 96. Nothing about the matrix decides it. Only how much of the matrix can be held at once decides it, which is a property of the machine.

The assertion allows a slack of two rather than requiring equality, because the integer effects above are real: at M = 144 the optimum is 10 at four of six sizes and 6 at the other two, and the two are 4% apart on a curve that is not monotone between them. Requiring exact equality would be asserting something about which integer divides into n.

What the floor is, and how far above it this sits #

There is a lower bound, and it is one of the genuinely surprising results in the area. Hong and Kung proved by a pebbling argument that any schedule of the classical n³ matrix multiply must move

words  ≥  Ω(n³ / √M)

between a fast memory of size M and a slow one, and Ballard, Demmel, Holtz and Schwartz extended it to essentially every direct dense linear algebra algorithm, LU included. It is a statement about every possible ordering, not about any particular one — no cleverness gets below it.

The constant depends on which proof is used, so what is published here is the ratio of the measured traffic to n³/√M rather than a claim that the bound is attained:

The ratio settles at about 2.05 and stays there. That is the useful number: it says how much of the available saving a given implementation has taken, it has a floor no ordering can pass, and an implementation reporting a ratio of a hundred has not started while one reporting two has finished.

And the exponent does not move #

This is the finding, and it is a correction to the phrase the area is usually described with.

Fitted across n = 16, 24, 32, 48 and 64:

blocked traffic    ∝  n^3.182
unblocked traffic  ∝  n^3.185

The same exponent, to within the fit’s own noise on five points. What blocking changes is the constant, and the constant is 2.19, 2.11, 1.98, 2.17 and 2.17 across those sizes — flat.

Which is exactly what the floor says it must be. The bound is n³/√M with M fixed, so both orderings are Ω(n³) and no reordering can make either of them anything else. Blocking moves the constant towards the floor’s; it cannot move the exponent, because the exponent is a property of the computation rather than of the schedule.

The word communication-avoiding suggests something stronger, and for a fixed M on a two-level memory it is not stronger. Where it does become stronger is where M is allowed to grow with the problem, or where the memory has more than two levels, or where “communication” means messages between processors rather than words between caches — and none of those is this measurement. The assertion here is written in both halves for that reason: the exponents are required to be near three and the ratio is required to stay above 1.2 at every size, so a reading in which blocking is a complexity improvement and a reading in which it is nothing both fail.

This is the same shape as nested dissection losing to minimum degree in the sparsity field, and the same shape as Szegő’s limit never being reached in the structure one. A statement about a limit and a statement about a computation are different statements, and this collection has a whole verdict category for the cases where they get confused.

The refusal, which is about the model rather than the algorithm #

assertTheCostAssertionsReject feeds the claim that blocking helps a matrix that already fits in fast memory, and requires it to fail.

That is worth having because it is the failure mode a memory model invites. Set M large enough and every miss is a compulsory one — the matrix is read once, whatever order it is read in — and both orderings move exactly the same number of words. Measured at n = 16 with M = 576: 256 against 256. A figure drawn there would show two coincident curves under a caption about a gap, so the generator filters its sizes to those with n² > 2M and asserts that some remain.

The general form of that mistake is a model with nothing to model in it, and it is not specific to caches. It is the same defect as a refusal that waits for a rare event, which this site hit in its expansion phase: an assertion that cannot fire is not evidence, and one that fires on everything is not either.

What this field will not do #

No durations. Worth restating, because the temptation is constant and the payoff would be immediate: a wall-clock comparison would produce a bigger and more impressive factor than 2.12, and it would be a fact about one machine on one afternoon. There is nothing to check it against, and this collection’s rule is that a number arrives with something it was measured against.

No claim about a real cache. Fully associative, LRU, one level, no prefetching, no line size. Every one of those is a simplification and every one of them makes the measured factor smaller than a real machine’s, which is the direction that keeps the argument honest.

And no tuning. Whatever a real library does to pick a block size — and ATLAS and its descendants do it by exhaustive search on the target machine, which is the strongest possible admission that the model above is incomplete — is engineering rather than numerics. What is worth having here is that the number is a property of the machine rather than of the algorithm, that a three-line count finds it, and that both facts are checkable.

What the count got wrong, and why that is the normal outcome #

It is worth being explicit about the shape of the failure, because “the model was out by a factor” covers two very different situations and this is the instructive one.

A model can be wrong because the mechanism is wrong. Then the scaling is wrong, the error grows with the parameter, and no constant rescues it — which is what happens if the block size is chosen by reasoning about the flop count, since the flop count does not contain M at all and would give a block size independent of the machine.

Or a model can have the right mechanism and mis-count the pieces. Then the scaling survives, the error is a constant factor, and the model remains useful for the thing scaling laws are for: saying what happens when the machine changes. That is this case. Double M and both √(M/3) and √M − 2 say the block should grow by about √2, and they agree about that at every M in the table while disagreeing about the value at every one of them.

The distinction matters because of what each is good for, and the usual next sentence — nobody should choose a block size from √(M/3), it is wrong by 1.56 — measures the wrong thing. A block size is chosen for what it costs, not for how close it is to the argmin, and those are different numbers.

Cache misses at every block from 1 to 11, relative to the best block at each fast memory, n = 48:

Using the model’s block costs 3%, 7%, 0%, 9% and 14%. Never more than a seventh of the traffic, and twice under a tenth. The model is wrong about the value by a factor of 1.56 and wrong about the cost by at most fourteen per cent, because the curve is nearly flat near its minimum. The harsh verdict is therefore harsher than the quantity being chosen supports: √(M/3) is a perfectly serviceable block size, and its real defect is that it does not know it is serviceable.

The asymmetry has a shape worth stating. Undershooting is expensive and gradual — a block of one costs two and a half to three and a half times the best, and the penalty climbs steadily as the block shrinks. Overshooting is free and then catastrophic: a plateau within twenty per cent of the best for several blocks past the optimum, then a step to between 1.9 and 2.9 times the moment the working set stops fitting.

And the plateau is bought by the machine. At M = 64 there is none — the optimum is at b = 6 and b = 7 already costs 1.88 times as much — while at M = 256 the cost stays within twenty per cent of the best from b = 6 all the way to b = 14. The room to be wrong on the high side is a function of the fast memory, and on a small cache there is no room at all.

The far end of the curve is worth reading across the five memories too, because it goes the way nobody expects. A block of 24 costs 2.02, 2.56, 2.91, 3.34 and 3.54 times the best as M runs from 64 to 256 — the penalty for a fixed too-large block grows with the fast memory. Nothing about the block changed; what changed is that the optimum moved away from it, and past the cliff the cost is set by how badly the working set misses rather than by how large it is.

Which gives the practical rule the constant never could: choose a block a little below the cliff, rather than at the argmin. The argmin is not sharply defined, and the only feature of this curve that costs real traffic is the edge.

Three things about that rule are worth separating, because they are what a scan of one machine cannot say and what a model of the wrong kind gets confidently wrong.

The edge is where the working set stops fitting, and that is the quantity the model is about. A block of b needs roughly 3b² words resident — the panel, the row block, and the piece of the trailing update being touched — so the step happens where 3b² passes M, which is exactly √(M/3) rearranged. The model being dismissed for its constant is a good model of the cliff and a poor one of the minimum. Those are different features of the same curve, and only one of them is worth predicting, because only one of them costs anything.

The penalty for being below the edge is bounded and the penalty for being above it is not. Below, the worst case is a block of one, which is the unblocked algorithm and costs about three times the best. Above, the cost does not settle: it is 1.88 times at the first block past the edge on the smallest memory and continues upward, because every further increase makes a working set that fits less well. A rule that errs low has a floor under it. A rule that errs high does not.

And the whole shape moves with the machine, in the one direction a scan cannot follow. The plateau at M = 64 is nothing and at M = 256 it is eight blocks wide, so a block tuned by scanning on a small cache is tuned against a curve that does not exist on a large one. The scan finds the argmin of the machine it ran on; the model finds the feature that survives the move.

That is the division of labour the two numbers are for, and it is why printing only the survivor would have concealed it. What √(M/3) is good for is still the sentence the block size scales as the square root of the fast memory — the thing a person needs when the machine changes, and the thing no amount of scanning on one machine will produce.

Three lines of counting that get the exponent and miss the constant is the usual return on three lines of counting, and printing both numbers rather than the survivor is the whole of what this collection does differently.

Two other places this site has fitted rather than quoted #

Worth naming, because the habit is more visible across three instances than in one.

The fill growth exponent in the sparsity field: 1.044 for the matrix and 1.488 for its factor, fitted across grid sizes 6 to 12 rather than quoted from the complexity class the ordering is famous for.

The convergence order of the QR algorithm in the spectra field: 2.87, fitted from the recorded subdiagonal history, against a cubic convergence that is quoted everywhere and is asymptotic.

And the multigrid work exponent at 0.079 against conjugate gradients’ 1.00, in the iterative field — the number that says a V-cycle’s cost does not contain the grid.

In each case the quoted value is a limit and the fitted one is what happens at the sizes drawn, and in each case the two differ by an amount that is worth a sentence. The block size is the fourth, and it is the first where the form of the relation turned out to be wrong rather than only its constant.

What is left #

More than two levels. Every real machine has three or four caches plus main memory, and the optimal blocking is recursive — a block for each level, which is what a cache-oblivious layout gets for free without knowing any of the sizes. That is algorithms-data-structures.com’s subject and it has measured the van Emde Boas layout against a B-tree across seven block sizes; what would belong here is the numerical version, and this field does not have it.

Parallel communication, where the quantity is messages between processors rather than words between caches, the lower bound has a different shape, and the algorithms that attain it are genuinely different algorithms rather than reorderings. A reduction that changes the order is the closest this field comes: a QR that touches the matrix once instead of n times, and whose real difficulty is exactly the part left out there.

And the panel factorisation itself, which is unblocked by construction and is therefore the part of a blocked LU that does not block. On a wide enough matrix it becomes the bottleneck, and the standard answer is to block it recursively — which is where a blocked algorithm stops being a loop transformation and becomes a different algorithm.