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

An iteration that walks out of the set

Every sweep of alternating least squares is the exact minimiser of its own subproblem, so the objective can only fall. What it cannot do is converge, when the target's nearest rank-r point is not in the rank-r set — and a plateau at a small residual looks identical to slow convergence unless the size of the terms is plotted beside it.
33 min read 26 figures Where the algebra stops being the arithmeticMeasured, not assertedAssertions that reject

Worth reading first: A nearest point that is not there · The road that squares the problem.

The model this essay fits is the one the first two essays of this field are about: a sum of r rank-one terms, whose least r is the tensor’s rank.

Fix every factor but one and the model is linear in what is left, so each of the d subproblems is an ordinary least-squares solve — the object this collection has had from the beginning. Sweep round the modes until nothing moves. Every step is the exact minimiser of its own subproblem, so the objective can only fall, and the method has been the default for four decades.

What it does not have is a limit.

Alternating least squares on a tensor with a rank-two answer and on one without, over 20,000 sweepsThe steeply falling curve is a rank-three fit to a tensor built from rank-three factors: it reaches 9.74·10⁻¹⁵ in 58 sweeps and its largest term does not move — 12.202 at sweep twenty and 12.202 at the end. The other is a rank-two fit to the border-rank tensor, whose error goes from 0.00677 to 0.00141 and has not finished, while the rising curve is its largest term: 4.698 to 10.15, a factor of 2.16. Fitted over 199 points, the error falls as the 2.03 power of the term size — a better path than the explicit sequence's 1/n, and still one that leaves the set rather than converging inside it.110¹10²10³10⁴10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹sweeprelative error, and the largest term's sizerising: the swamp's largest rank-one termfalling, slowly: its errorfalling, once: a fit with an answera plateau with a rising floorswamp error0.0014swamp term10term growth2.2benign error9.7·10⁻¹⁵benign term growth1the error alone cannot tellthe size of the terms can
Fig. 1 Two runs of the same code at the same rank. One finishes in fifty sweeps; the other has an error that keeps falling, terms that keep growing, and twenty thousand sweeps is not enough for either.

One sweep

The arithmetic is worth writing out, because the whole of the essay’s numerical caveat is in one line of it.

For mode k, the model’s unfolding is A_k · (⊙_{i≠k} A_i)ᵀ, where ⊙ is the Khatri–Rao product — the column-wise Kronecker product of the other factor matrices, which is n^{d−1} × r. The least-squares solution for A_k is therefore

A_k = M_k · KR · (⊛_{i≠k} A_iᵀA_i)⁻¹

where M_k is the tensor’s mode-k unfolding and ⊛ is the entrywise product of the small Gram matrices.

The inverse there is of an r × r matrix, and it is formed through the normal equations of that system rather than through a decomposition of the tall Khatri–Rao matrix. That is the road this collection has a whole essay warning about, and it is taken here deliberately: the matrix being squared is r × r with r at most five, and the alternative is a decomposition of an n^{d−1} × r matrix at every step of every sweep.

The consequence is stated rather than hidden. A ridge of the machine epsilon times the trace is added, and what it is for is precisely the case this essay is about — a swamp is a run in which that matrix becomes singular by design, and without the ridge the run would end with a division rather than with a measurement.

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. 2 The road being taken, from the least-squares field: two routes to the same answer, one of which squares the condition number. Here the squared object is 5 × 5 and the alternative is a decomposition of a million rows.
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. 3 And what squaring costs when the matrix is not small, which is why the trade above needs stating rather than assuming.

The objective can only fall

The first measurement is the reassuring one and it is worth having because it establishes what the failure is not.

A rank-three fit to a tensor built from random rank-three factors stops itself after fifty-eight sweeps: the error runs from 0.487 to 9.7·10⁻¹⁵ and every sweep is at least as good as the one before it, which the assertion checks sweep by sweep rather than end to end. There is nothing wrong with the method’s monotonicity, and there is nothing wrong with its arithmetic.

