The road that squares the problem

Every derivation of least squares arrives at the same place. The residual must be perpendicular to the column space, Aᵀ(b − Ax) = 0, therefore

AᵀAx = Aᵀb

and the problem is now a square system that can be solved by the methods of the previous field. It is three lines, it is correct, and it is the method that appears in every statistics textbook, in every introductory numerical methods course, and in a great deal of production code.

It is also the reason a fit that ought to give eight good digits gives none.

The squaring, measured #

The singular values of AᵀA are the squares of those of A. So

κ(AᵀA) = κ(A)²

which is checked here rather than quoted: at two values of ε, log₁₀ κ(AᵀA) is required to equal 2 log₁₀ κ(A) to within 0.02, by a route that computes both condition numbers with a one-sided Jacobi SVD of the assembled matrices and shares nothing with the algebra above.

The practical translation uses the arithmetic from the condition number is an amplifier: if κ = 10ᵏ then about 16 − k correct digits can be expected in double and 7 − k in single. Forming AᵀA takes k to 2k. Ten digits become four. Four digits become none.

A polynomial fit of degree six on equispaced points has κ(A) ≈ 10⁵. Through QR that is eleven digits in double, which is plenty. Through the normal equations it is κ = 10¹⁰ and six digits, which is usually still enough — and that is exactly why the method survives. It works, until the degree goes up by three.

The cliff #

The degradation is the expected part. What the figure shows that a bound does not is that below a computable ε the method does not degrade at all — it fails completely, and the failure is a cliff.

Läuchli’s matrix has columns of ones with an ε on the diagonal below, so the entries of AᵀA are 1 + ε² on the diagonal and 1 off it. When ε² falls below the gap next to 1, the diagonal entries round to 1, every entry of AᵀA becomes exactly 1, and the matrix is the all-ones matrix: rank one, exactly singular, with no rounding subtlety about it at all.

That happens at ε = √u, which is 2.4·10⁻⁴ in single precision and 1.5·10⁻⁸ in double. The prediction is made before the measurement and the assertion requires the measured cliff to sit within a factor of three of it, at every precision on the slider.

Two things about that are worth holding on to. The κ of the original problem at the cliff is only about 7·10³ in single precision — a thoroughly ordinary matrix, one that QR handles to seven digits. And the failure is not gradual: at ε just above the cliff the normal equations return an answer with a few digits, and at ε just below it they return nothing, because the solve refuses a singular matrix.

A method whose failure mode is “returns nothing” is at least honest. A slightly different matrix would have produced a nearly singular AᵀA rather than an exactly singular one, and the answer would have been noise with no indication at all.

Why QR is immune #

The QR route computes A = QR and then solves Rx = Qᵀb.

Every operation in that chain is either an orthogonal transformation or a triangular solve. An orthogonal transformation has κ = 1 and amplifies nothing; the triangular system has κ® = κ(A), because R and A differ by an orthogonal factor and orthogonal factors do not change singular values.

So the conditioning of the computation is κ(A). Not κ(A)². The problem’s difficulty is passed through unchanged, which is the most that can be asked of any method, and the guarantee comes from the structure rather than from care.

Why the normal equations survive anyway #

Being fair to the method, because it is not simply a mistake.

AᵀA is smaller. For a fit with a million observations and eight parameters, AᵀA is 8×8 and can be accumulated in one pass with no need to hold A in memory. That is a real advantage and it is why the method persists in statistics, where m is often enormous and n small.

AᵀA is symmetric positive definite, so it can be solved by Cholesky at half the cost of LU, with no pivoting required. The arithmetic is genuinely cheaper: about mn² + n³/3 against 2mn² for Householder.

When the fit is well conditioned it does not matter. With κ(A) = 100, squaring gives 10⁴ and double precision has twelve digits left. Most fits people actually run are in this regime, and for them the normal equations are fine and faster.

The failure is therefore not “this method is wrong” but “this method has a hidden precondition”. The precondition is κ(A) ≪ √(1/u), which is 10⁸ in double and 10³·⁵ in single, and nothing in the derivation mentions it.

Three ways to keep the cheapness #

Where the size argument is what recommends the normal equations, there are ways to keep it.

Accumulate AᵀA in higher precision than the data. The squaring happens at the accumulation, so doing that step in double when the data is single moves the cliff by four orders of magnitude. This is the classic use of extended-precision accumulation and it costs almost nothing.

Use a QR that streams. Householder needs the whole matrix, but a QR can be built by processing rows in blocks and updating R with Givens rotations, holding only R in memory. This is standard in recursive least squares and gives the QR conditioning at the normal equations’ memory cost.

Regularise, and say so. Adding a multiple of the identity to AᵀA — ridge regression — makes the smallest eigenvalue at least λ and the conditioning at most (σ₁² + λ)/λ. That is a genuine fix for the arithmetic, and it is a change of question: the answer is now the solution of a different, better-posed problem. The valley with no bottom is about why that substitution is often the honest thing to do rather than a fudge.

Where the ill-conditioning came from #

The Läuchli matrix is a construction, and it is fair to ask whether the situation is contrived.

It is not, and the reason is a coincidence worth knowing: the Hilbert matrix is the normal-equations matrix of polynomial fitting in the monomial basis. The inner product of xⁱ and xʲ over [0,1] is 1/(i+j+1), which is the Hilbert matrix exactly. So fitting a polynomial by the normal equations is solving a Hilbert system, and the most famously ill-conditioned matrix in the subject is what comes out from the most obvious modelling choice.

κ(H₁₀) = 1.6·10¹³. That is a tenth-degree fit — a routine request — and it means the coefficients are determined to three digits at best in double precision and to nothing at all in single.

What is asserted here #

κ(AᵀA) = κ(A)², to within 0.02 in the exponent, at two values of ε, by an independent route.

