Numerical Linear Algebra
where the algebra stops being the arithmetic
When the index is a tuple

A solve that is d decompositions

A Kronecker sum is closed under nothing useful — its inverse is not a Kronecker sum and no factorisation of it is one. What it has instead is eigenvectors that are Kronecker products, so a solve with 1,728 unknowns takes one decomposition of a 12 × 12 matrix and nothing else.
33 min read 27 figures Where the algebra stops being the arithmeticExact ground truthTwo routes to a number

Worth reading first: An index that is a pair · Orthogonal is a number.

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.

The Kronecker spectrum of the inverse of a Kronecker sum on 8 points a sideA Kronecker product B ⊗ C, read as a four-index array and cut between its two index pairs, is exactly rank one. The inverse of T ⊕ T is rank 8 at the same cut, so the format is not closed under inversion — which is why a solve in it goes through the eigenbasis rather than through an inverse. What the spectrum says is that it is nearly closed: the singular values are 1, 0.19, 0.0262, 0.00214 of the first, and the number of Kronecker terms needed runs 3, 5, 6, 7, 8, 8 at 10⁻², 10⁻⁴, 10⁻⁶, 10⁻⁸, 10⁻¹⁰ and 10⁻¹². That is 0.50 terms a decade — the same shape, and nearly the same number, as the 0.554 columns a decade the hierarchy field measures for a kernel block, arrived at from a different direction entirely.024681012141610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹index of the singular value at the cutσ ⁄ σ₁eight digits10^-4: 5 Kronecker terms10^-8: 7 Kronecker terms10^-12: 8 Kronecker termsnot closed, and nearly closedrank at the cut8a Kronecker product's1terms at 10⁻⁴5terms at 10⁻⁸7terms a decade0.5the inverse leaves the formatby half a term a decade
Fig. 1 The inverse of a Kronecker sum, read as a four-index array and cut between its two index pairs. A single Kronecker product is exactly rank one there. This is not, and its singular values say by how much it misses.

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, over ten decades, with no cliff.

The number is worth stopping on. 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. Both are the same statement in the end: a smooth function of a well-separated argument is nearly separable, and nearly costs a constant per digit.

How many columns a decade of accuracy costs, measured and predicted, at q = 0.500The number of singular values above ε, against the number of digits ε asks for, on a 128 × 128 block between two intervals separated by a gap of 1. The measured curve is a straight line at 0.55 columns a decade: a digit costs the same handful of columns wherever you buy it, which is why the accuracy is a knob and not a cliff. The upper line is what the geometry alone promises — q^(p+1)/(1 − q) below ε, solved for p, using four numbers and no entry of the matrix — at 3.32 columns a decade. Both are straight and they are not the same straight: the bound is right about the shape and loose by 5.3× about the constant, which is the safe direction for a quantity you have to allocate storage from.0369121501122334455digits asked for, −log₁₀ εcolumns keptwhat the geometry promiseswhat the matrix costsa rank is a number of digitscolumns a decade0.55bound, a decade3.3rank at 10⁻⁸5bound at 10⁻⁸28q0.5the shape is rightand the constant is not
Fig. 2 The hierarchy field’s version of the same count, from the essay that measured it: rank against accuracy for a kernel block, a straight line with the same slope this page finds for an inverse.
The Kronecker spectrum of the inverse of a Kronecker sum on 14 points a sideA Kronecker product B ⊗ C, read as a four-index array and cut between its two index pairs, is exactly rank one. The inverse of T ⊕ T is rank 11 at the same cut, so the format is not closed under inversion — which is why a solve in it goes through the eigenbasis rather than through an inverse. What the spectrum says is that it is nearly closed: the singular values are 1, 0.205, 0.045, 0.00694 of the first, and the number of Kronecker terms needed runs 3, 5, 7, 8, 10, 11 at 10⁻², 10⁻⁴, 10⁻⁶, 10⁻⁸, 10⁻¹⁰ and 10⁻¹². That is 0.80 terms a decade — the same shape, and nearly the same number, as the 0.554 columns a decade the hierarchy field measures for a kernel block, arrived at from a different direction entirely.024681012141610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹index of the singular value at the cutσ ⁄ σ₁eight digits10^-4: 5 Kronecker terms10^-8: 8 Kronecker terms10^-12: 11 Kronecker termsnot closed, and nearly closedrank at the cut11a Kronecker product's1terms at 10⁻⁴5terms at 10⁻⁸8terms a decade0.8the inverse leaves the formatby half a term a decade
Fig. 3 And the same measurement on a finer grid. The curve barely moves — the number of terms is set by the spectrum’s dynamic range, which grows like a logarithm, not by the number of unknowns.

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^{2d} 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.

