When the matrix is wrong too

The projection and the right angle set up least squares the way everyone does. min ‖Ax − b‖ finds the point in the column space of A closest to b, the residual is orthogonal to that space, and the picture is a right triangle.

There is an assumption in it that is never stated, and it is visible in the picture: b moves and A does not. The whole of the correction is applied to the right-hand side. The columns of A are drawn as a fixed plane.

That is the right model when A holds the design of an experiment — the times at which measurements were taken, the settings of a dial, a basis evaluated at known points — and b holds the readings.

It is the wrong model whenever A holds measurements too.

What total least squares asks for #

$$\min \bigl\| [\, \Delta A \;\; \Delta b \,] \bigr\|_F \quad \text{subject to} \quad (A + \Delta A)x = b + \Delta b$$

The smallest perturbation of the whole augmented matrix that makes the system consistent, rather than the smallest perturbation of b alone.

Its solution comes out of one SVD, and the derivation is three sentences. Consistency says [A + ΔA, b + Δb]·[x; −1] = 0, so the perturbed augmented matrix must have a null vector. The smallest Frobenius perturbation that creates one is the one that zeroes σₙ₊₁, by Eckart–Young — which this site checked as an equality in the best approximation there is. So the null vector is v, the last right singular vector of [A b], and

x  =  −v(1:n) / vₙ₊₁

That is the whole algorithm. One SVD of an m×(n+1) matrix, and a division.

The model, not the method #

The hero figure sweeps the location of a fixed amount of noise, which is the comparison worth making: both methods see problems of the same difficulty, and only where the error sits changes.

With all of it in b, total least squares is worse — it corrects a matrix that did not need correcting, and the correction is error. assertOrdinaryLeastSquaresWinsWhenOnlyBIsNoisy measures that at three noise levels and finds ordinary least squares ahead at every one.

With all of it in A, total least squares is ahead by a factor of 2.5, and the crossing is between three tenths and six tenths of the way across.

Neither is the better method. Each is the maximum-likelihood estimator under a different noise model, and choosing between them is choosing a model. A reader who takes away “total least squares is the more sophisticated one” has taken away the wrong thing, and the left-hand end of that curve is there to prevent it.

What ordinary least squares is doing wrong when A is noisy #

The failure has a name and a direction: attenuation bias.

If A is the true design and the observed one is A + E, then (A+E)ᵀ(A+E) = AᵀA + EᵀE + cross terms, and EᵀE has expectation mσ²I — a positive definite matrix added to the normal-equations matrix before it is inverted. The estimate is therefore shrunk towards zero, systematically, and the obvious factor to write down is 1/(1 + mσ²/σₘᵢₙ(A)²) — which is a bound rather than the shrinkage, for a reason measured two sections below.

It is a bias, not a variance: more data does not remove it, because EᵀE grows with m exactly as AᵀA does. That is the one thing separating this from every other error on this site, where more samples or more precision eventually help.

Total least squares removes the bias — the SVD of the augmented matrix is exactly the total-least- squares estimator, and it is consistent under the errors-in-variables model — and pays for it in variance. Hence the crossing.

Where the bias becomes visible #

The figure’s own assertions change at a noise level of 0.08, and the two frames either side of that line are worth putting next to each other, because the shape of the curve changes rather than its size.

Both frames report a crossing, and only the second one has found anything. The figure at 0.07 puts it “between 0 and 0.3” because the curve happens to pass 1.0 there on its way up; at 0.05 the same detector would put it in the same place, on a curve whose highest and lowest points differ by 27%. This is the reason the generator asserts in two halves rather than one — an assertion of separation that held at every stop would be asserting something false at the bottom of the slider, and the site had already written the sentence for it: an assertion that must hold at every drag position cannot be stated unconditionally when the slider crosses a threshold.

And above it the advantage is quadratic in the noise #

Once the sweep separates, the size of the separation follows the mechanism closely enough to check. The bias is second order in the noise and the shared variance is first, so the excess advantage at the right-hand end — the advantage minus one, which is what the bias removal actually buys — should go as the square.

Divide the excess by the square of the noise and the three stops agree:

Constant to within 6% across the first three rows, and then it is not. The quadratic law is the mechanism showing through, and the last two rows are where “second order” stops being an adequate description of a noise level that is no longer small — 0.4 is a fortieth of the way to the assertion’s own ceiling of 0.5, and the estimator is being asked to separate signal from an error of comparable size. A reader who extrapolates the first three rows to a noisy problem will over-promise by a factor of two.

The penalty is larger than the prize #

One more reading of the same five curves, and it is the one that should decide a default.

Total least squares used where the matrix really is exact costs a factor of 1/advantage at share 0. Ordinary least squares used where all the noise is in the matrix costs the advantage at share 1. Those are the two ways of being wrong, and they are not the same size:

Choosing the wrong model in the total-least-squares direction is worse than choosing it in the ordinary direction, at every noise level on the slider, and the gap widens with the noise. At 0.4 it is three to one.

That is not an argument that ordinary least squares is the better method — the section above said it is not, and it is not. It is an argument about which mistake to make under uncertainty, and it points the opposite way to the intuition that the more elaborate estimator is the safer default. The asymmetry is in the geometry: correcting a matrix that needed no correcting moves the answer by the full size of the correction, while failing to correct one that did leaves a bias that is second order and therefore smaller than the noise that caused it.

The shrinkage is by direction, and σₘᵢₙ is the worst of them #

That factor is worth evaluating rather than quoting, because a caller reading it will use it to decide whether the bias is worth removing. Against the measured projection of the ordinary estimate onto the true coefficients, medians over forty seeds, all the noise in the matrix:

The σ_min formula overstates the bias by a factor of two to three — 0.437 of shrink predicted against 0.200 measured at σ_A = 0.35 — and the factor is not a constant to divide out: it runs from 3.0 at the smallest noise to 2.1 at the largest.

The reason is that σ_min is the worst direction and the coefficients are not in it. Every direction shrinks by its own σᵢ²/(σᵢ² + mσ²), and the estimate’s shrinkage is the average of those weighted by how much of the true coefficient vector lies along each. Resolved that way the prediction lands within one per cent of the measurement at all sixteen cells of the grid. So the closed form is there and it is a sum rather than a single ratio, and the single ratio is what it collapses to when the coefficients happen to sit in the worst direction.

That is the same shape as several other bounds this collection has taken apart, and it has the same practical consequence: read as a value, the σ_min formula tells a caller their coefficients are shrunk by 44 per cent when the answer is 20, and the decision to switch estimators is being made on a number twice too large.

And the bias really does survive the data #

The claim that separates this failure from every other one on the site — more data does not remove it — is worth having as a measurement, because it is the reason a total-least-squares estimator exists at all. Over the same grid, at σ_A = 0.2:

Seventeen times the data, and the ordinary estimate’s error does not move — 1.10·10⁻¹ to 1.30·10⁻¹, with no trend in it. The total-least-squares error falls by a factor of 3.2 over the same range, which is √17 = 4.1 to the accuracy a median over forty seeds supports. One estimator is consistent and the other is not, and this is the picture of it: two curves against the number of observations, one of them flat.

Which makes the choice sharper than the hero figure alone does. At a fixed m the two methods trade a bias against a variance and the crossing is where the trade balances. Across m they do not trade at all — the variance goes away and the bias does not — so at any noise level with a matrix in it there is an m past which total least squares wins and stays winning. On this family, at σ_A = 0.2, that m is under sixty.

assertTheAttenuationIsByDirectionNotBySigmaMin holds all of it: the direction-resolved prediction to two per cent at every cell, the σ_min rule’s overstatement and its variation, the flatness of the ordinary estimate’s error in m, and the improvement of the total one.

The other end of the sweep, and what the residual says there #

Everything above compares the two answers by their distance from coefficients only a constructed problem has. A caller has no such thing, and the quantity a caller does have is about to order the two methods backwards — so it is worth looking first at the one case where it does not.

The residual orders them backwards #

Here is the part this site exists for.

‖Ax − b‖ is the quantity ordinary least squares minimises. No other vector can make it smaller — not the total-least-squares answer, and not the coefficients the problem was actually built from.

So wherever total least squares is more accurate, the residual prefers the other answer. Not sometimes. By construction.

assertTheBetterAnswerHasTheLargerResidual states it as a strict inequality on the residual ratio at every share, and then checks that at the shares where total least squares is more accurate, the two orderings disagree.

And it pushes the same statement one step further: the exact coefficients the problem was built from have a larger residual than the least-squares ones. That is not a defect of anything; it is what “minimises” means. It is also the sharpest possible form of a small residual is not a small error: here the residual is not merely uninformative about the error, it is anti-correlated with it, over a whole family of problems, provably.

What makes a total least-squares problem hard #

Not κ(A), which amplifies a perturbation and says nothing here. The conditioning is governed by the gap between σₙ(A) and σₙ₊₁([A b]), and those two numbers are interlaced — the augmented matrix has one more column, so its singular values interlace A’s from below.

As the noise rises, σₙ₊₁([A b]) climbs towards σₙ(A), the gap closes, and the estimator blows up. κ(A) does not notice.

assertTheTLSConditioningIsTheGapNotKappa asserts that the gap is monotone, that the error is monotone with it, that κ moves by less than a factor of eight, and that the error grows by more than twenty times whatever κ did.

The drag is over the number of observations rather than over the noise, and that is deliberate: more rows open the gap — σₙ(A) grows like √m while the noise the augmented matrix picks up does not keep pace — while κ, being a ratio, is nearly scale-free in m. Two quantities move together and the third does not move at all, which is the same statement the sweep makes along its own axis, reached by turning a different knob.

