Numerical Linear Algebra
where the algebra stops being the arithmetic
Neither sparse nor dense

Which pairs are allowed to be small

A hierarchical representation is a partition of the matrix into blocks, and the rule that produces it reads four numbers per pair of index clusters and not one entry of the matrix. On a 256-square it yields 112 blocks, 66 of them stored as two thin factors, none of rank above five.
34 min read 25 figures Where the algebra stops being the arithmeticMeasured, not assertedAssertions that reject

Worth reading first: A block nobody can call sparse · The order decides the memory · The factor is not sparse.

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.

A 256 × 256 kernel matrix partitioned by the strong rule, with each compressed block's rank112 blocks: 46 kept dense and 66 stored as two thin factors, whose ranks run from 4 to 5. The rule decides one thing — whether a pair of index clusters may be compressed — and it decides it from four numbers about where those clusters sit, before any entry of the matrix is read. The strong rule refuses any pair whose clusters touch, and subdivides instead, so the diagonal is fringed with small dense blocks and no rank on the picture exceeds 5. The whole thing stores 27,008 numbers against 65,536 entries, and reproduces the matrix to 3.41·10⁻¹⁰.545545545455554545545545545454555555454545545545545455554545545545rows and columns, in the order the points arrivedark: kept dense · light: two thin factors, rank printedthe strong partitionblocks112kept dense46largest rank5numbers stored2.7·10⁴‖A − A_H‖ ⁄ ‖A‖3.4·10⁻¹⁰the picture is decidedbefore a number is read
Fig. 1 A 256-square cut into 112 blocks by a rule that reads four numbers per pair. The filled blocks are kept dense; the rest carry the rank they cost.

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.

A 256 × 256 kernel matrix partitioned by the weak rule, with each compressed block's rank46 blocks: 16 kept dense and 30 stored as two thin factors, whose ranks run from 8 to 12. The rule decides one thing — whether a pair of index clusters may be compressed — and it decides it from four numbers about where those clusters sit, before any entry of the matrix is read. The weak rule compresses every off-diagonal block there is, including the two touching halves at the top level, whose rank is the largest number on the picture and is the one that climbs as the matrix grows. The whole thing stores 24,064 numbers against 65,536 entries, and reproduces the matrix to 1.08·10⁻⁹.121210109988889988881010998888998888rows and columns, in the order the points arrivedark: kept dense · light: two thin factors, rank printedthe weak partitionblocks46kept dense16largest rank12numbers stored2.4·10⁴‖A − A_H‖ ⁄ ‖A‖1.1·10⁻⁹the picture is decidedbefore a number is read
Fig. 2 And what the same matrix looks like under a rule that asks nothing: 46 blocks, every off-diagonal one compressed, including the two touching halves at the top level.

What comes out

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

blocks kept dense compressed largest rank numbers stored
the test applied 112 46 66 5 27,008
the test omitted 46 16 30 12 24,064

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.

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 the ones between the first and last quarters of the index set, and each of them is 128 × 64 stored in five columns.

And every compressed block carries the same small number. There are 66 of them and the largest rank anywhere is 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 rank of an admissible block and of a touching one, against how finely they are sampledBoth blocks are of the kernel 1/r; both are 32, 64, 128, 256 points a side; both are truncated at 10⁻⁸. The admissible pair — [0, 1] against [2, 3] — needs 5, 5, 5, 5 columns, which is one number. The touching pair — [0, 1] against [1, 2] — needs 9, 11, 12, 13, climbing by about one per doubling, which is a logarithm. Neither of them grows like the block, and only one of them stops. That difference is what the admissibility test in a partition is buying, and it is why the touching pair is kept dense rather than compressed at all.45678903691215log₂ of the points a sidecolumns above 10⁻⁸two intervals that touchtwo that do notthe measurement that is a null resultadmissible, n = 325admissible, n = 2565touching, n = 329touching, n = 25613stored ⁄ dense at largest0.039the rank belongs to the geometryand not to the sampling
Fig. 3 Where the uniform rank comes from, from the essay before this one: an admissible block’s cost is one number across a factor of eight in sampling, and a touching pair’s is not.

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.