The largest rank-one term over the same run does not move at all: 12.202 at sweep twenty and 12.202 at sweep two hundred, to five figures.

Alternating least squares on a tensor with a rank-two answer and on one without, over 200 sweepsThe steeply falling curve is a rank-three fit to a tensor built from rank-three factors: it reaches 9.74·10⁻¹⁵ in 58 sweeps and its largest term does not move — 12.202 at sweep twenty and 12.202 at the end. The other is a rank-two fit to the border-rank tensor, whose error goes from 0.00677 to 0.00616 and has not finished, while the rising curve is its largest term: 4.698 to 4.918, a factor of 1.05. Fitted over 149 points, the error falls as the 2.07 power of the term size — a better path than the explicit sequence's 1/n, and still one that leaves the set rather than converging inside it.110¹10²10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹sweeprelative error, and the largest term's sizerising: the swamp's largest rank-one termfalling, slowly: its errorfalling, once: a fit with an answera plateau with a rising floorswamp error0.0062swamp term4.9term growth1benign error9.7·10⁻¹⁵benign term growth1the error alone cannot tellthe size of the terms can
Fig. 4 The benign run alone, at the scale where it finishes. Both of its curves are flat by the fiftieth sweep — one at the rounding level and one at the size of the data.

And it can fall for ever

The same code, at rank two, on the tensor whose rank is three and whose distance to the rank-two set is zero.

Over twenty thousand sweeps the error goes from 0.614 to 1.4·10⁻³ and the largest rank-one term goes from 3.98 to 10.15. Neither has finished. The error is still falling at the last sweep and the term is still growing, and a run of two hundred thousand would produce the same picture with both axes extended.

That is a swamp, and it is not a numerical accident, a bad starting point or a precision problem. It is the geometry of the previous essay arriving as a plateau in an iteration count: the sequence of iterates is a sequence of rank-two tensors approaching a rank-three limit, and there is nothing in the set for it to converge to.

The measurement that separates the two readings

A long plateau at a small residual looks the same as slow convergence, and the objective alone cannot tell them apart. What can is the size of the terms.

The trace records, per sweep, the magnitude of each rank-one term — the product of its columns’ norms — rather than the largest column norm. That distinction matters: a single factor’s norm can be held at one by rescaling the others, without changing the model at all, so a plot of column norms can be made to show anything. The product is invariant under that rescaling and is the quantity that has to diverge.

Fitted over two hundred points, the relation between them is

error ∝ (term size)^{−2.03}

which is not the relation the explicit sequence of the previous essay gives. That sequence has error √3/n against a term size of n, so a product that is constant. The path an alternating fit takes is better — at a term size of 10.15 the explicit sequence’s error would be 0.17 and this fit’s is 0.0014 — and it is still walking out of the set rather than converging inside it.

Both statements are useful. The first says the diagnostic is real: a run whose terms are growing in proportion to its residual’s fall is not converging. The second says the rate is a property of the path rather than of the geometry, so quoting the explicit sequence’s constant as a prediction would be wrong.

A rank-two sequence approaching a rank-three tensor: the residual falls like 1/n and the terms grow like nA_n = n(e₁ + e₂/n)⊗³ − n·e₁⊗³ has rank two for every n and converges to a tensor of rank three. The falling curve is ‖A_n − A‖, which is √(3/n² + 1/n⁴) exactly and reaches 0.00169 at n = 1024; the rising one is the norm of the larger of its two rank-one terms, 1024. Their product runs 2.519, 1.917, 1.777 … 1.73205, descending onto √3 = 1.73205. So getting one digit closer costs a factor of ten in the size of the pieces, for ever, and the infimum of the distance is zero while no rank-two tensor attains it.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
Fig. 5 The explicit sequence, from the previous essay, whose product is √3 — the path the algebra hands over rather than the one an optimiser takes.
Alternating least squares on a tensor with a rank-two answer and on one without, over 1,000 sweepsThe steeply falling curve is a rank-three fit to a tensor built from rank-three factors: it reaches 9.74·10⁻¹⁵ in 58 sweeps and its largest term does not move — 12.202 at sweep twenty and 12.202 at the end. The other is a rank-two fit to the border-rank tensor, whose error goes from 0.00677 to 0.00469 and has not finished, while the rising curve is its largest term: 4.698 to 5.618, a factor of 1.20. Fitted over 189 points, the error falls as the 2.06 power of the term size — a better path than the explicit sequence's 1/n, and still one that leaves the set rather than converging inside it.110¹10²10³10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹sweeprelative error, and the largest term's sizerising: the swamp's largest rank-one termfalling, slowly: its errorfalling, once: a fit with an answera plateau with a rising floorswamp error0.0047swamp term5.6term growth1.2benign error9.7·10⁻¹⁵benign term growth1the error alone cannot tellthe size of the terms can
Fig. 6 A thousand sweeps, where the plateau is already indistinguishable from convergence on the error curve alone and the term curve has already moved.

