Numerical Linear Algebra
where the algebra stops being the arithmetic
Least squares, and the road not to take

When the matrix is wrong too

Every least-squares problem on this site has assumed A is exact and b is not, and moved b onto the column space of A. Where both were measured, the smallest correction that makes the system consistent moves the matrix as well — and on the problems where that answer is more accurate, it has the larger residual, by construction rather than by luck.
21 min read 13 figures Whose fault is itMeasured, not assertedTwo routes to a numberAssertions that reject

Worth reading first: The projection and the right angle · The best approximation there is.

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.

Which of the two least-squares methods is more accurate, as a fixed amount of noise moves from b into AOne curve: the ordinary least-squares error over the total least-squares error, at each share of the noise placed in the matrix, as the median over 40 seeds. It runs from 0.28 when all the noise is in b — where ordinary least squares is the more accurate — to 2.47 when all of it is in A. The line at one is where the two methods are equally right.00.250.50.75110⁻¹110¹share of the noise placed in the matrixleast-squares error ÷ total least-squares errorequally accuratetotal leastsquares aheadordinary leastsquares aheadthe model, not the methodadvantage, all noise in b0.28advantage, all noise in A2.5seeds at each share40the same total noise at every pointand only where it sits changes
Fig. 1 A fixed amount of noise moved from the right-hand side into the matrix, and which method is more accurate at each point. On the left, all the noise in b: ordinary least squares is three and a half times closer to the true coefficients. On the right, all of it in A: total least squares is 2.5 times closer. The line at one is where they are equally right. Drag the total noise.
The least-squares solution as a projection, with the right angle measuredThe column space drawn edge-on as a plane, the data vector above it, and the perpendicular dropped to the plane, with the residual marked at a right angle to it.everything Ax can reachb = (1.3, 0.5, 1.4)Ax, the closest reachable pointr = b − Ax‖Aᵀr‖ / (‖A‖‖r‖)3.2·10⁻¹⁶‖b‖² − ‖Ax‖² − ‖r‖²1.6·10⁻¹⁵‖r‖1.2200 random nearby points of the plane were tried; none is closer.a 3×2 system, Householder QRperpendicularity is checked
Fig. 2 The picture this site drew for ordinary least squares. The plane is fixed and b moves onto it, which is the assumption the whole of this essay is about.

What total least squares asks for 

min ‖[ΔA  Δb]‖_F   subject to   (A + ΔA)x = b + Δ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 σ_{n+1}, 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_{n+1}

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

Which of the two least-squares methods is more accurate, as a fixed amount of noise moves from b into AOne curve: the ordinary least-squares error over the total least-squares error, at each share of the noise placed in the matrix, as the median over 40 seeds. It runs from 0.17 when all the noise is in b — where ordinary least squares is the more accurate — to 3.14 when all of it is in A. The line at one is where the two methods are equally right.00.250.50.75110⁻²10⁻¹110¹share of the noise placed in the matrixleast-squares error ÷ total least-squares errorequally accuratetotal leastsquares aheadordinary leastsquares aheadthe model, not the methodadvantage, all noise in b0.17advantage, all noise in A3.1seeds at each share40the same total noise at every pointand only where it sits changes
Fig. 3 The same sweep at a higher noise level. The crossing moves left as the noise grows, because the bias ordinary least squares carries is second order in it and the variance both methods share is first.

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, by a factor of about 1/(1 + mσ²/σ_min(A)²).

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.

Which of the two least-squares methods is more accurate, as a fixed amount of noise moves from b into AOne curve: the ordinary least-squares error over the total least-squares error, at each share of the noise placed in the matrix, as the median over 40 seeds. It runs from 0.99 when all the noise is in b — where ordinary least squares is the more accurate — to 0.91 when all of it is in A. The line at one is where the two methods are equally right.00.250.50.75110⁻¹110¹share of the noise placed in the matrixleast-squares error ÷ total least-squares errorequally accuratetotal leastsquares aheadordinary leastsquares aheadthe model, not the methodadvantage, all noise in b0.99advantage, all noise in A0.91seeds at each share40the same total noise at every pointand only where it sits changes
Fig. 4 And at a low one, where the sweep is flat: neither method ahead by more than a few per cent at any share, and not monotone. There is nothing to remove, so removing it buys nothing.

And below a noise level there is nothing to remove

The drag has a stop at total noise 0.05, and the sweep there is flat: 0.99, 1.16, 1.14, 1.00, 0.91 — neither method ahead by more than a few per cent anywhere, and not monotone.

