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

The same matrix, numbered twice

One symmetric permutation. The condition number is 24.3948 either way to eight digits and the Frobenius norm is 6.13996414·10³ either way to twelve. The partition that stored 27,008 numbers now finds no admissible pair anywhere and stores all 65,536, and the format that compresses regardless stores 118,208.
33 min read 25 figures Where the algebra stops being the arithmeticMeasured, not assertedAssertions that reject

Worth reading first: Which pairs are allowed to be small · The order decides the memory · The factor is not sparse.

Everything in this field so far has been a statement about a matrix. This essay is the one that says none of it was.

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. 1 Five ways of storing one 256-square. Two of them are the same matrix in a different order, and one of those two is more expensive than storing every entry.

The experiment

Take the 256 × 256 kernel matrix this field has been using and a permutation P of its indices. Form PAPᵀ — the same rows and columns, renumbered together — and compress it exactly as before.

A symmetric permutation is about as harmless a transformation as linear algebra has. It has the same eigenvalues, the same singular values, the same determinant, the same rank, the same trace, the same norm in every norm the site uses, and the same condition number. Two of those are asserted in the code rather than argued:

clustered shuffled
κ₂ 24.3948 24.3948
‖A‖_F 6.13996414·10³ 6.13996414·10³

Eight digits and twelve digits respectively. Nothing a norm can see has moved.

What happens to the representation

numbers stored share of n² largest rank
clustered, the test applied 27,008 41% 5
clustered, no test 24,064 37% 12
the dense matrix 65,536 100% —
shuffled, the test applied 65,536 100% —
shuffled, no test 118,208 180% 119

The fourth row is the test working. Every node of the cluster tree over a shuffled numbering contains points scattered over the whole interval, so every node’s half-width is the whole domain’s and every pair of centres is nearly coincident. q is infinite for every pair, nothing is admissible, and the partition consists entirely of dense blocks. The representation is the matrix, entry by entry, and it took a decomposition of nothing to find that out.

The fifth row is the format with no test to fail. It compresses every off-diagonal block regardless, gets ranks of 119 out of a possible 128, and stores 118,208 numbers to represent 65,536 — because a rank-119 factorisation of a 128-column block costs 119 × 256 numbers where the block costs 128 × 128. The error of that representation is 9.1·10⁻¹⁰, so it is a perfectly accurate way of writing down the matrix at 1.8 times the price.

That fifth row is the useful one. A format that degenerates to dense storage is disappointing; a format that degenerates to worse than dense storage is a trap, because the code reports success, the accuracy is fine, and the memory is gone.

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 The blocks the fifth row compresses, in the numbering where it works. Renumber, and the same block structure is applied to a matrix in which each of those blocks mixes near interactions with far ones.

Why a permutation can do this

Nothing about A changed. What changed is the relationship between a contiguous run of indices and a cluster of points, and every construction in this field is built on that relationship.

The cluster tree splits the index set in half and reads the geometry of the resulting halves. If the indices are ordered along the line, a half of the index set is a half of the line and its half-width is half the domain’s. If the indices are shuffled, a half of the index set is 128 points scattered from one end to the other, and its half-width is the whole domain’s. The tree is still a tree; its nodes have stopped being clusters.

And the blocks follow. The (1, 2) block of the shuffled matrix is not the interaction between the left half and the right half; it is the interaction between one arbitrary half of the points and the other, which contains every near-field pair that used to sit on the diagonal along with every far-field pair that used to be compressible. A block containing a near-singular entry is not a block with a smooth kernel on it, and its singular values do not fall.

So the object with the low-rank structure was never the matrix. It was the matrix together with an ordering of its indices, and the ordering came from the geometry of the problem rather than from anything the matrix carries.

