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

The valley with no bottom

A degree-nine fit's coefficients can be moved by a third of their own size before the residual changes in the sixth significant figure. The arithmetic did not lose those digits. The data never contained them.
15 min read 6 figures Where the algebra stops being the arithmeticWhose fault is itMeasured, not asserted

Fit a degree-nine polynomial to twenty-four points of a smooth function. The fit is excellent — the residual is at the level of rounding, the curve passes through the data, and any plot of it looks perfect.

Now take the ten coefficients and move them along a particular direction, all together, by 39% of their own size. Plot the new curve.

It is the same curve. The residual has not changed in its sixth significant figure.

How far the coefficients can move without changing the fit, degree 9Relative 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 doubled39% change, fit unmoved in the sixth digit308×: the third digit movesκ(A) = 3.6·10⁶. Exact arithmetic would pick one point on this floor. It would not raise it.24 points, degree 9, monomial basisthe data leaves them free
Fig. 1 The increase in the residual against the size of the change in the coefficients, along the least determined direction — the right singular vector belonging to the smallest singular value. The coefficients can be moved by a third of themselves before the fit degrades in the sixth digit, and by orders of magnitude more before it moves in the third.

What this says about the arithmetic

Nothing. That is the point.

The usual reading of an ill-conditioned fit is that floating point has destroyed the coefficients. It has not. The coefficients were never determined by the data to the precision one might hope for, and an infinitely precise computer would face the same situation: it would return one point on that valley floor, exactly, and a different point would fit the data just as well.

The distinction is the one this whole site is arranged around, and this is the place where it is most often got wrong, because the symptom — wildly varying coefficients from nearly identical data — looks exactly like a numerical failure.

Test it directly. Solve the same fit twice, in double and in single precision. The two sets of coefficients differ enormously. Now solve it twice in double precision, with the data perturbed in its twelfth digit. The two sets of coefficients differ enormously in the same way. The second experiment has no arithmetic difference in it at all.

The two questions are different questions

A least-squares problem contains two questions and they have different answers.

Where is the projection? That is well conditioned. The point Ax̂ in the column space closest to b is determined by b and the space, and it moves by no more than the perturbation. This is the question the projection and the right angle answers, and it is computable to the last bit.

Which combination of columns produced it? That is the coefficients, and when the columns are nearly dependent it is nearly undetermined. Many very different x give nearly the same Ax.

The conditioning of the second question is κ(A), and κ(A) is large exactly when the columns are nearly dependent. The flatness of the valley and the largeness of the condition number are the same fact described from two directions: a small singular value means a direction in coefficient space that A barely moves, so travelling along it changes x a great deal and Ax hardly at all.

Where the direction comes from

The direction in the figure is not chosen by search. It is the right singular vector vₙ belonging to the smallest singular value σₙ, and it has the defining property that ‖Avₙ‖ = σₙ, which is the smallest that ‖Av‖ gets for any unit v.

So moving the coefficients by tvₙ changes the fitted values by tσₙ, and the residual by no more than that. With σₙ small, the residual barely notices. This is not an empirical finding about a particular matrix; it is what the singular value decomposition is for, and the figure is that statement plotted.

The same object appears in rank is a decision as the direction that becomes unresolvable when the noise reaches it, and in the condition number is an amplifier as the direction in which a perturbation is amplified worst. Small singular values, undetermined coefficients, amplified errors and near rank deficiency are four names for one thing.

The basis is part of the problem

Here is the most actionable fact in this field, and it is not about arithmetic.

“Fit a degree-eleven polynomial” sounds like one request. It is two, because the matrix depends entirely on which polynomials the coefficients multiply. The monomials 1, x, x², … on [0,1] give a design matrix with κ = 1.2·10⁸. The Chebyshev polynomials on [−1,1], on the same data, give κ = 2.5.

Eight orders of magnitude, for a change of basis that does not alter the space being fitted, the fitted curve, or the residual.

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. 2 The same degree-eleven fit to the same thirty points, in two bases. The two curves are drawn on top of one another and are indistinguishable. The condition numbers of the two design matrices differ by eight orders of magnitude, and so does the size of the largest fitted coefficient: 113 against 0.51.

The mechanism is that the monomials on [0,1] look increasingly alike as the degree rises — x⁹ and x¹⁰ agree to within a few percent over most of the interval — so the columns are nearly dependent by construction. Chebyshev polynomials are designed to be as unlike each other as possible on the interval, and are very nearly orthogonal with respect to the natural inner product.

