Five indices are cheaper than two

A Kronecker sum is the operator an earlier essay described and then declined to assemble: d matrices of size n, each applied along one index of an array with n^d entries, standing in for a matrix with n^d rows that nobody stores. The essay that solved with one established the route — transform along every index, divide by n^d numbers, transform back — and priced it at a single dimension, where the fitted exponent comes out d + 1 against a dense factorisation’s 3d.

This is about what happens to that price when the number of indices is the quantity that moves, and the answer runs against every other essay in this field. Everywhere else here, an index is the difficulty. The Tucker core is r^d numbers, which is why a format that does not notice the dimension had to be invented to get past five of them. Rank stops being a single number the moment a third index arrives, and with it the guarantee that a nearest low-rank point exists. The word for all of it is the curse of dimensionality, and it is earned.

For this solve it is a discount, and the discount is exact rather than asymptotic. The cost is 2dn multiplications per unknown. Hold the number of unknowns fixed at N = n^d and that is 2d·N^(1/d), which falls as indices are added, because n falls faster than d rises. The residual falls with it. Both statements are read off the same five sweeps, and the whole of the arithmetic behind them is one line: two passes of d mode products, each of which is n multiplications for every entry of the array.

The count per unknown is 2dn, and there is nothing else in it #

The five counts on that figure divided by the unknowns they solve for are 32, 48, 64, 96 and 128, against grids of 8, 12, 16, 24 and 32 points a side. Every one of them is 4n, which at d = 2 is 2dn. Nothing has been fitted: the transform along one index multiplies each of the n^d entries by an n × n matrix once, so it is n multiplications per entry, and there are d of them in each of two passes.

That count is worth holding still, because it is the whole essay. It does not contain N. It does not contain a constant that has to be measured on a machine. It contains the number of indices and the length of one of them, and the number of unknowns enters only through the fact that there are n^d of them to pay it for.

The comparison the count invites is with the operator itself. Applying a Kronecker sum without assembling it is d·n^(d+1) multiplications, which is the same n per entry per index and only one pass. So the solve is exactly two applications of the operator, plus n^d divisions. At two indices and n = 32 the product is 65,536 and the solve is 131,072. At five indices and n = 6 the product is 233,280 and the solve is 466,560. The divisions are not multiplications and are not in either count, and there are 1,024 of them in the first case and 7,776 in the second.

A direct solve that costs two matrix-vector products is an unusual object. An iterative method on the same operator is a matrix-vector product per step and takes as many steps as its conditioning demands; here the answer arrives after two, in a fixed number of operations known before anything runs. Nothing about that depends on the factors being discrete Laplacians. It depends on their eigenvectors being computable and on the operator being a sum rather than a product, which is the distinction the essay that established the spectrum had to make a refusal about.

The first row of that figure is the one to compare against the first row of the one above it. At three indices and n = 4 the problem has 64 unknowns and costs 1,536 multiplications. At two indices and n = 8 the problem has the same 64 unknowns and costs 2,048. The three-index version of the same count of unknowns is cheaper, and it is cheaper by exactly the ratio the formula predicts: 24 per unknown against 32, which is 6·4 against 4·8.

The exponent is d + 1 at every dimension, and it is measured at five of them #

The exponent against n is the statement the earlier essay made at one dimension and could not make at more, because the figure it drew refused any dimension above three. It now draws six, and the fitted slopes over the five sweeps are 3.00, 4.00, 5.00, 6.00 and 7.00 — d + 1 to two decimals every time, against 3d, which over the same range is 6, 9, 12, 15 and 18.

The fit is not doing much work, and saying why is worth a sentence. These are exact multiplication counts rather than timings, so the curve through them is a monomial and the fit recovers its exponent. That is a different kind of claim from a measured running time, which would carry a constant that depends on the machine and a slope that bends where the working set stops fitting — the distinction the cost field spends a run on, where two eliminations perform identical arithmetic at different prices.

The gap between the two exponents is 2d − 1, so it widens with every index added as well as with every point added to an axis. At four indices the largest grid drawn is seven points a side, which is 2,401 unknowns and a dense factorisation of 9.23·10⁹ — a number that already forbids the comparison being made by running both routes. It is computed rather than run, at (2/3)N³, which is what a dense factorisation of a matrix with N rows costs and not what any code would do to a matrix this sparse.

Two things about those sweeps are properties of the drawing rather than of the method. Each spans about a decade of unknowns, and the grid shrinks as indices are added so that it does — six points a side is 7,776 unknowns at five indices and 46,656 at six. And the dense comparison is a formula rather than a run at every size above the smallest, for the reason above.

Four thousand and ninety-six unknowns, indexed five ways #