The singular values of one off-diagonal block, for four kernels on the same 128 pointsTwo intervals that do not touch — [0, 1] and [2, 3] — and the 128 × 128 block between them, for four kernels. 1/r falls a factor of 52 a column and is below 10⁻⁸ after 5; log r behaves the same way, for the reason the expansion makes obvious. An independent draw per entry gives 128 singular values above 10⁻⁸ out of 128, on the same size and the same density, which is what makes the other curves a measurement rather than a property of sorted numbers. The block has full algebraic rank in every case; what differs is where the numbers stop mattering, and that is a decision rather than a fact about the matrix.051015202530354010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, in orderσ ⁄ σ₁eight digitsan independent draw per entry1 ⁄ ra cliff, and a control1/r rank at 10⁻⁸5log r rank at 10⁻⁸5noise rank at 10⁻⁸128σ₂ ⁄ σ₁0.024σ₆ ⁄ σ₁2.8·10⁻⁹the block has full rankand five useful columns
Fig. 3 What an off-diagonal block looks like when it is one, from the first essay in this field. The shuffled matrix’s off-diagonal blocks look like the flat curve on that figure.

Where the ordering comes from

In practice nobody shuffles. The reason this matters is not that a code might permute its matrix by accident; it is that a code has to choose the numbering, and the choice is made somewhere that has nothing to do with linear algebra.

The numbering that works here is the one that comes out of a spatial sort — a k-d tree, a Morton order, a recursive bisection of the geometry. It is computed from the points before the matrix exists, it costs O(n log n), and it is the whole of the preprocessing a hierarchical solver does.

That is a strong claim about where the difficulty of a problem lives and this collection has made it once before, in a different field and about a different resource. the-order-decides-the-memory measures two elimination orders on one sparse matrix and finds factors of a thousand entries and ten thousand — the same matrix, the same arithmetic, the same accuracy, one of them fitting in memory. The sentence there is that nothing numerical chooses between them.

The sentence here is the same one about a different resource, and it is worth putting the two side by side because the resemblance is not superficial. In both cases:

  • the quantity being optimised is a storage cost, not an accuracy;
  • the decision is made from the graph or the geometry, before any arithmetic;
  • and no norm, spectrum or condition number of the matrix contains any information about it at all.
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. 4 The sparsity field’s version, from the essay that measured it: two orderings of one matrix, and the fill each of them leaves.
The matrix and its LU factors at τ = 0.1Two sparsity patterns. The left is the matrix; the right is L and U together, with the entries elimination created drawn in a second colour.A156 entries, 36 unknownsL + U372 entries, 23 interchanges→‖PA − LU‖/‖A‖2.9·10⁻¹⁶growth factor38fill created216forward error1.5·10⁻¹⁵the 6×6 grid with a drift and one ruinous cornerred marks are entries that were zero
Fig. 5 And what the difference looks like as a picture rather than a number.

The rank-119 block, which is worth staring at

The fifth row of the table has a number in it that does not look like a numerical result, and it is the clearest single thing in this essay.

The block being compressed is 128 × 128. Its rank, at a tolerance of 10⁻⁸, is 119. Not 128 — the matrix is a kernel matrix and a shuffle does not destroy every trace of smoothness, so nine of its directions are still nearly redundant. But 119 out of 128 is, for every practical purpose, full.

Storing it as two factors costs 119 × (128 + 128) = 30,464 numbers. Storing it as a block costs 128 × 128 = 16,384. The factorisation is 1.86 times the size of the thing it factorises, and because that happens in both off-diagonal blocks at every level, the whole representation comes to 1.8 times the dense matrix.

The reason is arithmetic and it applies to every low-rank format there has ever been: two factors of an m × n block at rank k cost k(m + n), which is smaller than mn only when k < mn/(m + n) — for a square block, only when k is under half the side. A low-rank representation of a block whose rank is more than half its width is a way of making the block bigger. Nothing about it is wrong; it is simply the wrong container.

A code that checks for this is a code that has one more line in it than a code that does not, and the line is if (2k >= n) store it densely. Real implementations have that line. What they cannot check is whether the ranks they are getting are the ranks the geometry should have produced, and this essay’s fifth row is what that looks like from the inside: everything succeeded, the error is 9.1·10⁻¹⁰, and the memory is gone.

