A nearest point that is not there

Half of this collection rests on one theorem. Eckart and Young: the best rank-k approximation of a matrix in the Frobenius or spectral norm is its truncated singular value decomposition, the error is the next singular value, and — the part nobody states, because for matrices it needs no stating — the best one exists.

For an array with three or more indices it is false. Not approximately false, not false in a corner case: there are tensors whose distance to the set of rank-two tensors is exactly zero and which no rank-two tensor equals.

The example, which is eight numbers #

With e₁ = (1, 0) and e₂ = (0, 1) in two dimensions, let

A = e₁⊗e₁⊗e₂ + e₁⊗e₂⊗e₁ + e₂⊗e₁⊗e₁

which is a 2 × 2 × 2 array with three entries equal to one and five equal to zero. Its rank is three: three rank-one terms are written above, and no two suffice.

Now consider, for each n,

Aₙ = n · (e₁ + e₂/n)⊗³ − n · e₁⊗³

which is a difference of two rank-one terms and therefore has rank at most two for every n. Expanding the cube,

Aₙ = A + (1/n)·(the three terms with two e₂’s) + (1/n²)·e₂⊗e₂⊗e₂

so ‖Aₙ − A‖ = √(3/n² + 1/n⁴), which goes to zero.

The infimum of the distance from A to the rank-two set is therefore zero, and A has rank three, so nothing in the set attains it. The set of tensors of rank at most two is not closed.

What that costs, in numbers #

The hero figure plots both quantities. The residual falls like √3/n and the larger of the two rank-one terms grows like $n(1 + 1/n^2)^{3/2}$, so their product is

$$\sqrt{3 + 1/n^2}\;(1 + 1/n^2)^{3/2}$$

which descends onto √3 from above: 1.7774 at n = 8, 1.73205 at n = 2048, a fall of 2.6 per cent across two and a half decades. That expression is what the assertion checks, rather than a constant, because the algebra hands over a formula and asserting “constant” where a formula exists is describing the picture instead of testing it.

The practical reading of a constant product is severe. Getting one digit closer to A costs a factor of ten in the size of the pieces being subtracted from one another, for ever. At n = 10¹⁶ the two terms are 10¹⁶ and the answer is their difference — which is the arithmetic this collection’s first field is about, and where the cancellation takes the whole answer.

Both laws are exact over the four decades the slider covers, which is what makes the product constant rather than nearly constant.

Across n = 16, 64, 256, 1,024, 4,096, 16,384 and 65,536 the residual reads 0.108, 0.0271, 0.00677, 0.00169, 4.23·10⁻⁴, 1.06·10⁻⁴ and 2.64·10⁻⁵ — each one a quarter of the last, at every stop, which is 1/n with no correction term — while the larger of the two rank-one terms reads 16.09, 64.02, 256, 1,024, 4,096, 16,380 and 65,540, which is n. The product descends onto √3 and arrives: 1.74334, 1.73276, 1.73209, then 1.73205 at every remaining stop, to six digits.

Those two exact laws fix where the sequence stops being computable, and the answer is the square root of the unit roundoff. The residual is √3/n; the terms are of size n, so forming their difference in a format with unit roundoff u loses an absolute nu. The best attainable accuracy is where the two meet, √3/n = nu, at n = (√3/u)^(1/2) ≈ 1.3·10⁸ in double precision — and the residual there is √3/n ≈ 1.4·10⁻⁸, which is √u.

So a code chasing A along this sequence in double precision cannot get closer than about eight digits, however long it runs, and the n at which it does best is 10⁸ rather than anything larger. That is the same √u threshold this collection meets in the road that squares the problem and in the definiteness test, arriving here from geometry rather than from a Gram matrix: the infimum is zero in the algebra, and half the digits is what the arithmetic will pay for it.

The word that is missing, and the one that replaces it #

