Concept

Blocked algorithm — where it appears

A rearrangement of a factorisation that works on submatrices rather than on rows, performing identical arithmetic in a different order. It exists to turn memory traffic into cache reuse, and the best block size is a property of the machine rather than of the matrix.

Named by 9 essays across 3 fields — each of them below, with the objects they name alongside it.

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.

cost · Blocking
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.

cost · Blocking
10¹10³10⁵10⁷10⁹10¹¹10¹³10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classical once, Householdermodified once, Householderclassical twice, Householderclassical twice, Cholesky QRκ²uκuat κ = 10⁸classical once, Householder0.0042modified once, Householder5.7·10⁻⁹classical twice, Householder3.1·10⁻¹⁵classical twice, Cholesky QR10⁻¹⁵64×16 in blocks of 4, three seedsHouseholder inside does not help between

A stable block is not a stable basis

Block Gram–Schmidt orthogonalises twice over — between blocks, and inside each one. Householder inside the blocks does not stop the classical between-block step losing orthogonality like κ², 4.2·10⁻³ at κ = 4.3·10⁷, and a second pass does not stop Cholesky QR inside the blocks breaking down at κ = 10⁸. Each level fails only on ill-conditioning placed at its own level, and one variant holds 3·10⁻¹⁵ on every placement.

orthogonality · Gram–Schmidt
110¹10²10⁵block size bwords movedbest block: b = 10the recursion: no block size126,742 wordsbest block, scanned1.1·10⁵the recursion1.3·10⁵recursion ÷ best1.2M = 144 words; the recursion never reads it1.16× the best of 24 blocks

The recursion that was never told the memory

A blocked elimination has to be tuned to its fast memory, and tuned to one memory it costs up to 2.9 times the best at another. A recursive elimination splits the columns in half down to one and reads no memory size at all. On eight fast memories from 36 to 576 words it moves between 0.94 and 1.28 times the words of the best tuned block, with the same 585,200 operations and the same pivots — and on a machine with two caches it beats the block tuned to either cache on six machines of seven.

cost · Blocking
110¹123base case, columnswords ÷ the best block's√M − 2 = 10M = 144M = 144base of one1.2panels ≤ 10, worst1.2panels of 123the base case is a block size, rounded to a halvingand it has the block's cliff at √M − 2

The block size a recursion still has

A recursive elimination is sold as having no block size, and every real one switches to plain loops below some width. Swept over that width, the traffic is a staircase with its steps at the halvings of n, and its cliff sits where the blocked elimination's does — the first panel wider than √M − 2 moves 1.53 to 4.47 times the words, on six memories of six. On three caches the innermost decides, and a third cache costs every tuned block up to 14 per cent and the recursion nothing.

cost · Blocking
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

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.

machine · Algorithm selection
-12-10-8-6-4-2024681011.522.533.5leaf width minus (√M − 2)words ÷ pure recursionM = 64M = 144M = 256the edge is the same column at every memoryand the cheapest leaf sits on it

The leaf that sits on the edge

A recursive elimination's base case was found to be a block size in disguise, with a cliff where the blocked elimination's is, and the choice read as a trade: the processor wants wide leaves, the cache wants narrow ones. Measured at every width rather than at the halvings of 96, there is no trade inside the edge. Leaves exactly √M − 2 wide are the cheapest the recursion can have in words as well as calls — 0.81, 0.80 and 0.84 of the pure recursion's traffic at 64, 144 and 256 words — and one column wider moves 1.77 to 3.15 times it. And a matrix of 100 columns, which halves unevenly, meets the cliff in two steps rather than one.

cost · Blocking
M = 144, edge 10halving, worst size1.1aligned, worst size0.86961041121201280.70.80.911.1columns nwords ÷ pure recursionhalving to the edgesplit at the edgedashed: the pure recursion's wordsone base case, two ways to reach it

Leaves cut to the edge on purpose

A recursive elimination's cheapest leaf is exactly the square root of M, less 2, columns wide, and the rule drawn from it was to set the base case there and let the halvings put the leaves at or below it. On thirty-two sizes from 96 to 127 columns, halving to that base case moves more words than the pure recursion on sixteen of them at 144 words of fast memory, because the halvings stop at six and seven columns, not ten. Cut every dimension at a multiple of the edge instead and every size keeps the saving: 0.79 to 0.86 of the pure recursion's words, against halving's 0.91 to 1.07, with the same arithmetic and the same pivots. What it cannot make full is the one leftover leaf, and a leftover of one column is where it loses.

cost · Blocking
median departureone reflector at a time, formed7.5·10⁻¹⁵one reflector at a time, applied3.4·10⁻¹⁵blocks of 64, applied3.7·10⁻¹⁵10⁻¹⁵10⁻¹⁴block sizedeparture from orthogonality12481632642·10⁻¹⁵4·10⁻¹⁵6·10⁻¹⁵8·10⁻¹⁵Q formedapplied through the blocksopen dots: each of the five seedsblocking helps the formed matrix and not the operator

The factor nobody forms

A blocked Householder factorisation's orthogonal factor, multiplied out, departs from orthogonality half as far in blocks of sixteen as one reflector at a time, and that was read as blocking buying a factor of two. Libraries do not multiply it out. Applied to vectors through its stored blocks — which is how every caller uses it — the same factor departs by 3.1 to 4.0·10⁻¹⁵ at every block size from one to sixty-four, and stops growing after about twenty reflectors instead of adding them up. The factor of two was the price of forming the product, and a factor that is never formed never pays it.

orthogonality · Householder

Named alongside it

The objects these essays reach for when they reach for this one.

Memory hierarchyData movementFlop countLU factorisationRecursive factorisationCache obliviousLoop orderPartial pivotingCacheCommunication lower boundHouseholder reflectionLoss of orthogonality

All concepts