The exponent is a statement about refining a grid at a fixed number of indices. The claim in the title is the other one, and it needs the number of unknowns held still while the indexing changes.

Four thousand and ninety-six is 2¹², so it factors as n^d in five ways with n and d both integers beyond the trivial one, and every one of those is a model problem this collection can solve:

A dense factorisation of the assembled operator is 4.581·10¹⁰ on every row, because it sees 4,096 unknowns and nothing else. The eigenbasis route sees the indexing and charges a factor of 5.33 less for twelve indices than for two, and the ratio against the dense route accordingly runs from 4.37·10⁴ at the top of the table to 2.33·10⁵ at the bottom.

The middle rows of that table are rows of the figures above. Sixteen points at three indices is the fourth point of the three-index sweep; eight points at four indices is the second point of the four-index one. The first row is off the right-hand end of the two-index sweep by one size, and the last two rows need the six-index figure.

The smallest point of that sweep is worth reading separately. Two points a side at six indices is 64 unknowns for 1,536 multiplications — the same 64 unknowns and, to the multiplication, the same cost as four points a side at three indices, and three quarters of what the 8 × 8 grid on the first figure paid for them. Three ways of indexing the same sixty-four numbers, and the two-dimensional one is the most expensive.

The accuracy moves in the same direction, which is the part that is not arithmetic #

A saving that costs digits is not a saving, and this field is full of them. The measurement here runs the other way, and the residual column of the table is the whole of it: 7.82·10⁻¹³ at two indices against 1.69·10⁻¹⁵ at six, a factor of 462 across a set of problems with identical unknowns and identical right-hand sides.

The reason is not the number of indices. Over all twenty-three rows of the five sweeps, the log of the relative residual against the log of n fits a slope of 2.05 with R² = 0.913, and against the log of the number of unknowns it fits a slope of 0.09 with R² = 0.009 — which is to say that the size of the problem predicts the accuracy of the solve not at all, and the length of one axis predicts it almost exactly, as its square.

The square is the tell. The condition number of the one-dimensional factor grows like n², and the condition number of a Kronecker sum is the condition number of one factor, independent of d — a fact the essay that established the operator derived from the closed-form spectrum, where it reads as a favourable accident. Here it is the lever. Against the log of κ the same twenty-three residuals fit a slope of 1.14 with R² = 0.925, and the κ column of the table falls by a factor of 571 down the same five rows the cost falls down.

The direct check separates the two candidate causes completely. Hold n at eight and raise d from two to six: the unknowns go from 64 to 262,144, a factor of 4,096, and the residual goes from 5.26·10⁻¹⁵ to 6.20·10⁻¹⁵, a factor of 1.18. Adding indices at a fixed axis length costs almost no accuracy at all. Adding unknowns by lengthening the axis costs it as n². So re-indexing a fixed N into more indices improves the residual for exactly one reason: it shortens every axis, and every transform is then a multiplication by a better-conditioned matrix.

This is a rare shape on this site, and it is worth naming as such. The usual arrangement is that the cheap route is the inaccurate one — that is the entire argument of the inverse nobody forms, where the cheaper-sounding route is the one that is not backward stable, and of the product nobody had to form, where forming a convenient product costs half the significant digits. Here the two move together, and they do because they are governed by the same quantity.

The discount has a bottom, and it is at d = ln N #

Per unknown the cost is 2d·N^(1/d). Differentiating in d gives zero at d = ln N, exactly, and the value there is 2e·ln N. At N = 4,096 that is d = 8.32 and 45.22 multiplications per unknown.

The bottom is visible in the table without the calculus. Six indices and twelve give the same 48 per unknown, which is what a minimum between them looks like when only integers are available. The minimising d is not one of them: 4,096^(1/8) is 2.83, so eight indices of equal length do not exist at this size and the two achievable neighbours bracket the optimum at 6% above it.

Two consequences follow, and the second is the more useful.

The discount is bounded. Falling from 256 per unknown to 48 is a factor of 5.33, and no indexing of 4,096 unknowns does better than 5.66 against the two-index case, because that is 256 divided by the continuous bottom. The curse of dimensionality is not repealed here; it is charged to a different account, and the account is N.

And the accuracy bottoms out at the same place. The residual in the last row of the table is 2.00·10⁻¹⁵, which is slightly worse than the 1.69·10⁻¹⁵ two rows above it, and is eighteen times the unit roundoff of the arithmetic it was computed in. Twelve indices of two is a problem whose factors are 2 × 2 matrices with a condition number of 3; there is no conditioning left to improve, so what remains is the accumulation of a longer chain of transforms, and it is visible. Both curves flatten in the same neighbourhood, for two unrelated reasons that happen to arrive together.

What is actually being held fixed #