What A has and Aₙ does not is a rank; what Aₙ has and A does not is membership of a set. The quantity that reconciles them has a name.

The border rank of a tensor is the least r such that it is a limit of tensors of rank at most r. For a matrix, border rank and rank are the same number — that is what closedness means. For A they are two and three, and the gap between them is the whole subject of this page.

Border rank is the more useful of the two for anything an algorithm does, because it is what a truncation can approach. Rank is the more useful for anything a decomposition claims, because it is what a decomposition returns. A field that could only quote one of them would be describing either the algorithms or the objects, and this collection’s habit of putting a computed number beside a claimed one requires both.

There is a third quantity, and it is the one every practical code actually uses. A numerical rank at a tolerance ε is the least r for which some rank-r tensor is within ε — which is a border rank in disguise, since the set being searched is the closure. So the object a code reports is not the rank of the tensor and is not a mistake; it is a different quantity, and this collection has an essay about the matrix version of exactly that confusion.

The rank of a 2 × 2 × 2 tensor is a sign #

The example above is not an isolated construction. For arrays of this smallest interesting shape the rank is decided exactly, by a polynomial in the eight entries, with no iteration anywhere.

Take the two 2 × 2 slices A₀ and A₁ and form det(A₀ + λA₁), which is a quadratic c₀ + c₁λ + c₂λ². Its discriminant

Δ = c₁² − 4c₀c₂

is Cayley’s hyperdeterminant, and it decides the rank outright:

The example has A₀ = [[0, 1], [1, 0]] and A₁ = [[1, 0], [0, 0]], so det(A₀ + λA₁) = −1 identically: c₀, c₁ and c₂ are −1, 0 and 0, and Δ = 0. Its three unfoldings all have rank two, which is what separates the rank-three case from the rank-one one.

Both of those show the sequence from the inside: terms growing, residual falling, and no limit in the set to arrive at. The other way to see it is to stop following the sequence and look instead at the neighbourhood it is travelling through.

Why the matrix case has no analogue of this #

It is worth being precise about what the matrix case has that this does not, because the difference is not that matrices are simpler.

The set of matrices of rank at most k is closed: it is the zero set of all (k+1) × (k+1) minors, which is a system of polynomial equations, and the zero set of polynomials is closed by definition.

The set of tensors of rank at most r is not a zero set of anything. It is the image of a polynomial map — take r triples of vectors, form the sum of their outer products — and the image of a polynomial map need not be closed. Its closure is a variety, and the difference between the image and its closure is exactly the tensors whose border rank is r and whose rank is more.

So the failure is not a defect in the definition of tensor rank; it is the ordinary behaviour of images of maps, and the matrix case is the accident. What makes the matrix case work is that rank-≤k matrices happen to be describable by equations as well as by a parametrisation, and nothing hands that over for three indices.

One more difference is worth naming because it is the one that bites in practice. For matrices, the rank-≤k set is irreducible and its dimension is known: k(m + n − k). For tensors, the set of rank-≤r tensors can have several components of different dimensions over the reals, which is what the next essay’s two typical ranks are a symptom of. So even the question how many parameters does a rank-r tensor have does not have one answer.

What it does to an algorithm #

An optimisation over a set that does not contain its own limit points has a specific and recognisable failure, and it is not divergence.

The objective decreases, monotonically, for as long as anyone is prepared to wait. The parameters — here the rank-one terms — grow without bound. The two are in exact proportion, so the iterate keeps improving and the representation keeps degenerating, and neither an objective plot nor a residual test can tell that from slow convergence.

This collection has a name for that shape from its own hierarchy field: a quantity going to zero is one arithmetic event with two opposite meanings. Here it is a quantity going to zero that means the answer is being left behind rather than approached.

The field’s essay on alternating least squares measures it. What is worth carrying from here is that the phenomenon is geometric and not numerical: it is not a rounding effect, it does not go away in higher precision, and it is not a bad starting point. It is the set.