Multiplications in a hierarchical solve and in a dense factorisation, and where they crossBoth counted rather than estimated: the recursion carries a counter and reports what it actually did. The dense factorisation is n³/3 and rises at exactly three per doubling. The hierarchical solve rises at 2.13, 1.93, 1.74 — falling, because the cost is n log²n and its exponent is on its way to one. The two cross between n = 64 and n = 128: below it the format is the more expensive way to solve the system, at 1.48 times the dense count, and at n = 512 it is 6.2 times cheaper. Every point returns an answer at a backward error of about 1.4·10⁻¹⁰, so the comparison is between two ways of getting the same thing.567891010⁴10⁵10⁶10⁷10⁸10⁹log₂ nmultiplicationsdense factorisation, n³⁄3the recursion, countedwhere the format starts payingratio at n = 641.5ratio at n = 5120.16exponent, first doubling2.1exponent, last doubling1.7backward error1.4·10⁻¹⁰cheaper is a sizenot a property
Fig. 4 What the weakly admissible partition is for, from three essays ahead: a recursion that solves without assembling anything, and where it starts costing less than a dense factorisation.
Nonzeros in the Cholesky factor of the 12×12 grid Laplacian, by orderingA horizontal bar chart comparing the number of nonzero entries in the Cholesky factor under four elimination orderings, with reference marks for the matrix itself and for a dense factor.natural1739reverse Cuthill–McKee1354minimum degree1026nested dissection1413matrix: 408 entries · dense factor: 10440bandwidth 12 · 4.26× the matrixbandwidth 12 · 3.32× the matrixbandwidth 123 · 2.51× the matrixbandwidth 108 · 3.46× the matrixn = 144, five-point stencilevery ordering fills in; none avoids it
Fig. 5 And the sparsity field’s version of the same trade, from the essay that measured it: two orderings of one matrix, chosen for two different resources.

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 the reason the trade 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.

The largest rank in each partition, and what refusing to compress a touching pair costsThe strong rule's worst block is rank 5 at every size — one number across a factor of eight — because it never compresses a pair of clusters that touch. The weak rule compresses them and its worst rank climbs 9, 10, 12, 13, at about one per doubling, which is the touching block's logarithm arriving inside a whole partition. That is what the test buys. What it costs is on the badge: at every size measured, the finer partition stores MORE — 67,968 numbers against 61,440 at n = 512 — because it pays in blocks what it saves in rank, and the blocks near the diagonal are dense. The bounded rank is an asymptotic argument and the sizes here are not asymptotic.567891003691215log₂ nlargest rank in the partitionweak: touching pairs compressedstrong: touching pairs refusedthe test bounds a rank and costs storagestrong, blocks250weak, blocks94strong, numbers6.8·10⁴weak, numbers6.1·10⁴strong ⁄ weak1.1the better partitionis the more expensive one
Fig. 6 The trade, measured across four sizes: one rule whose worst rank is a constant, one whose worst rank climbs, and the storage that goes with each.
Products with A to return all 2 copies, against the block widthOne bar per block width. Widths below 2 return fewer copies than the eigenvalue has, at every step count up to a basis of half the problem's dimension, and are marked rather than drawn. The cheapest width that works is 2, at 32 products with A; the widest drawn costs 84.an eigenvalue repeated 2 times, in a 200-dimensional problemblock of 1never returns them allblock of 232 products, 16 stepsblock of 351 products, 17 stepsblock of 460 products, 15 stepsblock of 570 products, 14 stepsblock of 684 products, 14 stepsthe narrowest that workscheapest width2its products with A32the widest drawn84narrower than the multiplicity is not slowit is absent
Fig. 7 Why an optimum in a trade like this belongs to the machine, from the cost field: the best block size is a property of a memory hierarchy rather than of a matrix.

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.

