Built from products alone

An operator with no entries ends on a sentence this essay contradicts: an operator that can only be applied is an operator that cannot be factorised, and every solver in that essay is an iteration for that reason.

A hierarchical representation is the exception, and the reason is structural rather than clever. Every block of it is a submatrix, and a submatrix can be sampled.

The sampling #

Take a block A(τ,σ) — rows in cluster τ, columns in cluster σ. Build a random matrix Ω that is nonzero only on the rows σ, and apply A. Then

(A Ω)|_τ = A(τ,σ) Ω_σ ,

which is the block applied to a random matrix, which is exactly what the randomised field’s range finder turns into an orthonormal basis for the block’s range.

Once a basis Q for the range is in hand, the second factor needs one more product. Put Q’s columns into n-vectors supported on τ, apply A, and read the rows σ:

(A Q)|_σ = A(σ,τ) Q = A(τ,σ)ᵀ Q (for symmetric A),

so A(τ,σ) ≈ Q (A(τ,σ)ᵀQ)ᵀ. Two batches per block and no entry read — the same sketch that is spent after one use, drawn fresh for each.

The symmetry is used once and deliberately. On a non-symmetric operator the second step needs a product with Aᵀ instead, the count doubles, and nothing else changes — which is worth knowing and is why it is said rather than assumed.

The count #

The blocks at one level of the tree have disjoint column supports, so one batch serves all of them at once. Two batches a level for the two triangles, two more for the projections, and one batch of leaf products recovers every diagonal block in the matrix together, because the leaves are disjoint too. So

products = leaf + 2(2k + p)·levels ,

with k the rank asked of each block and p the oversampling. At leaf 16, k 10, p 4 and four levels that is 16 + 2(24)(4) = 208, and the counter inside the operator says 208.

Forty-eight more products per doubling, which is 2(2k + p) exactly, and the entries column squares. Products per unknown falls at every doubling: the construction gets cheaper per degree of freedom as the problem gets bigger, which is the same shape the storage claim has and for the same reason.

The subtraction that makes it correct #

There is one step in the construction whose omission produces an answer that looks like an answer, and it is worth the section because it cost this build a day.

One batch samples every block at a level at once, because their column supports are disjoint. But reading rows τ of the result picks up the contribution of every cluster in the batch’s support, not only the one paired with τ — and those extra contributions are the blocks at coarser levels, which have already been recovered.

So each batch has the coarser levels peeled off it before the range is taken. Without that, the representation comes out at a relative error of 0.25 with every intermediate quantity looking correct: the bases are orthonormal to 8·10⁻¹⁶, the projections are projections, the counts are right, and the answer is wrong.

There is a second version of the same mistake one level finer, and it is the one that survives a first fix. A block recovered by the other triangle of the same level is the transpose of the one being recovered now — and the projection step reads exactly it — so peeling against it subtracts the answer. The subtraction has to be against completed levels only, which is one line and is the difference between 0.25 and 4·10⁻¹³.

Both are instances of a shape this collection keeps meeting: an error whose symptom is a plausible number. Every gate this site has asks whether a quantity is finite, whether a label fits, whether an assertion held. None of them asks whether the thing being sampled is the thing intended, and the only check that caught this was a comparison against the representation the entries give.

The error that grows with n is the fixed rank #

The table above holds k = 10 and walks n, and the error climbs — 8.2·10⁻¹¹, 5.6·10⁻⁹, 6.7·10⁻⁸, 4.0·10⁻⁷ — while the products per unknown falls. Read as a property of the construction that is an unpleasant trade: cheaper per degree of freedom and less accurate at every doubling.

It is a property of the rank. Holding n at 512 and walking k instead, at a leaf of 16:

k = 6 gives 176 products and an error of 2.6·10⁻⁴. k = 10 gives 256 and 4.0·10⁻⁷ — the table’s row. k = 14 gives 336 and 4.0·10⁻¹⁰. k = 16 gives 376 and 8.7·10⁻¹².

And with a leaf of 32, which allows a larger rank: k = 10 gives 224 products and 3.3·10⁻⁷, k = 16 gives 320 and 7.2·10⁻¹², and k = 22 gives 416 products and 1.0·10⁻¹⁵.

The working precision, from products alone, at 416 products against 262,144 entries — still 0.79 products per unknown. So the construction is not less accurate at scale. A fixed rank is, and the size sweep varies the wrong thing to see it.

