Which pairs are allowed to be small

The first four essays in this field are about one block. A matrix is not one block, and turning a fact about a pair of clusters into a way of storing an operator needs a rule that says which pairs.

The tree #

Number the points and split them in half; split each half in half; keep going until a piece is small enough to store densely. That is a binary tree over the index set, and every node of it is a cluster — a set of indices, and the set of points those indices stand for.

For a matrix of 256 unknowns with a leaf size of 16 the tree has four levels below the root and sixteen leaves. Nothing in building it looks at the matrix.

Each node carries the only three numbers the rest of this essay needs: the interval its points span, its centre, and its half-width. On a one-dimensional problem those are two minima and two maxima. In three dimensions they are a bounding box and its diameter, and everything below is unchanged.

The test #

Take two nodes, τ and σ, with centres c and half-widths r. The quantity

q = (r_τ + r_σ) / dist(c_τ, c_σ)

is the same q the first essay in this field derived: the expansion of the kernel about the two centres converges when q < 1, and its truncation at order p has relative error at most q^(p+1)/(1 − q).

The rule is one comparison:

A pair of clusters is admissible when q ≤ η, for a constant η the code chooses, and the block between them may be stored as two thin factors. Otherwise it is inadmissible, and it is subdivided.

η = 0.8 throughout this field, which is enough below one that the expansion converges with room and enough near it that most pairs qualify.

Three things about that rule are worth saying out loud, because each of them is doing work.

It is decided from geometry. Four numbers per pair, two subtractions and a division. No entry of the matrix is read, no decomposition is taken, and nothing about the kernel appears — which is why the same partition serves 1/r and log r and, disastrously, cos(κr)/r.

It is decided before anything is computed. The whole block structure of the matrix exists as a data structure before a single entry is evaluated, which is what makes it possible to allocate for and to distribute across processors. That is the same property the sparsity field’s symbolic phase has, and it fails in the same way for the same reason: a bound computed from structure alone is a bound about the worst arithmetic could do.

And it is recursive. A pair that fails the test is not stored dense; it is replaced by the four pairs of its children, which are half the size, at the same centres’ distance — so their q is roughly halved and three of the four now pass. Only the pair that still touches is subdivided again, and the recursion bottoms out on leaves.

What comes out #

At 256 unknowns, a leaf of 16, η = 0.8 and ε = 10⁻⁸:

Both reproduce the matrix to about 10⁻⁹. The first has more than twice the blocks and less than half the worst rank; the second has fewer of everything except rank, and stores fewer numbers.

That last column is the finding the next essay is about and it goes the opposite way from what the theory suggests. It is on this table because a partition figure with only the flattering columns on it would be an advertisement.

And the column that matters most — the rank — behaves completely differently on the two kinds of pair.

The admissible rank is the same integer at every size, and the touching rank is not. Across tolerances of 10⁻², 10⁻⁴, 10⁻⁶, 10⁻¹⁰ and 10⁻¹² the admissible pairs need 2, 3, 4, 6 and 7 — one number per tolerance, unchanged across a factor of eight in n — while the touching pairs need 3–5, 6–8, 8–11, 11–16 and 13–19, growing with the size at every tolerance.

That is what the test is buying, stated as a property rather than as a count of blocks. An admissible pair’s rank is a function of the tolerance alone — it depends on how many digits are wanted and not at all on how large the matrix is — and that is the whole basis of an O(n log n) method, because a cost that is (blocks) × (rank) only works if the second factor stops growing. A touching pair’s rank depends on both, so a partition that keeps touching pairs has a cost with the size in it twice.

The gap widens with the tolerance, too, which decides how much the test is worth: the ratio of touching rank to admissible rank at n = 256 runs 2.5, 2.7, 2.8, 2.7 and 2.7 across the five tolerances — near enough constant — so the test buys a factor of about 2.7 in the worst rank at any accuracy, and it buys the n-independence at all of them.

The picture, which is worth reading #

The strong partition has a shape and the shape is the argument.

Along the diagonal there is a fringe of small filled squares. Those are the pairs the test refused all the way down: a cluster with itself, and a cluster with its immediate neighbour, whose q is 1 or worse at every level. They are kept dense, they are small, and there are O(n/leaf) of them — which is why the dense part of a hierarchical matrix costs O(n) and not O(n²).

Away from the diagonal the blocks get larger as they get further out, which is the geometry doing exactly what it should: two clusters far apart may be large and still separated. The largest blocks in the picture are 64 × 64 and there are six of them, which is a smaller number than the picture suggests and is worked out below.

