QR, Householder, and Numerical Least Squares

An overdetermined system $A x = b$ , with more equations than unknowns, usually has no exact solution: the data $b$ does not lie in the column space of $A$ . The least-squares problem asks for the next best thing, the $x$ making $A x$ as close to $b$ as possible. The classical answer is the normal equations, which are numerically the wrong thing to compute directly. The stable alternative is orthogonalization, built from the Householder reflector.

The least-squares problem and normal equations

Fitting a linear model to noisy measurements gives $b = A p + e$ , where $A \in R^{m \times n}$ is the design matrix ( $m$ measurements, $n$ parameters, $m > n$ ), $p$ the parameters, and $e$ unknown noise. Replacing the noise by a computable residual $r = b - A x$ , the least-squares estimator minimizes its length.

Minimizing $F (y) = \frac{1}{2} ∥ b - A y ∥_{2}^{2}$ sets its gradient $\nabla F (y) = A^{⊤} A y - A^{⊤} b$ to zero. This is the normal equation.

The geometry carries the method: the best approximation $A x$ is the orthogonal projection of $b$ onto $Col A$ , and what makes it best is that the leftover $r$ is perpendicular to that subspace, which is exactly $A^{⊤} r = 0$ .

The least-squares solution projects

b

onto the column space of

A

; the residual

r = b - A x

is orthogonal to that plane, which is the content of the normal equation

A^{⊤} r = 0

The direct algorithm forms $A^{⊤} A$ , then solves it by Cholesky, at a cost of $m n^{2} + \frac{1}{3} n^{3}$ flops. It is cheap, and for anything but a well-conditioned $A$ it is inaccurate.

Conditioning of the normal equations

The condition number of the least-squares problem grows with $κ_{2} (A) = σ_{1} / σ_{n}$ , the ratio of extreme singular values of $A$ . Forming $A^{⊤} A$ squares it:

κ_{2} (A^{⊤} A) = κ_{2} (A)^{2} .

Solving the normal equation therefore inherits a condition number of $κ_{2} (A)^{2}$ , even though the underlying least-squares problem is only $κ_{2} (A)$ -conditioned when the residual is small. By the rule of thumb, squaring the condition number doubles the number of significant digits lost. If $κ_{2} (A) \approx 1 0^{8}$ in double precision — a middling value — then $κ_{2} (A^{⊤} A) \approx 1 0^{16}$ , and the normal equations can lose every digit while the problem itself still permits eight. The information destroyed in computing $A^{⊤} A$ cannot be recovered afterward.

Forming

A^{⊤} A

squares the condition number: a moderately conditioned

A

becomes a severely conditioned normal-equation matrix, doubling the digits lost.

Orthogonalization: the stable route

The fix is Stiefel's principle of direct attack: operate on $A$ itself, never on $A^{⊤} A$ . A QR factorization supplies it.

Substituting $A = QR$ into the normal equation collapses it. Since $A^{⊤} A = R^{⊤} Q^{⊤} QR = R^{⊤} R$ and $A^{⊤} b = R^{⊤} Q^{⊤} b$ , cancelling the invertible $R^{⊤}$ leaves a triangular system:

R x = Q^{⊤} b .

This is solved by one back substitution. Its stability follows from the unitary factor: reading the solve as factorization followed by back substitution, the inverse of the factorization step $(Q, R) \mapsto QR = A$ is well-conditioned because $∥ Q ∥_{2} ∥ R ∥_{2} = ∥ R ∥_{2} = ∥ A ∥_{2}$ . The factors never grow, so the benign case of the previous lesson holds unconditionally.²

There is even a $Q$ -free variant: factor the augmented matrix $[A b]$ , and the extra column delivers $z = Q^{⊤} b$ and the residual norm $ρ = ∥ b - A x ∥_{2}$ directly, so $Q$ never needs to be formed. In MATLAB the entire least-squares solve is the single backslash x = A\b, which runs exactly this orthogonalization.

The stable least-squares pipeline: factor

A = QR

directly, form

Q^{⊤} b

, and back-substitute

R x = Q^{⊤} b

— never touching the ill-conditioned

A^{⊤} A

Householder reflectors

Three algorithms build the QR factorization: modified Gram–Schmidt (project each column against the earlier orthonormal ones), Givens rotations (zero one subdiagonal entry at a time), and Householder reflectors (zero a whole column below the diagonal at once). The last is the standard, used by LAPACK's qr, because it is the most stable and produces the most accurate $Q$ .³

A Householder reflector reflects space across the hyperplane orthogonal to a chosen vector $v$ :

Q = I - \frac{2}{v ^{⊤} v} v v^{⊤} .