That is not a failure of the measurement. The bias is second order in the noise — it goes like σ² — while both methods carry the same first-order variance. So below a noise level the bias is under the variance and there is nothing to remove; above it the bias dominates and the sweep separates cleanly.

Asserting the separation at every stop would be asserting something false at the lowest one, so the figure asserts it in two halves, with the boundary stated. The site has been here before, and the sentence it wrote then was an assertion that must hold at every drag position cannot be stated unconditionally when the slider crosses a threshold.

Distance from the truth and residual, for both methods, with 0% of the noise in the matrixFour bars, medians over 40 seeds. The upper pair is how far each answer is from the coefficients the problem was built from; the lower pair is ‖Ax − b‖ on the problem as given. Ordinary least squares minimises the lower quantity by definition, so its bar is the shorter of the two whatever happens above — and at this share it is the less accurate answer.the upper pair is distance from the truth; the lower pair is ‖Ax − b‖least squares · error0.04439total least squares · error0.1565least squares · ‖Ax − b‖4.734total least squares · ‖Ax − b‖4.971two orderingserror ratio (ls ÷ tls)0.28residual ratio (tls ÷ ls)1.1seeds40no vector makes the residual smallernot even the one the problem was built from
Fig. 5 The same four bars with all the noise in b. Here the two orderings agree — ordinary least squares wins both — which is the case where checking the residual reaches the right conclusion for the wrong reason.

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.

Distance from the truth and residual, for both methods, with 100% of the noise in the matrixFour bars, medians over 40 seeds. The upper pair is how far each answer is from the coefficients the problem was built from; the lower pair is ‖Ax − b‖ on the problem as given. Ordinary least squares minimises the lower quantity by definition, so its bar is the shorter of the two whatever happens above — and at this share it is the less accurate answer.the upper pair is distance from the truth; the lower pair is ‖Ax − b‖least squares · error0.1098total least squares · error0.04436least squares · ‖Ax − b‖4.536total least squares · ‖Ax − b‖4.698two orderingserror ratio (ls ÷ tls)2.5residual ratio (tls ÷ ls)1seeds40no vector makes the residual smallernot even the one the problem was built from
Fig. 6 Four bars, medians over forty seeds, with all the noise in the matrix. The upper pair is distance from the true coefficients; the lower pair is ‖Ax − b‖. Total least squares wins the first comparison by 2.48× and loses the second by 3.6%. Drag the share of the noise in the matrix and watch only the upper pair change order.

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.

Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u60×60, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 7 The site’s standing version of that gap. This essay is the case where the two quantities do not merely fail to track each other — they run in opposite directions, by construction.
Distance from the truth and residual, for both methods, with 60% of the noise in the matrixFour bars, medians over 40 seeds. The upper pair is how far each answer is from the coefficients the problem was built from; the lower pair is ‖Ax − b‖ on the problem as given. Ordinary least squares minimises the lower quantity by definition, so its bar is the shorter of the two whatever happens above — and at this share it is the less accurate answer.the upper pair is distance from the truth; the lower pair is ‖Ax − b‖least squares · error0.05224total least squares · error0.03071least squares · ‖Ax − b‖3.398total least squares · ‖Ax − b‖3.467two orderingserror ratio (ls ÷ tls)1.7residual ratio (tls ÷ ls)1seeds40no vector makes the residual smallernot even the one the problem was built from
Fig. 8 And at 60% of the noise in the matrix, just past the crossing. The accuracy ordering has flipped and the residual ordering has not, because it cannot.

What makes a total least-squares problem hard

