A decomposition made only of SVDs

The previous two essays are about what happens when the definition of a matrix decomposition is carried to three indices. The rank is the least number of rank-one terms; the best approximation of that rank need not exist; the number depends on the field; and computing it is NP-hard.

This one is about carrying the algorithm across instead, and everything comes back.

The algorithm, transcribed #

The singular value decomposition’s method, stated without the word rank, is: take the leading left singular subspace of the matrix and project onto it. That is one instruction, and a tensor has d matrices attached to it rather than one — its unfoldings.

So: for each index k, unfold the tensor with that index on the rows, take the leading rₖ left singular vectors, call them Uₖ. Then form

C = T ×₁ U₁ᵀ ×₂ U₂ᵀ … ×_d U_dᵀ

which is r₁ × r₂ × … × r_d, and the approximation is C put back through the same factors. That is the higher-order SVD, and everything the definition lost is restored by it:

• it exists, because a projection always does;
• it is computed by d matrix decompositions and nothing else, so it inherits their cost and their numerical behaviour exactly;
• it is quasi-optimal: its error is within √d of the best possible for the same ranks, with no iteration — which is as close as this field gets to the best approximation there is, where the matrix case has no factor at all.

The last of those is the one worth stating carefully, and the rest of this page is about what it does and does not say.

What the projection is, and what it is not #

Two things about the construction are worth separating, because conflating them is the source of most of the confusion this format attracts.

It is not a decomposition of the tensor into rank-one terms. The core is not diagonal, so writing the approximation out as a sum of outer products takes r^d of them rather than r. What the format gives is a subspace per index and a small dense object inside them; the rank-one reading is available and is not the point.

Its ranks are ranks of matrices. The triple (r₁, r₂, r₃) is called the multilinear rank, and every one of its entries is the ordinary rank of an ordinary matrix. That is why the set of tensors with multilinear rank at most a given triple is closed — it is an intersection of three sets each defined by vanishing minors — and why the two failures of the previous essays do not occur here. There is a best approximation of a given multilinear rank; it is attained; and the field the entries are read over does not change any of the three numbers.

The price of both properties is the same object: the core. It is dense, it is r^d, and it is where the difficulty of the field goes when the definition is repaired.

The two bounds #

The theorem gives two computable numbers either side of the best possible error, and neither of them is the best possible error.

Write tail_k for the energy in the discarded singular values of the k-th unfolding: the root of the sum of squares of everything past rₖ. Then

maxₖ tail_k ≤ best ≤ ‖T − T̂‖ ≤ √(Σₖ tail_k²)

The lower bound holds because truncating mode k alone already costs tail_k, and no approximation of that multilinear rank can do better in that mode. The upper bound holds because the projections in different modes are orthogonal to one another in the relevant sense, so their errors add in squares.

The two differ by at most √d — that is where the quasi-optimality factor comes from — and both are made entirely of singular values of matrices, so both are computed rather than estimated.

Which of the two it sits on #

The theorem says the error is between the bounds. It does not say where, and the measurement is the unflattering-sounding one: it sits on the upper one.

For the smooth family at n = 12 and d = 3, the ratio of the measured error to its own upper bound runs

0.688, 0.879, 0.986, 0.999, 0.99998

at ranks one, three, six, eight and ten. By rank six the bound is attained to two decimals and by rank ten to five. The ratio to the lower bound therefore climbs to √3 = 1.7320 and sits there.

Read alone, that says the projection is as bad as the theorem permits. Read alone, it is also the natural conclusion, and it is wrong.

The first thing to do with a constant like √d is to move d, and the generator will draw two and four as well as three.

So the ratio to the lower bound reads 1.0000, 1.7320 and 1.9994 at d = 2, 3 and 4, against a permitted √d of 1.4142, 1.7321 and 2.0000. At three and four indices it is not merely inside the factor, it is the factor, to four digits. At two it is 1.0000 and the factor is 1.4142.

That row is the control the rest of the page needed. At d = 2 the construction is the singular value decomposition and Eckart–Young applies: the truncation is optimal, so the measured error must equal the best possible error, and it does — the ratio is one and stays one. Meanwhile the upper bound is 1/0.70711 = √2 above it, because a matrix’s two unfoldings are one matrix and its transpose and their discarded tails are therefore equal, so √(Σ tail_k²) is √2·tail where the truth is tail.

Which settles what the √d is. It is not a defect of the projection that grows with the dimension — at the one dimension where the answer is known independently, the projection is exactly optimal and the ratio to the lower bound is still 1.4142 short of its permitted value. It is the weakness of the lower bound, present already at d = 2 where nothing else is, and at d = 3 and 4 the measured ratio climbing to √d says the projection has caught up with a bound that was always loose rather than that it has fallen behind.

What the second measurement says #

The way to find out whether a truncation is near the best is to try harder and see how much is available. Alternating over the modes — fix every factor but one, take the leading singular subspace of what is left, sweep — is a genuine improvement procedure: each step is the exact minimiser of its own subproblem, so the error can only fall.