CG iterations on the model problem, with the two preconditioners a matrix-free code can and cannot haveAt 256 unknowns, conjugate gradients takes 55 iterations unpreconditioned. Diagonal preconditioning — n probes, or the analytic diagonal, and therefore available without any entries — takes 55, which is the same number, because this operator's diagonal is constant and scaling by it is scaling by a scalar. An incomplete Cholesky takes 21, a factor of 2.62, and it is defined by the sparsity pattern: there is no way to ask a subroutine for it. That factor is the price of an operator with no entries.36861361862360102030405060unknownsCG iterations to 10⁻¹⁰no preconditionerdiagonal (available)incomplete Choleskyat 256 unknownsunpreconditioned55diagonal55incomplete Cholesky21what the pattern buys2.6a Krylov method needs only productsand a preconditioner needs the entries
Fig. 4 What the matrix-free field pays for giving up entries, from the essay that measures it: a factor in the iteration count, bought with a preconditioner that has to be built from something other than the matrix.

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 λ_i + μ_j, 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^{d+1} 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.

Multiplications in a 3-dimensional model solve, through the eigenbasis against a dense factorisationA Kronecker sum's eigenvectors are the Kronecker products of its factors' eigenvectors, so the change of basis that diagonalises an operator with 1,728 rows is 3 changes of basis along 3 indices. The whole solve is a transform, 1,728 divisions and a transform back: 1.24·10⁵ multiplications at n = 12, against a dense factorisation's 3.44·10⁹, a factor of 2.76·10⁴. The fitted exponent is 4.00 against 3d = 9. The only decompositions taken are of the 3 one-dimensional factors, which on the model problem are the same matrix — so there is exactly one, of size 12. Every solve reproduces its right-hand side to 8.36·10⁻¹⁵.10¹10²10⁴10⁶10⁸10¹⁰n, points along one axismultiplicationsa dense factorisationthrough the eigenbasis3 decompositions of an n × nunknowns1728multiplications1.2·10⁵dense factorisation3.4·10⁹fitted exponent4‖Ax − b‖ ⁄ ‖b‖8.4·10⁻¹⁵nothing of size n^dis ever factorised
Fig. 5 The count, against a dense factorisation of the assembled operator. At twelve points a side the solve is 1.24·10⁵ multiplications and the factorisation would be 3.44·10⁹.
Multiplications in a 3-dimensional model solve, through the eigenbasis against a dense factorisationA Kronecker sum's eigenvectors are the Kronecker products of its factors' eigenvectors, so the change of basis that diagonalises an operator with 1,728 rows is 3 changes of basis along 3 indices. The whole solve is a transform, 1,728 divisions and a transform back: 1.24·10⁵ multiplications at n = 12, against a dense factorisation's 3.44·10⁹, a factor of 2.76·10⁴. The fitted exponent is 4.00 against 3d = 9. The only decompositions taken are of the 3 one-dimensional factors, which on the model problem are the same matrix — so there is exactly one, of size 12. Every solve reproduces its right-hand side to 8.36·10⁻¹⁵.10¹10²10⁴10⁶10⁸10¹⁰n, points along one axismultiplicationsa dense factorisationthrough the eigenbasis3 decompositions of an n × nunknowns1728multiplications1.2·10⁵dense factorisation3.4·10⁹fitted exponent4‖Ax − b‖ ⁄ ‖b‖8.4·10⁻¹⁵nothing of size n^dis ever factorised
Fig. 6 The same picture in three dimensions, where the dense line is drawn only so the other has something to be compared against: at twelve points a side it is a factorisation nobody would attempt.

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^{3d}. The route above is 2d·n^{d+1} — 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.