The singular values of a sum of two rank-8 blocks, and the 8 a truncation has to discardA rank-8 block times a vector is a rank-8 block times a vector. A rank-8 block times a rank-8 block is a rank-8 block. A rank-8 block PLUS a rank-8 block is a rank-16 block, exactly, and the 16 bars here are why: the sum of two 8-dimensional spaces is generally 16-dimensional, and none of the 16 singular values is small. Truncating back to 8 costs 54.4 per cent of the block. Below the 16th the values are the unit roundoff, which is the check that the doubling is exact rather than approximate. Every product, every factorisation and every Schur complement inside this format is a chain of these, and there is nothing else to do: without the truncation the ranks double at every level and the format is dense by the bottom.σ ⁄ σ₁ of the sum, 64 × 64σ1, kept1σ2, kept0.914σ3, kept0.899σ4, kept0.842σ5, kept0.827σ6, kept0.773σ7, kept0.749σ8, kept0.681σ9, discarded0.655σ10, discarded0.62σ11, discarded0.589σ12, discarded0.562σ13, discarded0.543σ14, discarded0.499σ15, discarded0.461σ16, discarded0.389σ172.96·10⁻¹⁶σ182.77·10⁻¹⁶the operation that is not closedrank of each term8rank of the sum16truncated back to8cost of the truncation0.54the best there is0.54two planesmake a four-space
Fig. 6 Where a rank of that size comes from, from later in this field: the singular values of a block whose structure has been mixed, and how many of them a truncation is not allowed to discard.
Error of the best rank-k approximation to a 10×10 matrixApproximation error against k on a logarithmic axis, with the measured error and the next singular value drawn as separate curves lying exactly on top of one another.12345678910⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1rank k of the approximation‖A − Aₖ‖measured, 2-normσₖ₊₁, from theorymeasured, Frobeniusthe first two agreeto 4.3·10⁻⁹worst |‖A−Aₖ‖₂ − σₖ₊₁| / σₖ₊₁4.3·10⁻⁹worst Frobenius discrepancy4.3·10⁻⁹κ = 10⁹; 30 random rank-3 matrices, none closerthe error is σₖ₊₁
Fig. 7 And the spectrum that makes a factorisation pointless, from the spectra field: a decay slow enough that no truncation saves anything.

What a norm cannot see, which is a theme here

The site has now measured four quantities that decide a computation and are invisible in every norm of the matrix.

the-units-the-matrix-is-measured-in is the first: a diagonal scaling changes the condition number by orders of magnitude without changing the problem, so κ is a property of how the problem was written down. where-the-drift-lands is the second: two perturbations of identical relative Frobenius norm cost 19 preconditioned iterations and 5, because the norm cannot see which end of the spectrum they landed on. the-order-decides-the-memory is the third. This is the fourth.

The common shape is that a norm is a summary, and a summary discards exactly the structure a structural method exploits. That is not a defect of norms — a summary that kept the structure would not be a summary — and it is a reason to be suspicious of any claim about a method’s applicability stated in terms of one.

The practical form of the suspicion: if a method exploits structure, no scalar computed from the matrix says whether the method applies. The thing to look at is what the structure is about, which here is a set of points and in the sparsity field is a graph.

