Where the last guide left us
You arrive here holding a complaint. To fit a model to more data than unknowns you set up an overdetermined system A x = b and asked for the x that makes the residual ||A x - b||_2 smallest. The previous guide handed you the normal equations A^T A x = A^T b, a clean square system you can solve with Cholesky — but it also delivered the bad news: forming A^T A squares the condition number. If A has condition number 10^4, then A^T A has condition number about 10^8, and via accuracy = conditioning times stability you have silently thrown away about half of your roughly 16 double-precision digits before the solver even starts.
The natural question is whether that loss is fundamental or just self-inflicted. It is self-inflicted. The original least-squares problem usually has condition number around that of A, not A^2; the squaring is an artifact of the method, not the problem. This guide is about a different route to the same x — one that never forms A^T A, never squares anything, and is provably backward stable on the least-squares problem as it actually stands. The trick, as so often in numerical linear algebra, is orthogonality.
Why orthogonality saves the day
Recall the geometry from earlier in this rung: the least-squares solution is an orthogonal projection of b onto the column space of A, and the residual sticks out perpendicular to that space. Orthogonal directions are special for a reason that is also numerical: an orthogonal transformation Q preserves length, ||Q v||_2 = ||v||_2 for every vector v. It rotates and reflects but never stretches. So if you hit both sides of A x - b with such a Q, the residual norm you are minimizing does not change at all — you are looking at the same problem from a rotated camera angle.
That is the whole strategy. Choose the rotation cleverly so that, in the new view, the minimization becomes trivial to read off. Because the rotation changes no lengths, it changes no condition numbers either: it cannot make the problem harder or easier, only clearer. Compare this to the normal-equations route, which transforms the problem by A^T — an operation that emphatically does change lengths and squares the conditioning. Orthogonal first, everything else second: that instinct will serve you across the whole subject, from eigenvalue solvers to the SVD itself.
QR: the everyday workhorse
The QR factorization writes the tall, skinny matrix A (more rows than columns) as A = Q R, where Q has orthonormal columns and R is upper-triangular and square. Substitute it into the residual and use the length-preserving property of the orthogonal part: minimizing ||A x - b||_2 turns into minimizing ||R x - Q^T b||_2. But R is upper-triangular, so this is just a small triangular system R x = Q^T b — solved instantly by the back substitution you have done since the direct-solvers rung. No A^T A is ever formed, and the conditioning of R matches that of A, not its square.
Least squares two ways for min || A x - b ||_2 :
normal equations : A^T A x = A^T b cond ~ cond(A)^2 (risky)
QR route : A = Q R cond ~ cond(A) (stable)
then solve R x = Q^T b by back substitution
why the residual norm is unchanged:
|| A x - b || = || Q(R x) - b || = || Q^T(Q R x - b) || = || R x - Q^T b ||
(Q^T preserves length, so the thing we minimize never changes)How is the factorization itself computed? You met Gram-Schmidt in the eigenvalue rung — orthogonalize the columns one by one — and it does produce a QR factorization, but the plain version loses orthogonality badly in floating point and is not recommended for serious work (modified Gram-Schmidt is better but still second-tier). The professional method builds Q out of Householder reflections: each reflection is an orthogonal matrix that zeroes out everything below a diagonal entry in one column, and chaining n of them triangularizes A. Householder QR is backward stable and is what LAPACK actually runs when you ask for a least-squares solve. For a matrix that is already nearly triangular or arrives one row at a time, Givens rotations — orthogonal matrices that zero a single entry — are the tidier choice.
The SVD: the truth-teller for sick problems
QR is the right default, but it assumes A has full column rank — that its columns are genuinely independent. When they are nearly dependent (you will see exactly this with the Vandermonde basis of high-degree polynomial fitting in the next guide), R has a tiny diagonal entry, the back substitution divides by something near zero, and even a perfectly stable algorithm returns a wildly amplified, untrustworthy x. The problem itself is ill-conditioned, and no stable method can rescue an ill-conditioned problem — a stable algorithm gives the exact answer to a nearby problem, and here nearby problems have wildly different answers. What you need now is not a faster solve but a diagnosis.
The singular value decomposition is that diagnosis. From the eigenvalue rung you know it factors A = U S V^T, where U and V are orthogonal and S is diagonal carrying the singular values sigma_1 >= sigma_2 >= ... >= 0. For least squares the SVD does something QR cannot: it lays the conditioning bare. The ratio sigma_max / sigma_min is exactly the condition number, so a glance at the singular values tells you whether your fit is healthy or doomed. Tiny singular values are the precise fingerprint of near-dependent columns — the directions in which the data barely constrains the model.
The SVD also delivers the most honest single formula for the least-squares answer: the Moore-Penrose pseudoinverse A^+ = V S^+ U^T, where S^+ inverts each nonzero singular value (replacing sigma by 1/sigma) and leaves the zeros alone. Then x = A^+ b is the least-squares solution — and, when the columns are dependent so the answer is not unique, it is the specific solution of smallest norm. Be careful with the word 'inverse': the pseudoinverse is a clean way to describe the solution, but you almost never form A^+ as an explicit matrix to solve a real problem, exactly as you never form a true inverse to solve A x = b. You compute the SVD and apply its pieces.
Which route, and when
You now own three routes to the same least-squares x, and choosing between them is a real engineering decision rather than a matter of taste. Here is the honest ranking, cheapest first, safest last.
- Normal equations (A^T A x = A^T b, then Cholesky): cheapest, about half the flops of QR, and fine when A is well-conditioned and you only need a few digits — but it squares the condition number, so retire it the moment accuracy or conditioning is in question.
- QR via Householder (A = Q R, then solve R x = Q^T b): the everyday default for full-rank problems. Backward stable, conditioning matches A not A^2, costs about twice the normal equations — a price you should almost always pay.
- SVD (A = U S V^T, then the pseudoinverse): the most expensive but the most informative. Use it when the columns may be rank-deficient or near-dependent, when you need the minimum-norm solution, or when you want to read the condition number straight off the singular values before deciding to regularize.
A practical rule of thumb that ties the rung together: default to QR; reach for the SVD the instant you suspect ill-conditioning, because only it tells you how ill-conditioned you are and gives you the singular-value handle to do something about it. The next and final guide of this rung picks up exactly there — it meets the ill-conditioned Vandermonde matrix of high-degree polynomial fitting head-on, and shows how truncating small singular values or adding a Tikhonov penalty turns a doomed fit into a useful one by trading a little bias for a lot of stability.