Multiplications in a 2-dimensional model solve, through the eigenbasis against a dense factorisationA Kronecker sum's eigenvectors are the Kronecker products of its factors' eigenvectors, so the change of basis that diagonalises an operator with 1,024 rows is 2 changes of basis along 2 indices. The whole solve is a transform, 1,024 divisions and a transform back: 1.31·10⁵ multiplications at n = 32, against a dense factorisation's 7.16·10⁸, a factor of 5461. The fitted exponent is 3.00 against 3d = 6. The only decompositions taken are of the 2 one-dimensional factors, which on the model problem are the same matrix — so there is exactly one, of size 32. Every solve reproduces its right-hand side to 1.44·10⁻¹³.10¹10²10⁴10⁶10⁸n, points along one axismultiplicationsa dense factorisationthrough the eigenbasis2 decompositions of an n × nunknowns1024multiplications1.3·10⁵dense factorisation7.2·10⁸fitted exponent3‖Ax − b‖ ⁄ ‖b‖1.4·10⁻¹³nothing of size n^dis ever factorised
Fig. 7 The two-dimensional version over a wider range of grids, where the assembled matrix can be formed at every size drawn and the comparison is against a factorisation that could actually be attempted.

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.

The same preconditioner in one dimension and in two (ρ = 0.9)Iteration count against the number of unknowns, on a logarithmic size axis. In one dimension the preconditioned count is 7, 10, 10, 10 — flat across a factor of six in size. On square grids with the same unknown counts it is 10, 18, 20, 21, climbing, against an unpreconditioned 14, 37, 73, 115. Nothing about the construction changed.10²020406080100120unknownsiterations2D, no preconditioner2D, preconditioned1D, preconditionedone construction, two dimensions2D steps at 16 unknowns102D steps at 100 unknowns211D steps at 100 unknowns10the same averaging, the same transformand a count that no longer stops growing
Fig. 8 The structure field’s version of a transform doing the work in more than one dimension, from the essay that takes it there: the same reshaping argument applied to a doubly structured matrix.

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.

The 196 eigenvalues of the 2-dimensional Laplacian on 14 points a side, computed against their closed formThe matrix has 196 rows and 38,416 entries, and is a sum of 2 Kronecker products of one 14 × 14 matrix — 392 numbers. Its eigenvalues are every sum of 2 numbers drawn from 4sin²(kπ/2(n+1)), so the whole spectrum is written down before anything runs. The line is that closed form and the marks are what a Jacobi decomposition of the assembled matrix returns; the largest disagreement anywhere is 4.01·10⁻¹³. The smallest eigenvalue is 0.08741 and the largest 7.913, so the condition number is 90.52 — which is the part the structure does not help with.033669913216502468eigenvalues in orderλthe closed formmarks: the assembled matrix, decomposeda spectrum nobody computedrows of the matrix196numbers that describe it392λ smallest0.087λ largest7.9worst |computed − exact|4·10⁻¹³the matrix is never neededand neither is its decomposition
Fig. 9 The other check on the same machinery, from the previous essay: every eigenvalue of the assembled matrix against the sum of sines that gives it.
The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 10 And the collection’s oldest version of the habit, from the error field: a solution compared against one that is known rather than against a better approximation of itself.
Three routes to the error of a rank-k approximation of one block, at q = 0.5The upper line is q^(p+1)/(1 − q), which used four numbers about two intervals and no entry of the matrix. The middle line is the truncated expansion actually evaluated — a rank p + 1 matrix written down as a table of powers, with no decomposition anywhere in it — and it falls at -0.36 decades a column, close to the log₁₀ q = -0.30 the bound predicts. The lower line is the decomposition, which is optimal by construction, and it falls at -1.74 — 4.8 times as fast. Two of these curves share no arithmetic. What they agree on is that the error is geometric in the rank; what they disagree on is the base, by a factor in the exponent rather than a constant, and the disagreement is why a partition allocated from the bound is safe and wasteful at the same time.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹rank of the approximationrelative errorthe boundthe expansion, evaluatedthe decompositiontwo routes, one shapebound at rank 50.063expansion at rank 50.0048decomposition at rank 52.8·10⁻⁹written, decades a column0.36best, decades a column1.7one route used the matrixand the other used four numbers
Fig. 11 The hierarchy field’s version of the same discipline, from the essay that introduced it: a rank obtained by a decomposition and by arithmetic on two interval endpoints, made to agree.

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.

