Why solve a least-squares system by QR factorization rather than forming A transpose A?

## Two routes to the same mathematical answer Given a tall matrix `A` with m rows and n columns and a right-hand side `b`, the least-squares problem is to minimise `||b - Ax||`. Two standard computational routes exist. The **cross-product route** forms `A^T A` and `A^T b`, then solves the small n-by-n system - usually with Cholesky, since `A^T A` is symmetric and, with independent columns, positive definite. The **QR route** factors `A` directly. Both give the same answer in exact arithmetic. They differ sharply in floating point. ## Why forming the cross product hurts In the 2-norm, `cond(A^T A) = cond(A)^2`. Since you lose about `log10(cond)` significant digits when you solve a system, squaring the condition number roughly doubles the loss. A matrix with `cond(A) = 10^6` - not unusual once columns sit on different scales - hands the solver a system with conditioning `10^12`, leaving about 4 of your 16 double-precision digits. Worse, the loss happens during formation, not during the solve. Consider the 3-by-2 matrix with rows `(1, 1)`, `(e, 0)` and `(0, e)` for a small positive `e`. Its columns are plainly independent for any nonzero `e`. The cross-product matrix has diagonal entries `1 + e^2` and off-diagonal entries `1`. Take `e = 10^-9`; then `e^2 = 10^-18`, which is below double precision's roughly `2.2 * 10^-16` resolution, so `1 + e^2` rounds to exactly `1`. The computed cross-product matrix is `[[1, 1], [1, 1]]` - exactly singular. The information that distinguished the columns was destroyed by the multiplication, and no solver, however careful, can recover it. ## How QR sidesteps it QR factorization writes `A = QR`, where `Q` has orthonormal columns (each of unit length and mutually perpendicular) and `R` is upper triangular. Orthogonal transformations preserve the 2-norm, which is the whole point: applying `Q^T` to a vector does not change its length, so it cannot magnify error. Substituting the factorization, ``` ||Ax - b|| is minimised by solving R x = Q^T b ``` which is a triangular system dispatched by back substitution. The conditioning you actually face is `cond(A)`, not its square, and Householder QR - which builds `Q` implicitly out of reflections rather than forming it - is backward stable. The residual can be read off directly as the part of `Q^T b` that `R` cannot reach. ## The cost side of the tradeoff The cross-product route is cheaper. Forming `A^T A` costs about `m n^2` operations and factoring the small n-by-n result adds only about `n^3 / 3`. Householder QR runs at roughly `2 m n^2`. So you pay about a factor of two for the stability. When `m` is enormous relative to `n` and the columns are well behaved, that factor can be worth saving - and the cross-product route has a genuine operational advantage: `A^T A` and `A^T b` can be accumulated a chunk of rows at a time, so a data set that never fits in memory can still be reduced to a small n-by-n system in a single streaming pass. ## Choosing between them Use the cross-product route when `n` is small, the columns are on comparable scales and known to be well separated, and either speed or streaming accumulation matters. Use QR when the columns are close to dependent, when they sit on wildly different scales, when the condition number is unknown, or whenever the answer feeds a decision that must be defensible. Since the extra cost is a constant factor rather than a change in the growth rate, QR is the sensible default and the cross-product route is the considered optimisation. ## When the columns are actually dependent If `A` does not have full column rank, no minimiser is unique - many coefficient vectors achieve the same minimum residual. Plain QR does not resolve that ambiguity; the standard extension is **column-pivoted QR**, which reorders columns to expose the dependent ones so you can see which columns carry no independent information and drop or combine them. Detecting the deficiency is the useful part: it tells you the measurement design, not the arithmetic, needs fixing. ## What a strong answer contains Name the squaring of the condition number, give the mechanism - information lost while forming the product, not while solving - explain that orthogonal transformations preserve lengths and so cannot amplify error, quote the roughly two-times cost, and concede the genuine cases where the cross-product route still wins. Simply asserting that the cross-product route is bad, with no numbers and no counter-case, is the shallow version.