And the excess is under-ranking, not sampling #

The section above prices the construction against the representation the same partition gives when every block is decomposed from its entries, and calls the ratio the cost of sampling. Watching it while the rank rises says it is not.

At n = 512 and leaf 16 the excess runs 8.0, 7.3, 6.5, 4.8 at k = 6, 10, 14, 16. At leaf 32 it runs 6.0, 4.0, 1.0, 1.0 at k = 10, 16, 22, 28.

It goes to one. At an adequate rank the representation built from products is as accurate as the one built from every entry — not within a factor, not quasi-optimally, but equal to two digits. So the excess is a symptom of a starved rank rather than a price paid for randomness: where the range finder has enough columns to find the range, it finds it, and the sampling costs nothing in accuracy at all.

That is a stronger claim than the field usually gets to make about a randomised method, and it is worth being exact about why it is available here. The randomised field’s bounds are statements about a truncation at rank k of a matrix whose spectrum continues past k, and the excess they describe is the tail the sketch fails to capture. Here the blocks are numerically rank deficient — that is the whole premise of the partition — so once k passes the block’s numerical rank there is no tail left to miss, and the two routes are computing the same object.

Which relocates the decision the construction actually asks. It is not how much accuracy does sampling cost, because the answer is none. It is how large must k be, which is the rank that has to be guessed — and the cost of guessing high is 2(2k + p) products per level, linear in k, against an error that falls by decades. Guessing high is cheap and guessing low is not, and the sweep above is the price list.

What it costs #

The construction is not free and the price is on the same figure. Against the representation the same partition would give if every block were decomposed at the same rank:

A factor of three to seven, widening slowly with the depth. The widening is the peeling: every level’s data has the coarser levels’ approximations subtracted from it rather than the coarser levels themselves, so a coarse error propagates into every finer level.

That makes this the one accumulation in this whole phase that accumulates. The essay next door measures ninety-eight truncations inside a factorisation costing nothing measurable; this measures five levels of peeling costing a factor of seven, and the difference is that one of them is a projection onto a fixed target and the other is a subtraction of an approximation from a measurement.

Seven is a factor rather than an order of magnitude, which is what keeps the route usable. It is also the thing to watch on a deeper tree, and it is the argument for the variants of this construction that re-sample rather than peel.

The trade the compression route does not make #

compress is handed an accuracy and returns whatever rank it costs. This is handed a rank and returns whatever accuracy it buys, because the random matrix has to be drawn before a single entry of the block has been seen.

A straight line at 0.73 decades a column. So the guess is a guess about a number of digits, being two columns short costs about a decade and a half, and there is no cliff — which is the same flatness the hierarchy field’s second essay measures from the other direction and is what makes guessing tolerable.

And the trade improves as the rank rises, which is the opposite of what a cost table suggests.

Seven orders of accuracy for 2.4 times the products, and the penalty falls as they are bought. Across k = 4, 6, 8 and 14 the error at the largest size runs 7.8·10⁻³, 2.7·10⁻⁴, 10⁻⁵ and 4·10⁻¹⁰, the products per unknown run 0.27, 0.34, 0.42 and 0.66, and the excess over the compression route runs 10.7, 8.0, 7.1 and 6.5.

So the two costs move in opposite directions and the ratio is the one that matters. Absolute work rises — 0.27 to 0.66 products an unknown — while the relative penalty for never assembling the block falls by 40%, because the compression route has to form the block whatever rank it then finds, and that cost is fixed while the sampling cost is not. A code choosing between the two on a low-rank benchmark is choosing at the worst point for the sampling route, and the choice gets better for it at exactly the ranks where accuracy is actually wanted.

It is the same trade the randomised field already has an essay about. the-sketch-that-is-not-the-answer measures what a sketch buys and what it gives up, and the answer there is that the sample size is chosen against a target the sample cannot see. Here the quantity being guessed is a storage decision rather than an answer, which is a milder failure: guess low and the representation is inaccurate, which is visible; guess high and it is expensive, which is also visible. Neither is silent.

The oversampling, and why the figure is a band #

Every number in this construction is a draw, so a line would be a claim about one run.

Both columns matter and only one of them is visible in a single run. Oversampling buys a better answer and a more predictable one, and the second is the half a line would hide — at p = 0 the spread across five seeds is 44 per cent of the mean, which means a code that ran once and reported its error has reported a number that could have been half again as large.

