Least Squares with QR

In Least Squares you found the $\hat{\mathbf{x}}$ that makes $\lVert\mathbf{b} - A\mathbf{x}\rVert_2$ as small as possible by solving the normal equations $A^{T}A\hat{\mathbf{x}} = A^{T}\mathbf{b}$. In exact arithmetic that is fine. In floating point it can throw away half of the available digits, for a reason you met in The Condition Number of a Matrix.

Factor the $m \times n$ matrix $A$ (with $m > n$ and independent columns) as $A = QR$, as in Gram–Schmidt and QR: $Q$ is $m \times n$ with orthonormal columns ($Q^{T}Q = I$) and $R$ is $n \times n$, upper triangular and invertible.

Key idea

The least-squares solution solves the triangular system $$R\hat{\mathbf{x}} = Q^{T}\mathbf{b},$$ by back substitution. This avoids $A^{T}A$, whose condition number is $\kappa_2(A^{T}A) = \kappa_2(A)^{2}$. With QR the loss of digits is governed by $\kappa_2(A)$ (for problems that fit the data well); with the normal equations it is governed by $\kappa_2(A)^{2}$.

Why it works

Substitute $A = QR$ into the normal equations: $A^{T}A = R^{T}Q^{T}QR = R^{T}R$ and $A^{T}\mathbf{b} = R^{T}Q^{T}\mathbf{b}$, so $R^{T}R\hat{\mathbf{x}} = R^{T}Q^{T}\mathbf{b}$. Since $R^{T}$ is invertible, cancel it: $R\hat{\mathbf{x}} = Q^{T}\mathbf{b}$. For the conditioning, the singular values of $A^{T}A$ are the squares $\sigma_i^{2}$ of those of $A$, so $\kappa_2(A^{T}A) = \frac{\sigma_1^{2}}{\sigma_n^{2}} = \kappa_2(A)^{2}$. Multiplying by $Q^{T}$ never changes lengths of vectors in the column space, so it adds no amplification of its own.

Worked example: A line through four points

Fit $y = c_0 + c_1t$ to $(0, 1)$, $(1, 2)$, $(2, 2)$, $(3, 4)$ using QR.

$A$ has columns $(1, 1, 1, 1)$ and $(0, 1, 2, 3)$. Gram–Schmidt: $\mathbf{q}_1 = \frac{1}{2}(1, 1, 1, 1)$, $r_{11} = 2$, $r_{12} = \mathbf{q}_1 \cdot (0, 1, 2, 3) = 3$.
Subtract: $(0, 1, 2, 3) - 3\mathbf{q}_1 = (-1.5, -0.5, 0.5, 1.5)$, of length $\sqrt{5}$, so $r_{22} = \sqrt{5}$ and $\mathbf{q}_2 = \frac{1}{2\sqrt{5}}(-3, -1, 1, 3)$.
$Q^{T}\mathbf{b} = \left(\frac{1 + 2 + 2 + 4}{2},\ \frac{-3 - 2 + 2 + 12}{2\sqrt{5}}\right) = \left(4.5,\ \frac{9}{2\sqrt{5}}\right)$.
Back substitution: $\sqrt{5}\,c_1 = \frac{9}{2\sqrt{5}}$ gives $c_1 = 0.9$; then $2c_0 + 3(0.9) = 4.5$ gives $c_0 = 0.9$.

Answer: $y = 0.9 + 0.9t$

Worked example: Where the normal equations break

Let $\delta = 10^{-8}$, $A = \begin{bmatrix} 1 & 1 \\ \delta & 0 \\ 0 & \delta \end{bmatrix}$ and $\mathbf{b} = (2, \delta, \delta)$, so the exact solution is $(1, 1)$. Compare the two methods in double precision.

$A^{T}A = \begin{bmatrix} 1 + \delta^{2} & 1 \\ 1 & 1 + \delta^{2} \end{bmatrix}$, and $1 + 10^{-16}$ rounds to $1$ (it is less than half a gap above $1$).
So the computed $A^{T}A$ is $\begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}$, which is singular: the normal equations have no unique solution at all.
QR computed with Householder reflections (the way libraries do it) returns $(1, 1)$ to full accuracy here, even though $\kappa_2(A) \approx 1.4 \times 10^{8}$.

Answer: normal equations fail; QR succeeds

Worked example: Counting digits

A least-squares problem has $\kappa_2(A) = 10^{5}$. About how many correct digits can each method give in double precision (about 16 digits)?

QR: lose up to about $5$ digits, leaving about $11$.
Normal equations: $\kappa_2(A^{T}A) = 10^{10}$, so lose up to about $10$, leaving about $6$.

Answer: about $11$ versus about $6$

Fitting polynomials to 30 points in double precision. As the degree grows, κ(A) grows fast. The QR error (green) tracks κ(A)·ε; the normal-equations error (orange) tracks κ(A)²·ε.

The normal equations are not always wrong: they cost about half as much as QR, and when $A$ is well conditioned (say $\kappa_2(A) < 10^{3}$) the lost digits do not matter. But a library cannot know in advance, so it uses QR (or the SVD for the hardest cases). One more caution: classical Gram–Schmidt, done naively in floating point, can lose orthogonality in $Q$; Householder reflections do not.

Common mistake

Forming $A^{T}A$ for every least-squares problem. Since $\kappa(A^{T}A) = \kappa(A)^{2}$, a problem that QR solves to 8 digits can come out with none from the normal equations. Check: estimate $\kappa(A)$; if its square is anywhere near $10^{16}$, do not form $A^{T}A$.

Least Squares with QR

Practice

Builds on

Leads to

Sources and further reading