Unique, and not determined

A factorisation that is unique for once measured the property that makes the CP model worth its other troubles. A rank-$r$ matrix factorisation is never unique, since $AB = (AM)(M^{-1}B)$ for any invertible $M$; for three indices, Kruskal’s condition — the k-ranks of the three factors summing to at least $2r + 2$ — makes the decomposition unique up to permuting and scaling its terms, and it holds generically. Fits from different starts on a tensor that meets it agreed; on one that missed it by two they did not. The essay closed by calling it “the only place in this collection where a decomposition’s factors are the answer rather than a coordinate system for it”.

A k-rank is an integer. Two columns are either parallel or not, a set of $k$ columns is either independent or not, and the condition either holds or it does not. Every other uniqueness in this collection comes with a number for how unique: the gap decides the eigenvector found an eigenvector determined to the extent its eigenvalue is separated from the next. This essay asks what number goes with Kruskal’s yes, by walking a tensor up to the edge of the condition while the condition keeps saying yes.

A condition that holds all the way #

The tensors are $4 \times 4 \times 4$ and of rank three, built from three random $4 \times 3$ factors $A$, $B$ and $C$. Then the third column of $C$ is turned toward its first, keeping its length, to an angle $\theta$ from it: 1 radian, 0.3, 0.1, 0.03, 0.01, 0.003 and 0.001. Every pair of columns of $C$ stays independent at every $\theta$ above zero, so the k-ranks are 3, 3 and 3 and the sum is nine against the eight the condition needs. At $\theta = 0$ the turned column is parallel to the first, $C$’s k-rank drops to one, the sum is seven, and the condition fails. Three random tensors are built this way, from three seeds.

Each tensor is then perturbed by one part in $10^8$ of its norm, in a random direction, and two things are computed. The exact decomposition of the perturbed tensor, found by Gauss–Newton from the original factors, which converges in a few steps because the starting point is a perturbation away; and what alternating least squares returns from the same start, run until its own stopping test fires or to a cap of 3,000 sweeps. Distances between decompositions are measured as the largest relative change of a rank-one term, which ignores the scaling within a term as uniqueness does, and are reported in units of the perturbation.

The decomposition’s own condition is the angle #

The first number is the Jacobian’s. The map from the 36 numbers in $A$, $B$ and $C$ to the 64 entries of the tensor has a $64 \times 36$ Jacobian, and six of its singular values are zero to rounding — the six directions that rescale one column of a term and unscale another, which are exactly the freedom uniqueness allows. The smallest of the other thirty is the decomposition’s sensitivity: a perturbation of the tensor along its singular vector moves the factors by the perturbation’s size divided by it.

Over the largest singular value it is 0.122, 0.040 and 0.013 at $\theta$ = 1, 0.3 and 0.1 on the first tensor, and from there it is exactly proportional to the angle: $0.1331\,\theta$ to four digits at every $\theta$ from 0.1 down to $10^{-3}$. The other two tensors give $0.0968\,\theta$ and $0.0867\,\theta$, each constant to ten per cent over the same range. The condition number of the decomposition, the inverse of that ratio, is about ten over $\theta$: 75 at $\theta = 0.1$, 7,500 at $10^{-3}$. Kruskal’s verdict at those two angles is the same.

One number is not the whole of the Jacobian, and the rest of its spectrum says how much of the decomposition the angle reaches.

Of the thirty non-trivial singular values, exactly two fall with the angle. On the first tensor they are 0.13 and 0.12 of the largest at $\theta = 1$, in line with the 0.16 to 0.21 of the three above them, and at $10^{-3}$ they are $2.1 \cdot 10^{-4}$ and $1.3 \cdot 10^{-4}$, each a hundred times smaller than at 0.1. The third smallest sits at 0.157 at every angle from 0.1 down, and none of the twenty-eight above the falling pair moves by more than 1.9 per cent over that range. The other two tensors draw the same picture: two singular values falling by factors of 99 to 102 between 0.1 and $10^{-3}$, and twenty-eight moving by at most seven per cent.

Two is the count the limit predicts. At $\theta = 0$ the two terms share a third-factor direction $c$, and their sum is a rank-two matrix, $a_1 b_1^T + a_3 b_3^T$ up to the columns’ lengths, times that one vector. A rank-two matrix splits into two rank-one pieces in a four-parameter family of ways: any invertible two-by-two change of basis applied to the $a$’s and undone on the $b$’s. Two of those four parameters only rescale the pieces, and they are among the six directions already set aside; the other two mix the pieces, and they are the pair that goes soft. Above zero the mixing is no longer free, but it is cheap, and its price is proportional to the angle that separates the columns.

The answer moves #