And every compressed block carries one of two small numbers. There are 66 of them, twenty-two at rank four and forty-four at rank five, because the test guarantees q ≤ 0.8 and the rank a block needs is a function of q and ε. That uniformity is the thing the test buys, and it is what the previous essay’s flat row becomes when applied to a whole matrix rather than to one pair.

The other rule, and why anybody uses it #

The second row of the table is not a straw man. It is the format most codes that solve with a hierarchical matrix actually use, and there are two reasons.

The first is that it is trivial. At every level the matrix is [[A₁₁, A₁₂], [A₂₁, A₂₂]]; compress the two off-diagonal blocks; recurse on the two diagonal ones. There is no test, no bounding box and no η. On a one-dimensional index set it produces a correct representation with no geometric information at all beyond the ordering.

The second is that it is what a solve needs. At each level the off-diagonal part of a weakly admissible matrix is a single low-rank correction to a block-diagonal matrix, so the Woodbury identity turns a solve into two solves one level down plus a small dense one. A strongly admissible partition has no such structure — its off-diagonal part at a given level is a mosaic of blocks at several sizes — and inverting it is a substantially harder piece of software.

That is a recurring shape in this collection and it has a name here: a format is chosen for what is going to be done to it. The sparsity field’s version is that an ordering good for fill is not necessarily an ordering good for parallelism; the cost field’s is that a factorisation with the right flop count can have the wrong memory traffic.

What η does #

The constant is a knob and it is worth knowing which direction each way turns.

Lower η is a stricter test. Fewer pairs qualify, more are subdivided, the partition has more blocks and each compressed one has a smaller rank because its q is smaller. In the limit η → 0 nothing is admissible and the matrix is stored densely.

Higher η is a looser test. More pairs qualify, at larger q, so the ranks grow. At η = 1 the touching pairs qualify — and the expansion at q = 1 does not converge, so their ranks are the ones that climb with the sampling. Past η = 1 there are pairs whose clusters overlap and the whole justification is gone; the code refuses those outright rather than producing a representation nobody can bound.

So η trades block count against rank, which is the same trade the two rows of the table make in a cruder way, and what that trade costs against what it saves is the next essay; the reason it has no obvious optimum is that the two resources are not commensurable. More blocks cost bookkeeping and worse locality; larger ranks cost arithmetic that scales like the cube of the rank in any formatted operation. Which dominates is a property of the machine and not of the matrix, which is the cost field’s standing observation arriving here.

What the test cannot do #

Two limits, and the second is the one that produces the next essay but one.

It cannot see the kernel. q is computed from points. A pair of clusters at q = 0.5 is admissible whether the kernel between them is 1/r, log r or cos(40r)/r, and the previous essay is what happens in the third case: the test passes, the block is compressed, and the rank it takes is fifty-three instead of six. The test is a necessary condition for the expansion to converge and it is not a sufficient one for the expansion to converge quickly.

And it cannot see the numbering. The tree is built by splitting the index set in half, and the geometry attached to each node is read from the points those indices stand for. If the indices are numbered so that a contiguous run of them is a scattered set of points, every node of the tree spans the whole domain, every half-width is the whole domain’s, every distance between centres is nearly zero, and no pair is admissible anywhere. The test does not fail in that case — it correctly reports that there is nothing to compress — and the representation degenerates to dense storage exactly.

The second limit is a feature dressed as a bug and the essay after next takes it seriously, because what it means is that the compressibility of a matrix lives in the numbering rather than in the matrix, and a numbering is a thing somebody chose.

Counting the blocks, which is where n log n comes from #

The storage claim in this field is usually stated as an asymptote and it is worth deriving once, since the derivation is four lines and explains the shape of the picture.

At level ℓ of the tree there are 2^ℓ clusters, each of n/2^ℓ indices. The pairs that are admissible at level ℓ are the ones whose parents were not, and for a one-dimensional index set the number of those is bounded by a constant times 2^ℓ — each cluster has a bounded number of neighbours-of-neighbours, independent of how deep the tree is. Each such block is (n/2^ℓ) × (n/2^ℓ) stored in rank k, which costs 2k·n/2^ℓ numbers.

So each level costs a constant times k·n, there are log₂(n/leaf) levels, and the dense fringe on the diagonal costs n·leaf. Total:

O(k · n · log n) + O(n · leaf) .

