Where the format starts paying

Every asymptotic claim in numerical linear algebra is a claim about sizes nobody has. This essay is about the size somebody does have, which is a different question with a different answer — and it is the question a block size is a property of the machine asks one field along, where the constant that decides it is tuned rather than derived.

What the recursion is #

At the top level, a hierarchical matrix in the weakly admissible partition is

A = [ A₁₁ U₁V₂ᵀ ] = D + Ũ Ṽᵀ , [ U₂V₁ᵀ A₂₂ ]

with D block diagonal and the off-diagonal part of rank k₁ + k₂. The Woodbury identity says

A⁻¹ = D⁻¹ − D⁻¹Ũ (I + ṼᵀD⁻¹Ũ)⁻¹ ṼᵀD⁻¹ ,

so a solve with A is a solve with D against k₁ + k₂ + 1 right-hand sides, plus one dense solve of size k₁ + k₂ — and a solve with D is two solves one level down. The recursion bottoms out on dense leaf blocks. Nothing anywhere in it forms an n × n matrix — which is the inverse that is never formed’s discipline applied to a whole factorisation — and nothing in it is approximate: the answer is the exact solution of A_H x = b, and how far that is from the answer to Ax = b is the previous field’s measurement.

Counting what it costs is not a modelling exercise. The recursion carries a counter, every multiplication it performs increments it, and the numbers below are what the code did.

The measurement #

The last column is what makes the comparison legitimate. Every one of those runs returns an answer at the same backward error, so the two methods are being asked for the same thing.

The ratio column has a crossover in it, between 64 and 128 unknowns. Below it the format is the more expensive way to solve the system. Above it, and increasingly, it is not.

The exponent, which is the interesting number #

A cost of n log²n has no exponent. What it has is an exponent that falls towards one and never arrives, and the way to see that is to take it between consecutive sizes rather than fitting it over a sweep:

Against the dense factorisation’s 3.00 at every step, exactly, because it is a cube.

Three things are worth reading off that column. The first is that it is well under three throughout, which is why the ratio improves. The second is that it is falling, which is what a logarithm looks like from inside — at n = 512 there are five levels, so log₂ n contributes about a fifth of the answer, and it will contribute a tenth at n = 2¹⁰ and a twentieth much later. The third is that at n = 64 it is above two: at the small end this is not merely slower than the dense factorisation, it is growing faster than a quadratic, because the recursion’s overhead is concentrated in the small dense solves at the leaves and there are more of them per unknown when the tree is shallow.

Reporting a fitted single exponent over the whole sweep would give 1.43 and would say none of that. It would also be a number with no meaning, since the quantity it purports to describe has no exponent.

What the accuracy does to the crossover #

The dense line is fixed. A factorisation has no tolerance in it, so n³/3 is n³/3 at every accuracy, and the whole trade lives in the other curve.

Tightening ε raises the ranks, raises the cost of every level of the recursion, and walks the crossover right. That is the economics of the format in one sentence: what the accuracy knob buys is a size at which the format starts paying, and asking for more digits than the problem needs pushes that size out of reach.

The sizes are worth writing down rather than describing:

Each pair of decades in the tolerance multiplies every ratio by about 1.5, and does it uniformly across the four sizes — 0.61 to 0.92 at the smallest, 0.06 to 0.10 at the largest. The tolerance moves the whole curve rather than tilting it, so it changes where the crossover sits without changing how steeply the format wins once past it.

The dense column is the same four numbers at every tolerance, which is the fixed line the section opens with, measured: 8.74·10⁴, 6.99·10⁵, 5.59·10⁶ and 4.47·10⁷ at ε = 10⁻⁴, 10⁻⁶ and 10⁻¹⁰ alike. So every change in the ratio column is the hierarchical curve moving, and it moves by a factor of 3.2 across those six decades — 5.3·10⁴ to 1.72·10⁵ at the smallest size — while the ranks move by a factor of a little over two.

The rank grows like log(1/ε) — four to sixteen across ten decades — which is the format’s defining property, measured here rather than assumed. And the crossover grows faster than the rank: a rank three and a half times larger moves it eight times further out, so it goes roughly as k^1.7. Tightening the tolerance costs the crossover more than it costs the rank, which is the sharper form of the sentence above.

The leaf, and a caution about the counter #

The recursion’s leaf size is fixed at sixteen throughout and never swept, and the section above attributes the small-n overhead to it. Sweeping it, the counter says bigger is always better — at n = 512 the count falls monotonically from 9.97·10⁶ at a leaf of 4 to 3.92·10⁶ at 256 — and that is an artefact of the counter rather than a property of the method.

