From projection to a square system
The previous guide left us with a clean picture: in the least-squares problem we cannot solve the tall system A x = b exactly, because b almost never lives in the column space of A. The best we can do is land on the closest point in that column space, and the closest point is the orthogonal projection of b onto it. The defining feature of that projection is that the leftover, the residual r = b - A x, points straight out of the column space — it is perpendicular to every column of A.
Now turn that geometric sentence into algebra. 'The residual is perpendicular to every column of A' means every column, dotted with r, gives zero. Stacking those dot products is exactly the matrix product A^T r = 0. Substitute r = b - A x and you get A^T (b - A x) = 0, which rearranges into the famous normal equations A^T A x = A^T b. We started with m equations in n unknowns that had no solution; we have ended with n equations in n unknowns that does.
Why the recipe looks irresistible
On paper the normal equations are a dream. The matrix A^T A is square, symmetric, and (when A's columns are independent) positive definite — exactly the class of matrix that the cleanest, fastest direct solver loves. You do not even need a general-purpose method: a Cholesky factorization, which writes A^T A = L L^T using only the lower triangle, solves the system in about half the work of standard elimination. For a small number of parameters n, the cost is trivial and the code is three lines.
- Form the small square matrix M = A^T A and the vector c = A^T b — for n parameters this is just n-by-n and length n.
- Cholesky-factor it, M = L L^T, using only the lower triangle of the symmetric positive-definite M.
- Solve L y = c by forward substitution, then L^T x = y by back substitution; the resulting x minimizes ||A x - b||_2.
So the normal equations are the first thing every course teaches, and for understanding the geometry they are perfect. The trouble is purely numerical, and it hides entirely inside that innocent first line, M = A^T A. The act of forming that product — before any solver even runs — quietly damages the problem. To see how, we have to talk about conditioning.
The hidden danger: squaring the condition number
Recall the master rule from the conditioning rung: accuracy = conditioning x stability. Even a flawless, backward-stable solver can only deliver digits that the problem's condition number leaves intact; a condition number near 1e8 already costs you about 8 of your roughly 16 double-precision digits. The conditioning of a least-squares fit is governed by the condition number of A, written cond(A). The brutal fact is this: the matrix you actually hand to the solver, A^T A, has condition number cond(A^T A) = cond(A)^2.
Read that exponent slowly, because it is the whole guide. Squaring a condition number doubles the number of digits it eats. If A is mildly ill-conditioned with cond(A) = 1e4 — a perfectly ordinary value for a real fit — then your original problem loses about 4 digits, which is fine; you keep 12. But A^T A has condition number 1e8, so the normal-equations route loses about 8 digits. You have thrown away 4 good digits of accuracy not because of the data, not because of the solver, but solely because of the algebraic step of forming A^T A. The damage is done before Cholesky touches a single number.
A tiny example you can almost do by hand
Here is the classic miniature that makes the danger concrete — the Lauchli matrix. Let A have columns that are nearly parallel: take A = [[1, 1], [eps, 0], [0, eps]] where eps is small, say eps = 1e-8. The columns are independent, so the fit is perfectly well posed in principle, and cond(A) is roughly 1/eps = 1e8 — large but survivable in double precision, which carries about 16 digits.
A = [ 1 1 ] eps = 1e-8
[ eps 0 ]
[ 0 eps]
A^T A = [ 1 + eps^2 1 ]
[ 1 1 + eps^2 ]
In double precision, 1 + eps^2 = 1 + 1e-16 -> ROUNDS to exactly 1.
so the stored matrix becomes [ 1 1 ] <- SINGULAR (rank 1!)
[ 1 1 ]
cond(A) ~ 1e8 (the honest problem: solvable)
cond(A^T A) ~ 1e16 (what you actually solve: numerically singular)Watch what the first line of the recipe does. The diagonal of A^T A is 1 + eps^2 = 1 + 1e-16, and in double precision that addition rounds straight back to 1 — recall that floating point cannot represent 1 + 1e-16, and that addition is not exact. The off-diagonals are exactly 1. So the matrix you store is [[1,1],[1,1]], which is singular: its two rows are identical. Cholesky will fail outright or return garbage, and yet the underlying fit was solvable. The information that distinguished the two columns lived in the eps^2 terms, and forming A^T A rounded it out of existence.
This is the danger in one frame. The problem you wanted to solve had condition number 1e8 and was comfortably inside double precision's reach. The matrix you chose to solve had condition number 1e16, right at the edge of total breakdown. Nothing about the data was bad; the loss was entirely self-inflicted by the algorithm. A stable solver on an unstable formulation still gives you a bad answer — because accuracy is conditioning times stability, and you doubled the conditioning yourself.
When it is fine, and where to go instead
None of this means the normal equations are forbidden. When A is well-conditioned — cond(A) of order 1 up to maybe 1e3 or so — squaring it still leaves you plenty of digits, and the normal equations are fast, simple, and entirely respectable. Their compactness genuinely matters when you have millions of data rows and only a handful of parameters, since A^T A is tiny and can be accumulated row by row without ever storing all of A. The rule is not 'never use them'; it is 'know the condition number before you do.'
But this danger is no edge case, because one of the most common fitting tasks walks straight into it. Polynomial regression — fitting a degree-d curve through data — builds its matrix A from columns 1, x, x^2, x^3, and so on. Those power columns look more and more alike as the degree rises, so this Vandermonde matrix becomes spectacularly ill-conditioned; even degree 10 on equispaced points can push cond(A) past 1e10. Run that through the normal equations and you square it to 1e20, far beyond any double-precision rescue. The next guide is devoted to the cure.
The escape route is to never form A^T A at all. Instead of squaring the problem and then solving, the stable methods work on A directly: a QR factorization of A, or in the hardest cases the singular value decomposition and the pseudoinverse. These pay for the least-squares solution at condition number cond(A), not cond(A)^2 — they keep all the digits the problem genuinely allows. That is exactly the road the next two guides build, so hold this danger in mind as the reason they exist.