Two condition numbers of one 12×12 system, as its rows are put into different unitsFour curves against the spread of the row units, in decades. κ_∞ of the scaled matrix rises from 15 to 3.3·10⁶ while the componentwise condition number stays at 13.2 throughout — the same system, the same solution, and one of the two numbers is a fact about the units. Hilbert's two numbers are drawn flat beside them at 3.8·10¹⁶ and 1.1·10¹⁶: a matrix whose sensitivity no scaling repairs.0123456110³10⁶10⁹10¹²10¹⁵10¹⁸spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread15κ_∞ at 6 decades3.3·10⁶cond, either end13Hilbert, equilibrated1.3·10¹⁶the solution is the same at every spreadand one of these curves knows it
Fig. 8 The first of the four, from the error field: a condition number moved by six decades by a choice of units, on a problem that did not change.
Preconditioned iterations against the size of the drift, for a drift on each half of the spectrumThe preconditioner is the exact Cholesky factorisation of the first matrix, so it solves that matrix in one step and every iteration counted here is bought by the drift. Both curves are drifts of the same relative Frobenius norm; the upper one sits on the half of the spectrum with the small eigenvalues and the lower on the half with the large. At a relative drift of 10⁻² they cost 38 and 9 iterations, because the drift is divided by the eigenvalue it lands on: the preconditioned matrix has eigenvalues 1 + w/λ, and λ runs over 10⁴.10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10112233445566relative drift ‖E‖ / ‖A₀‖preconditioned conjugate gradient iterationson the small eigenvalueson the large onessame drift, two placessmall end at 10⁻²38large end at 10⁻²9κ(M⁻¹A), small end36κ(M⁻¹A), large end1.3‖E‖/‖A₀‖, both0.01how much the matrix changedis not what the preconditioner cares about
Fig. 9 And the second, from the sequence field, where two perturbations a norm calls identical cost four times as much work as each other.

The partial cases, which are the realistic ones

A shuffle is an extreme and no code produces one. The interesting question is what a partly wrong ordering costs, and the answer is that the degradation is smooth rather than sudden.

The mechanism says why. A cluster’s half-width is the largest distance between any two of its points, so a cluster that is ninety-nine per cent local and one per cent scattered has the half-width of the scattered part. One misplaced point in a cluster of a hundred is enough to make its bounding box the whole domain, and the pair fails the test.

That is a genuinely awkward property and it is why real implementations use bounding boxes computed from a spatial tree rather than from a numbering somebody supplied. It also explains a failure mode worth naming: a matrix whose points are mostly in a nice order with a few stragglers — the boundary nodes of a mesh, say, appended at the end of the numbering — will produce a partition with a few enormous ranks and no obvious reason for them. The diagnosis is not in the linear algebra.

The reverse is also true and more cheerful. A numbering does not have to be optimal; it has to be local, in the sense that a run of indices is a set of nearby points. Any spatial sort achieves that, and the difference between a good one and a mediocre one is a constant factor in the ranks rather than the difference between working and not.

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. 10 What a working ordering buys, from two essays back: storage per unknown that adds a constant per doubling. None of that line exists without the numbering.
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, 4, 5, 6 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 236renumbered, largest11the block, largest11share of the square stored0.52the fill is totaland it is not independent
Fig. 11 And the same comparison in the sparsity field, from the last essay of this phase: a Schur complement whose block is rank six in the separator’s own order and full rank in a shuffled one.

What this makes of the first six essays

It is worth restating the field’s results with the qualification attached, because every one of them needs it and none of them said so.

The off-diagonal blocks of a kernel matrix are numerically low rank — of a kernel matrix whose indices are ordered so that a contiguous run of them is a cluster of nearby points.

The rank does not grow with the block — the rank of the block between two clusters does not grow, and a block of a shuffled matrix is not the block between two clusters.

Storage is n log n — under a numbering produced by a spatial sort, which costs its own n log n and is not part of any of the counts.

The admissibility test decides from four numbers per pair — four numbers computed from bounding boxes, which are meaningful only if the index sets they bound are spatially coherent.

None of that weakens the results. What it does is move them: they are results about a problem, in which a matrix and a set of points arrive together, rather than results about a matrix. On the problems this format exists for the points always do arrive with the matrix — a boundary element code knows where its panels are, a covariance matrix knows where its observations were taken — so the qualification is satisfied by construction and is invisible until somebody hands the format a matrix with no geometry attached.

That last case is worth a sentence because it is the common one in practice. A matrix that arrives without its points can still be compressed, by finding an ordering: cluster the rows by their similarity, or run a spectral partitioner on the graph of the large entries. That is a real technique and it is a different subject, and what it is looking for is exactly the thing this essay shows a norm cannot see.