And what it does to a bound #

The other consequence is that a whole class of statement stops being available.

For matrices, every low-rank approximation on this site is priced against the best one: the hierarchy field’s blocks are compared with bestError, the randomised field’s sketches are compared with the truncated decomposition, and the word quasi-optimal means within a factor of the best. All of that needs the best to exist.

For tensor rank it does not, so the phrase has nothing to attach to. What replaces it is a different notion of rank — the multilinear rank, and the train ranks, which are ranks of ordinary matrices — for which the best does exist and the projection onto it is computable. Two essays on, that is exactly the repair: carry the algorithm across rather than the definition.

What survives #

It is worth stating plainly what is not damaged, because the list is longer than the damage.

The Frobenius norm is still a norm and distances still mean what they mean. The unfoldings are still matrices and their singular values are still computable, exactly and cheaply. Every statement in the previous two essays about Kronecker structure is untouched, because a Kronecker sum is a matrix. And the tensor’s own entries are as available as they ever were.

What is damaged is one specific inference: there is a best approximation of this shape, therefore what an algorithm returns can be compared against it. That inference is used, unstated, in every low-rank argument this collection makes, and the rest of the field is about the two ways of restoring it — change the notion of rank until the set is closed, which is what the Tucker and train formats do, or accept the non-closedness and measure what an iteration does inside it, which is what the alternating-least-squares essay does.

What it means for a claim about data #

One more consequence, and it is the one that reaches outside numerical linear algebra.

A CP decomposition is used because its factors are interpretable: the columns are supposed to be the components a measurement is a mixture of, and the whole reason for preferring it to a matrix factorisation is that they are recoverable rather than arbitrary. The last essay of this field measures that and finds it true.

What this page adds is the condition under which the sentence is meaningful at all. If the best rank-r approximation of the data does not exist, then the components a rank-r model would recover is not a well-formed phrase — there is no rank-r model that is nearest, so there is nothing for the recovered factors to be an estimate of. A fit that has entered a swamp is returning terms of size 10⁶ that cancel to within 10⁻³, and reading those terms as components is reading an artefact of where the iteration stopped.

The diagnostic is on the hero figure and costs nothing: report the size of the terms beside the residual. A model whose terms are the size of its data is a model; a model whose terms are a thousand times its data is a subtraction.

In double precision the infimum is not zero #

Everything above is exact arithmetic. This site’s whole premise is that the exact statement and the computed one are different statements, so it is worth asking what the distance from A to the rank-two set actually is in a working floating-point format.

The construction subtracts two rank-one terms of size n to leave something of size √3/n. Forming them costs about n·u in rounding, so the computed distance is roughly √3/n + c·n·u, which does not go to zero: it bottoms out where the two terms are equal, at n ≈ √(√3/u) ≈ 10⁸, at a distance of a few times √u ≈ 1.05·10⁻⁸. Past that, pushing n further makes the approximation worse.

Measured, that is what happens — and only in the right basis. Writing A and Aₙ in a basis rotated by 30°, the computed distance follows √3/n down to 5.5·10⁻⁸ at n = 10⁷·⁵, turns round, and is 1.5·10⁻⁷ at n = 10⁹ and 2.5·10⁻⁵ at n = 10¹². The floor is 7.5·10⁻⁹, at n = 2.6·10⁸, which is √u to within a factor of one and a half.

In the standard basis it does not happen at all. The computed distance tracks √(3/n² + 1/n⁴) to eleven digits at n = 10¹⁰ and keeps falling. The reason is that the entry which cancels holds exactly n in both terms — n·1·1·1 and n — so the subtraction n − n is exact, and the other entries are n·(1/n), n·(1/n²) and n·(1/n³), every one of them formed without loss. The example that every account of non-closedness uses is arranged, accidentally, so that the arithmetic is free.

