Numerical methods for linear least squares

Remarks

Example: for fitting a line $y = c_0 + c_1 t$ to the points $(0, 0)$, $(1, 1)$, $(2, 1)$, take $A$ with rows $(1, 0)$, $(1, 1)$, $(1, 2)$ and $b = (0, 1, 1)$. Then $A^T A = \begin{pmatrix} 3 & 3 \\ 3 & 5 \end{pmatrix}$ and $A^T b = (2, 3)$, so $c_0 = 1/6$ and $c_1 = 1/2$; the residual $(-1/6, 1/3, -1/6)$ is orthogonal to both columns of $A$, as it must be. Non-example: solving the square system formed by the first $n$ rows of $Ax = b$ fits those rows exactly, but it does not give the minimizer in general. Here it gives $c_0 = 0$, $c_1 = 1$, with residual norm $1$ instead of $1/\sqrt{6}$. The normal equations are the cheapest method but square the condition number, $\kappa(A^T A) = \kappa(A)^2$. Householder QR is backward stable and is the standard choice for dense full-rank problems; the SVD is the most robust for rank-deficient problems, at a higher cost. For large sparse problems one uses iterative methods such as LSQR or the conjugate gradient method applied to the normal equations.