And no draw, at any p, beats the flat line. That is the honesty check: the flat line is the best representation of that rank, a sample cannot beat the decomposition it approximates, and an experiment in which one did would be an experiment with a bug in it.

Why the blocks can be sampled at all #

The step that makes this work is easy to skip past and it is the one that distinguishes a hierarchical matrix from every other structured format on this site.

Every block of the partition is a submatrix of A. Not a transform of one, not a projection of one — the literal set of entries Aᵢⱼ for i in τ and j in σ. So a product with A, restricted to the right rows and columns, is a product with the block, and the randomised field’s machinery applies without modification.

Nothing else in this collection has that property. A circulant’s structure lives in a transform, and it cannot sample the diagonal of F A F* by applying A to a vector supported on part of the index set, because the transform mixes everything. A sparse factorisation’s structure lives in an elimination order, and the factors are not submatrices of anything. A Krylov method’s structure lives in a subspace that does not exist until the products have been done.

That is why this essay exists in this field and not in one of those. The randomised range finder is a general tool for finding the column space of a matrix that can be applied; what a hierarchical partition supplies is a set of matrices that can be applied by restricting products with A, which is the one ingredient the tool needs and almost nothing has.

There is one place the property is nearly lost and it is worth flagging. The peeling above works because a block at a coarse level is also a submatrix, so its action can be subtracted exactly. In the strongly admissible partition — where a level’s off-diagonal part is a mosaic of blocks at several sizes rather than one clean pair — the same subtraction is possible and considerably more bookkeeping, which is a third reason, after storage and after the solve, that the weakly admissible partition is the one codes build.

What this makes possible #

The construction is the piece that turns everything else in this phase from an observation into a method, and it is worth being explicit about why.

The hierarchy field’s storage and solve results are all measured on representations built by decomposing every off-diagonal block. That reads every entry, costs O(n³) in the worst case, and on the problems this format exists for is the one thing nobody can afford. A field whose central object can only be built that way has described a compression scheme for matrices that have already been formed.

With this, the input is a routine. A boundary element code that can apply its operator by any means — a fast multipole method, a coarse quadrature, another hierarchical matrix — can hand that routine over and get back a representation it can factorise. That is the step from an iteration to a direct method, and it is why this construction is the one part of the phase with a whole essay in a different field.

The three quantities a code has to choose #

Nothing in this construction is adaptive, and a reader deciding whether to use it needs the three numbers it asks for and what each of them costs to get wrong.

The rank k. Chosen before anything is known, at 0.73 decades a column. Too low and the representation is inaccurate by a knowable amount; too high and it is expensive by a knowable amount. Neither failure is silent, which is the property that makes the guess acceptable, and an adaptive version — sample a few columns, check whether the residual has fallen enough, sample more — exists and costs an extra batch per block per round.

The oversampling p. Chosen for the spread rather than for the mean. Four is the conventional answer and this measurement supports it: p = 4 costs 24 more products at the largest size and takes the spread across seeds from 44 per cent to 13.

The leaf size. It sets the levels, and the levels appear in the product count linearly — but it also appears additively, as one batch of leaf products for the diagonal blocks. At 512 unknowns with a leaf of 16 that additive term is 16 of the 256 products; at a leaf of 64 it would be 64 of 208. So there is an optimum in the count, it is shallow, and it sits near the leaf sizes a code would choose for other reasons anyway.

The thing none of the three controls is the factor of seven against the decomposition. That is the peeling, it grows with the depth, and no setting of any of the three parameters removes it — which is the honest form of what never seeing the matrix costs.

The refusal #

The claim under test is the sentence an-operator-with-no-entries ends on and this essay contradicts: that a matrix that can only be applied is a matrix that cannot be factorised.

The refusal fed to the assertions is the opposite mistake, and it is the one worth guarding against here. A construction that produced a better representation than the decomposition of the same rank would be a construction with a bug in it — a sample cannot beat the optimum it is sampling — so the assertion that the sampled error is below the decomposition’s is fed 6.65·10⁻⁸ against 1.01·10⁻⁸, and it fails.

The claim above is refuted by the table rather than by an assertion, which is worth saying plainly. Two hundred and fifty-six products, a representation at 4·10⁻⁷, and no entry of the matrix ever formed. What the assertions guard is the direction of the comparison, which is where an experiment of this kind goes wrong.