How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 12 What the untouched half means, from the error field: an amplification from backward error to forward error that no storage format changes.
Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 13 And the pair of errors the identity separates, which is the reading every cost figure on this site needs beside it.

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^{d−1} 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 condition number a sketched preconditioner leaves, against the condition number it was givenTwo curves against κ(A), both axes logarithmic. The matrix's own condition number climbs the diagonal from 100 to 10¹⁰; κ(AR⁻¹), where R comes from a QR of a 4n-row sketch, is 2.2284 at every one of them — the same number to ten digits, not a similar one. The reason is two lines of algebra: with G = SU the preconditioned singular values are those of (GᵀG)⁻¹, which has no spectrum of A in it at all.10²10⁴10⁶10⁸10¹⁰110²10⁴10⁶10⁸10¹⁰condition number of the matrixcondition number seen by the iterationκ(A), unchangedκ(AR⁻¹)a bound with no κ(A) in itκ(AR⁻¹), every κ(A)2.2κ(SU), the other route2.2κ(A), across the sweep10⁸κ(AR⁻¹), across the sweep1the sketch never sees the spectrumand the spectrum cancels out of the answer
Fig. 14 What a preconditioner is measured by, from the field that owns it: the condition number of the preconditioned operator rather than the cost of building it.
The 32 eigenvalues of a circulant, two waysA circulant matrix of size 32 has its eigenvalues in closed form: they are the discrete Fourier transform of its first column. Plotted against an eigensolver's answer for the same matrix, the two curves lie on top of each other to 4.4·10⁻¹⁵ relative. The eigenvectors are the same for every circulant of this size and are known before any entry is looked at.04812162024283210⁻¹1index kλthe transformthe eigensolverC = F* Λ F is a factorisationworst relative disagreement4.4·10⁻¹⁵‖Cx − b‖/‖b‖ from the transform solve9.6·10⁻¹⁶imaginary part of a real spectrum1.9·10⁻¹⁶n = 32, and the whole matrix is 32 numberseigenvectors known in advance
Fig. 15 The structure field’s version of a basis that is fixed in advance, from the essay that measures it: eigenvalues that are a transform of one row, with no decomposition anywhere.
The spectrum of C⁻¹T against T's own, n = 64, ρ = 0.8Two sorted spectra on a linear vertical axis. The unpreconditioned Toeplitz matrix spreads its eigenvalues across its whole range; the preconditioned operator puts 97% of them within a tenth of 1, with 0 below zero.08162432404856640123456789index, sortedeigenvalueT aloneC⁻¹Teigenvalues, not singular valuesfraction within 0.1 of 10.97eigenvalues below zero0λ_min of the preconditioner0.11the dashed line is 1, where the cluster formsthe solid line is zero
Fig. 16 And what happens to the same idea when the matrix is only nearly circulant, which is the honest boundary of the approach.

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^{d+1} 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 three unfoldings of a 10 × 11 × 12 hilbert tensor, and the singular values of eachA tensor has one matrix per index — put that index on the rows and every other index down the columns — and each of those matrices has an ordinary rank. Here they are 8, 8, 8 at a relative tolerance of 10⁻⁸, from a tensor of 1320 entries whose modes are of different lengths. Nothing requires the three numbers to agree, and nothing requires any of them to be the tensor's own rank: they are three different matrices built from one array. The leading singular values are 3.2, 3.2, 3.2, each normalised to its own mode below.02468101210⁻¹⁶10⁻¹²10⁻⁸10⁻⁴1index of the singular valueσ ⁄ σ₁the tolerance the ranks are read atmode 1 · rank 8mode 2 · rank 8mode 3 · rank 8hilbert: three matrices, one arrayentries1320mode-1 rank8mode-2 rank8mode-3 rank8‖T‖3.3three ranksand none of them is the tensor's
Fig. 17 The first sign of it, from the previous essay: one array, three matrices, three ranks, and no reason for any two of them to agree.

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^{d−1} 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.