A 30° rotation changes nothing mathematically. The Frobenius norm is orthogonally invariant, rank is basis-independent, and border rank is too; A has rank three and border rank two in both bases and the infimum is zero in both. What the rotation changes is whether the reader can see that the infimum is not attainable twice over — once because the set is not closed, and once because the format is not the reals.

So the honest statement of what a fitting code is chasing is narrower than the geometry suggests. The best rank-two approximation of A does not exist; the best computable in double precision rank-two approximation does exist, it is about 10⁻⁸ away, and it is reached at terms of size 10⁸. That number is not a fact about tensors — it is √u, so it is a fact about the format, and in single precision it would be 3·10⁻⁴ and in double-double 10⁻¹⁶. The swamp has a bottom, and the bottom is where the format puts it.

What that changes about the diagnostic #

The section above proposes reporting the size of the terms beside the residual, and the floor sharpens that from a habit into a stopping rule.

The two quantities are not independent: on this family their product is √3 to within three per cent across two and a half decades, so a fit reporting a residual of 10⁻⁶ and terms of size 10⁶ has not found two facts, it has found one. What the floor adds is the third number that makes the pair readable. Terms of size s carry an arithmetic error of about s·u, so a residual below s·u is not a measurement of anything — the iterate is inside the noise of its own representation, and further decrease is the optimiser fitting rounding error.

That is a test a code can apply at every iteration and it needs nothing the fit does not already have: stop when the residual falls below the largest term times the unit roundoff.

One limit on it, measured after this page was written and recorded here rather than in a footnote. The rule is exact on the family above, where the product is √3 and the test fires at n ≈ 10⁸. It is not a swamp detector, because a real fit does not follow this family: an alternating least squares walking into the same geometry reaches an error of 1.4·10⁻³ with terms of 10.15, where this sequence needs 1,230 — a hundred and twenty times more — so the test sits twelve orders from firing after twenty thousand sweeps. What separates the regimes there is the exchange rate between size and accuracy rather than the arithmetic floor, and the alternating-least-squares essay measures it at 1.00 on this sequence, 0.49 in a swamp and zero on a converging fit. On the family here that fires at n ≈ 10⁸ and a residual of 10⁻⁸, which is exactly the floor measured above, and it fires for the right reason rather than on a tolerance somebody chose.

It also distinguishes the two failures that look identical on a convergence plot. A fit converging slowly towards an attained minimum has terms that stay bounded, so the product falls and the test never fires. A fit descending into a swamp has terms growing in exact proportion to its residual, so the product is constant and the test fires as soon as the arithmetic runs out. The plot is the same monotone decrease in both cases, which is the difficulty the alternating-least-squares essay measures, and the product is what separates them.

The limitation is worth stating too. The test says the iterate has stopped being informative; it does not say the tensor has border rank below its rank. A well-conditioned rank-two tensor never approaches the floor at all, and a badly scaled one might approach it for reasons that have nothing to do with the geometry — which is the same caution the numerical-rank essay attaches to every threshold on a small quantity.

The refusal #

The claim under test is the assumption every low-rank argument on this site has made without stating it, and the refusal has to be the other direction: not that the phenomenon exists, but that it does not happen everywhere.

If the classification reported rank three for every 2 × 2 × 2 tensor, the whole page would be a bug. So it is fed an ordinary rank-two tensor — a sum of two rank-one terms with disjoint supports, which is as far from the boundary as this shape allows — and required to refuse the claim that it has rank three.

That matters more than it looks. A phenomenon that applied to every tensor would be a statement about the definition rather than about the geometry, and the measured fact that the light region of the plane figure occupies most of the window is what makes the dark boundary a boundary.

The second refusal in the same file is aimed at the reading that most nearly works. The unfoldings of a tensor are matrices and their ranks are ordinary ranks, so it is tempting to take the largest of them as the tensor’s rank. It is fed the Δ < 0 tensor, whose unfoldings are all rank two and whose real rank is three, and required to fail.