Two things fall out of that. The first is that the leaf size is a knob on the dense part only, so a larger leaf costs more numbers and fewer levels — and this field’s later measurement is that it costs nothing in accuracy either way, which is not obvious in advance. The second is that the whole derivation needs the rank to be a constant, which is the third essay’s null result and the fourth’s counterexample, and needs the number of admissible blocks per level to be bounded, which is a property of the geometry rather than of the matrix.

That second condition is the one that gets harder in higher dimensions. In three dimensions a cluster has more neighbours, the constant is larger, and the practical crossover moves — which is why the sizes at which a hierarchical format starts to pay are different for a boundary element code and a volume one.

The counts are exact, which the derivation does not claim #

The four lines above are an order argument: a bounded number of blocks per level, times a rank, times a size, summed over levels. Counted rather than bounded, on this ordering with leaf 16 and η = 0.8, the bound is not merely satisfied — it is attained, exactly, at every size.

and the dense fringe comes out at 10, 22, 46 and 94 blocks for n = 64, 128, 256 and 512, which is 3n/leaf − 2 in each case. Neither of those is an asymptotic statement or a fit. They are the counts, at four sizes spanning a factor of eight, and they are integers that agree.

So the constant the derivation leaves unnamed is 3: three admissible partners per cluster, at every level below the second, forever. The storage figure can therefore be written with its constant in place — about 6·k·n·log₂(n/leaf) for the compressed part and 3·n·leaf for the fringe — which is worth having, because an order statement cannot tell anybody whether a hierarchical format is cheaper than a dense one at the size they actually have.

The reason it is exactly three is geometric and one line. Clusters at a level sit on a line at equal spacing; two of them at k cluster-widths apart have q = 1/k; the test admits k ≥ 2. A cluster’s partners at k = 2 and k = 3 are new at this level, its partners at k ≥ 4 were already admissible as part of a larger block one level up, and its partner at k = 1 is the neighbour that never qualifies. Two new partners on each side, minus the ends, is 3·2^ℓ − 6.

Two ranks, in a fixed proportion #

The same count explains why the ranks are not one number. The k = 2 partners sit at q = 0.492 and the k = 3 partners at q = 0.328, so they are not equally compressible, and at ε = 10⁻⁸ that is the difference between five columns and four. Across the whole partition:

Exactly one to two, at every size. That is a stronger statement than the uniformity claimed above and a weaker one than “the same small number”: the ranks are two values rather than one, and which value a block takes is decided by its q, which is decided by the test, which reads four numbers. Nothing in the partition is a spread; it is two populations.

It also settles what the largest block is. The two halves of the index set have q = 1 exactly — the refusal this essay ends on — so no 128-row block is ever admissible, and the widest compressed block in a 256-square is 64 × 64. There are six of them, being the pairs among the four quarters at two or three quarter-widths apart. The pair between the first and last quarters is at q = 0.328 and costs four columns, which makes it the cheapest block in the picture rather than the dearest; the five-column blocks are the nearer pairs. assertTheBlockCountsAreExactAndTheRanksAreTwo holds all of that, including the two integer identities, so a change to the tree or to η that broke the count would be reported rather than absorbed.

The relative in the sparsity field #

The sparsity field has a phase that decides a structure before any arithmetic happens and it is called symbolic factorisation. It reads the graph of the matrix, works out which entries the elimination will fill in, and allocates for them — and this collection has an essay about the fact that what it produces is a bound, because an entry it predicts may turn out to be numerically zero and a pivot choice may change the graph underneath it.

The admissibility test is the same kind of object with the same kind of caveat. It decides a structure from geometry, in advance, cheaply, and what it decides is where compression is allowed, not what it will cost. The rank of each block is discovered afterwards, and a code that allocated from the bound rather than from the outcome has over-allocated by the factor of six the second essay in this field measures.

The two phases are worth reading together because they fail differently. A symbolic phase over-predicts fill and wastes memory. An admissibility test over-predicts rank and wastes memory. Neither ever under-predicts, which is the property that makes both usable, and neither is tight.

The refusal #

The claim under test is the one the whole test exists to deny: that two clusters of points which are next to each other interact weakly enough to compress.

For [0, 1] and [1, 2] the half-widths are ½ and ½ and the centres are 1 apart, so q is exactly 1. The assertion that q < 1 for that pair is fed that number and it fails.

Exactly one, and not more than one, is the whole point. Nothing about the pair is pathological; it is the most ordinary pair of adjacent intervals there is, and it sits precisely on the boundary where the expansion stops converging. Every partition in this field spends most of its blocks dealing with that one case, and a rule that admitted it would be a rule with a slowly growing rank in the one place a partition has the most blocks.