What a practitioner sees

It is worth being concrete about what the failure looks like from outside, because none of the usual stopping rules catch it.

A relative change test stops early and reports success. The error changes by less than 10⁻⁶ between sweeps for tens of thousands of sweeps, so any rule of the form stop when the objective stops moving fires almost immediately, at a residual that looks respectable.

A residual test reports success too. The fit reproduces the data to a fraction of a per cent, which on real data is well inside the noise.

The factors are the output, and they are the artefact. The returned terms have magnitude 10.15 against data of magnitude 1.7, and they cancel. Reading them as components — which is the entire reason for preferring this model to a matrix factorisation — reads the cancellation as structure.

So the diagnostic has to be the term size, it costs one line, and it is the number this page exists to put on a figure.

Alternating least squares on a tensor with a rank-two answer and on one without, over 5,000 sweepsThe steeply falling curve is a rank-three fit to a tensor built from rank-three factors: it reaches 9.74·10⁻¹⁵ in 58 sweeps and its largest term does not move — 12.202 at sweep twenty and 12.202 at the end. The other is a rank-two fit to the border-rank tensor, whose error goes from 0.00677 to 0.00266 and has not finished, while the rising curve is its largest term: 4.698 to 7.42, a factor of 1.58. Fitted over 197 points, the error falls as the 2.04 power of the term size — a better path than the explicit sequence's 1/n, and still one that leaves the set rather than converging inside it.110¹10²10³10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹sweeprelative error, and the largest term's sizerising: the swamp's largest rank-one termfalling, slowly: its errorfalling, once: a fit with an answera plateau with a rising floorswamp error0.0027swamp term7.4term growth1.6benign error9.7·10⁻¹⁵benign term growth1the error alone cannot tellthe size of the terms can
Fig. 7 Five thousand sweeps, which is a generous budget by any practical standard and is a third of the way through the picture.
A plane of 2 × 2 × 2 tensors, coloured by the sign of the hyperdeterminant that decides their rankTwo of the eight entries are varied over ±2 and the other six are held at the values that make the centre of the picture the rank-three tensor A = e₁⊗e₁⊗e₂ + e₁⊗e₂⊗e₁ + e₂⊗e₁⊗e₁. The light region is Δ > 0, where the tensor has real rank two and two real rank-one terms exist; the dark region is Δ < 0, where the real rank is three and the complex rank is two. 54.2 per cent of this window is rank two. The boundary between them is the curve Δ = 0, on which the rank is three and the distance to the rank-two set is zero — every point of it is a tensor with no nearest rank-two approximation.rank three, Δ = 0a₁₁₁ across, a₀₀₀ up, both over ±2light: two real roots, rank two · dark: a conjugate pair, rank threeone sign, two rankscells sampled4096rank two2222rank three1874Δ at the centre0rank at the centre3the rank is a signcomputed from eight numbers
Fig. 8 Where the target has to be for this to happen: on the boundary curve, which has measure zero and is therefore never sampled and always constructed.

Why the objective is the wrong thing to watch

The general point is worth separating from the tensor that occasions it, because it applies to every descent this collection runs.