What one symmetric permutation does to the storage, on a matrix it does not changeThe same 256 × 256 matrix, its rows and columns renumbered by one permutation and its columns by the same one. The condition number is 24.3948 either way, to eight digits; the Frobenius norm is 6139.964 either way, to twelve. Nothing a norm can see has moved. Clustered, the partition stores 27,008 numbers. Shuffled, the admissibility test finds no admissible pair anywhere — every cluster of a shuffled numbering spans the whole interval, so every q is infinite — and the format degenerates to dense storage exactly. The rule with no test to fail does worse than that: it compresses every off-diagonal block regardless, gets ranks up to 119 out of 128, and stores 118,208 numbers — 1.80 times the matrix it was compressing. A rank-119 factorisation of a 128-column block is a more expensive way to write down the block than the block.numbers stored, 256 × 256clustered, strong27,008clustered, weak24,064the dense matrix65,536shuffled, strong65,536shuffled, weak118,208the same matrix, twiceκ, clustered24κ, shuffled24‖A‖_F, clustered6140‖A‖_F, shuffled6140shuffled weak ⁄ dense1.8the compressibility is in the numberingand the numbering is not in the matrix
Fig. 8 The second limit, measured: the same matrix under one permutation, where the test finds no admissible pair anywhere and the format stores every entry.
Two kernels on identical geometry, as the picture is scaled up at a fixed wavelengthTwo parallel segments of length L a distance 2L apart, so q = 0.50 at every point on this sweep and nothing about the arrangement changes. 1/r needs 6, 6, 6, 6, 6 columns — one number. cos(κr)/r at κ = 40 needs 12, 16, 22, 33, 53, growing like a power of the size with no size at which it stops, because what decides it is the number of wavelengths across the pair — 3.2, 6.4, 12.7, 25.5, 50.9 — and no ratio of lengths can see that. The segments are parallel rather than collinear on purpose: in one dimension cos(κ|x − y|) obeys an addition formula and is exactly rank two, which makes the collinear version of this look like a confirmation of the smooth case and is an artefact of the arrangement.-1012301020304050log₂ of the segment lengthcolumns above 10⁻⁸cos(κr) ⁄ r at κ = 401 ⁄ r, same pointsthe case with no answerq, held fixed0.51/r at L = 0.561/r at L = 86cos(κr)/r at L = 0.512cos(κr)/r at L = 853the geometry did not moveand the rank did
Fig. 9 And the first, from the previous essay: a test that passes on a pair whose block needs nine times the columns.

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.

Numbers stored per unknown, against the size of the matrix, at 10⁻⁶The dense matrix stores n numbers per unknown, which is the straight line through the origin and doubles whenever n does. The hierarchical representations store 51, 72, 93, 114 and 42, 58, 76, 96, adding about 21 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 22 and 19 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.49, 1.38, 1.30, falling towards one and never arriving, which is what a logarithm looks like from inside.56789100108216324432log₂ nnumbers stored per unknowndenseweak: every off-diagonal blockstrong: only the admissible onesa constant per doublingstrong, n = 5125.9·10⁴weak, n = 5124.9·10⁴dense, n = 5122.6·10⁵per doubling21‖A − A_H‖ ⁄ ‖A‖4.1·10⁻⁸the dense line doublesand the other two add a constant
Fig. 10 The derivation, measured, at six digits: a straight line against log n where a dense matrix is a straight line against n.
Stored entries per row at each level, 400 unknowns and 1192 edgesA rising curve of stored entries per row against level in the hierarchy, beside a flat one for a geometric hierarchy on a grid.01234567110¹10²levelstored entries per rowalgebraicgeometrictwo complexities, one hierarchygrid complexity3.1operator complexity18geometric, for comparison1.5a ring plus random chords — nothing is anywherelevel 6 is 100% dense
Fig. 11 And the multigrid field’s version of the same arithmetic, from the essay that counts it: a hierarchy whose levels each cost a fraction of the one below, and a total that is a constant times the finest.

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.

