An index that is a pair

Three fields of this collection are answers to one sentence: the matrix is too large to store. Sparsity says most of its entries are zero. A matrix-free operator says never form it, and let it be a subroutine. A hierarchical representation says its off-diagonal blocks are numerically low rank, and where they are cut is a decision.

This field is the fourth answer, and it is the only one that changes what an index is.

The size nobody quotes #

A discretisation on a grid of n points along each of d axes has n^d unknowns. That number is quoted constantly and it is not the difficult one. The difficult one is the size of the matrix, which is the square of it: n²ᵈ.

At n = 100 and d = 2 that is 10⁸ entries, which is 800 megabytes and merely annoying. At d = 3 it is 10¹², which is eight terabytes. At d = 5 it is 10²⁰, which is not a large number of entries — it is a number of entries with no physical meaning, exceeding the count of atoms in a gram of anything by enough that no improvement in storage will ever matter.

The sparsity field’s answer applies and is the right one for a finite-difference stencil: the matrix has at most 2d + 1 non-zeros a row, so it is 5n² numbers in two dimensions rather than n⁴. That is the answer this collection has used since the two-dimensional model problem arrived, and it is not what this field is about.

What this field is about is a stronger statement, available for the same matrices and for many that are not sparse at all:

The matrix is a sum of d Kronecker products of n × n matrices.

That is d·n² numbers. Nothing has been approximated, no tolerance has been chosen, and no property of the entries has been exploited. The matrix was never anything else. Writing it out as n²ᵈ numbers was the mistake, and it is a mistake that is made because the notation for the alternative is not taught alongside the notation for the matrix.

What the product is #

For A of size m × n and B of size p × q, the Kronecker product A ⊗ B is the mp × nq matrix built by replacing every entry of A with that entry times the whole of B. It has one row for each pair (i, p) and one column for each pair (j, q), which is the sentence this field turns on: its index is a pair.

The two-dimensional Laplacian on an interior grid is the standard example. Number the unknowns so that uᵢ,ⱼ is at position i·n + j, let T = tridiag(−1, 2, −1) be the one-dimensional operator, and the five-point stencil is exactly

A = T ⊗ I + I ⊗ T

The first term differences along the first index and leaves the second alone; the second does the reverse. That is what a Kronecker sum is, and the whole of the two-dimensional operator is in it.

The two exponents are the whole field in two numbers. Assembling costs n²ᵈ per product because the matrix has that many entries; the factored form costs d·nᵈ⁺¹, because each of its d terms applies an n × n matrix along one axis of an n^d array. 2d against d + 1, measured as fitted slopes rather than quoted: 4.00 and 3.00 at two axes, 6.00 and 4.00 at three, 8.00 and 5.00 at four, 10.00 and 6.00 at five, 12.00 and 7.00 at six.

Across d = 2, 3, 4, 5 and 6 the ratio at n = 256 reads 128, 2.18·10⁴, 4.19·10⁶, 8.59·10⁸ and 1.83·10¹¹ — each about 190 times the last, which is n^(d−1) gaining a factor of n at every axis added. The saving is not a constant factor a bigger machine catches up with; it is a power of the grid whose exponent is the dimension.

For d axes there are d terms and the same construction, with the identity in every position but one. The storage is d·n² whatever d is, which is the property that makes the format worth having: the matrix’s size is exponential in d and its description is linear in it.

The same unknowns, spread over more axes #

The figure holds the problem size roughly fixed at two hundred unknowns and moves d, so what changes along the slider is only how those unknowns are arranged.

At d = 1, 2, 3 and 4 the same construction on 216, 196, 216 and 81 unknowns has condition numbers of 19,080, 90.52, 19.2 and 5.828. Three and a half orders of magnitude, for the same number of unknowns, changed by nothing but how many directions they are laid out along.

The reason is one line and it is worth having: κ is about 4/(π²h²) with h the spacing along an axis, and holding N = n^d fixed makes n = N^(1/d), so h grows with d and κ falls like N^(−2/d). Spreading a fixed budget of unknowns over more axes makes each axis coarser, and it is the coarseness of an axis that a Laplacian’s conditioning is about.

That is worth stating because the field this essay opens is about the curse of dimensionality, and the first thing the measurement says is that the curse is in the storage and not in the conditioning. The spectrum’s lower end climbs from 2.1·10⁻⁴ to 2.343 as the axes multiply — the small eigenvalues, which are the ones every iterative method’s rate is decided by, simply are not there in four dimensions at this budget.

And the two routes agree at every dimension, which is what makes the rest of it a measurement: the worst disagreement between the decomposition of the assembled matrix and the closed form is 2.62·10⁻¹³, 4.01·10⁻¹³, 7.08·10⁻¹³ and 3.5·10⁻¹³ at d = 1, 2, 3 and 4 — no trend, and at the level a decomposition of a two-hundred-row matrix is entitled to.

The spectrum, which is written down #

The one-dimensional Laplacian has its whole spectrum in closed form. Its eigenvalues are

