Solving Least Squares Reliably: Cholesky, QR, and SVD
The formula $\hat\beta=(X^\top X)^{-1}X^\top y$ describes the OLS solution, but it is usually a poor recipe for computing it. Forming $X^\top X$ squares the condition number, and explicitly taking an inverse adds more numerical error. A model can therefore have a tiny residual and still have unreliable coefficients.
What the Condition Number Means
For full-rank $X$, the condition number is $\kappa_2(X)=\sigma_{\max}/\sigma_{\min}$. A large value means some predictor direction is weakly determined. Small changes in the data can then cause large changes in the coefficients.
Ill-conditioning belongs to the data and feature design. Numerical instability belongs to the algorithm. A stable solver prevents avoidable error, but it cannot create information that nearly duplicate columns do not contain.
Four Solver Families
| Method | Work for dense $n\times p$ | Numerical behavior | Best fit |
|---|---|---|---|
| Cholesky on $X^\top X$ | $O(np^2+p^3/3)$ | Fast, but squares the condition number | Well-scaled, clearly full-rank data |
| Householder QR | $O(2np^2-2p^3/3)$ | Stable and avoids forming $X^\top X$ | Default dense least squares |
| Pivoted QR | QR plus pivot search | Reorders columns to reveal weak directions | Possible rank problems |
| SVD | Typically $O(np^2)$ with a larger constant | Shows rank and returns a minimum-norm solution | Rank deficiency and diagnostics |
QR writes $X=QR$ and solves a triangular system. Cholesky factors $X^\top X$, which is faster but less stable. SVD writes $X=U\Sigma V^\top$ and solves one singular direction at a time:
$$\hat\beta=X^+y=\sum_{\sigma_j>t}\frac{u_j^\top y}{\sigma_j}v_j,$$The threshold $t$ decides which directions count as numerically nonzero. Removing a direction changes the model and acts like regularization, so the threshold should be reported rather than hidden.
Rank, Scaling, and Thresholds
A common automatic threshold is $t=\epsilon\max(n,p)\sigma_{\max}$, although noisy applications may need a larger value. Feature scaling changes singular values, so rank decisions should be made after units are handled sensibly.
Iterative Solvers for Large or Sparse Systems
Dense QR and SVD require the matrix to fit in memory. Very large or sparse problems often use iterative methods that need only products with $X$ and $X^\top$.
CGNR applies conjugate gradient to $X^\top X\beta=X^\top y$. It avoids storing $X^\top X$, but it does not avoid the squared condition number. It can therefore stop with a small residual while the coefficients are still inaccurate.
LSQR and LSMR work directly with $X$ through bidiagonalization and usually behave better in finite precision. LSQR focuses on the data residual $\|X\beta-y\|$, while LSMR focuses on the normal-equation residual $\|X^\top(X\beta-y)\|$. Both can include ridge-style damping.
Preconditioning and Stopping Rules
A preconditioner rescales or transforms the problem so iterative steps converge faster. Simple column scaling already removes large unit differences. Stronger preconditioners approximate part of a QR factorization.
Stopping should use more than an iteration count. Useful checks include a relative residual, a maximum budget, lack of recent improvement, and finite-value guards.
Small Singular Values Amplify Noise
▶ Dividing by Smaller Singular Values
The same signal, $u_j^\top y=0.04$, is divided by successively smaller singular values. The resulting coefficient grows rapidly.
Practical Solver Choice
Moderate dense problems usually favor QR. Very large or sparse problems favor sparse QR, LSQR, or LSMR with suitable preconditioning. SVD is the clearest fallback when rank is uncertain. Solver logs should include numerical rank, a condition estimate, residual norm, stopping tolerance, and any discarded directions.
For $n<p$, many exact fits may exist. SVD returns one minimum-norm solution, but Ridge may express a more useful stability preference. Duplicate columns, invalid values, zero weights, and extreme unit differences need explicit checks.
The C++ QR and SVD guide implements these choices. The normal-equation baseline appears in the C++ OLS guide.