Strong couplings and the points the matrix kept, ε = 1A grid of points with lines drawn between the pairs the matrix couples strongly, and the points the coarsening kept drawn larger. At small anisotropy only the vertical lines remain and whole rows are kept.the grid the operator came froman edge is a coupling the matrix calls strongkeptinterpolatedwhat the entries decidedcoupling ratio, x against y1strong couplings across x132strong couplings along y132rows kept or dropped whole012×12 grid, θ = 0.25the coarse grid, from the matrix alone
Fig. 12 Finding structure in a matrix with no geometry attached, from the algebraic multigrid field: a hierarchy built from the entries alone, which is the same manoeuvre for a different purpose.
One entry 64× the rest, at 6-bit block significandsTwo bar charts: how many of 320 entries were rounded to zero, and the median entry's relative error, for the block format in two orderings and for a per-element format.entries rounded to zero, of 320block, as given310block, sorted by size22E4M3, either order0median entry's relative errorblock, as given1block, sorted0.013E4M30.022the same numbers, three waysdeleted, as given310deleted, sorted22deleted, per-element06-bit significands, blocks of 32sorting is free and changes no value
Fig. 13 And an ordering chosen from a graph rather than from a space, in the sparsity field.

The refusal

The claim under test is the one the first six essays of this field would leave standing: that whether a matrix can be compressed is a property of the matrix.

The assertion that the shuffled representation stores under half of n² is fed 65,536 out of 65,536, and it fails. The two conditioning figures are asserted equal first, to eight and twelve digits, so that the failure cannot be read as the shuffled matrix is a different matrix.

It is worth saying which reading this refusal is guarding against, because it is not the obvious one. Nobody believes a permutation changes a matrix. What the essays before this one invite is the belief that a number — a decay rate, a numerical rank, a fraction of n² — is something a matrix has, and that a program could go and measure it. The number in the fourth row of the table is what that belief costs: the same matrix, measured twice, 41 per cent and 100 per cent.

What one symmetric permutation does to the storage, on a matrix it does not changeThe same 128 × 128 matrix, its rows and columns renumbered by one permutation and its columns by the same one. The condition number is 20.8802 either way, to eight digits; the Frobenius norm is 2155.267 either way, to twelve. Nothing a norm can see has moved. Clustered, the partition stores 10,112 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 57 out of 64, and stores 28,096 numbers — 1.71 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, 128 × 128clustered, strong10,112clustered, weak8,960the dense matrix16,384shuffled, strong16,384shuffled, weak28,096the same matrix, twiceκ, clustered21κ, shuffled21‖A‖_F, clustered2155‖A‖_F, shuffled2155shuffled weak ⁄ dense1.7the compressibility is in the numberingand the numbering is not in the matrix
Fig. 14 The same five bars at half the size, where the clustered representation is a smaller saving and the shuffled one is exactly as expensive.
What one symmetric permutation does to the storage, on a matrix it does not changeThe same 64 × 64 matrix, its rows and columns renumbered by one permutation and its columns by the same one. The condition number is 17.4280 either way, to eight digits; the Frobenius norm is 752.7504 either way, to twelve. Nothing a norm can see has moved. Clustered, the partition stores 3,456 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 29 out of 32, and stores 6,656 numbers — 1.63 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, 64 × 64clustered, strong3,456clustered, weak3,200the dense matrix4,096shuffled, strong4,096shuffled, weak6,656the same matrix, twiceκ, clustered17κ, shuffled17‖A‖_F, clustered753‖A‖_F, shuffled753shuffled weak ⁄ dense1.6the compressibility is in the numberingand the numbering is not in the matrix
Fig. 15 And at a quarter, where the format is barely worth having in either order and the fifth bar is still larger than the third.
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. 16 The partition the numbering makes possible, which under a shuffle has no light blocks in it at all.