The counter charges a leaf L² per right-hand side, which is the triangular solve, and never charges the leaf’s L³/3 factorisation, because the build takes no counter. At a leaf of 16 the omission is 4.4·10⁴ against 7.2·10⁶ — 0.6%, so the table at the top of this essay is unaffected. At a leaf of 128 it is 2.8·10⁶ against 4.7·10⁶, which is 59%.

Adding the leaf factorisations back, at n = 512 and ε = 10⁻⁸:

An interior optimum at a leaf of 64, twenty-one per cent below the sixteen used throughout, and a sharp rise past it. So the leaf is a real knob worth about a fifth, and the reason it looks like a free lunch is a counter that omits half of what a large leaf costs.

One thing that does not settle. The compressions — an SVD of each off-diagonal block — are uncounted too, and they are the other half of the build. So the table at the top of this essay compares a hierarchical solve against a dense factorisation: the right comparison when the same matrix is solved against many right-hand sides, and one that flatters the format for a single solve. That asymmetry is the same one the same arithmetic at a different price has to be careful about, where two eliminations performing identical arithmetic are separated by what is counted rather than by what is done. And the distinction between counting words and counting rounds, which the message and the word is entirely about, is the same lesson a third time: a cost model reports what it was built to report, and the quantity it leaves out is invisible rather than small.

The practical form of all three is one sentence. Before believing a cost comparison, ask what each side counts — and the answer here is that one side counts a factorisation and the other counts a solve, which is exactly right for a code that reuses the matrix and exactly wrong for one that does not.

An accuracy that is a backward error is what closes this. The accuracy is a backward error and the digits it should be set to are the digits the answer needs divided by the condition number — so a well-conditioned problem needs a loose tolerance, gets small ranks, and crosses over early. An ill-conditioned one needs a tight tolerance, gets larger ranks, and crosses over late, which is a second and quite separate way in which conditioning costs money.

Where the multiplications go #

The counter reports one number and the number is a sum of four quite different things, which is worth separating because three of them behave well and one does not.

The leaves. The fill they avoid is what the factor is not sparse prices for a sparse elimination. At the bottom of the recursion sit n/leaf dense solves of size leaf, each performed against every right-hand side that reaches it. That is where the small-n overhead lives: at n = 64 with a leaf of 16 there are four of them and the tree is two levels deep, so the fixed cost of the recursion is spread over very little work.

The block products. The rank they are applied at is the quantity a rank that is a number of digits is about. Each level applies the thin factors Ũ and Ṽᵀ, which costs O(k·n) per right-hand side per level. Linear in the rank, linear in the size, and there are log n levels — this is the term the whole format is built around.

The small dense solves. At each level the Woodbury correction requires factorising (I + ṼᵀD⁻¹Ũ), which is (k₁ + k₂) square. That is a cube in the rank and a constant in n, which is why the rank matters more here than the storage tables suggest.

And the right-hand sides multiply. This is the term that surprises. Solving with A against one vector requires solving with D against k₁ + k₂ + 1 vectors, and solving with those requires solving one level down against (k₁ + k₂ + 1) times the next level’s rank plus one. The number of right-hand sides at depth ℓ grows like the product of the ranks above it. For small ranks that is fine; for the ranks a touching block acquires at 512 unknowns it is the reason the exponent is 1.74 rather than closer to one.

That last term is the strongest argument in this whole phase for the admissibility test the hierarchy field measures as more expensive. A rank that climbs from 9 to 13 costs 44 per cent in storage and considerably more than that here, because it appears in a product along the depth of the tree.

What a multiplication count is not #

This is the cost field, and the cost field’s standing result is that an operation count stopped predicting the time some decades ago. It is worth applying that to the figure above rather than letting the figure stand as a runtime claim.

The dense factorisation is the best-behaved kernel in scientific computing. It is three nested loops over contiguous memory, it blocks perfectly, and a tuned implementation reaches a large fraction of a machine’s peak rate. Its 4.47·10⁷ multiplications at n = 512 are cheap multiplications.

The hierarchical solve is not that. It is a recursion over a tree, its leaves are small dense solves that do not block, its intermediate objects are thin factors with poor locality, and it spends a significant fraction of its time on bookkeeping that this counter does not count at all. Its 7.23·10⁶ multiplications are expensive multiplications.

So the ratio of 0.16 at n = 512 is a ratio of two counts and not a ratio of two runtimes, and the runtime ratio is certainly worse. How much worse is a property of the machine, which is this field’s oldest observation: the-same-arithmetic-at-a-different-price measures two orderings of one factorisation with identical flop counts and a memory traffic differing by a factor that decides which is usable.

The honest form of the claim is therefore about the exponent rather than the constant. A method whose exponent is falling towards one beats a method whose exponent is three, at some size, whatever the constants are — and where that size is, in seconds rather than in multiplications, is a question no count answers.

