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.

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. Stated in the previous essay and 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. Nobody should choose a block size from √(M/3): it is wrong by 1.56 and the right-hand cliff makes being wrong in that direction cheap and being wrong in the other direction expensive, so a model that under-predicts is at least failing safely. What √(M/3) is good for is the sentence the block size scales as the square root of the fast memory, which is the thing a person needs to know when the machine changes and which 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. The next essay 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.