Not κ(A). The conditioning is governed by the gap between σ_n(A) and σ_{n+1}([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, σ_{n+1}([A b]) climbs towards σ_n(A), the gap closes, and the estimator blows up. κ(A) does not notice.

What governs a total least-squares problem: the gap σ_n(A) − σ_{n+1}([A b]), not κ(A)Three curves against the noise level, both axes logarithmic, as medians over 12 seeds of a 60×3 fit. The gap between the smallest singular value of A and the smallest of the augmented matrix closes from 2.1 to 0.108, and the total least-squares error rises with it from 0.0067 to 1.69 — a factor of 252. κ(A) moves from 3.66 to 2.48 across the same sweep and predicts none of it.10⁻¹10⁻³10⁻²10⁻¹110¹noise level, relative to the datagap, error, and condition numberκ(A)TLS errorthe gapone of these predicts the errorgap, at the least noise2.1gap, at the most0.11error, ratio across the sweep252κ(A), ratio across the sweep1.5the condition number is nearly flatand the error moves by two orders
Fig. 9 Three curves against the noise level. The gap closes by a factor of nineteen and the total least- squares error rises by 252; κ(A) moves by a factor of 1.5 across the same sweep and predicts none of it. Drag the number of observations, which opens the gap and changes nothing about κ.

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 — σ_n(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.

What governs a total least-squares problem: the gap σ_n(A) − σ_{n+1}([A b]), not κ(A)Three curves against the noise level, both axes logarithmic, as medians over 12 seeds of a 240×3 fit. The gap between the smallest singular value of A and the smallest of the augmented matrix closes from 4.1 to 0.373, and the total least-squares error rises with it from 0.00434 to 1.48 — a factor of 342. κ(A) moves from 3.73 to 2.45 across the same sweep and predicts none of it.10⁻¹10⁻³10⁻²10⁻¹110¹noise level, relative to the datagap, error, and condition numberκ(A)TLS errorthe gapone of these predicts the errorgap, at the least noise4.1gap, at the most0.37error, ratio across the sweep342κ(A), ratio across the sweep1.5the condition number is nearly flatand the error moves by two orders
Fig. 10 The gap sweep with four times the observations. More rows open the gap — σ_n(A) grows like √m — and the error at every noise level falls with it, while κ(A) does not move.

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.

Least squares by QR and by the normal equations in 32-bitRelative 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 = 1.5·10⁻⁵
Fig. 11 The site’s other essay about a least-squares route chosen for its cost. There the cheap route squares a condition number; here the expensive route is skipped and the model changes. The two decisions are made in the same place by the same person.
What governs a total least-squares problem: the gap σ_n(A) − σ_{n+1}([A b]), not κ(A)Three curves against the noise level, both axes logarithmic, as medians over 12 seeds of a 24×3 fit. The gap between the smallest singular value of A and the smallest of the augmented matrix closes from 1.4 to 0.239, and the total least-squares error rises with it from 0.0117 to 0.722 — a factor of 62. κ(A) moves from 3.54 to 2.46 across the same sweep and predicts none of it.10⁻¹10⁻³10⁻²10⁻¹110¹noise level, relative to the datagap, error, and condition numberκ(A)TLS errorthe gapone of these predicts the errorgap, at the least noise1.4gap, at the most0.24error, ratio across the sweep62κ(A), ratio across the sweep1.4the condition number is nearly flatand the error moves by two orders
Fig. 12 And with a quarter of them, where the gap starts closer and the collision arrives at a lower noise level. Two knobs move the same pair of quantities together and neither moves κ.

When there is no answer at all

If v_{n+1} = 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 σ_min(A) ≤ σ_{n+1}([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_{n+1} 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 diagonal of a pivoted R against the singular values, on a matrix of planted rank 4Two curves down the index, lying on top of each other. Both fall by more than five orders of magnitude between index 4 and index 5, so both report rank 4 at any threshold placed in the gap.1611162126313610⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹indexmagnitude|r_kk|σ_kone factorisation, two verdicts‖AP − QR‖/‖A‖1.1·10⁻¹⁵|r_nn|10⁻¹²σ_min10·10⁻¹³column interchanges37|r_nn| is never below σ_minso the cheap verdict errs one way only
Fig. 13 The other thing an SVD of the augmented matrix is asked for. A rank verdict and a total least-squares fit come from the same factorisation, and the cheap substitute for the first has nothing to offer the second.

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.

Where this goes next

everything Ax can reachb = (1.1, 0.4, 1.5)Ax, the closest reachable pointr = b − Ax‖Aᵀr‖ / (‖A‖‖r‖)1.7·10⁻¹⁶‖b‖² − ‖Ax‖² − ‖r‖²1.3·10⁻¹⁵‖r‖1.3200 random nearby points of the plane were tried; none is closer.a 3×2 system, Householder QRperpendicularity is checked The projection and the right angle The least-squares solution is the one whose residual is perpendicular to everything the columns can reach. That is not a mnemonic — it is an equation, Aᵀr = 0, and the computed answer satisfies it to 10⁻¹⁶. 110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method A small residual is not a small error Substituting the answer back and finding that it fits is the most natural check there is, and it verifies the wrong thing. A residual of 10⁻¹⁷ is entirely compatible with an answer whose second digit is wrong. 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.

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.

Condition numberEckart–YoungErrors-in-variablesForward errorLeast squaresNormal equationsOrthogonal projectionResidualSingular value decompositionTotal least squares