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

A decomposition made only of SVDs

Everything the definition of tensor rank loses comes back if the SVD's algorithm is carried across instead of its definition — take the leading left singular subspace of every unfolding and project onto all of them. It exists, it costs d matrix decompositions, and its error is within √d of the best there is.
35 min read 26 figures Where the algebra stops being the arithmeticMeasured, not assertedTwo routes to a number

Worth reading first: A nearest point that is not there · When the answer is a choice · A bound that holds with probability.

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 truncation error of a smooth tensor against the rank kept, between the two bounds the theorem givesThe middle curve is the measured error of the projection; the upper dashed one is √(Σ_k tail_k²), which the theorem says it cannot exceed, and the lower one is max_k tail_k, which the best possible error cannot fall below. They are a factor of √3 apart. The measurement is that the projection sits on the upper one, and not between them: the ratio of error to bound runs 0.688, 0.717, 0.879, 0.933 … 0.999982, so by rank 10 the bound is attained to five decimals and the ratio to the lower bound is 1.7320 against √3 = 1.7321. That reads as a bad result and is not one — what it says is that the lower bound is weak, which only a second measurement can establish.024681010⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹rank kept in every moderelative errordashes above: √(Σ tail²), the upper bounddashes below: max tail, a floor under the bestsolid: what the projection returnssmooth: pinned to the upper boundrank 10 error1.1·10⁻¹¹its upper bound1.1·10⁻¹¹the lower bound6.3·10⁻¹²error ⁄ bound1error ⁄ lower1.7inside the boundand sitting on it
Fig. 1 A truncation between the two bounds the theorem gives, which are a factor of √d apart. Which of the two it lands on is not in the theorem.

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_k left singular vectors, call them U_k. 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.

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.

The three unfoldings of a 10 × 11 × 12 smooth 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 2.37, 2.37, 2.37, 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 8smooth: three matrices, one arrayentries1320mode-1 rank8mode-2 rank8mode-3 rank8‖T‖2.4three ranksand none of them is the tensor's
Fig. 2 The d matrices the whole construction is made of, and their spectra — three ordinary decompositions of three ordinary matrices.

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 share of a HOSVD core's energy on its superdiagonal, by family, at rank 4A matrix SVD hands over orthonormal factors and a diagonal middle at the same time. For three indices they come apart, and this is the half that does not survive. Every core here is all-orthogonal — the largest inner product between two slices perpendicular to a mode, relative to the core's own energy, is 3.5·10⁻¹⁶ — and none of them is diagonal. The bars are the fraction of the squared norm carried by the 4 entries on the superdiagonal: smooth 97.7%, hilbert 94.5%, wave 26.6%, noise 0.4%. Keeping only those entries costs 0.152, 0.235, 0.857, 0.998 in relative error against the full core's 8.81·10⁻⁴, 4.67·10⁻⁴, 1.33·10⁻¹⁵, 0.925.smooth97.7%hilbert94.5%wave26.6%noise0.4%share of the core's energy on its 4 superdiagonal entrieskeeping only them: 0.152 against 8.81·10⁻⁴keeping only them: 0.235 against 4.67·10⁻⁴keeping only them: 0.857 against 1.33·10⁻¹⁵keeping only them: 0.998 against 0.925orthogonal, and not diagonalsmooth on-diagonal0.98hilbert on-diagonal0.94wave on-diagonal0.27noise on-diagonal0.0039worst slice pair3.5·10⁻¹⁶the slices are orthogonalthe core is not diagonal
Fig. 3 The core, measured: how much of its energy sits on the superdiagonal a rank-one reading would need, and how much does not.

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_k. Then