λₖ = 4 sin²(kπ / 2(n + 1)), k = 1 … n

with eigenvectors the discrete sines, and this collection has used them since the iterative field opened, to give it a rate that is known rather than measured.

The Kronecker sum inherits it exactly. If A has eigenvalues α and B has eigenvalues β, then A ⊗ I + I ⊗ B has eigenvalues αᵢ + βⱼ — every sum, one for each pair of indices — and its eigenvectors are the Kronecker products of theirs. So the spectrum of the d-dimensional model problem is every sum of d numbers drawn from that list of sines, and it is known before anything runs.

The hero figure is that statement, checked. The matrix is assembled at a size where a Jacobi decomposition of it can be afforded, its eigenvalues are computed and sorted, and they are compared against the closed form entry by entry. The largest disagreement anywhere is at the rounding level.

That is exact-ground-truth in a place it is rarely available. A collection that measures error has to have something to measure against, and its usual recourse is a better approximation — a solve in higher precision, a decomposition of a smaller problem — which is a comparison against a computation rather than against an answer. Here there is an answer, for an operator nobody can store, at every d.

The product, which never needs the matrix #

The practical half is a product. Given x with n^d entries, the multiplication (A₁ ⊗ A₂ ⊗ … ⊗ A_d)x is computed by reading x as a d-way array and applying each factor along its own index.

In two dimensions that is one line, and the solve built on it is the same reshape used d times. Reshape x into an n × n matrix X and

(A ⊗ B) x = vec(A X Bᵀ)

with vec the same row-major reading that turned the grid into a vector in the first place. Two dense matrix multiplications of n × n matrices — 2n³ multiplications — against the assembled route’s n⁴.

The count is the argument, and it is returned by the code beside the answer rather than described. At n = 10 and d = 3 the assembled route is 10⁶ multiplications and the reshaped one is 3·10⁴, and the two answers agree to 4·10⁻¹⁶ relative. At n = 256 and d = 3 the assembled route is 2.8·10¹⁴ and the reshaped one is 1.3·10¹⁰.

The general form of the saving is worth stating because it is the one place in this collection where a difficulty and its remedy scale together. The assembled product costs n²ᵈ and the reshaped one costs d·nᵈ⁺¹, so the ratio is nᵈ⁻¹/d. Every index added multiplies the saving by another factor of n. Everywhere else on this site a remedy is a constant factor and the difficulty is not.

What it does not buy #

A format that made a difficulty vanish would be worth suspecting, and this one does not.

The condition number is untouched. κ(A ⊗ B) = κ(A)·κ(B) exactly, which is easy to see from the singular values — those of a Kronecker product are all the products of the two sets — and is checked here on a pair of small matrices rather than asserted. Two well-conditioned factors give a product whose condition number is their product, and the d-fold version is a d-th power.

For the model problem the arithmetic is worse than that. The one-dimensional Laplacian’s condition number is about (2(n+1)/π)², so it grows like n². The Kronecker sum’s is the ratio of the largest sum to the smallest, and both scale with d — the largest eigenvalue is d times the largest one-dimensional eigenvalue and the smallest is d times the smallest — so the sum has the same condition number as one of its terms, independent of d.

That is a genuinely favourable accident and it is an accident of the sum. A Kronecker product of d copies of T has a condition number that is the d-th power of T’s, and T’s is already 4,135 at n = 100, so the cube of it is 7·10¹⁰. The distinction between the two constructions is not decoration.

The conditioning goes the other way #

The section above notes that a Kronecker sum of d copies of T has the same condition number as T, because the largest eigenvalue is d·λ_max and the smallest is d·λ_min and the d divides out. That is stated there as an accident, and it is worth following, because the consequence is the opposite of what the phrase curse of dimensionality leads anyone to expect.

Hold the number of unknowns fixed rather than the number of points a side. A problem with N unknowns in d dimensions has n = N^(1/d) points along each axis, so its condition number is κ(T) at that n — and κ(T) grows like n². At N = 10⁶:

in one dimension n is 10⁶ and κ is 4.05·10¹¹. In two, n is 1,000 and κ is 4.06·10⁵. In three, n is 100 and κ is 4,134. In six, n is 10 and κ is 48.4.

Ten orders of magnitude, in the direction nobody expects. The same number of unknowns, the same operator, the same discretisation, and the six-dimensional problem is the well-conditioned one — because resolution is what costs conditioning, and spreading a fixed budget of unknowns over more axes buys less of it along each.

The practical form of that is worth stating plainly, because it decides how a solver behaves. Conjugate gradients converges at a rate governed by √κ, so on the six-dimensional problem it needs a few dozen iterations where the one-dimensional problem of the same size needs hundreds of thousands. The difficulty of high dimensions is entirely in the storage, and the storage is what this field’s representation removes. What is left after the removal is a better-conditioned problem than the one that was easier to store.

