A solve that is d decompositions

The previous essay established what a Kronecker sum is and what it stores. This one is about what can be done with it, and the answer starts badly: almost nothing survives the operation that matters.

What the format is not closed under #

A representation is worth having when the operations a code performs on it stay inside it. The sparsity field’s is not — the whole point of that field is that the factors of a sparse matrix are dense — and this collection has an essay measuring exactly that.

A Kronecker sum is worse. It is closed under addition of two sums on the same factors, and under multiplication by a scalar, and that is the end of the list.

• Its square is not a Kronecker sum. (T ⊗ I + I ⊗ T)² expands to T² ⊗ I + 2T ⊗ T + I ⊗ T², which is a sum of three Kronecker products rather than two.
• Its inverse is not a Kronecker sum, and is not a Kronecker product either.
• Its factors are not. There is no LU of a Kronecker sum whose factors are Kronecker sums, because the elimination couples the two indices at the first step.

The hero figure is the second of those, measured rather than cited. The inverse of T ⊕ T on eight points a side is formed column by column, read as a four-index array, and cut between the index pairs. A Kronecker product is exactly rank one at that cut — the same fact from a different direction, since B ⊗ C is one outer product of the two matrices read as vectors. The inverse’s rank there is eight.

And the sense in which it nearly is #

The singular values at that cut are the interesting half. They fall fast: 1, 0.19, 0.026, 0.0021 and on down, so the inverse is within any accuracy anyone asks for of a short sum of Kronecker products.

Counting the terms an accuracy buys gives 3, 5, 6, 7, 8, 8 at 10⁻², 10⁻⁴, 10⁻⁶, 10⁻⁸, 10⁻¹⁰ and 10⁻¹². That is half a term a decade at eight points a side, over ten decades, with no cliff.

The half is a property of eight points a side, and the slider says so.

Across n = 4, 6, 8, 10, 12 and 14 the slope reads 0.20, 0.30, 0.50, 0.60, 0.70 and 0.80 terms a decade and the rank at the cut reads 4, 6, 8, 9, 10 and 11. So the cost of a digit is not a constant of the format; it rises with the grid, at about 0.06 of a term a decade for each point added to a side.

Reading the two ends of each column separately says where the growth is. At 10⁻² the term counts across the six sizes are 2, 3, 3, 3, 3 and 3 — flat — and at 10⁻¹² they are 4, 6, 8, 9, 10 and 11. Two digits of the inverse of a Kronecker sum cost three terms whatever the grid; twelve digits cost a number that grows with it. The separable approximation is cheap and size-independent where it is coarse, and neither where it is fine.

That changes what the comparison below is worth, and it is worth making the comparison anyway. The hierarchy field measures 0.554 columns a decade for the singular values of a kernel block, on an entirely different object, arrived at from an entirely different argument — a multipole expansion of 1/r rather than an exponential-sum approximation of 1/x over a spectrum. The agreement with the half a term a decade above is an agreement at one grid size: this slope passes 0.554 somewhere between eight and ten points a side and keeps going.

What the two do share is the shape rather than the constant — a straight line of terms against digits, with no cliff, on both objects. A smooth function of a well-separated argument is nearly separable, and nearly costs a fixed number of terms per digit for a fixed problem. Which of the two problems is being held fixed is the part the coincidence of constants obscured.

Why any of that matters #

A representation that is not closed under an operation does not forbid the operation; it forbids performing it in the representation. That distinction is the whole practical content of this section and it is worth being explicit about, because the three failures above sound fatal and only one of them is even inconvenient.

A code that wanted to precondition a Kronecker sum by an incomplete factorisation of it would have to give up the format at the first step, allocate n²ᵈ numbers, and stop. That is the fatal reading, and it is what happens if the format is treated as a storage scheme — a way of writing down a matrix that some other algorithm then consumes.

The alternative reading is that the format is an operator, and that the algorithms which use it are the ones expressible in the operations it does support. That is the same move the matrix-free field makes, and this collection has an essay about what survives it. What is different here is how much survives: a matrix-free operator gives up its entries and keeps only a product, while a Kronecker sum gives up its factorisations and keeps a product, a transpose, a diagonal, and an exact spectrum.

The one thing it is closed under, which is enough #

A Kronecker sum has no useful factorisations and it does have an eigendecomposition, of the only shape that matters here.

If A = QΛQᵀ and B = RMRᵀ are the two factors, then

A ⊗ I + I ⊗ B = (Q ⊗ R)(Λ ⊗ I + I ⊗ M)(Q ⊗ R)ᵀ