max_k tail_k ≤ best ≤ ‖T − T̂‖ ≤ √(Σ_k 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.

The truncation error of a hilbert tensor against the rank kept, between the two bounds the theorem givesThe middle curve is the measured error of the projection; the upper dashed one is √(Σ_k tail_k²), which the theorem says it cannot exceed, and the lower one is max_k tail_k, which the best possible error cannot fall below. They are a factor of √3 apart. The measurement is that the projection sits on the upper one, and not between them: the ratio of error to bound runs 0.689, 0.829, 0.906, 0.950 … 0.999812, so by rank 10 the bound is attained to five decimals and the ratio to the lower bound is 1.7317 against √3 = 1.7321. That reads as a bad result and is not one — what it says is that the lower bound is weak, which only a second measurement can establish.024681010⁻¹²10⁻⁹10⁻⁶10⁻³1rank kept in every moderelative errordashes above: √(Σ tail²), the upper bounddashes below: max tail, a floor under the bestsolid: what the projection returnshilbert: pinned to the upper boundrank 10 error1.4·10⁻¹²its upper bound1.4·10⁻¹²the lower bound8.3·10⁻¹³error ⁄ bound1error ⁄ lower1.7inside the boundand sitting on it
Fig. 4 The same three curves on the d-way version of the matrix this site has used for exact ground truth since it began.

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.

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:

family projection after ten sweeps gain
smooth 8.753·10⁻³ 8.752·10⁻³ 1.0001
hilbert 3.9613·10⁻³ 3.9608·10⁻³ 1.0001
noise 0.9412 0.8897 1.058

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.

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.

The truncation error of a noise tensor against the rank kept, between the two bounds the theorem givesThe middle curve is the measured error of the projection; the upper dashed one is √(Σ_k tail_k²), which the theorem says it cannot exceed, and the lower one is max_k tail_k, which the best possible error cannot fall below. They are a factor of √3 apart. The measurement is that the projection sits on the upper one, and not between them: the ratio of error to bound runs 0.618, 0.648, 0.687, 0.725 … 0.927736, so by rank 10 the bound is attained to five decimals and the ratio to the lower bound is 1.5839 against √3 = 1.7321. That reads as a bad result and is not one — what it says is that the lower bound is weak, which only a second measurement can establish.024681010⁻¹rank kept in every moderelative errordashes above: √(Σ tail²), the upper bounddashes below: max tail, a floor under the bestsolid: what the projection returnsnoise: pinned to the upper boundrank 10 error0.5its upper bound0.54the lower bound0.31error ⁄ bound0.93error ⁄ lower1.6inside the boundand sitting on it
Fig. 5 The one family where iterating is worth anything, and where the error is 0.94 — the counterweight this whole field needs.
The truncation error of a wave tensor against the rank kept, between the two bounds the theorem givesThe middle curve is the measured error of the projection; the upper dashed one is √(Σ_k tail_k²), which the theorem says it cannot exceed, and the lower one is max_k tail_k, which the best possible error cannot fall below. They are a factor of √3 apart. The measurement is that the projection sits on the upper one, and not between them: the ratio of error to bound runs 0.780, 1.103, 1.292, 1.530 … 0.779766, so by rank 10 the bound is attained to five decimals and the ratio to the lower bound is 7.1766 against √3 = 1.7321. That reads as a bad result and is not one — what it says is that the lower bound is weak, which only a second measurement can establish.024681010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹rank kept in every moderelative errordashes above: √(Σ tail²), the upper bounddashes below: max tail, a floor under the bestsolid: what the projection returnswave: pinned to the upper boundrank 1 error0.89its upper bound1.1the lower bound0.69error ⁄ bound0.78error ⁄ lower1.3inside the boundand sitting on it
Fig. 6 And a family that is reproduced exactly at rank two, where both bounds are at the rounding level and the comparison between them is a comparison of two roundings.

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, and the measurement above says it will be very close to exactly ε, because the bound is attained.

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.

The backward error of a hierarchical solve, against the error of the representation it usedSolve A_H x = b exactly and the residual against the matrix that was wanted is b − Ax = (A_H − A)x, so ‖b − Ax‖ ⁄ ‖A‖‖x‖ cannot be anything but the representation's own error. Measured across six accuracies spanning ten decades, the two track at a slope of 0.994 and sit a constant 6.7× apart, which is the difference between a norm of a matrix and a norm of that matrix applied to one vector. The reading is the one this field is built on: the accuracy knob is not an accuracy, it is a backward error chosen in advance. Everywhere else on this site a backward error is something an algorithm produced and somebody then measured; here it is a line in the program, and its size is known before the solve starts.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.99representation at 10⁻⁸1.2·10⁻⁹backward error there1.3·10⁻¹⁰representation at 10⁻¹²1.2·10⁻¹³backward error there2.5·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem
Fig. 7 The hierarchy field’s version of the same sentence, from the essay that establishes it: an accuracy that is a backward error, tracking it at a slope of 1.000 over ten decades.
The accuracy asked for, against the accuracy obtained, for two kernels on one partitionε is applied to each block against that block's own largest singular value; the error is then reported against the whole matrix's norm. The two are not the same number and the dashed diagonal is where they would be. For 1/r the obtained error is 57 times smaller than the tolerance asked for; for log r on the identical partition it is 46 times smaller. Both curves are straight and parallel to the diagonal, so the knob does what a knob should — a decade in buys a decade out — and neither of them sits on it. The factor is how much of the matrix's mass lives off the diagonal, which is a property of the kernel; it is not the number of blocks, which is 16 and would put the curves on the other side of the diagonal.10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³ε asked for, per block‖A − A_H‖ ⁄ ‖A‖ obtainedasked = obtainedlog r1 ⁄ rtwo numbers, not one1/r, obtained at 10⁻⁸1.8·10⁻¹⁰log r, obtained at 10⁻⁸2.2·10⁻¹⁰1/r, obtained ⁄ asked0.018log r, obtained ⁄ asked0.022blocks in the partition16the tolerance is per blockand the error is per matrix
Fig. 8 And what a per-block budget actually delivers there, which is the caveat this page inherits with the rule.

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. 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.

Numbers stored against the number of indices, at n = 20 and rank 4: the tensor, its core, and a trainThe tensor is n^d, which at d = 20 is 1.05·10²⁶. A Tucker representation of it is r^d + d·n·r — the core is still exponential in d, so at rank 4 it is 1.1·10¹², smaller than the tensor by 9.54·10¹³ and still unstorable. A train is (d − 2)nr² + 2nr, which is 5,920 — linear in d. The three lines are the field's whole argument: fixing the definition of the decomposition does not fix the size of what it returns, and the second fix is the same projection cut in a different place.15913172110¹10⁶10¹¹10¹⁶10²¹10²⁶number of indicesnumbers storedthe tensor: n^dthe core: r^d + dnrthe train: (d − 2)nr² + 2nrthe curse, movedentries at d = 2010²⁶core1.1·10¹²train5920core ⁄ train1.9·10⁸entries ⁄ core9.5·10¹³the definition is repairedthe size is not
Fig. 9 The three lines: the tensor, exponential in d; the core, exponential in d with a smaller base; and the train, linear in it, which is the next essay.
Numbers stored against the number of indices, at n = 20 and rank 8: the tensor, its core, and a trainThe tensor is n^d, which at d = 20 is 1.05·10²⁶. A Tucker representation of it is r^d + d·n·r — the core is still exponential in d, so at rank 8 it is 1.15·10¹⁸, smaller than the tensor by 9.09·10⁷ and still unstorable. A train is (d − 2)nr² + 2nr, which is 23,360 — linear in d. The three lines are the field's whole argument: fixing the definition of the decomposition does not fix the size of what it returns, and the second fix is the same projection cut in a different place.15913172110¹10⁶10¹¹10¹⁶10²¹10²⁶number of indicesnumbers storedthe tensor: n^dthe core: r^d + dnrthe train: (d − 2)nr² + 2nrthe curse, movedentries at d = 2010²⁶core1.2·10¹⁸train2.3·10⁴core ⁄ train4.9·10¹³entries ⁄ core9.1·10⁷the definition is repairedthe size is not
Fig. 10 At twice the rank, where the core’s line lifts by 2^d and the train’s by four.
Numbers stored against the number of indices, at n = 20 and rank 3: the tensor, its core, and a trainThe tensor is n^d, which at d = 20 is 1.05·10²⁶. A Tucker representation of it is r^d + d·n·r — the core is still exponential in d, so at rank 3 it is 3.49·10⁹, smaller than the tensor by 3.01·10¹⁶ and still unstorable. A train is (d − 2)nr² + 2nr, which is 3,360 — linear in d. The three lines are the field's whole argument: fixing the definition of the decomposition does not fix the size of what it returns, and the second fix is the same projection cut in a different place.15913172110¹10⁶10¹¹10¹⁶10²¹10²⁶number of indicesnumbers storedthe tensor: n^dthe core: r^d + dnrthe train: (d − 2)nr² + 2nrthe curse, movedentries at d = 2010²⁶core3.5·10⁹train3360core ⁄ train10⁶entries ⁄ core3·10¹⁶the definition is repairedthe size is not
Fig. 11 And at rank three, which is where the core and the train cross for these shapes.

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_k × n^{d−1}: 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_k × n_k, 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 in the least-squares field.

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_k 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.

Least squares by QR and by the normal equations in binary32Relative error of the computed coefficients against epsilon, on log axes. The QR route is a flat line near the bottom; the normal-equations route climbs and then stops, at the epsilon where the cross-product matrix becomes exactly singular.10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴
Fig. 12 The road not taken, from the least-squares field: the forward error of two routes to the same answer, on a problem whose condition number the first of them squares.
Fitting the same degree-11 polynomial in two basesOn the left, the data and two fitted curves that lie on top of each other. On the right, the condition numbers of the two design matrices, ten orders of magnitude apart.00.10.20.30.40.50.60.70.80.910.20.40.60.81xyboth fits, drawn on top of one anothercondition number of the design matrixmonomial1.2·10⁸Chebyshev2.5largest fitted coefficientmonomial113Chebyshev0.51rms residual: 3.9·10⁻⁶ and 2.9·10⁻⁷ — the data is fitted either way.30 points, degree 11, single precisionthe basis is part of the problem
Fig. 13 And what squaring a condition number costs before any arithmetic happens.

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.

What makes that 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.

The error of a rank-4 Tucker representation built from a sketch, against the columns of oversamplingThe flat line is the deterministic decomposition — d matrix SVDs — at 8.81·10⁻⁴, and nothing random can go below it, since every curve here is a projection onto a subspace of the same size. The upper two are medians over five seeds: a dense Gaussian sketch, which the randomised field's bounds cover, and a Khatri–Rao sketch, whose columns are outer products of small random vectors and which no bound in that field applies to. With no oversampling the dense sketch is 3.65 times the deterministic answer; at twelve extra columns it is 1.15. The structured one costs 6.8 per cent more at four extra columns and 9.5 per cent at twelve — measured, because there is nothing else to say about it.-113579111310⁻³10⁻²extra columns in the sketchrelative errora Khatri–Rao sketch: no bound covers ita dense Gaussian sketchthe decomposition, which nothing beatsrank 4, five seedsdeterministic8.8·10⁻⁴dense, p = 00.0032dense, p = 120.001structured, p = 120.0011structured ⁄ dense1.1a bound and a measurementand only one of them is available
Fig. 14 The repaired version, from the essay that owns it: three curves, none of which rises.

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.

A 256 × 256 kernel matrix partitioned by the strong rule, with each compressed block's rank112 blocks: 46 kept dense and 66 stored as two thin factors, whose ranks run from 4 to 5. The rule decides one thing — whether a pair of index clusters may be compressed — and it decides it from four numbers about where those clusters sit, before any entry of the matrix is read. The strong rule refuses any pair whose clusters touch, and subdivides instead, so the diagonal is fringed with small dense blocks and no rank on the picture exceeds 5. The whole thing stores 27,008 numbers against 65,536 entries, and reproduces the matrix to 3.41·10⁻¹⁰.545545545455554545545545545454555555454545545545545455554545545545rows and columns, in the order the points arrivedark: kept dense · light: two thin factors, rank printedthe strong partitionblocks112kept dense46largest rank5numbers stored2.7·10⁴‖A − A_H‖ ⁄ ‖A‖3.4·10⁻¹⁰the picture is decidedbefore a number is read
Fig. 15 The other format’s structure, decided from geometry before any entry is evaluated — the thing this one has no analogue of.
The rank of an admissible block and of a touching one, against how finely they are sampledBoth blocks are of the kernel 1/r; both are 32, 64, 128, 256 points a side; both are truncated at 10⁻⁸. The admissible pair — [0, 1] against [2, 3] — needs 5, 5, 5, 5 columns, which is one number. The touching pair — [0, 1] against [1, 2] — needs 9, 11, 12, 13, climbing by about one per doubling, which is a logarithm. Neither of them grows like the block, and only one of them stops. That difference is what the admissibility test in a partition is buying, and it is why the touching pair is kept dense rather than compressed at all.45678903691215log₂ of the points a sidecolumns above 10⁻⁸two intervals that touchtwo that do notthe measurement that is a null resultadmissible, n = 325admissible, n = 2565touching, n = 329touching, n = 25613stored ⁄ dense at largest0.039the rank belongs to the geometryand not to the sampling
Fig. 16 And its central measurement, for reading against the mode spectra above: a rank that does not grow with the size of the block.

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.

At other settings

The share of a HOSVD core's energy on its superdiagonal, by family, at rank 6A matrix SVD hands over orthonormal factors and a diagonal middle at the same time. For three indices they come apart, and this is the half that does not survive. Every core here is all-orthogonal — the largest inner product between two slices perpendicular to a mode, relative to the core's own energy, is 3.5·10⁻¹⁶ — and none of them is diagonal. The bars are the fraction of the squared norm carried by the 6 entries on the superdiagonal: smooth 97.7%, hilbert 94.5%, wave 26.6%, noise 0.4%. Keeping only those entries costs 0.152, 0.235, 0.857, 0.998 in relative error against the full core's 7.94·10⁻⁶, 1.77·10⁻⁶, 1.34·10⁻¹⁵, 0.842.smooth97.7%hilbert94.5%wave26.6%noise0.4%share of the core's energy on its 6 superdiagonal entrieskeeping only them: 0.152 against 7.94·10⁻⁶keeping only them: 0.235 against 1.77·10⁻⁶keeping only them: 0.857 against 1.34·10⁻¹⁵keeping only them: 0.998 against 0.842orthogonal, and not diagonalsmooth on-diagonal0.98hilbert on-diagonal0.94wave on-diagonal0.27noise on-diagonal0.0044worst slice pair3.5·10⁻¹⁶the slices are orthogonalthe core is not diagonal
Fig. 17 The property the format keeps and the one it cannot, which is the next essay: slices that are mutually orthogonal and a core that is not diagonal.
The share of a HOSVD core's energy on its superdiagonal, by family, at rank 2A matrix SVD hands over orthonormal factors and a diagonal middle at the same time. For three indices they come apart, and this is the half that does not survive. Every core here is all-orthogonal — the largest inner product between two slices perpendicular to a mode, relative to the core's own energy, is 3.5·10⁻¹⁶ — and none of them is diagonal. The bars are the fraction of the squared norm carried by the 2 entries on the superdiagonal: smooth 97.6%, hilbert 94.5%, wave 26.6%, noise 0.2%. Keeping only those entries costs 0.154, 0.235, 0.857, 0.999 in relative error against the full core's 0.0671, 0.0447, 1.32·10⁻¹⁵, 0.971.smooth97.6%hilbert94.5%wave26.6%noise0.2%share of the core's energy on its 2 superdiagonal entrieskeeping only them: 0.154 against 0.0671keeping only them: 0.235 against 0.0447keeping only them: 0.857 against 1.32·10⁻¹⁵keeping only them: 0.999 against 0.971orthogonal, and not diagonalsmooth on-diagonal0.98hilbert on-diagonal0.94wave on-diagonal0.27noise on-diagonal0.0017worst slice pair3.5·10⁻¹⁶the slices are orthogonalthe core is not diagonal
Fig. 18 The same measurement at rank two, where the compressible families put nearly all their energy on the diagonal and the incompressible one puts none.
Numbers stored against the number of indices, at n = 20 and rank 2: the tensor, its core, and a trainThe tensor is n^d, which at d = 20 is 1.05·10²⁶. A Tucker representation of it is r^d + d·n·r — the core is still exponential in d, so at rank 2 it is 1.05·10⁶, smaller than the tensor by 9.99·10¹⁹ and still unstorable. A train is (d − 2)nr² + 2nr, which is 1,520 — linear in d. The three lines are the field's whole argument: fixing the definition of the decomposition does not fix the size of what it returns, and the second fix is the same projection cut in a different place.15913172110¹10⁶10¹¹10¹⁶10²¹10²⁶number of indicesnumbers storedthe tensor: n^dthe core: r^d + dnrthe train: (d − 2)nr² + 2nrthe curse, movedentries at d = 2010²⁶core10⁶train1520core ⁄ train690entries ⁄ core10·10¹⁹the definition is repairedthe size is not
Fig. 19 At the smallest rank drawn, where the core is small enough that the factors are most of the storage.
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. 20 The three matrices for the Hilbert family, whose spectra decay fast enough that the ranks are close together.
The three unfoldings of a 10 × 11 × 12 noise 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 10, 11, 12 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 14.2, 13.6, 13.4, 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 10mode 2 · rank 11mode 3 · rank 12noise: three matrices, one arrayentries1320mode-1 rank10mode-2 rank11mode-3 rank12‖T‖36three ranksand none of them is the tensor's
Fig. 21 And for a tensor with no structure, where all three are full and the format has nothing to find.
Train ranks at each of the 4 cuts of a 5-index tensor on 6 points a side, four familiesCut the index list after position k, put the first k indices on the rows and the rest on the columns, and take the rank of the matrix that results. There are 4 such cuts and each rank is an ordinary matrix rank. sin(x₁ + … + x_d) has rank exactly two at every one of them, for every d, because the addition formula separates it into two terms at every cut — a rank written down rather than measured, and the computed values are 2, 2, 2, 2. The reciprocal family climbs to 9, the product family is one everywhere, and independent normal entries reach 36, which is the largest rank the cut allows. Storage runs sinsum 96, reciprocal 1,206, product 30, noise 10,440 against 7,776 entries.012345061218243036cut after index krank of the reshapesinsum · 2 2 2 2reciprocal · 6 9 9 6product · 1 1 1 1noise · 6 36 36 64 cuts, 4 rankssinsum stored96reciprocal stored1206product stored30noise stored10⁴entries7776one rank per cutand one of them is a theorem
Fig. 22 The next essay’s version of the same picture: ranks of cuts rather than of modes, which is what makes the storage linear in d.
Error of the best rank-k approximation to a 12×12 matrixApproximation error against k on a logarithmic axis, with the measured error and the next singular value drawn as separate curves lying exactly on top of one another.123456789101110⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1rank k of the approximation‖A − Aₖ‖measured, 2-normσₖ₊₁, from theorymeasured, Frobeniusthe first two agreeto 2·10⁻⁸worst |‖A−Aₖ‖₂ − σₖ₊₁| / σₖ₊₁2·10⁻⁸worst Frobenius discrepancy2·10⁻⁸κ = 10⁹; 30 random rank-3 matrices, none closerthe error is σₖ₊₁
Fig. 23 The matrix theorem this is the generalisation of, from the spectra field: an error that is the next singular value and is attained.
Singular values of a rank-4 matrix with noise of relative size 10⁻⁸Ten singular values on a logarithmic axis. The first four sit near one; the rest sit at the noise level, and the vertical distance between the two groups is the evidence for the rank.1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 24 And the decision behind every rank on this page, which for a multilinear rank is three decisions rather than one.
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. 25 The hierarchy field’s storage curve, for reading against the r^d above: a logarithm, on a format whose cost is also chosen.
Random numbers a sketch needs, against the number of indices, at n = 8A dense Gaussian sketch of a mode-k unfolding multiplies an n × n^{d−1} matrix by a random one with n^{d−1} rows, so the random object is the size of the tensor divided by n — the line that crosses the tensor's own entry count at d = 2 and is 1.51·10⁸ by d = 8. A Khatri–Rao sketch replaces it with d − 1 small Gaussians per column, costing 4,032 numbers at the same point — a factor of 3.74·10⁴ — and is applied as mode products, so nothing the size of the tensor is ever formed. What it is not is Gaussian, which is what the previous figure has to measure rather than bound.12345678910¹10³10⁵10⁷number of indicesrandom numbers drawna dense Gaussian sketchdashes: the tensor's own entriesa Khatri–Rao sketcha random matrix nobody can afforddense at d = 81.5·10⁸structured4032the tensor's entries1.7·10⁷dense ⁄ structured3.7·10⁴crossing at d2the sketch outgrows its tensorand the structured one does not
Fig. 26 And what it costs to build this decomposition from a sketch instead, from two essays on.

Where this goes next

110¹10²10³10⁻³10⁻²10⁻¹110¹10²10³n‖A_n − A‖ and the largest term's norm‖A_n − A‖the larger of its two termsan infimum that is not attainedn1024‖A_n − A‖0.0017largest term1024their product1.7√31.7the distance goes to zeroand nothing reaches it A nearest point that is not there Eckart and Young guarantee that a matrix has a best rank-k approximation and that the truncated SVD is it. For three indices the guarantee is false in the strongest available way — there are tensors whose distance to the rank-two set is zero and which no rank-two tensor equals. 12345678910⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1rank k of the approximation‖A − Aₖ‖measured, 2-normσₖ₊₁, from theorymeasured, Frobeniusthe first two agreeto 4.3·10⁻⁹worst |‖A−Aₖ‖₂ − σₖ₊₁| / σₖ₊₁4.3·10⁻⁹worst Frobenius discrepancy4.3·10⁻⁹κ = 10⁹; 30 random rank-3 matrices, none closerthe error is σₖ₊₁ The best approximation there is The error of the best rank-k approximation is not bounded by the next singular value. It is equal to it. That is an unusually sharp theorem, and it makes the theorem itself usable as an independent check on the computation. 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 The format that does not notice the dimension A Tucker core is r^d numbers, so the format that repaired the definition still cannot go past five indices. Cutting between the indices rather than across them gives d − 1 ranks instead of d, storage linear in the number of indices, and a family whose ranks are two everywhere by an addition formula.

Reads more easily once this is understood

Essays that name this one as worth reading first.

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.

Alternating least squaresCurse of dimensionalityEckart–YoungHigher-order SVDLow-rank approximationMultilinear rankOrthogonalityTruncated svdTucker decompositionUnfolding