The figure at the top of the page is what that condition does to an actual perturbation. On the first tensor the exact decomposition of the perturbed tensor moves by 1.61, 1.59 and 1.60 perturbations at $\theta$ = 1, 0.3 and 0.1 — the random perturbation happens to lie mostly along well-determined directions while the angle is large — and then by 6.1, 19.5, 66 and 200 at 0.03, 0.01, 0.003 and $10^{-3}$, proportional to $1/\theta$ to within a third from 0.03 down. On the other two tensors the moves at $10^{-3}$ are 422 and 311. A tensor known to eight digits determines its unique decomposition to between five and six.

That is the number Kruskal’s condition does not carry. The condition is a statement about the factors’ independence, which is a yes or a no; the conditioning is a statement about how independent, which is a distance — here the angle, directly. A rank that is not a property of the tensor found the same contrast for rank itself: the integer that a definition produces and the continuous quantity a computation meets are different objects, and the second decides what can be computed.

Alternating least squares slows, then stops #

The method nearly every CP code uses is alternating least squares, and it gets every advantage here: it starts at the exact factors of the unperturbed tensor, a perturbation away from the answer. At $\theta = 1$ it reaches the moved decomposition in 28, 25 and 47 sweeps on the three tensors and lands within $10^{-3}$ of a perturbation of it. At $\theta = 0.1$ it needs 284, 734 and 532 sweeps, ten to twenty-nine times as many, and still lands within a tenth of a perturbation. At $\theta = 0.03$ the first tensor takes 2,446 sweeps and the other two hit the cap. At 0.01 and 0.003 all three hit the cap.

At $10^{-3}$ something different happens. ALS stops after 11, 17 and 65 sweeps, by its own test: the change in its error between sweeps has fallen below $10^{-16}$ of the error, which is the test’s definition of convergence. Its terms at that point are 200, 423 and 311 perturbations from the decomposition — exactly as far as the decomposition moved, because ALS has not moved at all. A still error is not a settled one found the same failure from the other side, a stopping test fooled by an error that had stopped moving; here the error stops moving because a sweep of ALS, solving for one factor with the other two fixed, can only move along directions each factor controls alone, and the poorly determined direction is a coordinated motion of all three.

The residual cannot tell #

The residual does nothing to warn of it. At $\theta = 0.01$ on the first tensor ALS’s relative residual falls within half a per cent of the exact decomposition’s in its first sweeps and spends the remaining three thousand closing the last two tenths of a per cent, while its terms close from twenty perturbations away to twelve and a half. At $\theta = 1$ the residual and the terms arrive together, in 28 sweeps. On the dial the residual’s curve looks almost the same at every angle; the terms are what differ, and the residual does not show them.

Across all 21 fits the residual is within 0.6 per cent of the exact decomposition’s. The fits at $\theta \ge 0.1$ have found the decomposition to between $5 \cdot 10^{-4}$ and $2 \cdot 10^{-2}$ of its move; the fits at $\theta \le 0.03$ are between 0.08 of its move and all of it away. Every one of them has a residual at the noise level, which is the only number a fit on real data can report, since there the decomposition is not known. A small residual is not a small error is this collection’s oldest warning and it has seldom been sharper: a residual matching the true decomposition’s to a per cent, from factors that are the true decomposition’s to between zero and six digits.

Which terms are uncertain #

The move is not spread over the decomposition. Split term by term, at $\theta = 10^{-3}$ on the first tensor, the two terms whose third-factor columns are nearly parallel move by 200 and 100 perturbations and the third term by 1.6. On the second tensor, 420 and 250 against 0.47; on the third, 310 and 110 against 4.8. The well-separated term is determined about as well at $\theta = 10^{-3}$ as at $\theta = 1$, and the two nearly parallel terms carry the whole of the decomposition’s sensitivity between them. Those are the two soft directions of the Jacobian’s spectrum, the mixings counted above: a little of one term traded for a little of the other, in a way that changes the tensor only at second order when their columns nearly coincide.

Alternating least squares follows the same split. At $\theta = 10^{-3}$ it recovers the well-separated term to within 0.07, 0.09 and 0.47 of a perturbation on the three tensors, as well as it does at any angle, and misses the two nearly parallel terms by their full move. So a practitioner reading the factors of such a fit is right about one component and wrong about two, and nothing in the fit says which: the residual, as the section above showed, is the same either way. The gap decides the eigenvector found exactly this shape in a symmetric matrix, an invariant subspace determined while the vectors inside it are not, and here the subspace is the span of the two nearly parallel terms.

Two steps of Gauss–Newton #