and the middle term is diagonal. Its entries are the sums λᵢ + μⱼ, which is the statement the previous essay checked against a closed form; and Q ⊗ R is orthogonal, because a Kronecker product of orthogonal matrices is orthogonal.

So the change of basis that diagonalises an operator with n^d rows is d changes of basis along d indices, each of them an n × n matrix applied along one axis. And the whole solve is:

1. transform b along each index — d mode products, d·nᵈ⁺¹ multiplications;
2. divide by n^d numbers, each a sum of d eigenvalues;
3. transform back — the same d mode products again.

Nothing of size n^d is factorised anywhere. The only decompositions taken are of the d one-dimensional factors, and on the model problem those are the same matrix, so there is exactly one.

What it costs, measured #

The exponent is the argument. A dense factorisation of the assembled operator is (2/3)N³ with N = n^d, so it is n³ᵈ. The route above is 2d·nᵈ⁺¹ — two passes of d mode products, each of which is n·n^d multiplications. Fitted over the sizes the sweep runs, the measured exponent is 4.0 at d = 3, against 3d = 9.

At n = 12 and d = 3 that is 1,728 unknowns, 124,416 multiplications, and a dense factorisation that would be 3.44·10⁹ — a factor of 27,648. At n = 32 and d = 2 the ratio is 5,461. And the gap widens with every refinement of the grid, because the two exponents differ by 2d − 1.

Every solve reproduces its right-hand side to within 8.4·10⁻¹⁵ relative in three dimensions and 1.4·10⁻¹³ in two, which is the badge on the figure and is what the rule this site is named for asks of a routine that decomposes anything.

There is one further saving that this count does not include and a real code would take. The transform along each index is a multiplication by the eigenvector matrix of T, and for the model problem those eigenvectors are discrete sines — so the transform is a discrete sine transform, which costs n log n rather than n². That takes the exponent from d + 1 to d, plus a logarithm, and it takes the number of decompositions from one to zero.

The figures here do not use it, deliberately. A fast transform is a property of this factor and this field is about a structure that holds for any factors at all; measuring the general route and then noting that the model problem admits a faster one keeps the two claims separate, which is the same reason the structure field measures a general Toeplitz solve rather than only a circulant one.

Two routes, and one of them is the point #

The residual above is checked against the operator applied without assembling it, which makes the whole page circular unless something independent has confirmed that the operator is what it claims to be. It has: the previous essay’s two-route check multiplies a vector by the assembled Kronecker product and by the reshaped route and requires the answers to agree, which they do to 4·10⁻¹⁶.

That is the habit this collection runs on, and it matters more here than usual, because both the operator and the solve are written in the same reshaping vocabulary. A transposition error in unfold would produce a consistent pair — an operator and a solver that agree with each other and describe a different matrix — and nothing about the residual would notice. The only thing that catches it is the comparison against the assembled object at a size where assembling is affordable.

Where the method stops #

Three restrictions, and only the first is usually stated.

Every factor must be diagonalisable, and the transform must be affordable. For a symmetric factor that is free — an orthogonal similarity, and for the model problem the transform is a discrete sine transform which costs n log n rather than n². For a non-symmetric factor the eigenvector matrix can be ill-conditioned, and a change of basis by a matrix with a condition number of 10⁸ throws away half the digits before any division happens.

The right-hand side is unrestricted and the operator is not. Any b works; what has to separate is A. A variable coefficient that is not itself a product of one-variable functions destroys the structure, and so does a domain that is not a box.

The condition number is untouched. The previous essay’s arithmetic applies: for the Kronecker sum the condition number is the same as one factor’s, independent of d, which is a favourable accident; for a Kronecker product it is the d-th power. Neither has anything to do with the storage.

The relative in the structure field #

A circulant matrix is diagonalised by the Fourier matrix, which is the same matrix for every circulant of that size, so its solve is a transform, n divisions and a transform back. That sentence and the one three sections above are the same sentence.

The difference is where the fixed basis comes from. A circulant’s is fixed by the structure — every circulant of size n has the same eigenvectors, and nothing about the particular matrix enters. A Kronecker sum’s is fixed by the factors, so it is one decomposition per problem rather than none, and the decomposition is of an object nᵈ⁻¹ times smaller than the operator.

Both are the same move: find the basis in which the operator is diagonal, and pay for the change of basis rather than for an elimination. The structure field’s version is cheaper and applies to fewer matrices.