Two cautions, both of which the essay’s own arithmetic supplies. This is a property of the Kronecker sum with equal factors, and a Kronecker product of d copies of T has κ(T)^d — 7·10¹⁰ for three copies at n = 100, which is the same n and eight orders worse. And a fixed N in six dimensions is ten points an axis, which is a resolution nobody would accept for a physical problem: the comparison is honest about conditioning and says nothing about whether either discretisation is accurate enough to use.

Which vec, and why the model problem cannot tell #

The identity (A ⊗ B)x = vec(A X Bᵀ) is written above with vec the row-major reading, and it is written the other way round — vec(B X Aᵀ) — in most of the literature, which stacks columns. Both are correct. Neither notation records which convention is in force, and the two differ by a transpose in the middle of the only line of code the whole field turns on.

That is an ordinary hazard and it has an unusual property here: the standard example cannot detect it. The two-dimensional model problem is T ⊗ I + I ⊗ T. Under the row-major reading its action is vec(T X) + vec(X Tᵀ); under the column-major one it is vec(X Tᵀ) + vec(T X). Those are the same two terms in the other order, so a routine written to the wrong convention returns the right answer, to the last bit, on the matrix every account of this field uses to introduce it.

It stays right for as long as the factors are equal and the operator is a sum. It goes wrong the first time somebody uses a different one-dimensional operator along each axis — an anisotropic diffusion, a different mesh spacing, a convection term along one direction only — which is exactly the point at which the format stops being a demonstration and starts being useful. The failure is then silent in the ordinary way: the output has the right shape, the right norm and no NaN in it, and it is the action of the operator with its axes swapped.

So the test that catches it has to be built deliberately, and the site’s own habit says what it looks like: feed the routine two different factors of different sizes, so that a swap changes the shape and not merely the values. That is the second refusal this essay’s file carries, and it is the reason it is written with a 3 × 3 × 2 array rather than a cubic one — a cube is the same kind of blind test as the model problem, one dimension further up.

Where it sits beside the other three #

The four answers to the matrix is too large are not alternatives, and the difference between them is worth naming precisely because they are usually all available at once.

Sparsity is a property of the entries. An entry is zero or it is not; the problem decides, and a code either stores the non-zeros or it does not. Nothing is chosen.

A matrix-free operator is a property of the code. The entries may be anything; what is available is the action, and this collection has an essay about what survives when a matrix is only a subroutine. Again nothing is chosen — the problem hands over a subroutine or it does not.

A hierarchical representation is the first one whose cost is chosen. The off-diagonal blocks are numerically low rank, the place their singular values are cut is a number somebody types, and the storage is a function of that number rather than of the matrix.

A Kronecker representation is none of these. It is not an approximation at any accuracy, so there is nothing to choose; it is not a property of the entries, since a Kronecker product of two dense matrices is dense; and it is not a property of the code, since the entries are perfectly available. It is a property of how the problem was written down — of the fact that the domain is a product of intervals and the operator acts along one axis at a time.

The failure of the format is a failure of the geometry #

The condition under which the representation exists is that the operator separates. A finite-difference Laplacian on a rectangle separates. The same operator on an L-shaped domain does not, because the index set is not a product; a variable coefficient a(x, y) that is not itself a product does not; and an unstructured mesh does not, because there are no axes to speak of.

That is a much sharper restriction than sparsity’s, which is why this field does not replace the others. What it replaces is the class of problems where it applies — tensor-product discretisations, which is most of what a finite-difference or spectral code produces, and separable kernels.

The refusal #

The claim under test is the easiest mistake this field offers, and it is easy because both statements are true of something.

The eigenvalues of a Kronecker product are the products of the eigenvalues. The eigenvalues of a Kronecker sum are the sums. Confusing the two gives a spectrum with the right number of entries, the right sign, and a plausible spread, and it is wrong at every entry.

So the assertion is fed exactly that: the whole spectrum of T ⊗ I + I ⊗ T is computed from the assembled matrix, the set of pairwise products of T’s spectrum is formed, the two sorted lists are compared, and the claim that they agree is required to fail. It does, at the first entry — the smallest sum is 2λ₁ and the smallest product is λ₁², and for the five-point Laplacian at n = 5 those are 0.536 and 0.0719.

The second refusal in the same file is narrower and is about the code rather than the algebra. A mode-k unfolding is a specific matrix, and a routine that folded a mode-1 unfolding back along mode 2 would return an array of exactly the right shape whenever two of the dimensions happen to agree. That is fed a 3 × 3 × 2 array and required to fail.

What the rest of the field is #

Everything above is about a matrix. The word tensor has not appeared, because none of it needs one: a Kronecker sum is an ordinary matrix with an index that happens to be a tuple, and its arithmetic is ordinary matrix arithmetic in a different order.

The rest of this field takes the same step one level further, to arrays with three or more indices treated as objects in their own right, and the news there is worse. Every theorem this collection has relied on about rank stops being true: the best low-rank approximation need not exist, the rank depends on the field the entries are read over, and there is no decomposition that is orthogonal and diagonal at once. What survives is not the definition of a decomposition but the algorithm — take the singular value decompositions of the reshapes — and existence, computability and a quasi-optimality factor all come back with it.