The table is five different operators, not five ways of solving one. A 64 × 64 grid and a twelve-dimensional lattice with two points on each axis are two problems with 4,096 unknowns each and nothing else in common — different discretisations, different physics, different answers. No code re-indexes its unknowns to take this discount, and nothing here suggests one should.

What the table measures is the price of a dimension, and the honest way to state it is as a comparison between problems a modeller might actually be choosing between. It says that if a model with 4,096 unknowns must be solved, the number of indices those unknowns carry is not what makes it expensive, and having more of them makes it cheaper. That is the opposite of the received statement, and the received statement is about the same operator seen with n held fixed instead.

That figure is the received statement drawn. Both lines steepen as indices are added, because the horizontal axis is n and every point on it is a different number of unknowns — 256 points a side is 6.55·10⁴ unknowns at two indices and 1.10·10¹² at five. Reading it as a picture of the dimension being expensive is correct, and reading it as a picture of the solve being expensive in the dimension is the mistake the table above corrects.

Both readings come off the same formula. The reshaped product is d·n^(d+1), the solve is twice it, and which way the number moves as d rises depends entirely on whether n or n^d is being held still. That is the same trap the contraction-order essay sets in a different variable, where one arithmetic expression has one value and evaluation orders that differ by a factor of two million.

Where the discount stops being a discount #

Four restrictions, and the first is the one that decides whether any of this is reachable.

Every index has to be a genuine separable direction. The operator must be a Kronecker sum in the indexing being priced, and that is a statement about the problem rather than about the array. A variable coefficient that is not a product of one-variable functions destroys it, as does a domain that is not a box. This is the structure an operator known only by what it does can be built to have and cannot be given.

The factors must be diagonalisable and their transforms affordable. For symmetric factors this costs nothing; the eigenvector matrix is orthogonal and the change of basis is exact to rounding. For a non-symmetric factor the eigenvector matrix can be badly conditioned, and the accuracy argument above inverts — a change of basis by a matrix with a condition number of 10⁸ throws away half the digits before the divisions happen.

The counted quantity is multiplications. At twelve indices of two, an operator with 4,096 unknowns is applied as twenty-four passes over the array, each of them a multiplication by a 2 × 2 matrix along a different axis, and eleven of those axes have a stride the memory system will not enjoy. That is precisely the regime where a block size is a property of the machine and a multiplication count stops predicting a running time. The count is exact and the conclusion drawn from it is about arithmetic, not about seconds.

And the model problem is the easy case for the factors. All d of them are the same matrix here, so there is exactly one decomposition rather than d, and the same discrete Laplacian appears in the two-dimensional coarse problem for the same reason: it is the object whose spectrum is known in closed form, which is what makes the residual a check on the solve rather than on the eigensolver.

What follows for anything that indexes its unknowns #

A dimension count is not a difficulty estimate. It is one factor in N = n^d, and the difficulty belongs to N. Two codes with the same number of unknowns and different numbers of indices are not the same size of problem, and the one with more indices is the smaller of the two by every count on this page.

A separable operator should be told it is separable. The whole discount lives in the structure being visible to the solver. Assembling the operator and handing it to a general factorisation spends 4.581·10¹⁰ where 1.966·10⁵ would do, and the assembled matrix is also the one thing here that cannot be stored.

A tolerance is not the price of the saving. Nothing on this page is an approximation. Every count is exact, every residual is at the rounding level of the arithmetic, and no tolerance appears anywhere in the method — which distinguishes it sharply from the decomposition made only of SVDs, where d matrix decompositions also do the work and the error is bounded rather than absent.

And the shape of the saving is not the shape of the storage saving. The operator is described by d·n² numbers rather than n^(2d), which is a saving that grows with d without limit. The solve’s saving is bounded by 5.66 at this size and turns at d = ln N. The two are usually quoted together and only one of them keeps going.

The refusal #

The claim under test is the one this essay is most likely to invite, because the essay spends its length treating an index as a bookkeeping decision.

Every transform in the solve is applied to the same stored array, and the array is a flat run of n^d numbers with a shape attached. It is a short step from there to the idea that the shape is annotation — that unfolding along one index and folding back along another is a relabelling that returns what it was given.

It is not. A 3 × 3 × 2 tensor unfolded along its first mode and folded back along its second returns an array of exactly the right shape and the wrong contents: the difference has norm 44.09 against the tensor’s own 336.44. The assertion that folding along the wrong mode returns the tensor is fed that case and required to reject it, and it does.

The shape is the only thing that says which n entries a given transform multiplies together, and two of the three dimensions agreeing is enough for the error to survive every check that looks at dimensions. Which is the sense in which the indexing on this page is real: it is not a way of writing 4,096 numbers down, it is what decides which of them the arithmetic touches at once.