The rules that follow are short and are worth applying to every fit anyone writes:

  • Shift and scale the independent variable to [−1,1] before anything else. This alone recovers several orders of magnitude and costs two arithmetic operations.
  • Use an orthogonal polynomial family rather than monomials for anything above about degree four.
  • Do not read the coefficients of a monomial fit as if they meant something individually. They are large, they alternate in sign, and they cancel; the number 113 in that figure has no interpretation.

What the residual cannot say

The trap this essay exists for: both fits above have excellent residuals. The rms residual of the monomial fit is 3.9·10⁻⁶ and of the Chebyshev fit 2.9·10⁻⁷. Both are tiny; both curves pass through the data.

So a diagnostic based on the residual — which is what R², adjusted R², and every “goodness of fit” statistic is — cannot distinguish a fit whose coefficients are meaningful from one whose coefficients are noise. The residual is the well-conditioned question. It answers itself and says nothing about the other one.

The diagnostics that do say something: the condition number of the design matrix; the singular value spectrum, which shows how many directions are determined; and the sensitivity of the coefficients to a perturbation of the data, which can be measured directly by re-fitting with the data jittered.

Singular values of a rank-4 matrix with noise of relative size 10⁻⁸Ten singular values on a logarithmic axis. The first four sit near one; the rest sit at the noise level, and the vertical distance between the two groups is the evidence for the rank.1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 3 The spectrum as a diagnostic. Ten singular values of a matrix built with four nonzero ones plus noise. The gap says how many directions the data determines; the ones below it are directions the fit will happily assign coefficients to, on the strength of nothing.

Regularisation is a change of question

If the data does not determine the coefficients, no amount of computation will make it. What can be done is to ask a different question, and to say so.

Truncated SVD. Keep the k largest singular values and discard the rest. This is exactly deciding that the directions below the cutoff are not determined and setting them to zero — the most explicit form the choice can take, and the one that makes the number of retained directions visible.

Ridge, or Tikhonov. Minimise ‖b − Ax‖² + λ‖x‖², which adds a multiple of the identity to AᵀA and gives the smallest singular value a floor of √λ. The coefficients become determined, at the cost of biasing them towards zero, and λ is a dial between fidelity and stability.

A smaller model. Frequently the honest answer. If a degree-nine polynomial’s coefficients are undetermined, a degree-four one’s may not be, and the residual will be barely worse — because the extra degrees of freedom were fitting directions the data does not resolve.

All three are legitimate and all three are only legitimate when declared. The failure mode is running a regularised fit and reporting the coefficients as though they came from the unregularised problem, which is a statement about a question that was not asked.

What is asserted here

The flatness is measured, not asserted. The coefficients must be movable by more than 30% of their own size along the worst direction with the residual unchanged in the sixth digit — and that threshold has to be exceeded by orders of magnitude at the third digit, so the shape of the valley is checked and not only one point on it.

κ must exceed 10⁶, which is the claim that a condition number of this size is what a flat valley means.

The two bases must differ by at least a factor of a million in condition number.

And both must fit the data. The Chebyshev fit’s rms residual must be below 10⁻⁴, and the monomial one’s must be within a factor of a hundred of it. That assertion is the essay’s central point in its checkable form: if the badly conditioned fit did not fit the data, the story would be the ordinary one about a numerical failure, and it is not.

One thing had to be corrected while writing this. The first version of the figure measured the residual increase relative to the residual at the optimum, which for a near-perfect fit is about 10⁻¹⁵ — so any perturbation multiplied it by an enormous factor while changing the curve by nothing visible. The measurement is now relative to ‖y‖, which is what “the fit degrades” has to mean if it is to mean anything. The first version was not wrong arithmetically; it was measuring a quantity with no interpretation, and only a plot with a nonsensical axis revealed it.

The condition number of the fit, which is not κ(A)

A technical point that matters as soon as anyone quotes a number.

For a square system, the sensitivity of the answer is κ(A). For least squares it is not, and the correct expression involves the residual as well. Roughly: the coefficients are sensitive like κ(A) when the data lies close to the column space, and like κ(A)² when it does not, with the crossover governed by the angle between b and the space it is being projected onto.

The practical translation. A fit that goes nearly through the data — a small residual — behaves like κ(A), and QR delivers every digit available. A fit to noisy data that the model cannot capture — a large residual — behaves worse, and there is a κ² term that no algorithm removes, because it belongs to the problem rather than to the route.