Entries in U, against the bound the symbolic phase can compute, on the 5×5 gridA row of horizontal bars. The topmost is longest and is labelled as the bound; every measured bar below it is shorter.the bound211no pivoting127τ = 0.00198τ = 0.003100τ = 0.01106τ = 0.03111τ = 0.1116τ = 0.3117τ = 1138entries in Ua bound, and its slackthe bound, from the graph alone211the worst that occurs138loose by1.5no arithmetic was done to compute the bound — only the pattern of AᵀA and its elimination graphGeorge and Ng: U fits inside chol(AᵀA), whatever the swapsprovable, cheap, and loose
Fig. 12 The sparsity field’s version, from the essay that measured the gap: a structure decided in advance, and the difference between what it allocates and what the arithmetic uses.
How many columns a decade of accuracy costs, measured and predicted, at q = 0.500The number of singular values above ε, against the number of digits ε asks for, on a 128 × 128 block between two intervals separated by a gap of 1. The measured curve is a straight line at 0.55 columns a decade: a digit costs the same handful of columns wherever you buy it, which is why the accuracy is a knob and not a cliff. The upper line is what the geometry alone promises — q^(p+1)/(1 − q) below ε, solved for p, using four numbers and no entry of the matrix — at 3.32 columns a decade. Both are straight and they are not the same straight: the bound is right about the shape and loose by 5.3× about the constant, which is the safe direction for a quantity you have to allocate storage from.0369121501122334455digits asked for, −log₁₀ εcolumns keptwhat the geometry promiseswhat the matrix costsa rank is a number of digitscolumns a decade0.55bound, a decade3.3rank at 10⁻⁸5bound at 10⁻⁸28q0.5the shape is rightand the constant is not
Fig. 13 And this field’s, from two essays back: the bound and the outcome, both straight, six apart.

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.

A 256 × 256 kernel matrix partitioned by the strong rule, with each compressed block's rank112 blocks: 46 kept dense and 66 stored as two thin factors, whose ranks run from 4 to 5. The rule decides one thing — whether a pair of index clusters may be compressed — and it decides it from four numbers about where those clusters sit, before any entry of the matrix is read. The strong rule refuses any pair whose clusters touch, and subdivides instead, so the diagonal is fringed with small dense blocks and no rank on the picture exceeds 5. The whole thing stores 27,008 numbers against 65,536 entries, and reproduces the matrix to 3.41·10⁻¹⁰.545545545455554545545545545454555555454545545545545455554545545545rows and columns, in the order the points arrivedark: kept dense · light: two thin factors, rank printedthe strong partitionblocks112kept dense46largest rank5numbers stored2.7·10⁴‖A − A_H‖ ⁄ ‖A‖3.4·10⁻¹⁰the picture is decidedbefore a number is read
Fig. 14 The consequence, drawn: the diagonal fringed with the pairs the test refused, each of them small, each of them stored entry by entry.
How many columns a decade of accuracy costs, measured and predicted, at q = 0.800The number of singular values above ε, against the number of digits ε asks for, on a 128 × 128 block between two intervals separated by a gap of 0.25. The measured curve is a straight line at 0.84 columns a decade: a digit costs the same handful of columns wherever you buy it, which is why the accuracy is a knob and not a cliff. The upper line is what the geometry alone promises — q^(p+1)/(1 − q) below ε, solved for p, using four numbers and no entry of the matrix — at 10.34 columns a decade. Both are straight and they are not the same straight: the bound is right about the shape and loose by 12.7× about the constant, which is the safe direction for a quantity you have to allocate storage from.0369121503570105140digits asked for, −log₁₀ εcolumns keptwhat the geometry promiseswhat the matrix costsa rank is a number of digitscolumns a decade0.84bound, a decade10rank at 10⁻⁸7bound at 10⁻⁸90q0.8the shape is rightand the constant is not
Fig. 15 And what happens just outside the boundary, where q is 0.8 and the block is compressible at seven columns against a bound of ninety.
Numbers stored per unknown, against the size of the matrix, at 10⁻⁸The dense matrix stores n numbers per unknown, which is the straight line through the origin and doubles whenever n does. The hierarchical representations store 54, 79, 106, 133 and 50, 70, 94, 120, adding about 26 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 26 and 23 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.55, 1.42, 1.33, falling towards one and never arriving, which is what a logarithm looks like from inside.56789100108216324432log₂ nnumbers stored per unknowndenseweak: every off-diagonal blockstrong: only the admissible onesa constant per doublingstrong, n = 5126.8·10⁴weak, n = 5126.1·10⁴dense, n = 5122.6·10⁵per doubling26‖A − A_H‖ ⁄ ‖A‖3.8·10⁻¹⁰the dense line doublesand the other two add a constant
Fig. 16 What all of it adds up to: numbers stored per unknown under both rules, against a dense matrix that stores n.