Why the exponent is not one #

The literature’s headline for this family of methods is linear complexity, and the exponent column above says 1.74 at the largest size measured. Both are correct and the gap between them is worth closing.

The n log n storage argument in this field’s fifth essay needs the rank of every block to be a constant. The solve needs more than that: it needs the number of right-hand sides propagating down the recursion to stay bounded, and that number is a product of ranks along a path from the root to a leaf. With a bounded rank k and L levels it is O(k^L) in the worst case and O(kL) in practice, because the correction at each level acts on a subspace the level below has already reduced.

What the measurements here have is a weakly admissible partition whose ranks are not constant — 9, 10, 12, 13 across the sweep, as this phase’s third essay establishes — so the product along the path is growing at both ends. The exponent is therefore an honest 1.74 for this partition at these sizes, and it would be lower for the strongly admissible one whose ranks are five throughout, if that partition had a comparably simple solve, which it does not.

So linear is a statement about a limit, about a constant rank, and about a formulation with recompression at every level of the recursion rather than the plain Woodbury one used here. Each of those is a real thing and none of them is what the figure measures. What the figure measures is the exponent of a specific, simple, exactly-arithmetic recursion, and 1.74 falling is what that is.

Against the other things available #

The comparison on the figure is against a dense factorisation, which is the fair one — both are direct methods returning an answer at a stated backward error. Two other options are worth pricing.

An iterative solve on the dense matrix. Conjugate gradients at n = 512 on this matrix takes about 50 steps at n² multiplications each, which is 1.3·10⁷ — twice the hierarchical solve and three times cheaper than the dense factorisation. That is a reminder that a well-conditioned matrix rarely needs a factorisation at all, and the whole comparison above is between two methods that earn their keep when there are many right-hand sides or when the conditioning is bad.

A hierarchical solve used as a preconditioner. The previous essay measures that: at 128 unknowns the best preconditioned run costs 4.81·10⁵ multiplications against the direct hierarchical solve’s 5.66·10⁵, so the preconditioner wins slightly at that size. At 512 it does not, because the direct solve’s exponent is 1.74 and the preconditioned iteration’s is 2 — each of its steps is a dense product with A.

Which is a pleasing symmetry to end on. The preconditioned iteration is cheaper on small problems and the direct hierarchical solve is cheaper on large ones, for the same reason that a dense factorisation is cheaper than either on very small ones: at every size the winner is the method whose overhead matches the problem, and the overheads sort in the opposite order from the exponents.

The comparison the figure cannot make #

There is a size argument missing from all of this and naming it is more useful than pretending otherwise.

The dense factorisation at n = 512 needs 512² = 262,144 numbers of memory and it needs them all at once. The hierarchical representation at the same size needs 67,968. On this machine both are trivial; on the problems this format exists for, neither is, and the memory is what decides.

That is why the storage figure in the hierarchy field crosses over at a much smaller size than this one. Storage is already in the format’s favour at n = 64 — 3,456 numbers against 4,096 — while the arithmetic is still 48 per cent against it. So there is a band of sizes, roughly 64 to 128 here, where the format costs more arithmetic and less memory, and which of those a code should prefer is not a question about the matrix at all.

A second thing the figure cannot show is what happens with many right-hand sides. Both methods amortise a setup: the dense one factorises once and back-substitutes per vector at n² each; the hierarchical one builds a representation once and applies the recursion per vector at roughly k·n·log n each. So for many vectors the comparison is between n² and k n log n per solve, which crosses over far earlier than the figure’s crossover and in the format’s favour. Everything on the figure is the one-right-hand-side column.

And the third is the setup itself. The hierarchical representation on the figure was built by decomposing every off-diagonal block, which costs more than the dense factorisation it is being compared with — a fact this phase’s essay on building it from products alone exists to fix, and which the figure quietly assumes away by pricing only the solve.

The refusal #

The claim under test is the one every abstract about this subject makes: that a hierarchical solver is faster than a dense factorisation.

The verdict is asymptotic rather than false, because the claim is true and the word missing from it is eventually. The refusal is fed the ratio column: the assertion that the hierarchical solve costs less at every size on the sweep is handed 1.48, 0.81, 0.39, 0.16, and it fails on the first entry.

At 64 unknowns the format performs half again as many multiplications as the thing it replaces, for an answer that is no better, using a representation that took a decomposition of every off-diagonal block to build. It is comprehensively the wrong choice, and it is the wrong choice at a size that is perfectly ordinary — a 64 × 64 dense system is the sort of thing that appears inside larger algorithms by the thousand.

That is why the crossover is worth measuring rather than asserting, and why a figure of this kind belongs in the cost field rather than in the field that owns the method. A method’s own field is the last place a comparison it loses is likely to be drawn.