A monotone objective proves that an algorithm is behaving. It proves nothing about whether the problem has an answer, and the two are usually conflated because for most problems the second is not in doubt. Least squares on a matrix has an answer; a linear system with a non-singular matrix has an answer; conjugate gradients on a positive definite operator has an answer. In every one of those the objective falling is the whole story.

Here it is half the story, and the missing half is a property of the parameters rather than of the objective. That is a shape this collection has met before under other names. The sequence field measures a Newton iteration whose residual falls onto a plateau that is the linearisation’s own error rather than the answer’s. The iterative field measures a conjugate gradient run whose reported residual goes below the unit roundoff while the true residual does not. In both, a quantity that looks like progress is measuring something other than the distance to the answer.

The rule that comes out of all three is the same: report the size of the thing being returned beside the size of what is left over. It costs one line, it is scale-free once the terms are measured as products, and it is the only number on this page that distinguishes the two runs.

The residual conjugate gradients reports and the residual of the vector it holds, along one runA 50×50 diagonal matrix with one eigenvalue at 10⁻¹² and the rest between 1 and 50, so every product with it is exact to a rounding and nothing below can be blamed on the operator. The lower curve is r ← r − αAp, updated by the recurrence, which is what a stopping test reads. The upper curve is ‖b − Ax‖/‖b‖, recomputed from the iterate. They start as the same vector and end 1.6·10¹⁰ apart, with the reported one at 4.01·10⁻²¹ — below the unit roundoff of 1.11·10⁻¹⁶, which is the shortest proof available that it is not the residual of anything.01326395265789110410⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best4·10⁻²¹the answer's, at its best6.2·10⁻¹¹unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10¹¹iterations drawn105the recurrence remembers every roundingand the stopping test is written in it
Fig. 9 The iterative field’s version, from the essay that measures it: a reported residual and a true one, separating.
How far the two residuals drift apart, against how large the iterates got on the waySix runs, differing only in the smallest eigenvalue of the matrix, which sets how far the iterates travel: from 9309 to 9.31·10¹³, ten decades of it. The upper line is the standard bound, u·‖A‖·max‖x‖, which every run stays inside. The lower series is what the drift actually does, and its fitted slope is 0.507 where the bound's is 1 — so the roundings are a random walk rather than a sum, the bound is correct, and it is loose by a square root that reaches 7.5·10⁴ across the sweep.10³10⁶10⁹10¹²10¹⁵10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1largest iterate on the way, ‖x‖difference between the two residualsthe bound, linear in ‖x‖what it does, slope 0.51a bound of one, a walk of a halffitted slope0.51the bound's slope1smallest gap measured5.3·10⁻¹⁵largest gap measured7·10⁻¹⁰share of the bound, near end7·10⁻⁴share of the bound, far end9.3·10⁻⁹the bound is rightand loose by a square root
Fig. 10 And the mechanism behind that separation, which is the same mechanism as the one here — a quantity that grows while the thing it is subtracted from shrinks.

The near-boundary case, which is the common one

A target exactly on the boundary is a measure-zero event. What happens near it is the practical version and it is worse in one respect and better in another.

Near the boundary the best rank-two approximation exists — the set is closed except at the boundary itself, and a target at distance δ from a rank-three tensor with δ > 0 has a nearest rank-two point. So the iteration converges. What it converges at is a rate that degrades as the target approaches the boundary, and to factors whose size grows as the distance shrinks.

The error field’s essay on the conditioning of a decomposition measures exactly that: the condition number of the map from tensor to factors, along a sequence approaching the boundary, is 2n² + 1 — so a target ten times closer has a decomposition a hundred times more sensitive. And the practical consequence measured there is severe: every tensor on that sequence is exactly rank two, with its two terms written down in closed form, and a three-hundred-sweep fit from a random start does not find them.

