Least squares, three ways

Linear least squares asks for the coefficients $\hat\beta$ that minimize $\lVert X\beta - y \rVert_2$ for a tall matrix $X \in \mathbb{R}^{m \times p}$. Setting the gradient to zero gives the normal equations,

and in NumPy they take one line:

beta = np.linalg.solve(X.T @ X, X.T @ y)

That line is short, fast, and the least accurate of the three standard methods. This note measures how much accuracy it gives up, explains why, and says when the loss does not matter.

Polynomial fitting is the classic way to make a least-squares problem hard. On $m = 100$ equally spaced points $t_i \in [0, 1]$, the columns $1, t, t^2, \dots, t^{p-1}$ of the monomial basis look more and more alike as $p$ grows, and the matrix drifts towards rank deficiency.

To isolate the solver, I chose the right-hand side so that the exact solution is known: $\beta = (1, \dots, 1)$ and $y = X\beta$, with a residual of zero. Any error in $\hat\beta$ then comes from floating-point arithmetic, not from the data.

import numpy as np

def errors(p, m=100):
    t = np.linspace(0.0, 1.0, m)
    X = np.vander(t, p, increasing=True)
    beta = np.ones(p)
    y = X @ beta
    rel = lambda b: np.linalg.norm(b - beta) / np.linalg.norm(beta)

    normal = np.linalg.solve(X.T @ X, X.T @ y)
    Q, R = np.linalg.qr(X)
    qr = np.linalg.solve(R, Q.T @ y)
    svd = np.linalg.lstsq(X, y, rcond=None)[0]
    return np.linalg.cond(X), rel(normal), rel(qr), rel(svd)

np.linalg.qr computes a Householder factorization, and np.linalg.lstsq calls LAPACK’s SVD-based solver. The relative errors in the coefficients:

With twelve columns the normal equations return coefficients that are wrong in the first digit, while QR and the SVD still agree with the truth to eight or nine digits. Plotting every $p$ from 2 to 14 against the condition number shows two straight lines with different slopes.

The condition number $\kappa_2(A) = \sigma_{\max}(A) / \sigma_{\min}(A)$ measures how much a matrix can amplify relative perturbations. Forming $X^\top X$ squares it.

A backward-stable solver for a square system with matrix $A$ returns a solution with relative error of order $\kappa(A)\,u$, where $u = 2^{-53} \approx 1.1 \times 10^{-16}$ is the unit roundoff of double precision. Solving the normal equations is such a solve with $A = X^\top X$. Householder QR and the SVD never form that product: they work on $X$ itself, and for a problem with a small residual their error is of order $\kappa(X)\,u$.¹

The estimates are upper bounds up to modest constants, and the measured errors sit at or below the dashed guides in the figure. Each factor of ten in $\kappa(X)$ costs QR about one decimal digit and the normal equations about two.

When the residual is large, even the best algorithm inherits a term proportional to $\kappa(X)^2$ times the relative size of the residual. That term belongs to the problem, not to the method, so no choice of solver removes it.²

None of this makes the normal equations wrong. They are the right tool more often than the table suggests.

• The matrix is well conditioned. With $\kappa(X) = 100$, the normal equations lose about four digits and still deliver twelve. Standardized predictors with little collinearity are often in this range.
• The data do not fit in memory. $X^\top X$ and $X^\top y$ can be accumulated in one pass over the rows, in $p \times p$ storage, and solved at the end.
• Speed matters more than the last digits. For $m \gg p$, forming $X^\top X$ costs about $m p^2$ operations; Householder QR costs about $2 m p^2$.

The cheapest fix is often in the model rather than the algorithm. The monomials $t^k$ on $[0, 1]$ are a notoriously poor basis. Moving the same points to $[-1, 1]$ brings the condition number of the twelve-column matrix from $1.2 \times 10^{8}$ down to $6.9 \times 10^{3}$, and a Chebyshev basis on $[-1, 1]$ brings it to $2.9$. At that point all three methods agree to the last digit, and the normal equations are as good as any.

The general lesson is that accuracy is decided twice: once when the problem is posed, by the condition number, and once when it is solved, by whether the algorithm squares it.