The breakdown is at √u. The measured cliff must sit between √u/30 and 3√u, at every precision on the slider. Since the slider covers sixteen to forty bits, that is nine independent confirmations of a prediction made from the arithmetic rather than fitted to the data.

QR is accurate across the whole range: worst forward error below 10⁻⁴, over seven decades of ε — and the normal equations’ worst is at least fifty times larger at every precision on the slider. That comparison is between the two worst cases rather than between the values at the last surviving ε, because at that point both curves are sitting on the plot’s floor and their ratio says nothing. The first version of this check compared them there and failed at 19 bits, where the normal equations were momentarily the better of the two — an artefact of a clamped axis, caught because the assertion runs at every position of the slider rather than at the default.

AᵀA really is the all-ones matrix at the cliff. Not nearly singular — every entry is required to be exactly 1, which is a much stronger and much more checkable statement than “the matrix becomes ill-conditioned”.

And the loss is graded at milder ε, so the essay’s claim about degradation is separately supported: at ε = 10⁻², the normal-equations error is 146 times QR’s, which is a loss rather than a collapse.

The last one was added after the first version of this check tried to demonstrate degradation at an ε where the method had already broken down, and asserted a factor between an error and infinity. That is the kind of thing writing the assertions catches and reading the prose does not.

What Cholesky adds, and what it does not #

The normal-equations route is usually completed with Cholesky rather than LU, since AᵀA is symmetric positive definite, and it is worth being precise about what that buys.

Cholesky costs half of LU, needs no pivoting, and has a growth factor of 1 — the bound that is never attained does not apply, because for this class the worst case is provably benign. So the solve stage of the normal-equations route is as stable as anything in the subject.

None of which helps, because the damage was done in the previous line. The conditioning was destroyed when AᵀA was formed, and a perfectly stable solve of a badly conditioned system returns exactly the error the conditioning permits. Cholesky is backward stable for the system it is given; the system it is given is the wrong one.

There is a diagnostic hidden in it, though, and it is the one useful thing this route provides. Cholesky fails — hits a nonpositive number under a square root — precisely when AᵀA is not positive definite, which in exact arithmetic it always is. So a Cholesky failure on a normal-equations matrix is a definitive signal that the conditioning has exceeded what the arithmetic can carry. It is a cliff detector, and it is free.

The seminormal equations, and why they are not the answer either #

Between the two roads sits a third that is occasionally proposed and worth dismissing carefully.

If R has already been computed from a QR factorisation, then RᵀR = AᵀA, so the normal equations can be solved as two triangular solves without ever forming AᵀA explicitly. This is the seminormal equations method, and it looks like it should inherit QR’s conditioning while costing less.

It does not. The accuracy is governed by κ(A)² regardless, because the two triangular solves with R and Rᵀ compose to the same badly conditioned operator. What is true is that one step of iterative refinement on the seminormal equations restores a κ(A) error bound — the corrected seminormal equations — which is a genuine result and a genuinely useful method when Q was discarded.

The lesson generalises: avoiding the explicit formation of AᵀA is not the point. The point is avoiding the operator, and any factorisation of it inherits its conditioning.

What to do instead, concretely #

For a dense fit, the recommendation is short and has no interesting exceptions.

Use a QR-based least-squares routine. dgels in LAPACK, numpy.linalg.lstsq, scipy’s default. They are all Householder QR or SVD underneath, and they cost at most twice the normal equations.

Scale the columns. Dividing each column by its norm costs nothing and often removes an order of magnitude or two of conditioning that came from unit choices rather than from the problem.

Choose the basis deliberately. For polynomial fitting this is worth more than everything else combined: eight orders of magnitude for a change of basis that does not alter the fit. The valley with no bottom has the measurement.

If accumulation is unavoidable, because the data does not fit in memory, accumulate AᵀA in a wider precision than the data. The squaring happens at the accumulation, and doing it in double when the data is single moves the failure point by four orders of magnitude in κ.

What this field establishes #

The projection and the right angle is the geometry, with the perpendicularity condition measured to 10⁻¹⁶ rather than asserted. This essay is about the one modelling decision that destroys it. And the valley with no bottom is about the harder thing underneath: that even the correct route cannot recover coefficients the data does not determine.

Read in order, they separate three questions that arrive together in every fitting problem. Is the projection computed correctly? Is the route to the coefficients throwing away conditioning? And are the coefficients determined by the data at all? Only the first two are about arithmetic, and only the first two have fixes.

Two neighbouring essays are worth reading beside this one. The condition number is an amplifier is where the 16 − log₁₀κ arithmetic comes from, and an answer that is known is the case where the loss can be measured against an exactly known truth rather than inferred from a bound.

The number to check before choosing #

If only one diagnostic survives from this essay, it should be this one, because it takes a line and settles the question in advance.

Estimate κ(A) — from the R factor of a QR, from a condition estimator, or from the singular values if the matrix is small. Then compare it against √(1/u): about 10⁸ in double, about 3·10³ in single.

Below that, the normal equations lose at most half the available digits and often fewer, and if the remaining digits are enough then the memory advantage is real and worth having.

Above it, the normal equations lose everything, and somewhere not much further above it they stop returning an answer at all.

The comparison is a single line and it converts a methodological argument into an arithmetic one. It also explains why the debate persists: below the threshold both methods work, the normal equations are cheaper, and anyone whose problems have all been in that regime has correct experience that the method is fine. The disagreement is not about the mathematics; it is about which problems the participants have met.

Which is the argument for the threshold rather than for the rule. “Never use the normal equations” is a slogan that sounds like superstition to someone for whom they have always worked. “Compare κ against 10⁸” is a check anyone can run on their own data, and it produces the slogan when the slogan is right.