Along a sequence of exactly rank-two tensors: the step's condition number, and what a fit from a random start achievesEvery A_n on this sequence *is* a rank-two tensor and its two rank-one terms are written down, so nothing here is about existence. The rising curve is the condition number of the r × r system each alternating sweep solves, which on this sequence has a closed form in n whose asymptote is 2n² — the marks are measured and the dashed line is that closed form, agreeing to 4.5·10⁻¹³. The lower marks are what a three-hundred-sweep fit from a random start returns: 2.1·10⁻⁴ at n = 2 rising to 0.0118 at n = 128. Started at the answer instead, the same code stays within 1.4·10⁻¹⁰ of it at every n — not the rounding level, because the drift from an exact start is itself about κ times the unit roundoff, but nine orders below what a random start reaches. That is the control that says the failure is the conditioning and not the implementation. A tensor away from the boundary conditions its step at 6.03.110¹10²10⁻¹²10⁻⁸10⁻⁴110⁴ncondition number, and residual reachedmarks above: κ of the step · dashes: its closed formmiddle: a fit from a random startbelow: the same fit started at the answera decomposition that is ill-conditionedκ at n = 1283.3·10⁴its closed form3.3·10⁴cosine of the terms1from a random start0.012from the answer1.4·10⁻¹⁰the answer existsand cannot be found
Fig. 11 The near-boundary case, from the essay that prices it: a condition number that is a square, and a fit that fails to find an answer that exists.
How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 12 And the identity that governs it, from the error field, which is the same identity here as everywhere else on this site.

The seed, and what it decides

Every run above starts from random factors, so the natural next question is how much of the answer is the starting point.

On a target with a rank-r answer, every seed reaches it and every seed reaches the same factors — which is the next essay’s subject and is the one direction in which a tensor is better behaved than a matrix. On a target without one, the seeds do not agree about anything except that the residual falls, because there is no answer for them to agree about.

That is a cleaner separation than the randomised field’s, where the seed moves an answer that exists. Here it decides which member of a non-convergent family the run happens to be walking along, and the spread between seeds is not an error bar on an estimate — it is the absence of a thing to estimate.

How much two runs of the same fit agree about the factors, from 6 starting points eachEvery run here reaches its target to the rounding level — 3·10⁻¹³, 1.16·10⁻¹¹ and 2.39·10⁻¹⁵ at worst — so all three are successful factorisations. The bar is the worst agreement between any two of them about the *factors*, matched over permutations and scalings, which is exactly the freedom the uniqueness theorem allows. Kruskal's condition k_A + k_B + k_C ≥ 2r + 2 holds for the first (9 ≥ 8) and fails for the second (6 < 8), and the bars are 1.0000 and 0.0207. The matrix is the comparison the whole thing rests on: AB = (AM)(M⁻¹B) for every invertible M, so its runs agree to 0.1089 and its factors mean nothing on their own.6 × 6 × 6, rank 3 · 9 ≥ 81.00002 × 2 × 2, rank 3 · 6 < 80.02076 × 6 matrix, rank 3 · no condition0.1089worst agreement between two runs about the factors, 0 to 1every run fits to 3·10⁻¹³every run fits to 1.2·10⁻¹¹every run factorises to 2.4·10⁻¹⁵ and none agrees with anotherthe only thing that improvestensor, Kruskal holds1tensor, Kruskal fails0.021matrix0.11worst residual1.2·10⁻¹¹three successful fitsone recoverable answer
Fig. 13 The well-posed case, from the next essay: six starting points, six fits, one answer.
The randomised SVD against the optimum it cannot beat, with 1 power iterationA semi-logarithmic plot of approximation error against target rank. A shaded band shows the spread across seeds, a solid line the optimal error from the exact singular values, and a dashed line the published probabilistic bound well above both.04812162010⁻¹10⁻⁰.⁵1target rank k‖A − A_k‖₂published boundrandomisedσ_{k+1}, optimalhow far apart the three areworst seed spread1.2bound / median at k = 1211median / optimum at k = 12160×60, 6 seeds, oversampling p = 5band is best to worst
Fig. 14 And the randomised field’s version of a seed moving an answer that exists, for contrast.

Where the model is still the right one