At other settings

The Kronecker spectrum of the inverse of a Kronecker sum on 4 points a sideA Kronecker product B ⊗ C, read as a four-index array and cut between its two index pairs, is exactly rank one. The inverse of T ⊕ T is rank 4 at the same cut, so the format is not closed under inversion — which is why a solve in it goes through the eigenbasis rather than through an inverse. What the spectrum says is that it is nearly closed: the singular values are 1, 0.131, 0.00641, 9·10⁻⁵ of the first, and the number of Kronecker terms needed runs 2, 3, 4, 4, 4, 4 at 10⁻², 10⁻⁴, 10⁻⁶, 10⁻⁸, 10⁻¹⁰ and 10⁻¹². That is 0.20 terms a decade — the same shape, and nearly the same number, as the 0.554 columns a decade the hierarchy field measures for a kernel block, arrived at from a different direction entirely.024681012141610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹index of the singular value at the cutσ ⁄ σ₁eight digits10^-4: 3 Kronecker terms10^-8: 4 Kronecker terms10^-12: 4 Kronecker termsnot closed, and nearly closedrank at the cut4a Kronecker product's1terms at 10⁻⁴3terms at 10⁻⁸4terms a decade0.2the inverse leaves the formatby half a term a decade
Fig. 18 The coarsest grid the measurement runs at, where the cut is 16 × 16 and the whole spectrum fits on the axis.
The Kronecker spectrum of the inverse of a Kronecker sum on 10 points a sideA Kronecker product B ⊗ C, read as a four-index array and cut between its two index pairs, is exactly rank one. The inverse of T ⊕ T is rank 9 at the same cut, so the format is not closed under inversion — which is why a solve in it goes through the eigenbasis rather than through an inverse. What the spectrum says is that it is nearly closed: the singular values are 1, 0.199, 0.034, 0.00373 of the first, and the number of Kronecker terms needed runs 3, 5, 6, 7, 8, 9 at 10⁻², 10⁻⁴, 10⁻⁶, 10⁻⁸, 10⁻¹⁰ and 10⁻¹². That is 0.60 terms a decade — the same shape, and nearly the same number, as the 0.554 columns a decade the hierarchy field measures for a kernel block, arrived at from a different direction entirely.024681012141610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹index of the singular value at the cutσ ⁄ σ₁eight digits10^-4: 5 Kronecker terms10^-8: 7 Kronecker terms10^-12: 9 Kronecker termsnot closed, and nearly closedrank at the cut9a Kronecker product's1terms at 10⁻⁴5terms at 10⁻⁸7terms a decade0.6the inverse leaves the formatby half a term a decade
Fig. 19 And the middle of the range, for reading against the two either side of it.
Multiplications in one product with the 3-dimensional Laplacian, assembled against reshapedBoth curves are exact counts rather than estimates. The assembled matrix has n^3 rows, so multiplying by it is n^6 multiplications; applying the 3 one-dimensional factors along their own indices is 3·n^4. The fitted slopes are 6.00 and 4.00 against 6 and 4 exactly. At n = 256 that is 2.81·10¹⁴ against 1.29·10¹⁰, a factor of 2.18·10⁴ — and the ratio is n^2, so it grows with every size rather than settling.10¹10²10²10⁵10⁸10¹¹10¹⁴n, points along one axismultiplicationsassembled: n^6reshaped: 3n^4two exponentsfitted dense slope62d6fitted factored slope4d + 14ratio at n = 2562.2·10⁴nothing is approximatedthe matrix was never anything else
Fig. 20 The product this solve is built out of, counted: two exponents four apart, and no factorisation anywhere.
The 216 eigenvalues of the 3-dimensional Laplacian on 6 points a side, computed against their closed formThe matrix has 216 rows and 46,656 entries, and is a sum of 3 Kronecker products of one 6 × 6 matrix — 108 numbers. Its eigenvalues are every sum of 3 numbers drawn from 4sin²(kπ/2(n+1)), so the whole spectrum is written down before anything runs. The line is that closed form and the marks are what a Jacobi decomposition of the assembled matrix returns; the largest disagreement anywhere is 7.08·10⁻¹³. The smallest eigenvalue is 0.5942 and the largest 11.41, so the condition number is 19.2 — which is the part the structure does not help with.03672108144180216024681012eigenvalues in orderλthe closed formmarks: the assembled matrix, decomposeda spectrum nobody computedrows of the matrix216numbers that describe it108λ smallest0.59λ largest11worst |computed − exact|7.1·10⁻¹³the matrix is never neededand neither is its decomposition
Fig. 21 The diagonal the solve divides by, drawn: 216 numbers, each a sum of three sines, none of them computed.
The 81 eigenvalues of the 4-dimensional Laplacian on 3 points a side, computed against their closed formThe matrix has 81 rows and 6,561 entries, and is a sum of 4 Kronecker products of one 3 × 3 matrix — 36 numbers. Its eigenvalues are every sum of 4 numbers drawn from 4sin²(kπ/2(n+1)), so the whole spectrum is written down before anything runs. The line is that closed form and the marks are what a Jacobi decomposition of the assembled matrix returns; the largest disagreement anywhere is 3.5·10⁻¹³. The smallest eigenvalue is 2.343 and the largest 13.66, so the condition number is 5.828 — which is the part the structure does not help with.01428425670036912eigenvalues in orderλthe closed formmarks: the assembled matrix, decomposeda spectrum nobody computedrows of the matrix81numbers that describe it36λ smallest2.3λ largest14worst |computed − exact|3.5·10⁻¹³the matrix is never neededand neither is its decomposition
Fig. 22 And at four indices, where the grid is three points a side because that is where the assembled matrix stops being affordable — the difficulty this whole field is about, in one slider position.
Backward and forward error of two routes to x, on 24×24 systems across eight decades of κBoth routes start from the same LU factorisation. The LU solve's backward error is 1.7·10⁻¹⁷ at every conditioning — a flat line at the unit roundoff. Multiplying by the explicitly formed inverse gives 2.6·10⁻⁵ at κ = 10¹⁴, a slope of 0.96 against κ. The two forward errors, drawn above them, are 3.4·10⁻⁴ and 0.019 — within a factor of 57, which is why the difference between the two methods is invisible to anyone measuring how wrong the answer is.10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A)relative errorforward, A⁻¹bforward, A\bbackward, A⁻¹bbackward, A\bat κ = 10¹⁴η, LU solve1.7·10⁻¹⁷η, via the inverse2.6·10⁻⁵forward, LU solve3.4·10⁻⁴forward, via the inverse0.019one factorisation, two ways to use itand one of them forfeits the backward error
Fig. 23 What forming an inverse costs when it is done for its own sake, from the essay that measures it — the comparison that says this page’s inverse is a measurement rather than a method.
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. 24 The hierarchy field’s storage curve, for reading against the d·n² this format stores: a logarithm against a constant.
Stored entries per row at each level, 400 unknowns and 1192 edgesA rising curve of stored entries per row against level in the hierarchy, beside a flat one for a geometric hierarchy on a grid.01234567110¹10²levelstored entries per rowalgebraicgeometrictwo complexities, one hierarchygrid complexity3.1operator complexity18geometric, for comparison1.5a ring plus random chords — nothing is anywherelevel 6 is 100% dense
Fig. 25 And the multigrid field’s arithmetic on the same operator, which is the other way of not factorising it.
Convergence factor against grid size, with 2 pre-smoothing sweepsThree curves against grid size on a logarithmic horizontal axis. Two rise towards one; the third is a flat line near a tenth.10¹10¹.³10¹.⁶10¹.⁹⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰¹10².²00.250.50.751grid points nresidual reduction per stepJacobiGauss–SeidelV-cycleV-cycle spread, 8× in size0.003Jacobi at n = 1270.99work exponent, fitted0.079the dashed curve is cos(πh), Jacobi's closed formthe flat line is the whole method
Fig. 26 The rate that method achieves on the matrix this page solves directly, for scale.
Entries against numbers stored, for sin of a sum on 6 points a side, as indices are addedThe upper line is the tensor: 6^d entries, which is a straight line on a logarithmic axis and reaches 46,656 at d = 6. The lower one is the train, which for this family is 4n(d − 1) exactly — 24, 48, 72, 96, 120 — a straight line on a *linear* axis and therefore a logarithm on this one. Its fitted slope against d is 24.0, which is 4n. The two are the same object to within 1.39·10⁻¹⁵, so nothing has been given up: the ratio at d = 6 is 389, and it grows by a factor of n with every index added.123456710¹10²10³10⁴10⁵number of indicesnumbersentries: 6^dstored: 4n(d − 1)exponential against linearentries at d = 64.7·10⁴numbers stored120ratio389slope against d24‖T − T_tt‖ ⁄ ‖T‖1.4·10⁻¹⁵one line is n^dthe other is a constant per index
Fig. 27 And where the field goes next: a format whose storage is linear in the number of indices, for arrays that are not matrices at all.