It is symmetric, orthogonal, and its own inverse ( $Q^{2} = I$ ), with $det Q = - 1$ . Being unitary, it preserves lengths, which is what makes it numerically safe.

A Householder reflector across the hyperplane

⊥ v

sends the column

a

ρ e_{1}

, zeroing every entry below the first while preserving

∥ a ∥_{2} = ∣ ρ ∣

To zero the entries of a column $a$ below the first, choose $v$ so that $Q a = ρ e_{1}$ . Because a reflection is an isometry, $∣ ρ ∣ = ∥ a ∥_{2}$ , and the reflecting vector is

ρ = - sign (α_{1}) ∥ a ∥_{2}, v = a - ρ e_{1} .

The sign of $ρ$ is chosen opposite to the first entry $α_{1}$ on purpose: it makes $ω_{1} = α_{1} - ρ = sign (α_{1}) (∣ α_{1} ∣ + ∣ ρ ∣)$ an addition of like-signed numbers rather than a subtraction, avoiding the cancellation the opposite sign choice would cause. This is the cancellation-avoidance lesson applied to the reflector's own construction.

Algorithm:

\textsc{HouseholderQR}(A)

— full QR of

A \in \mathbb{R}^{m\times n}

1
for $k = 1$ to $n$ do
2
$a \gets A_{k:m,\,k}$
column below the diagonal
3
$\rho \gets -\operatorname{sign}(a_1)\,\lVert a\rVert_2$
4
$v \gets a - \rho\, e_1$
reflector vector, cancellation-free
5
$A_{k:m,\,k:n} \gets A_{k:m,\,k:n} - \tfrac{2}{v^\top v}\, v\,(v^\top A_{k:m,\,k:n})$
apply reflector
6
store $v$
Q kept implicitly as its reflectors
7
end for
8
return $R = \text{upper part of } A,\; \{v_k\}$

The reflector is applied without ever forming $Q$ : each application is a rank-one update $Q x = x - (w^{⊤} x) w$ with $w = v / ∣ ρ ∣ (∣ ρ ∣ + ∣ α_{1} ∣)$ , costing $O (m)$ flops. Storing the vectors $v_{k}$ instead of the matrix keeps $Q$ implicit, and products $Q x$ or $Q^{⊤} y$ cost only $O (mn)$ . A full Householder QR of $A$ costs

# flop \approx 2 m n^{2} - \frac{2}{3} n^{3},

roughly twice the normal-equation cost when $m ≫ n$ .

Choosing a least-squares method

Method	Solve	Cost ( $m ≫ n$ )	Conditioning	Use when
Normal equations	$A^{⊤} A x = A^{⊤} b$ (Cholesky)	$m n^{2} + \frac{1}{3} n^{3}$	$κ_{2} (A)^{2}$	$m ≫ n$ , well-conditioned $A$ , large residual
Orthogonalization	$R x = Q^{⊤} b$ (Householder QR)	$2 m n^{2} - \frac{2}{3} n^{3}$	$κ_{2} (A)$	ill-conditioned $A$ , small residual — the default

The normal equations keep a niche: when $A$ is well-conditioned and the residual is large, their lower cost is worth the squared condition number. Everywhere else, orthogonalization is the right choice: about twice the flops, with the condition number unsquared and the solve backward stable. Iterated, the same QR factorization becomes the standard algorithm for eigenvalues.

Bornemann, Numerical Linear Algebra, §16 — Normal Equation: the least-squares setup $b = A p + e$ , the normal equation $A^{⊤} A x = A^{⊤} b$ with $A^{⊤} A$ s.p.d., its solution via Cholesky, and the condition relation $κ_{2} (A) = κ_{2} (A^{⊤} A)$ . ↩
Bornemann, §17 — Orthogonalization: the stable reduction $R x = Q^{⊤} b$ , backward stability from $∥ Q ∥_{2} ∥ R ∥_{2} = ∥ A ∥_{2}$ , the $Q$ -free augmented-matrix variant, and the cost comparison with the normal equations. ↩
Bornemann, §9 — QR Decomposition (modified Gram–Schmidt, Givens rotations) and Appendix D — The Householder Method: the reflector $Q = I - 2 v v^{⊤} / (v^{⊤} v)$ , the cancellation-free sign choice $ρ = - sign (α_{1}) ∥ a ∥_{2}$ , and the $2 m n^{2} - \frac{2}{3} n^{3}$ cost. ↩

The least-squares problem and normal equations

Conditioning of the normal equations

Orthogonalization: the stable route

Householder reflectors

Choosing a least-squares method

Footnotes