Ten sweeps of it, on the same tensors, buy:

A hundredth of a per cent on the two families the format is for. The one place iterating is worth anything is the tensor with independent normal entries — 5.8 per cent — which is the tensor no low-rank format should be applied to.

Those two readings belong together. The families a compression format is for sit at 0.9998 and above of their upper bound and gain a hundredth of a per cent from ten sweeps of alternation; the family it is not for sits at 0.928 and gains 5.8 per cent. The slack against the bound and the gain from iterating are the same quantity seen twice, and on a tensor with structure there is neither.

So the pair of measurements says something neither half says alone. The truncation is essentially optimal; the lower bound is weak; and the √d in the theorem is the gap between two bounds rather than the gap between the projection and the best.

That is a general lesson about quoting a quantity against one of its bounds. A number pinned against its upper bound says nothing about where the optimum is until something independent has looked, and the independent look here is cheap.

Choosing the three ranks #

The construction takes r₁, …, r_d as given and says nothing about where they come from, which in practice is the whole of the decision.

For a matrix the corresponding question has a clean answer, and this collection has an essay about how clean it is not: a rank is a decision about a gap in a spectrum, taken at a tolerance. Here it is that decision d times, on d different spectra, with no requirement that they agree — and the interaction between them is not free, because the error is the root of the sum of the discarded energies rather than the largest of them.

The practical rule that follows from the upper bound is worth stating, because it is the same rule the train format uses and the same rule the hierarchy field arrives at from a different direction. To reach a total relative error of ε, give each mode a budget of ε/√d and truncate each spectrum where its tail falls below that. The total is then at most ε by construction.

What the rule actually delivers #

The natural next sentence is that it will deliver very close to exactly ε, since the upper bound is attained. That sentence is wrong, and the reason is worth having because it is not about the bound.

Applying the rule across two families, two dimensions and six tolerances, and measuring what came back:

family     d    ε        ranks chosen    achieved ⁄ ε
smooth     3   10⁻²        4, 4, 4          0.088
smooth     3   10⁻³        4, 4, 4          0.881
smooth     3   10⁻⁶        7, 7, 7          0.244
smooth     4   10⁻²        3, 3, 3, 3       0.930
smooth     4   10⁻⁶        7, 7, 7, 7       0.0079
hilbert    3   10⁻⁴        5, 5, 5          0.323
hilbert    4   10⁻⁵        6, 6, 6, 6       0.037
hilbert    4   10⁻⁸        7, 7, 7, 7       0.746

Never above ε, in all twenty-four cases, so the guarantee is sound and the rule does what it is for. And never close to it either: the achieved error runs from 0.0079·ε to 0.930·ε, a spread of a hundred and seventeen, with no pattern in the family or the dimension.

What decides where in that range it lands is that a rank is an integer. The rule stops at the first rank whose tail is under budget, and the tail falls in steps rather than continuously. On a spectrum decaying by a decade a singular value, the rank that first gets under the budget gets a decade under it, and the achieved error is a tenth of ε; on a spectrum decaying slowly it lands just under, and the achieved error is 0.9·ε. The two rows at 10⁻² and 10⁻³ for the smooth family choose the same ranks — 4, 4, 4 — because no rank sits between the two budgets, so one of them overshoots by ten and the other by nothing.

And the overshoot is not slack #

That looks like a rule that pays for accuracy nobody asked for, and it is worth checking whether the extra accuracy can be traded back for storage. Take the rule’s ranks and reduce any mode’s rank by one, repeatedly, for as long as the achieved error stays under ε:

In twenty of the twenty-four cases nothing can be reduced at all, and the largest storage saving available anywhere in the sweep is 1.62×. The rule is choosing ranks that are already minimal for the tolerance in almost every case; the accuracy overshoot is the quantisation of an integer rank, not headroom that a cleverer rule could spend.

So the honest form of the sentence has two clauses. The ε/√d rule is a guarantee and not a prediction — plan for ε, expect somewhere between a hundredth of it and all of it — and it is nevertheless close to the best rank choice available at that tolerance, which is why it stays the rule.

That is a backward error chosen in advance, which is the sentence the hierarchy field’s second finding turns on. The accuracy is not measured afterwards; it is typed, the ranks follow from it, and both factors of forward ⪅ κ × backward are known before the computation runs.

What the repair costs #

The format that fixed the definition has a difficulty of its own, and it is the difficulty this whole field is about, moved from the tensor to the core.

The core is r₁ × … × r_d, so at equal ranks it is r^d numbers — the curse this field is named for, arriving in the repair rather than in the thing repaired. At r = 4 and d = 3 that is 64, against a tensor of 512 — an eightfold saving. At d = 5 it is 1,024 against 32,768. At d = 20 it is 10¹² against 10²⁶, which is a compression scheme whose compressed form outgrows the machine before the tensor does.