What the numbering was holding up

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. 17 The line the numbering makes possible, from two essays back. Under a shuffle it is the dense line and nothing else.
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. 18 The picture the test produces on a clustered numbering. On a shuffled one every block in it is dark.
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. 19 The measurement the whole format rests on, which is a measurement about a pair of clusters rather than about a pair of index ranges.
The singular values of one off-diagonal block, for four kernels on the same 64 pointsTwo intervals that do not touch — [0, 1] and [2, 3] — and the 64 × 64 block between them, for four kernels. 1/r falls a factor of 52 a column and is below 10⁻⁸ after 5; log r behaves the same way, for the reason the expansion makes obvious. An independent draw per entry gives 64 singular values above 10⁻⁸ out of 64, on the same size and the same density, which is what makes the other curves a measurement rather than a property of sorted numbers. The block has full algebraic rank in every case; what differs is where the numbers stop mattering, and that is a decision rather than a fact about the matrix.051015202530354010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, in orderσ ⁄ σ₁eight digitsan independent draw per entry1 ⁄ ra cliff, and a control1/r rank at 10⁻⁸5log r rank at 10⁻⁸5noise rank at 10⁻⁸64σ₂ ⁄ σ₁0.024σ₆ ⁄ σ₁2.8·10⁻⁹the block has full rankand five useful columns
Fig. 20 And what an off-diagonal block’s spectrum looks like when it is one. A shuffled matrix’s blocks look like the flat curve.
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. 21 The solve that a compressible numbering buys, and that a shuffled one leaves as a dense factorisation with extra bookkeeping.
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 112, 160, 208, 256 at n = 64, 128, 256, 512 — 48 more per doubling, which is a logarithm. The other line is n², the entries the compression route reads. At n = 512 that is 256 products against 262,144 entries, and the representation it produces is within 7.3× 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 = 512256entries at n = 5122.6·10⁵per doubling48‖A − A_H‖ ⁄ ‖A‖4·10⁻⁷excess over the compression7.3no entry of the matrixwas ever read
Fig. 22 The construction that assumes it too: every batch it draws is supported on a cluster, which under a shuffle is the whole domain.
The singular values of the fill between a separator's two halves, and of the same block renumberedThe block is 11 × 12, every entry of it nonzero, and its singular values fall by a factor of 16.2 at the first step and keep falling: 1, 0.062, 0.0039, 1.3·10⁻⁴, 2·10⁻⁶. That is the same cliff a block of a kernel matrix between two intervals has, and it is the same cliff for the same reason — the Schur complement of a discrete Laplacian is a discrete Green's function, so away from the diagonal it is an integral operator with a smooth kernel. The other curve is the same block after one symmetric permutation of the separator's unknowns: 1, 0.84, 0.64, 0.63, 0.57, which is a spectrum with no cliff in it at all. Both matrices have exactly the same entries.02468101210⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹singular value, in orderσ ⁄ σ₁eight digitsrenumberedin the separator's orderthe same cliff, from the other endσ₂ ⁄ σ₁0.062σ₄ ⁄ σ₁1.3·10⁻⁴σ₆ ⁄ σ₁1.5·10⁻⁸renumbered σ₄ ⁄ σ₁0.63rank at 10⁻⁸6the same entriesin two orders
Fig. 23 The same comparison in the sparsity field, from the last essay of this run: one spectrum with a cliff and the same entries renumbered with none.
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, 4, 5, 6 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 236renumbered, largest11the block, largest11share of the square stored0.52the fill is totaland it is not independent
Fig. 24 And the growth it produces, where the ordering is handed over by the separator’s geometry rather than chosen.
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. 25 The sparsity field’s structure decided before any arithmetic, which is the same kind of object as an admissibility test and inherits the same dependence on a numbering.

Where this goes next

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. 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 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. 0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it The units the matrix is measured in One linear system, written twice. The rows of the second are the rows of the first in different units, the solution is identical to the last bit, and the condition number has moved by eight orders of magnitude. One of those two numbers is a fact about the problem and the other is a fact about the notation.

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.

AdmissibilityCluster treeCondition numberElimination orderFill-reducing orderingHierarchical matrixMatrix normOff diagonal rankOrthogonal invariantPermutation