The exact decomposition was found by Gauss–Newton, and how fast it got there is the measure of what alternating least squares is missing. From the unperturbed factors, one Gauss–Newton step lands within $3 \cdot 10^{-4}$ of a perturbation of the final answer at $\theta = 10^{-3}$ on the first tensor, and the second step reaches the rounding floor, a few millionths of a perturbation. At $\theta = 1$ the first step is already at the floor. On all three tensors and every angle, two steps are enough. A Gauss–Newton step solves one least-squares problem with the full $64 \times 36$ Jacobian, so it moves all three factors at once along whatever direction reduces the residual most, including the coordinated one; an ALS sweep solves three small problems, each with two factors frozen, and the coordinated direction is not available to any of them. The cost per step is larger and the count is smaller by three orders of magnitude.

What uniqueness is still worth #

None of this weakens the theorem; it locates it. Kruskal’s condition says the decomposition of an exact tensor is unique, and on every tensor here it is. What a practitioner needs to read factors as chemical species or sources is that the decomposition of a measured tensor is close to the decomposition of the true one, and that is a statement about conditioning, which the condition does not make. At $\theta = 10^{-3}$, eight-digit data give five- to six-digit factors; at $10^{-5}$, by the same proportion, they would give three to four.

An iteration that walks out of the set found alternating least squares plateauing where the nearest rank-$r$ point does not exist, and the test that is a deadline found the quantity that detects that plateau firing on badly conditioned fits that do converge, at sweep 758. This is that badly conditioned case, built on purpose, and a swamp’s milder cousin: no terms diverge, the answer exists and is unique, and the iteration slows anyway because the answer is poorly determined in one coordinated direction. A Gauss–Newton step, which moves all three factors at once along the Jacobian’s directions, crossed the same distance in a handful of steps; that is the difference between an algorithm that sees the coordinated direction and one that cannot.

A number to report beside the factors #

Everything that decides whether these factors can be read is available from the fit, at a price. The Jacobian’s smallest non-trivial singular value took one singular value decomposition of a $64 \times 36$ matrix at this size; for a tensor of $n^3$ entries and rank $r$ the Jacobian is $n^3 \times 3nr$, and its smallest singular values can be estimated by a few iterations of an alternating method on the normal equations without forming it. Reported beside the factors, as $0.13\,\theta$ was here, it says to how many digits the factors are determined by data of a stated accuracy: the data’s relative error divided by it. At $\theta = 10^{-3}$ that is eight digits of data and four of factors at worst — five to six in the random directions measured above — an honest summary of what the measurement found, and one Kruskal’s condition, which said only yes, could not give.

The singular vector that goes with it says which terms the uncertainty lives in, since its components are the coordinated change of the factors that moves the tensor least. Here that would have named the two nearly parallel terms and cleared the third, which is the split the per-term measurement found. A CP code that returned its factors with the Jacobian’s smallest singular value and its singular vector would return, for the price of one more decomposition, both the answer and the extent to which it is one.

None of this is new to the theory; the Jacobian’s role in CP identifiability, local and global, is standard in the literature on tensor decompositions. What the measurement adds is the size of the gap between the two statements on an ordinary small tensor: Kruskal’s yes at every angle, and between two and four digits of the factors lost by the time the angle is $10^{-3}$, with the residual of the method everyone uses indistinguishable throughout.

What three tensors do not show #

Three random tensors of one small size and rank, with the near-collinearity placed in one pair of columns of one factor. Near-collinearity in two factors at once would multiply the conditioning, and is the case an iteration that walks out of the set met at the boundary of the set itself. One perturbation direction per tensor, random; a perturbation aligned with the Jacobian’s worst singular vector would move the decomposition further than the random direction measured here does, by the condition number rather than by the share of it a random direction happens to pick up. ALS is started at the exact factors, the most favourable start there is; from random starts it would have the swamps to cross as well. And the cap of 3,000 sweeps is a choice: given enough sweeps ALS may eventually arrive at the small angles, and how many it needs is the question the measurement above stops short of.

Still open: how long the iteration needs, and a condition with a number in it #

The sweeps it needs. The sweeps to reach the moved decomposition rose from about thirty at $\theta = 1$ to several hundred at 0.1. The prediction with a sign is that they grow like $1/\theta^2$ — the square of the decomposition’s condition, as for any iteration whose contraction is set by the ratio of a Gram matrix’s extreme eigenvalues — so that at $\theta = 0.01$ ALS from the exact start needs between twenty and a hundred thousand sweeps, and a cap of three thousand misses by more than a factor of five.

A Kruskal test with a tolerance. The condition can be made numerical by asking whether every set of $k$ columns is independent to a tolerance $\tau$ rather than exactly. The prediction is that the tolerance at which the numerical k-rank of the third factor drops from three to one is within a factor of two of $\theta$ on every tensor here, and that setting $\tau$ to the data’s relative noise gives a test whose failures are exactly the fits whose factors move by more than a hundred perturbations — a version of the condition that carries its own conditioning.