The measured version of that is three sizes and one arithmetic statement: the core is exactly r^d numbers, so two more indices multiply it by r². The factors are only d·n·r and never matter.

That figure is the separation worth having. The approximation quality of this format does not deteriorate with the dimension: the ratio to the lower bound goes to √d, which is a factor of two at d = 4 rather than a catastrophe, and the tightness against the upper bound is 0.999 at both dimensions. What deteriorates is the size of the object the quality is a statement about, and that is arithmetic rather than analysis.

The arithmetic of the crossing is worth having in one line. A core is r^d + d·n·r and a train is (d − 2)·n·r² + 2·n·r, so at r = 4 and n = 20 the core is 576 against the train’s 800 at four indices and 1,424 against 1,120 at five. The crossing is between four and five, it moves left as the rank rises, and past it the core is never competitive again. That is the whole of the choice between the two formats at fixed rank, and it is a comparison of two counted numbers rather than of two asymptotic classes.

Two routes to the same subspace #

The implementation has a decision in it that is worth writing down, because taking it the obvious way would have put this site on the road it has an essay warning about.

A mode-k unfolding is nₖ × nᵈ⁻¹: a few rows and up to a million columns. The leading left singular vectors of it are wanted, and there are two ways to get them.

The Gram route forms M Mᵀ, which is nₖ × nₖ, and eigendecomposes it. It is cheap — one pass over the tensor — and it squares the condition number, so singular values below about 10⁻⁸ of the largest are lost. That is the normal-equations trap, and this collection measures it on the road that squares the problem.

The transpose route takes the decomposition of Mᵀ, which is tall and thin, and reads the left singular vectors of M off its right factor. It is the same decomposition read the other way round, it is not an approximation, and one-sided Jacobi over its nₖ columns costs a fraction of what running over the unfolding’s columns would.

The second is what the code does. The difference is invisible in every measurement on this page and would not be invisible on a tensor with a spectrum spanning more than eight decades — which is exactly the tensor a compression format is for.

The failure that was found in the drawing #

One defect on this page cost more than the rest of it and is worth recording, because its symptom was that the answer got worse as the method was given more.

The randomised version of this decomposition, two essays on, sketches each unfolding with a random matrix and takes an orthonormal basis for the result. The first implementation took a QR of the sketch and kept its first r columns. Every gate passed. The measurement said the error rose from 2.9 to 6.9 times the deterministic answer as the oversampling went from zero to ten.

The reason is that a QR’s columns are an orthonormal basis for the sketch’s column space in an arbitrary order, so truncating them keeps an arbitrary r-dimensional slice of it — and a wider sketch means a larger space to pick badly from. The repair is one line: take the decomposition of the sketch itself, which is small, and keep the directions that carry the most of it.

It is the same lesson as randomisation does not create structure: the sketch is a smaller thing to decompose and never a way past deciding which directions matter. What makes it worth a paragraph here rather than a comment there is the shape of the symptom. A monotone quantity going the wrong way is the only kind of defect this field’s assertions catch for free, because every other number involved was finite, ordered and plausible.

Where it sits beside the other decompositions here #

The hierarchy field’s format and this one are close relatives and the difference between them is one word.

Both take a large object, find that some part of it is numerically low rank, and store two thin factors instead. Both have a cost that is chosen rather than handed over: a tolerance goes in and a storage comes out. Both are built out of ordinary matrix decompositions and inherit their numerical behaviour.

What differs is what plays the role of “some part”. For a hierarchical matrix it is a block — a submatrix picked out by a geometric admissibility test, and the test reads four numbers per pair of clusters and no entry of the matrix. Here it is a mode — an entire index, with no geometry involved and no test to fail, and the low-rank property is discovered by decomposing rather than predicted by a rule.

That is a real advantage in robustness and a real disadvantage in cost. The hierarchy field’s partition is decided before any entry is read; this format’s ranks require d decompositions of matrices as large as the tensor. The randomised essay two on is the answer to that, and its answer is that the decompositions can be replaced by a sketch whose cost is d small Gaussians.

The refusal #

The claim under test is the flattering reading of the word quasi-optimal, which is that the factor in the bound is a formality and the truncation is optimal in practice.

The assertion that the measured error equals the lower bound is fed a rank-six truncation of the smooth family, where the ratio is 1.708 against a permitted 1.7321. It fails.

The refusal is doing something specific. It is not saying the truncation is bad — the alternating measurement above says it is within a hundredth of a per cent of the best available. It is saying that the lower bound is not the best, so the ratio to it is not the gap to the best, and a claim built on reading the ratio as that gap has to be refused whether the conclusion happens to be right or not.

The file’s other two refusals guard the ends. One is fed a rank-two truncation of a noise tensor and required to refuse the claim that it is exact, which is what keeps the untruncated exactness check from being read as a statement about truncations. The other is fed the core of a smooth tensor and required to refuse the claim that it is diagonal — which is the next essay’s whole subject.