None of this is an argument against the model, and it is worth saying so plainly because the previous sections read as one.

The reason to fit a sum of rank-one terms rather than a subspace format is that the terms mean something. A chemical spectrum measured at several excitation wavelengths and several times is a mixture of components, and a rank-r fit is supposed to return the components. A subspace format returns a basis for the space they span, which is not the same information and is not what was asked.

And the model delivers that, under a condition that is checkable and generically true — which is the next essay. What this page establishes is the other side of the same coin: the delivery is conditional, the condition can fail, and when it fails the failure is silent under every stopping rule anyone normally uses.

The pair is the honest description. A CP fit is the only decomposition in this field whose output is interpretable, and the only one that can return an artefact while reporting success.

The recovered signal at three λ, 0.10% noiseThree recovered signals against the truth. At the best λ the two bumps are recovered and the step is rounded off. Fifty times smaller and the answer is dominated by amplified noise; fifty times larger and the step has gone entirely. The step is the feature that costs the most to represent and it is the first thing any filter loses.0816243240485600.511.5positionvaluethe truthbest λtoo littletoo mucheither side of the bestbest λ, relative error0.11fifty times smaller0.73fifty times larger0.15the step is the first thing every filter losesand it does not come back
Fig. 15 The regularisation field’s version of the same pairing, from the essay that measures it: an answer that is a choice, and a measurement of what the choice cost.

What to do instead

Three answers, and this collection has essays about two of them.

Fit a different rank. The target has rank three; a rank-three fit reaches 10⁻¹⁵ in fifty sweeps and its terms do not move. If the rank is what is being estimated, a swamp is evidence that the rank is wrong, and the term size is the evidence.

Fit a different model. The multilinear rank and the train ranks are ranks of matrices; their sets are closed; and a projection onto them exists, is computable, and is quasi-optimal. Three essays of this field are about that trade, and what it gives up is the interpretability of the terms.

Constrain the factors. Requiring non-negativity, or orthogonality in one mode, or bounded norms, makes the feasible set closed and removes the phenomenon outright. What it does not do is make the original question well posed; it answers a different one, and the answer is a projection onto a set the data was not asked about.

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. 16 The second answer, priced: a projection with a bound either side of it, on a set that contains its own limit points.
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. 17 And its long-chain version, whose ranks are ordinary matrix ranks at every cut.

The refusal

The claim under test is the one every stopping rule implicitly makes: that a plateau is slow convergence, and the remedy is patience.

The assertion that the largest term stays within half again of its value at sweep ten is fed the twenty-thousand-sweep trace of the border-rank run. It fails, at a growth of 2.16.

The refusal is worth having in this direction rather than the other. It would be easy to assert that the terms do grow, and that assertion would pass on a run that had diverged for any reason at all — including a bug. What is asserted instead is the flattering reading, on the run where it is false, so what the gate is protecting is the diagnostic rather than the phenomenon.

The file’s other two refusals guard the neighbours. One is fed the 2 × 2 × 2 rank-three case and required to refuse the claim that every CP decomposition is unique, which is the next essay’s boundary. The other is fed a benign run and required to refuse the claim that alternating least squares can go uphill, which is what keeps the monotonicity measurement from being read as a claim about the method’s convergence.

At other settings