That is worth knowing because it inverts a natural expectation. The badly fitting model, the one where the residual is large and the modeller is already unhappy, is also the one where the coefficients are least determined. The two problems compound rather than trading off.

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: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 4 The general shape, from the square case. The lower line is what the algorithm contributes and does not move; the upper is the error in the answer. For a least-squares problem with a large residual the upper line rises faster than κ, and the extra factor is a property of the data’s distance from the model.

Cross-validation sees it, and R² does not

If the residual cannot distinguish a determined fit from an undetermined one, something else has to, and the standard tool happens to be the right one for a reason worth stating.

Cross-validation holds out part of the data, fits on the rest, and measures prediction error on what was held out. A fit whose coefficients are determined predicts held-out points well. A fit whose coefficients are floating on the valley floor predicts them badly, because the particular point on the floor that the fit landed on was chosen by the noise in the training data and the held-out data has different noise.

So cross-validation is, among other things, a conditioning diagnostic. The gap between training and held-out error is a measurement of how much of the fit was determined by the signal.

That is why it detects overfitting when R² cannot. R² is a function of the training residual, and the training residual is the well-conditioned question — it answers itself and says nothing about the other one. The same argument explains why adding predictors always increases R² and frequently worsens prediction: each new column reduces the residual a little and can reduce the smallest singular value a great deal.

What “the data does not determine it” means

A closing clarification, because the phrase can sound like an evasion.

It does not mean the answer is unknowable. It means the answer is determined only to a stated precision by the data at hand, and the precision is computable in advance: about 16 − log₁₀κ digits in double precision. That is a quantitative statement, and it can be checked against what is needed.

If four digits are needed and the data determines two, the options are more data, better data, a different parameterisation, or a smaller model — and which of those is available is a question about the experiment rather than about the arithmetic. What is not available is a better solver.

This is the same conclusion as the condition number is an amplifier and an answer that is known reach by other routes, and it is the most useful thing the framework provides: it turns “the numbers look wrong” into “the data supports two digits and four are needed”, which is a statement somebody can act on.

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. 5 The quantitative version. A perturbation of known size applied in two hundred directions, and the amplification measured. That distribution is what “the data does not determine it” means, and its position on the axis is the number of digits lost.
Loss of orthogonality against condition number, in binary64A log–log plot of the norm of Q-transpose-Q minus the identity against condition number. Classical Gram–Schmidt rises steeply, modified Gram–Schmidt rises gently, and Householder is flat.110²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classicalmodifiedHouseholderκ²uκu8×8, eight seeds per κ, binary64all three reconstruct A
Fig. 6 Why the route still matters even when the answer is undetermined. A flat valley means the coefficients are poorly determined by the data; a badly conditioned factorisation means they are additionally poorly determined by the arithmetic. The first cannot be fixed and the second is free to avoid, which is the argument of the road that squares the problem and of two Gram–Schmidts.

The distinction this field turns on is stated most plainly in the exact answer to a nearby problem: a wrong answer has two possible authors. For a fit, the algorithm’s contribution is measurable and usually tiny, and everything else belongs to the data.

A note on how this figure was nearly wrong

The measurement in this essay had to be rebuilt once, and the reason is worth recording because it is a failure mode of measurement rather than of arithmetic.

The first version plotted the increase in the residual relative to the residual at the optimum. That is the natural choice and it is meaningless here: a degree-nine fit to this data has a residual near 10⁻¹⁵, so any perturbation multiplies it by an enormous factor. The curve rose steeply, the conclusion would have been the opposite of the true one, and nothing about the plot looked wrong — it was smooth, monotone, and on sensible axes.

What revealed it was the assertion. The claim being checked was that the coefficients could move a long way with the fit unchanged, and the check failed, because by that measure they could not move at all. Following the failure rather than loosening the tolerance produced the fix: the residual increase has to be measured against the size of the data, which is what “the fit degrades” means to anyone looking at it.

The general lesson is that a ratio needs a denominator with an interpretation. Dividing by a quantity that happens to be near zero produces a number that is large, is correctly computed, and answers no question. That is not a rounding error and no amount of precision would have caught it — assertions that reject caught it, which is the argument for writing them before the figure looks right rather than after.

Where this goes next

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. 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⁻¹⁶. 1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision Rank is a decision A floating-point matrix does not have a rank. It has a spectrum of singular values, and somewhere in that spectrum is a place where the values stop being signal and start being noise. Deciding where is a judgement, and the evidence for it is a gap.