The same matrix, numbered twice

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

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:

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

What happens to the representation #

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.

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.

How much of the ordering has to be right #

The table compares the two extremes — a perfect spatial ordering against a complete shuffle — and neither is what a code has. Nobody shuffles, and no spatial sort is perfect. So the question that decides how fragile this is: what does a partially coherent ordering cost?

Shuffle within consecutive blocks of size b and leave the coarse order alone. At n = 256, with the format’s leaf at 16:

b = 1, 2, 4, 8 and 16 all store 27,008 numbers — the same 41 per cent, the same maximum rank of 5, byte for byte the same representation as the perfect ordering. At b = 32 it is 31,360 (48%). At 64, 43,776 (67%). At 128, 64,448 (98%). At 256, the full shuffle, 65,536.

Up to the leaf size, destroying the ordering costs exactly nothing. That is not approximately nothing or nothing to three digits — it is the identical representation, because a leaf is stored densely and the order of the indices inside one is not a quantity the format ever reads.

Past the leaf it rises, and it rises smoothly: half at twice the leaf, two thirds at four times, almost all of it at eight. So the failure is neither the cliff the two-row table suggests nor a decay that starts immediately. It is a plateau with a threshold, and the threshold is a parameter of the format — the leaf size — rather than anything about the geometry or the kernel.

One more reading of the plateau, because it says something about the leaf that the field has not had to say before. The leaf size is usually presented as a tuning parameter for arithmetic efficiency — small leaves give more blocks and more overhead, large leaves give more dense storage — and it is chosen from a cost curve. The plateau gives it a second role that is not about cost at all: the leaf is the scale below which the numbering stops being read. Raising it buys tolerance to a sloppier sort and costs dense storage on the diagonal; lowering it does the reverse. That trade is invisible in any measurement made on a perfectly ordered matrix, which is every measurement in this field before this one.

Which is why the preprocessing is robust #

That changes what the previous section’s claim demands of a spatial sort, and it demands much less.

The numbering has to get the coarse scales right and may be as sloppy as it likes below the leaf. Which is exactly what a k-d tree or a Morton order delivers, and not by design: both stop subdividing when a cell holds fewer than some number of points, and leave whatever order those points arrived in. The part of the sort that is arbitrary is precisely the part that is free.

It also says what to check if a hierarchical solver is storing more than it should. The diagnostic is not is the numbering spatially sorted, which is hard to test and mostly true; it is at what scale does the numbering stop being spatial, and the answer only matters if it is above the leaf. A numbering coherent down to groups of sixty-four with a leaf of sixteen is a numbering that will store about two thirds of the matrix and look, to every other check, entirely healthy.

And it locates the trap in the fifth row of the table above one level more precisely. The format with no admissibility test does not merely degrade when the ordering is wrong; it degrades on any block whose points are not clustered, which the plateau says begins at the leaf. So a code with no test and a numbering coherent only down to b = 64 is compressing blocks that are two thirds arbitrary, getting ranks that are two thirds of full, and paying 1.8 times dense storage — while reporting an error of 10⁻¹⁰ and no fault at all.

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 same place a sparse ordering comes from.

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.

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.

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.

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.

The elimination has a partial case of its own, and it is a trade rather than a degradation. Threshold pivoting accepts a diagonal pivot whenever it is at least τ times the largest candidate, so τ decides how much stability is traded for how much fill. Sweeping it on the same matrix:

The fill and the growth move in opposite directions, and by very different amounts. Across τ = 0.02, 0.05, 0.2, 0.3 and 0.5 the entry count runs 347, 359, 378, 381, 405 — a spread of 17% — while the growth factor runs 107, 37, 4.8, 4.2 and 1.0, a spread of a hundred.

Between 0.2 and 0.3 nothing decides differently — the same twenty-one interchanges, three more entries, a slightly lower growth. The threshold is not a continuous dial: it is a set of decision boundaries with flat ground between them, and a tenth of movement either buys nothing or buys everything depending on where it lands.

So seventeen per cent more fill buys two orders of magnitude of growth, and the forward error follows the growth rather than the fill: 7.2·10⁻¹⁵ at τ = 0.02 against 3.4·10⁻¹⁶ at τ = 0.5, a factor of twenty-one. That is a lopsided trade in one direction, which makes the usual advice — relax τ to save fill — a poor bargain on this matrix. It also puts this essay’s own subject in its place: the ordering decides the ranks and the block structure, and τ decides the stability, and the two are separate knobs that a reader might reasonably have expected to be one. A good numbering does not make a weak pivot safe, and a strict threshold does not make a scattered cluster local.

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.

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.