Alternating least squares on a tensor with a rank-two answer and on one without, over 40,000 sweepsThe steeply falling curve is a rank-three fit to a tensor built from rank-three factors: it reaches 9.74·10⁻¹⁵ in 58 sweeps and its largest term does not move — 12.202 at sweep twenty and 12.202 at the end. The other is a rank-two fit to the border-rank tensor, whose error goes from 0.00677 to 0.00101 and has not finished, while the rising curve is its largest term: 4.698 to 11.99, a factor of 2.55. Fitted over 199 points, the error falls as the 2.02 power of the term size — a better path than the explicit sequence's 1/n, and still one that leaves the set rather than converging inside it.110¹10²10³10⁴10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹sweeprelative error, and the largest term's sizerising: the swamp's largest rank-one termfalling, slowly: its errorfalling, once: a fit with an answera plateau with a rising floorswamp error0.001swamp term12term growth2.6benign error9.7·10⁻¹⁵benign term growth1the error alone cannot tellthe size of the terms can
Fig. 18 Forty thousand sweeps, twice the hero, for reading whether anything turns. Nothing does.
A rank-two sequence approaching a rank-three tensor: the residual falls like 1/n and the terms grow like nA_n = n(e₁ + e₂/n)⊗³ − n·e₁⊗³ has rank two for every n and converges to a tensor of rank three. The falling curve is ‖A_n − A‖, which is √(3/n² + 1/n⁴) exactly and reaches 2.64·10⁻⁵ at n = 65536; the rising one is the norm of the larger of its two rank-one terms, 65540. Their product runs 2.519, 1.917, 1.777 … 1.73205, descending onto √3 = 1.73205. So getting one digit closer costs a factor of ten in the size of the pieces, for ever, and the infimum of the distance is zero while no rank-two tensor attains it.110¹10²10³10⁴10⁻⁴10⁻³10⁻²10⁻¹110¹10²10³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 attainedn6.6·10⁴‖A_n − A‖2.6·10⁻⁵largest term6.6·10⁴their product1.7√31.7the distance goes to zeroand nothing reaches it
Fig. 19 The explicit sequence extended as far as it is drawn, which is the path this iteration is not taking.
A plane of 2 × 2 × 2 tensors, coloured by the sign of the hyperdeterminant that decides their rankTwo of the eight entries are varied over ±0.5 and the other six are held at the values that make the centre of the picture the rank-three tensor A = e₁⊗e₁⊗e₂ + e₁⊗e₂⊗e₁ + e₂⊗e₁⊗e₁. The light region is Δ > 0, where the tensor has real rank two and two real rank-one terms exist; the dark region is Δ < 0, where the real rank is three and the complex rank is two. 50.0 per cent of this window is rank two. The boundary between them is the curve Δ = 0, on which the rank is three and the distance to the rank-two set is zero — every point of it is a tensor with no nearest rank-two approximation.rank three, Δ = 0a₁₁₁ across, a₀₀₀ up, both over ±0.5light: two real roots, rank two · dark: a conjugate pair, rank threeone sign, two rankscells sampled4096rank two2048rank three2048Δ at the centre0rank at the centre3the rank is a signcomputed from eight numbers
Fig. 20 The boundary at four times the magnification, where it is still a boundary.
The share of random 2 × 2 × 2 tensors with real rank two, against the number drawn, and π/4Each draw is eight independent standard normal entries and its rank is decided exactly, by the sign of the hyperdeterminant, with no iteration involved. Over 19,953 draws 15,705 have rank two — a share of 0.7871 against the exact value π/4 = 0.7854, inside 0.6 standard errors. The band is ±2 of them. A random matrix has one typical rank; this is the picture of a random object having two, each with a probability that is a number rather than an experiment.10¹10²10³10⁴00.050.10.150.20.250.30.350.40.450.50.550.60.650.70.750.80.850.90.951drawsshare with real rank twoπ/4 = 0.78539815,705 of 19,953 have rank twoa probability with a closed formdraws2·10⁴rank two1.6·10⁴share0.79π/40.79standard errors out0.59two typical ranksand the split is π/4
Fig. 21 How often a random target of this shape has no rank-two answer at all: 21.5 per cent of the time.
Along a sequence of exactly rank-two tensors: the step's condition number, and what a fit from a random start achievesEvery A_n on this sequence *is* a rank-two tensor and its two rank-one terms are written down, so nothing here is about existence. The rising curve is the condition number of the r × r system each alternating sweep solves, which on this sequence has a closed form in n whose asymptote is 2n² — the marks are measured and the dashed line is that closed form, agreeing to 1.1·10⁻¹¹. The lower marks are what a three-hundred-sweep fit from a random start returns: 2.1·10⁻⁴ at n = 2 rising to 0.0118 at n = 1024. Started at the answer instead, the same code stays within 9.2·10⁻⁹ of it at every n — not the rounding level, because the drift from an exact start is itself about κ times the unit roundoff, but nine orders below what a random start reaches. That is the control that says the failure is the conditioning and not the implementation. A tensor away from the boundary conditions its step at 6.03.110¹10²10³10⁻¹²10⁻⁸10⁻⁴110⁴ncondition number, and residual reachedmarks above: κ of the step · dashes: its closed formmiddle: a fit from a random startbelow: the same fit started at the answera decomposition that is ill-conditionedκ at n = 10242.1·10⁶its closed form2.1·10⁶cosine of the terms1from a random start0.012from the answer9.2·10⁻⁹the answer existsand cannot be found
Fig. 22 The conditioning followed a decade further along the sequence.
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. 23 The other half of the trade, from two essays back: what a diagonal core costs, which is what this model insists on.
How far the coefficients can move without changing the fit, degree 11Relative increase in the residual against relative change in the coefficients, along the least determined direction. The residual does not move measurably until the coefficients have changed by more than a factor of one.10⁻⁴10⁻³10⁻²10⁻¹110¹10²10³10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹relative change in the coefficients, along the worst directionrelative increase in the residualcoefficients doubled1544% change, fit unmoved in the sixth digit308×: the third digit movesκ(A) = 1.4·10⁸. Exact arithmetic would pick one point on this floor. It would not raise it.24 points, degree 11, monomial basisthe data leaves them free
Fig. 24 The least-squares field’s picture of a fit with directions that cost nothing, for comparison — a valley whose floor is flat, on a problem whose set has a dimension.
Adding 0.1 to 256, a thousand times, at 8 significand bitsThree lines against the number of additions: the exact total rising steadily, a staircase following it, and a flat line at the starting value.02505007501000250275300325350375additionsrunning totalround to nearest: nothing arrivesexactstochasticnearesta thousand additionshalf an ulp at 2561moves, round to nearest0moves, stochastic46relative error, nearest0.28relative error, stochastic0.0228 significand bits, unbounded exponenta flat line is not a small error
Fig. 25 And the arithmetic field’s version of an iteration that stops moving without arriving, which is a different mechanism with the same trace.
A single Newton step solved to eleven inner tolerances, against how far the resulting point is from the rootThe iterate is 1.03 from a root that is known exactly by construction. The same linear system is solved to relative residuals from 0.3 down to 10⁻¹⁴, at 53 and 849 conjugate gradient iterations, and the resulting point is 0.5725 and 0.5071 from the root. The curve is flat below about 10⁻³: the linearisation is wrong at second order, so the step cannot land closer than the square of the distance it started at — 1.055 — however exactly it is computed.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110⁻¹110¹inner tolerance η, relative residual of the linear solvedistance from the root after the stepd² = 1.05distance before the step, 1.0353219543849the number by each point is the iterations it costeleven decades, one landing placedistance before the step1its square1.1where η = 10⁻³ lands0.51where η = 10⁻¹⁴ lands0.51iterations for the first219iterations for the second849the accuracy that is thrown awaymeasured against a root that is known
Fig. 26 The sequence field’s plateau, for reading against this one: a floor that is a linearisation rather than a boundary.

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. 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⁻⁴ The road that squares the problem The normal equations are the first method every course teaches and the method no library uses. Forming AᵀA squares the condition number, and below ε = √u it does not degrade — it produces a matrix that is exactly singular, from data that was perfectly usable. 0481216202428321rank keptrelative errormedianthe answer movesspread at rank 81.8spread at rank 241best median error0.14widest where the method is worstand the bound does not say so An answer that changes with the seed A randomised rank-k solve is a truncation computed in a random subspace, and it reaches the same floor as the deterministic ones. What it does not do is return the same answer twice — a factor of 1.84 across four seeds at rank 8, and 1.02 at the rank where the method is best.

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 squaresBorder rankCondition squaringCP decompositionIll posed problemKhatri–Rao productLow-rank approximationNormal equationsSwampTensor rank