There is a third relative, in the iterative field, and it is the one a large code actually uses. A Kronecker sum makes an excellent preconditioner for an operator that is nearly one — a variable coefficient close to a product, or a domain close to a box — because a preconditioner is applied and never inverted, and applying this one is the solve above. That is the standard construction behind fast Poisson preconditioning, and it is a use of the format that needs none of the closure properties the first section says it lacks.

The preconditioner, priced #

The paragraph above names the use a large code actually makes of this — a Kronecker sum is an excellent preconditioner for an operator that nearly is one — and then leaves it as a name. It is worth measuring, because the measurement changes the conclusion of the section before it.

A 24 × 24 grid, a five-point operator whose coefficient is 1 + s·sin(3πx)cos(2πy), preconditioned by the constant-coefficient T ⊕ T applied by the solve above. Conjugate gradient iterations to 10⁻¹⁰:

at s = 0 the plain run takes 86 and the preconditioned one takes 1, which is the solve and not a surprise. At s = 0.2, 107 against 11. At s = 0.5, 118 against 19. At s = 0.9 — a coefficient running from 0.1 to 1.9, a factor of nineteen across the domain — 170 against 55.

The first thing to read off is not the ratio. It is that the preconditioned count is the only quantity on the page that notices the variable coefficient. The plain run moves by a factor of two across the whole sweep; the preconditioned one moves by fifty-five. A coefficient that varies by a factor of nineteen is nearly invisible to the operator’s own conditioning and completely visible to a preconditioner built as though it were not there — so the iteration count of this preconditioner is a measurement of distance from separability, and nothing else here measures that.

And then the transform is the whole decision #

The iteration counts are the half everybody quotes. The other half is what an application costs, and on this construction it is large.

One application is two passes of two dense mode products: 55,296 multiplications at n = 24 and d = 2. One five-point matvec is 2,880. So an application is 19.2 matvecs, and the same runs priced in matvecs read

s = 0: 86 plain against 20. s = 0.2: 107 against 222. s = 0.5: 118 against 384. s = 0.9: 170 against 1,111.

It loses at every non-zero setting. The construction that cuts the iteration count by a factor of ten costs three times more than doing nothing, because it pays a full transform pair on every iteration to save nine matvecs.

Put the fast sine transform back and an application falls to 3.7 matvecs. The same runs then cost 5, 51, 89 and 257 against 86, 107, 118 and 170 — winning comfortably up to about s = 0.5 and losing past it.

That reverses the reason the section above gives for leaving the fast transform out. There the argument is that a discrete sine transform is a property of this factor while the structure holds for any factors, so measuring the general route keeps the two claims separate. For the solve that is right, and the reason is that the transform happens once: at n = 12 and d = 3 the whole solve is 124,416 multiplications against a dense factorisation’s 3.44·10⁹, and no transform is fast enough for that ratio to matter.

For the preconditioner it is not right, because the transform is inside the iteration and is paid once per step. So the fast transform is not an optimisation of this construction — it is the condition under which the construction is worth using at all, and a Kronecker-sum preconditioner built on factors with no fast transform is a preconditioner that improves the iteration count and increases the run time. That is the same distinction the cost field spends a run on: a saving counted in one currency and paid in another.

What the field has established so far #

Two essays, and between them one object: a matrix whose index is a tuple, described by d·n² numbers, with a spectrum in closed form, a product that costs d·nᵈ⁺¹ and a solve that costs the same. Nothing in either of them has been an approximation, and nothing has needed a tolerance.

That is the whole of the well-behaved part of this field. Everything after it is about arrays with three or more indices treated as objects in their own right rather than as matrices with paired indices — and there the news is that every theorem this collection has relied on about rank stops holding. The break is not gradual. A Kronecker sum is an ordinary matrix; a three-index array is not one, and the word rank stops naming a single quantity the moment the second index pair is dropped.

The refusal #

The claim under test is a habit of argument rather than a mathematical mistake, which makes it the easier one to commit.

A structural saving is usually stated as a ratio — nᵈ⁻¹ times fewer multiplications — and a ratio is true at every size, including sizes where it is one. The two-route assertion that checks the reshaped product against the assembled one also asserts that the reshaped route costs an order of magnitude less, and it is fed d = 1 and n = 2, where both routes are four multiplications.

It fails, and it is meant to. What that refusal records is that the saving on this page is a statement about an asymptote and about the sizes it was measured at, and not about the identity underneath it: the identity holds at n = 2 and buys nothing there.

The same file’s other refusals are about the algebra rather than the argument. One is fed the claim that a Kronecker sum’s spectrum is the set of products of its factors’ spectra — right shape, right count, wrong at every entry. The other is fed a fold along the wrong mode, which returns an array of exactly the right dimensions whenever two of them agree.