Where this goes next

03672108144180216024681012eigenvalues in orderλthe closed formmarks: the assembled matrix, decomposeda spectrum nobody computedrows of the matrix216numbers that describe it108λ smallest0.59λ largest11worst |computed − exact|7.1·10⁻¹³the matrix is never neededand neither is its decomposition An index that is a pair A discretisation on a two-dimensional grid of n points a side has n² unknowns and a matrix with n⁴ entries — 10⁸ at n = 100. What that matrix is instead is two Kronecker products of an n × n matrix, which is 2n² numbers, and nothing has been approximated: assembling it was the mistake. 10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A)relative errorforward, A⁻¹bforward, A\bbackward, A⁻¹bbackward, A\bat κ = 10¹⁴η, LU solve2.2·10⁻¹⁷η, via the inverse4.5·10⁻⁵forward, LU solve2.8·10⁻⁴forward, via the inverse0.015one factorisation, two ways to use itand one of them forfeits the backward error The inverse that is never formed x = A⁻¹b is how the solution of a linear system is written and it is not how it is computed. The usual reason given is cost — three times the arithmetic. The real reason is that one of the two routes is backward stable and the other is not, and at κ = 10¹⁴ they differ by twelve orders of magnitude in the number that says whose fault a wrong answer is. 024681012141610⁻¹1index kλthe transformthe eigensolverC = F* Λ F is a factorisationworst relative disagreement2.6·10⁻¹⁵‖Cx − b‖/‖b‖ from the transform solve4.7·10⁻¹⁶imaginary part of a real spectrum1.2·10⁻¹⁶n = 16, and the whole matrix is 16 numberseigenvectors known in advance The matrix that is one row A circulant of size 16 is sixteen numbers, has no zero entry anywhere, and hands over its entire spectrum in closed form — the discrete Fourier transform of its first column, exactly. An eigensolver spends a sweep of Jacobi rotations over 256 entries arriving at the same answer, and agrees to 1.2·10⁻¹⁵.

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.

Condition numberDiscrete laplacianExact ground truthFlop countKronecker productKronecker sumModel problemOff diagonal rankPreconditioningSeparability