The partition, at other settings

A 256 × 256 kernel matrix partitioned by the weak rule, with each compressed block's rank46 blocks: 16 kept dense and 30 stored as two thin factors, whose ranks run from 8 to 12. The rule decides one thing — whether a pair of index clusters may be compressed — and it decides it from four numbers about where those clusters sit, before any entry of the matrix is read. The weak rule compresses every off-diagonal block there is, including the two touching halves at the top level, whose rank is the largest number on the picture and is the one that climbs as the matrix grows. The whole thing stores 24,064 numbers against 65,536 entries, and reproduces the matrix to 1.08·10⁻⁹.121210109988889988881010998888998888rows and columns, in the order the points arrivedark: kept dense · light: two thin factors, rank printedthe weak partitionblocks46kept dense16largest rank12numbers stored2.4·10⁴‖A − A_H‖ ⁄ ‖A‖1.1·10⁻⁹the picture is decidedbefore a number is read
Fig. 17 The other rule’s picture, for reading against the hero: forty-six blocks where the test produces a hundred and twelve, and the two touching halves compressed like everything else.
Numbers stored per unknown, against the size of the matrix, at 10⁻¹²The dense matrix stores n numbers per unknown, which is the straight line through the origin and doubles whenever n does. The hierarchical representations store 60, 94, 131, 170 and 60, 88, 120, 156, adding about 37 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 33 and 30 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.65, 1.48, 1.37, falling towards one and never arriving, which is what a logarithm looks like from inside.56789100108216324432log₂ nnumbers stored per unknowndenseweak: every off-diagonal blockstrong: only the admissible onesa constant per doublingstrong, n = 5128.7·10⁴weak, n = 5128·10⁴dense, n = 5122.6·10⁵per doubling37‖A − A_H‖ ⁄ ‖A‖7.4·10⁻¹⁴the dense line doublesand the other two add a constant
Fig. 18 What both rules cost at twelve digits, where the compressed blocks dominate and the dense fringe does not.
Numbers stored per unknown, against the size of the matrix, at 0.01The dense matrix stores n numbers per unknown, which is the straight line through the origin and doubles whenever n does. The hierarchical representations store 45, 57, 67, 78 and 28, 34, 42, 50, adding about 11 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 15 and 10 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.33, 1.25, 1.21, falling towards one and never arriving, which is what a logarithm looks like from inside.56789100108216324432log₂ nnumbers stored per unknowndenseweak: every off-diagonal blockstrong: only the admissible onesa constant per doublingstrong, n = 5124·10⁴weak, n = 5122.6·10⁴dense, n = 5122.6·10⁵per doubling11‖A − A_H‖ ⁄ ‖A‖6.8·10⁻⁴the dense line doublesand the other two add a constant
Fig. 19 And at two, where almost all of both totals is the diagonal fringe and the test’s cost is at its largest.
How many columns a decade of accuracy costs, measured and predicted, at q = 0.200The number of singular values above ε, against the number of digits ε asks for, on a 128 × 128 block between two intervals separated by a gap of 4. The measured curve is a straight line at 0.39 columns a decade: a digit costs the same handful of columns wherever you buy it, which is why the accuracy is a knob and not a cliff. The upper line is what the geometry alone promises — q^(p+1)/(1 − q) below ε, solved for p, using four numbers and no entry of the matrix — at 1.45 columns a decade. Both are straight and they are not the same straight: the bound is right about the shape and loose by 3.5× about the constant, which is the safe direction for a quantity you have to allocate storage from.0369121505101520digits asked for, −log₁₀ εcolumns keptwhat the geometry promiseswhat the matrix costsa rank is a number of digitscolumns a decade0.39bound, a decade1.4rank at 10⁻⁸4bound at 10⁻⁸12q0.2the shape is rightand the constant is not
Fig. 20 What a generously separated pair costs, which is the easy end of what a partition contains and not where its blocks are.
The accuracy asked for, against the accuracy obtained, for two kernels on one partitionε is applied to each block against that block's own largest singular value; the error is then reported against the whole matrix's norm. The two are not the same number and the dashed diagonal is where they would be. For 1/r the obtained error is 57 times smaller than the tolerance asked for; for log r on the identical partition it is 46 times smaller. Both curves are straight and parallel to the diagonal, so the knob does what a knob should — a decade in buys a decade out — and neither of them sits on it. The factor is how much of the matrix's mass lives off the diagonal, which is a property of the kernel; it is not the number of blocks, which is 16 and would put the curves on the other side of the diagonal.10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³ε asked for, per block‖A − A_H‖ ⁄ ‖A‖ obtainedasked = obtainedlog r1 ⁄ rtwo numbers, not one1/r, obtained at 10⁻⁸1.8·10⁻¹⁰log r, obtained at 10⁻⁸2.2·10⁻¹⁰1/r, obtained ⁄ asked0.018log r, obtained ⁄ asked0.022blocks in the partition16the tolerance is per blockand the error is per matrix
Fig. 21 The accuracy a partition of this shape actually delivers against the accuracy asked of each of its blocks, at the smallest size the sweep runs.
‖A − LLᵀ‖ ⁄ ‖A‖ for a Cholesky computed in the format, against how many truncations it performedThe recursion is the textbook one and every step stays in the format: the off-diagonal factor is exact, because a triangular solve against one factor of a rank-k block leaves a rank-k block, and the Schur complement is a rank-k addition to a hierarchical matrix, which is the only approximate step there is. So the number of truncations is a count — one per admissible block in the subtree being updated, at every level — and it runs from 0 at a leaf of 128 to 98 at a leaf of 8. The residual goes from 3.856·10⁻¹² to 1.035·10⁻¹¹, and all of that movement is the representation, whose own error is the upper curve and moves by the same factor: the ratio between them is 1.00, 0.93, 0.93, 0.91, 0.86, never above one. A hundred approximate operations contributed nothing measurable to the answer.02040608010010⁻¹²10⁻¹¹10⁻¹⁰truncations performed by the factorisationrelative errorthe representation's own error‖A − LLᵀ‖ ⁄ ‖A‖no decomposition without its residualresidual, 0 truncations3.9·10⁻¹²residual, 98 truncations10⁻¹¹representation, deepest1.2·10⁻¹¹residual ⁄ representation0.86levels, deepest5a hundred approximate stepsand an exact-looking factorisation
Fig. 22 What a formatted factorisation of the partition costs in accuracy, which is the measurement that licenses choosing the leaf size from the machine.
Applications of an operator that is never assembled, against the entries of the matrix it stands forThe blocks at one level of the tree have disjoint column supports, so one batch of random vectors samples all of them at once: two batches a level for the ranges, two more for the projections, and one batch of leaf products for every diagonal block in the matrix together. That is leaf + 2(2k + p)·levels products, and the counter inside the operator says 96, 136, 176, 216 at n = 64, 128, 256, 512 — 40 more per doubling, which is a logarithm. The other line is n², the entries the compression route reads. At n = 512 that is 216 products against 262,144 entries, and the representation it produces is within 7.1× of the one that read them all.567891010¹10²10³10⁴10⁵10⁶log₂ nnumbers the construction touchedentries, n²products with the operatoran operator, applied a few hundred timesproducts at n = 512216entries at n = 5122.6·10⁵per doubling40‖A − A_H‖ ⁄ ‖A‖10⁻⁵excess over the compression7.1no entry of the matrixwas ever read
Fig. 23 And how a partition of this shape can be filled in without reading the matrix, from the randomised field: one batch of products a level, because the blocks at a level have disjoint supports.
The 12×12 grid Laplacian and its Cholesky factor, ordered by naturalTwo square sparsity plots side by side. The left shows the nonzeros of the matrix; the right shows the nonzeros of its Cholesky factor, with the entries created by elimination marked in a second colour.the matrix, lower triangle408 entriesits Cholesky factor1739 entries · 1331 created→‖A − LLᵀ‖/‖A‖1.4·10⁻¹⁶fill, symbolic1331fill, numeric1331n = 144 · density 3.2% · bandwidth 12same matrix, renumberedthe answer is identical to rounding
Fig. 24 The sparsity field’s structure decided in advance, for comparison: which entries an elimination will fill in, computed from the graph before any arithmetic.
Columns needed by the block between a separator's two halves, against the length of the separatorThe block is 3×4, 5×6, 7×8, 11×12 across the sweep, so its size more than triples. In the separator's own numbering — the unknowns in the order they sit on the line — it needs 3, 5, 6, 7 columns at 10⁻¹²: about one more per doubling, which is a logarithm. Renumber the same Schur complement by one symmetric permutation and the same block needs 3, 5, 7, 11 — all of them at eight digits and tighter, and a multiple of the geometric count at every tolerance drawn. Nothing about the matrix changed. What changed is that a cluster of a shuffled numbering is scattered along the whole separator, so every block mixes near interactions with far ones and there is no smooth kernel left to compress. The ordering that works here is not chosen: it is the one the separator's geometry already has, which is why nested dissection leaves a compressible fill and a fill-reducing ordering chosen on the graph alone need not.0481216202402468101214unknowns on the separatorcolumns above 10⁻¹²the same matrix, renumberedin the separator's own orderthe ordering the geometry hands overseparator 73separator 237renumbered, largest11the block, largest11share of the square stored0.61the fill is totaland it is not independent
Fig. 25 And what that structure turns out to contain, from the last essay of this run.

