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.

Worth reading first: The projection and the right angle · The condition number is an amplifier · The road that squares the problem.

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.

Which makes “eight orders” a number about degree eleven rather than about the two bases, and the sweep prices the difference.

Fitting the same degree-9 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 matrixmonomial3.5·10⁶Chebyshev2.3largest fitted coefficientmonomial107Chebyshev0.51rms residual: 1.2·10⁻⁵ and 1.2·10⁻⁵ — the data is fitted either way.30 points, degree 9, single precisionthe basis is part of the problem
Fig. 3 Degree nine. The two fitted curves are already indistinguishable; the monomial design matrix has κ = 3.5·10⁶ against Chebyshev’s 2.26 — 6.2 orders.
Fitting the same degree-13 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 matrixmonomial4.6·10⁹Chebyshev3.2largest fitted coefficientmonomial85Chebyshev0.51rms residual: 3·10⁻⁶ and 2.8·10⁻⁷ — the data is fitted either way.30 points, degree 13, single precisionthe basis is part of the problem
Fig. 4 Degree thirteen: κ = 4.6·10⁹ against 3.16, 9.2 orders, and the curves still lie on top of each other.

The gap opens by about 1.4 orders every two degrees, and it opens from one side. Across degrees 9, 11, 13, 15 and 17 the monomial condition number runs 3.5·10⁶, 1.2·10⁸, 4.6·10⁹, 1.9·10¹¹ and 8.5·10¹² — a factor of 2.4 million — while Chebyshev’s runs 2.26, 2.53, 3.16, 4.85 and 10.18, a factor of 4.5. So the orthogonal basis degrades too, and by an amount nobody would notice.

Fitting the same degree-15 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.9·10¹¹Chebyshev4.9largest fitted coefficientmonomial83Chebyshev0.51rms residual: 2.4·10⁻⁶ and 2.7·10⁻⁷ — the data is fitted either way.30 points, degree 15, single precisionthe basis is part of the problem
Fig. 5 Degree fifteen: 1.9·10¹¹ against 4.85, 10.6 orders.
Fitting the same degree-17 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 matrixmonomial8.5·10¹²Chebyshev10largest fitted coefficientmonomial88Chebyshev0.51rms residual: 3.7·10⁻⁶ and 3·10⁻⁷ — the data is fitted either way.30 points, degree 17, single precisionthe basis is part of the problem
Fig. 6 And degree seventeen: 8.5·10¹² against 10.18 — 11.9 orders, with the fits as indistinguishable as they were at nine.

The part that should be uncomfortable is the left-hand panel at every one of those degrees. The two curves lie on top of one another at degree nine and still lie on top of one another at seventeen, by which point one of the two design matrices has lost thirteen digits and the other has lost one. Nothing a reader can see distinguishes the fit that is fine from the fit that is being computed on a matrix a double cannot invert. The picture is identical and the arithmetic is not, which is the whole argument for choosing the basis rather than inspecting the result — by the time the curve looks wrong the conditioning has been unusable for several degrees. 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.

Which of the first two rules is doing the work

The figure changes two things at once — the interval and the family — and reports eight orders for the pair. The rules name them separately without saying what each is worth, so the third matrix, monomials on [−1,1], is what settles it. Thirty equispaced nodes:

degree monomials [0,1] monomials [−1,1] Chebyshev the shift buys the family buys
3 1.1·10² 7.7 2.1 15× 3.7×
5 3.4·10³ 4.0·10¹ 2.2 85× 18×
7 1.1·10⁵ 2.2·10² 2.2 495× 99×
9 3.5·10⁶ 1.2·10³ 2.3 2,890× 538×
11 1.2·10⁸ 7.2·10³ 2.5 17,100× 2,840×
15 1.9·10¹¹ 3.0·10⁵ 4.9 633,000× 61,400×

The shift is the larger half at every degree, and at the figure’s own degree of eleven it is four of the eight orders against the family’s three. Two arithmetic operations buy more than the change of basis they are conventionally mentioned before, which inverts the order the rules are usually given in and is the practical content of this section.

Neither gain is a constant. Both grow by roughly a factor of two per degree, because a Vandermonde on [0,1] becomes ill conditioned exponentially faster than one on [−1,1] rather than by a fixed amount — so the advice gets more valuable with degree and there is no degree at which either rule stops paying.

And the second rule’s threshold turns out to be a statement about unshifted monomials. Once the interval is fixed, monomials reach degree 7 at κ = 215 and degree 9 at κ = 1.2·10³, both of which are perfectly ordinary numbers to solve with. The degree at which the family starts to matter is therefore nearer nine than four — which does not make the orthogonal family a bad idea, since it costs nothing either, but does mean that a code that shifts and uses monomials is in a very different position from one that does neither, and the rule as stated does not distinguish them.

It is worth saying what the shift is doing, because it looks like a triviality and is not. On [0,1] every monomial is a non-negative function that rises to one at the right-hand end, so all eleven columns point into the same corner of the space and differ only in how fast they get there. On [−1,1] the odd powers are odd and the even powers are even, so half the columns are orthogonal to the other half before any argument about approximation theory — the separation is a symmetry rather than a construction. That is the whole of the four orders, and it is why the gain is exponential in the degree: each additional degree adds a column to whichever parity class it belongs to, and on [0,1] there is only one class.

That argument predicts that the shift should weaken when the sample points are not placed symmetrically about the centre, since then there is no parity to exploit. Measured with the thirty points cubed towards one end or the other, at degree eleven, it is the opposite that happens: the shift still buys 6.2·10³ clustered at the left and 2.2·10⁴ clustered at the right, against 1.7·10⁴ equispaced — the same order in every case.

What clustering erodes is the family’s advantage. Chebyshev’s own condition number rises from 2.5 to 37 under either clustering, and its gain over shifted monomials falls from 2,840 to 780, because Chebyshev polynomials are orthogonal with respect to their own measure and clustered sampling is not it. So the ranking survives a departure from the ideal sampling and widens: the two arithmetic operations are the robust half of the advice, and the orthogonal family is the half that assumes something about where the data was taken.

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.

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.

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.

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.

What a perturbation of 10⁻⁶ moves, against the gap it is applied acrossThree curves against the gap between two eigenvalues, both axes logarithmic. Two are flat; the third rises steeply as the gap closes and then levels off.10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻⁸10⁻⁵10⁻²gap between the two eigenvalueshow far it movedthe eigenvectorsthe eigenvaluestheir plane‖E‖ / gapone perturbation, three answerseigenvalue shift, spread over the sweep1plane angle, spread over the sweep1eigenvector angle, spread1.6·10⁵the dashed line is Davis–Kahan's ‖E‖/gaptwo of the three never noticed
Fig. 7 The same flat valley in the eigenvalue setting. When two eigenvalues are close, the direction within their plane is undetermined for exactly the reason a fit’s direction is here — and the plane, like the fitted value, is determined perfectly well.

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.

An observation that is holding a direction up

A flat valley is a direction the data does not determine. Its localised version is a single observation that determines a direction by itself, and there is a number for it.

What links here

Computed from the collection, not written here: the essays that point at this one.

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.

Basis choiceCondition numberIll-posednessRegularisationResidualSingular valuesSingular vectorsThe Vandermonde matrix