The geometry, and why “orthogonal regression” is the wrong summary #

The picture usually drawn for total least squares is a line fitted to a scatter of points by minimising perpendicular distances rather than vertical ones, and it is a correct picture of one special case that generalises badly.

It is correct when there is one predictor, both coordinates carry the same noise, and the model has an intercept — a straight line through a cloud. Then the total-least-squares line does minimise the sum of squared perpendicular distances, and the ordinary one minimises vertical distances, and the two-line summary is honest.

It stops being correct as soon as the noise levels differ between columns, and the reason is worth having: the objective is ‖[ΔA Δb]‖_F, a single Frobenius norm over the whole augmented matrix, which treats a unit of error in the first column as interchangeable with a unit in the last. If the first column is a temperature in kelvin and the last a pressure in pascals, that is a claim nobody made.

The repair is a scaling — divide each column by its own noise level and solve the scaled problem — and once that is said, this essay has met the units the matrix is measured in from the other side. There the scaling was a nuisance that moved a condition number and not the answer. Here it moves the answer, because the objective is a norm over columns that the scaling reweights, and it has to be chosen rather than removed.

The special case where all columns are equally noisy is exactly the case where no scaling is needed, which is why the perpendicular-distance picture is the one that survives into textbooks. It is also the case where the intercept column of ones is exactly known — noise-free — and therefore should not be scaled with the rest, which the picture quietly ignores.

What it costs, and the version that is affordable #

An SVD of an m×(n+1) matrix, against a QR of an m×n one for ordinary least squares. On a tall thin problem — m in the thousands, n in the tens, which is the shape most regressions have — that is a factor of two or three, and it is affordable.

On a large sparse problem it is not. There is no sparse total least squares in the sense there is a sparse least squares: the estimator needs the smallest singular value of the augmented matrix, and that is the hardest end of the spectrum to reach iteratively. What is used instead is a Rayleigh-quotient iteration on the augmented matrix, or a regularised variant, and both are solving a different problem than the one the SVD solves exactly.

That gap is worth noting because it is where the choice usually gets made in practice: total least squares is used on small problems where the SVD is free and not used on large ones where it is not, and the decision is made by the cost rather than by the noise model. The correct reason and the operative reason are different, which is a thing worth knowing about a method.

When there is no answer at all #

If vₙ₊₁ = 0 — the last right singular vector has no component along b — then no finite x satisfies the constraint and the total least-squares problem has no solution. It happens exactly when σₘᵢₙ(A) ≤ σₙ₊₁([A b]), which the interlacing makes possible only in the degenerate case.

tlsSolve reports it rather than dividing. A routine that computed −v(1:n)/vₙ₊₁ unconditionally would return a vector of enormous numbers, and a figure would draw it.

The constructed case is a rank-deficient A: its own smallest singular value is zero, so it can never exceed the augmented matrix’s, and no solution exists. Ordinary least squares on the same problem returns a perfectly finite answer, which is the point — the two methods do not fail on the same problems, and a code that switches between them has to handle that.

The regularised version, and why it exists #

Total least squares is less stable than ordinary least squares, not more, and the gap measurement is why: the estimator’s conditioning is a difference of two singular values, and a difference of two nearly equal quantities is the oldest hazard on this site.

The standard repair is regularised total least squares — a constraint ‖x‖ ≤ δ, or a Tikhonov penalty applied inside the augmented problem — which bounds the estimator away from the collision at the cost of a parameter nobody knows. That places it in the same position as every other rule choosing without knowing measured: a knob that has to be set from the data, by a heuristic, scored against a truth that exists only because the problem was constructed.

So the honest summary of the method is a pair. It removes a bias that no amount of data removes, and it introduces a sensitivity that no amount of data removes either, and which of those dominates is a property of the gap rather than of the noise.

Where the model is decided, and by whom #

A last observation, because the essay’s conclusion is a choice and choices have owners.

Nothing in the numbers says which model is right. The crossing in the hero figure is a fact about where the noise is, and where the noise is is a fact about the experiment rather than about the matrix. A solver cannot see it: A and b arrive as arrays, and no property of those arrays reports which of their entries were measured and which were set.

So the decision belongs to whoever knows how the data was collected, and it is made — implicitly and almost always — by calling the routine whose name has no adjective in front of it.

What is worth carrying #

Ordinary least squares assumes the matrix is exact, and that assumption is in the geometry of the picture everybody draws rather than in any sentence.

When it is wrong the error is a bias rather than a variance, so more data does not remove it — the one place on this site where that is true.

The more accurate answer has the larger residual, by construction, because the other method minimises the residual by definition. Where two methods disagree, the quantity a solver can see prefers the one that optimised it, and that is not evidence.

And the conditioning is a gap between two singular values of two different matrices. κ(A) is computable, familiar, and about a different question.