Where this goes next

45678903691215log₂ of the points a sidecolumns above 10⁻⁸two intervals that touchtwo that do notthe measurement that is a null resultadmissible, n = 325admissible, n = 2565touching, n = 329touching, n = 25613stored ⁄ dense at largest0.039the rank belongs to the geometryand not to the sampling The size the rank does not notice Sample a kernel block at 32, 64, 128 and 256 points a side and it needs five columns, five, five and five. Sample the touching block next to it at the same four sizes and it needs nine, eleven, twelve and thirteen. Same kernel, same accuracy, one number and a logarithm. natural1739reverse Cuthill–McKee1354minimum degree1026nested dissection1413matrix: 408 entries · dense factor: 10440bandwidth 12 · 4.26× the matrixbandwidth 12 · 3.32× the matrixbandwidth 123 · 2.51× the matrixbandwidth 108 · 3.46× the matrixn = 144, five-point stencilevery ordering fills in; none avoids it The order decides the memory Four elimination orderings on one matrix give factors of 1,739, 1,354, 1,413 and 1,026 entries. All four factorisations are exact, all four return the same answer, and the one with the better asymptotics is not the one that wins. the bound211no pivoting127τ = 0.00198τ = 0.003100τ = 0.01106τ = 0.03111τ = 0.1116τ = 0.3117τ = 1138entries in Ua bound, and its slackthe bound, from the graph alone211the worst that occurs138loose by1.5no arithmetic was done to compute the bound — only the pattern of AᵀA and its elimination graphGeorge and Ng: U fits inside chol(AᵀA), whatever the swapsprovable, cheap, and loose What the symbolic phase can only bound Without pivoting, the fill can be computed from the graph and the count is exact — 233 predicted, 233 measured. With pivoting it is 233 predicted and 242 measured, and what survives is a bound that is right at every threshold and loose by 1.7 times at the largest grid drawn.

Reads more easily once this is understood

Essays that name this one as worth reading first.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

AdmissibilityBlock methodsCluster treeHierarchical matrixKernel matrixMatrix structureOff diagonal rankSubdivisionSymbolic factorisationWeak admissibility