---
title: QR, Householder, and Numerical Least Squares
module: Numerical Linear Algebra
moduleNumber: 8
lessonNumber: 5
order: 805
summary: >
  The least-squares problem reduces to the normal equations, but forming AᵀA
  squares the condition number and can wreck accuracy. The stable route computes
  a QR factorization directly on A and solves Rx = Qᵀb. Householder reflectors
  build that QR one column at a time using length-preserving reflections, the
  unconditionally backward-stable building block behind every serious
  least-squares solver.
topics: [Numerical Linear Algebra]
sources:
  - book: Bornemann
    ref: "§9 QR Decomposition, §16 Normal Equation, §17 Orthogonalization"
  - book: Bornemann
    ref: "Appendix D — The Householder Method for QR Decomposition"
draft: false
---

An overdetermined system $Ax = 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 $Ax$ 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 = Ap + e$, where $A \in
\mathbb{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 - Ax$, the least-squares estimator minimizes its length.

> **Definition (Least-squares problem).** For $A \in \mathbb{R}^{m\times n}$ with
> full column rank and $b \in \mathbb{R}^m$,
>
> $$
> x = \arg\min_{y \in \mathbb{R}^n} \lVert b - Ay\rVert_2.
> $$

Minimizing $F(y) = \tfrac{1}{2}\lVert b - Ay\rVert_2^2$ sets its gradient
$\nabla F(y) = A^\top A y - A^\top b$ to zero. This is the normal equation.

> **Theorem (Normal equation).** The least-squares problem is equivalent to the
> uniquely solvable
>
> $$
> A^\top A\,x = A^\top b,
> $$
>
> where $A^\top A$ is s.p.d. because $A$ has full column rank. Geometrically, the
> residual $r = b - Ax$ is orthogonal to the column space of $A$.[^bor-ne]

The geometry carries the method: the best approximation $Ax$ is the
[orthogonal
projection](/linear-algebra/orthogonality-least-squares/orthogonal-sets-and-projections)
of $b$ onto $\mathrm{Col}\,A$, and what makes it best is that the leftover $r$ is
perpendicular to that subspace, which is exactly $A^\top r = 0$.

$$
% caption: The least-squares solution projects $b$ onto the column space of $A$;
% the residual $r = b - Ax$ is orthogonal to that plane, which is the content of
% the normal equation $A^\top r = 0$.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
% the plane Col A
\draw[black, fill=acc!8] (-2.2,-0.7) -- (2.6,-0.7) -- (3.4,0.7) -- (-1.4,0.7) -- cycle;
\node[anchor=west] at (2.5,0.35) {$\mathrm{Col}\,A$};
% projection point
\coordinate (O) at (0,0);
\coordinate (Ax) at (1.4,0.1);
\coordinate (b) at (1.7,2.0);
\fill[black] (O) circle (1.2pt);
\draw[->, acc, thick] (O) -- (b) node[anchor=south] {$b$};
\draw[->, acc, thick] (O) -- (Ax) node[anchor=north west] {$Ax$};
\draw[->, black, thick] (Ax) -- (b) node[midway, anchor=west] {$r = b - Ax$};
% right angle
\draw[black] ($(Ax)+(-0.18,0.02)$) -- ($(Ax)+(-0.14,0.22)$) -- ($(Ax)+(0.04,0.24)$);
\node[anchor=north] at (1.5,-0.75) {$r \perp \mathrm{Col}\,A$};
\end{tikzpicture}
$$

The direct algorithm forms $A^\top A$, then solves it by
[Cholesky](/linear-algebra/numerical-linear-algebra/lu-and-cholesky), at a cost of
$mn^2 + \tfrac{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 $\kappa_2(A) =
\sigma_1/\sigma_n$, the ratio of extreme singular values of $A$. Forming $A^\top A$
squares it:

$$
\kappa_2(A^\top A) = \kappa_2(A)^2.
$$

Solving the normal equation therefore inherits a condition number of
$\kappa_2(A)^2$, even though the underlying least-squares problem is only
$\kappa_2(A)$-conditioned when the residual is small. By the [rule of
thumb](/linear-algebra/numerical-linear-algebra/stability-and-error-analysis),
squaring the condition number doubles the number of significant digits lost. If
$\kappa_2(A) \approx 10^8$ in double precision — a middling value — then
$\kappa_2(A^\top A) \approx 10^{16}$, and the normal equations can lose _every_
digit while the problem itself still permits eight. The information destroyed in
computing $A^\top A$ cannot be recovered afterward.

> **Worked example.** Let
>
> $$
> A = \begin{bmatrix} 1 & 0 \\ 0 & 10^{-4} \\ 0 & 0 \end{bmatrix},
> $$
>
> whose nonzero singular values are $\sigma_1 = 1$ and $\sigma_2 = 10^{-4}$, so
> $\kappa_2(A) = 10^{4}$. The normal-equation matrix is
>
> $$
> A^\top A = \begin{bmatrix} 1 & 0 \\ 0 & 10^{-8} \end{bmatrix},
> \qquad
> \kappa_2(A^\top A) = 10^{8} = \kappa_2(A)^2.
> $$
>
> In double precision the QR route solves $Rx = Q^\top b$ at the true conditioning
> $10^4$, about twelve correct digits, while the normal equations inherit $10^8$,
> about eight, though both describe the same fit.

$$
% caption: Forming $A^\top A$ squares the condition number: a moderately
% conditioned $A$ becomes a severely conditioned normal-equation matrix, doubling
% the digits lost.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\draw[->, black] (0,0) -- (9,0);
\foreach \x/\l in {0/1, 2/{10^4}, 4/{10^8}, 6/{10^{12}}, 8/{10^{16}}}
  {\draw[black] (\x,-0.1)--(\x,0.1); \node[anchor=north] at (\x,-0.15) {$\l$};}
\node[anchor=south] at (4.5,0.15) {cond};
% A
\fill[acc] (4,0.6) circle (2.5pt);
\node[acc, anchor=south] at (4,0.75) {$\mathrm{cond}(A) = 10^8$};
% A'A
\fill[black] (8,0.6) circle (2.5pt);
\node[anchor=south] at (8,0.75) {$\mathrm{cond}(A^\top A) = 10^{16}$};
\draw[->, black, thick] (4.1,0.55) .. controls (6,1.3) .. (7.9,0.65);
\node[anchor=south] at (6,1.2) {square};
\end{tikzpicture}
$$

## Orthogonalization: the stable route

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

> **Definition (QR factorization).** For $A \in \mathbb{K}^{m\times n}$ with full
> column rank, a **QR factorization** is $A = QR$ with $Q \in \mathbb{K}^{m\times
> n}$ column-orthonormal ($Q^\ast Q = I$) and $R \in \mathrm{GL}(n;\mathbb{K})$
> upper triangular; it is **normalized** when $R$ has a positive diagonal. The
> columns of $Q$ are an orthonormal basis for $\mathrm{Col}\,A$.

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

$$
R x = Q^\top 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 $\lVert Q\rVert_2\lVert R\rVert_2 = \lVert R\rVert_2 = \lVert A\rVert_2$.
The factors never grow, so the benign case of the [previous
lesson](/linear-algebra/numerical-linear-algebra/stability-and-error-analysis) holds
unconditionally.[^bor-orth]

> **Backward stability.** Solving least squares via a QR factorization (by
> Householder, Givens, or modified Gram–Schmidt) is backward stable, because
> $\lVert Q\rVert_2\lVert R\rVert_2 = \lVert A\rVert_2$ puts it in the benign
> regime with no condition on $A$.

There is even a $Q$-free variant: factor the augmented matrix $[A\ \ b]$, and the
extra column delivers $z = Q^\top b$ and the residual norm $\rho = \lVert b -
Ax\rVert_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.

$$
% caption: The stable least-squares pipeline: factor $A = QR$ directly, form $Q^\top b$,
% and back-substitute $Rx = Q^\top b$ — never touching the ill-conditioned $A^\top A$.
\begin{tikzpicture}[scale=1.0, font=\footnotesize,
  box/.style={draw, minimum width=22mm, minimum height=10mm, align=center}]
\definecolor{acc}{HTML}{4A6FA5}
\node[box] (a) at (0,0) {$A$};
\node[box, draw=acc, text=acc] (qr) at (3.2,0) {$A = QR$};
\node[box] (qb) at (6.6,0) {$Q^\top b$};
\node[box, draw=acc, text=acc] (x) at (9.9,0) {$Rx = Q^\top b$};
\draw[->, acc, thick] (a) -- (qr) node[midway, above] {factor};
\draw[->, acc, thick] (qr) -- (qb) node[midway, above] {apply $Q^\top$};
\draw[->, acc, thick] (qb) -- (x) node[midway, above] {back-sub};
\end{tikzpicture}
$$

## 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$.[^bor-qr]

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

$$
Q = I - \frac{2}{v^\top v}\,v v^\top.
$$

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.

$$
% caption: A Householder reflector across the hyperplane $\perp v$ sends the
% column $a$ to $\rho e_1$, zeroing every entry below the first while preserving
% $\lVert a\rVert_2 = \lvert\rho\rvert$.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
% hyperplane (a line here)
\draw[black, thick] (-2.4,-1.2) -- (2.4,1.2);
\node[black, anchor=west] at (2.4,1.2) {hyperplane $\perp v$};
% a and its reflection
\coordinate (O) at (0,0);
\draw[->, acc, thick] (O) -- (1.2,1.8) node[anchor=south] {$a$};
\draw[->, acc, thick] (O) -- (2.15,-0.05) node[anchor=north west] {$\lVert a\rVert_2\, e_1$};
% v direction (normal)
\draw[->, black, dashed] (O) -- (-0.7,1.4) node[anchor=south east] {$v$};
% equal-length arc
\draw[black] (1.2,1.8) arc (56:-1:2.16);
\node[anchor=west] at (2.3,0.9) {length preserved};
\end{tikzpicture}
$$

To zero the entries of a column $a$ below the first, choose $v$ so that
$Q a = \rho e_1$. Because a reflection is an isometry, $\lvert\rho\rvert = \lVert
a\rVert_2$, and the reflecting vector is

$$
\rho = -\operatorname{sign}(\alpha_1)\lVert a\rVert_2,
\qquad
v = a - \rho e_1.
$$

The sign of $\rho$ is chosen _opposite_ to the first entry $\alpha_1$ on purpose:
it makes $\omega_1 = \alpha_1 - \rho = \operatorname{sign}(\alpha_1)(\lvert\alpha_1
\rvert + \lvert\rho\rvert)$ 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](/linear-algebra/numerical-linear-algebra/stability-and-error-analysis)
applied to the reflector's own construction.

```algorithm
caption: $\textsc{HouseholderQR}(A)$ — full QR of $A \in \mathbb{R}^{m\times n}$
for $k = 1$ to $n$ do
  $a \gets A_{k:m,\,k}$ // column below the diagonal
  $\rho \gets -\operatorname{sign}(a_1)\,\lVert a\rVert_2$
  $v \gets a - \rho\, e_1$ // reflector vector, cancellation-free
  $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
  store $v$ // Q kept implicitly as its reflectors
end for
return $R = \text{upper part of } A,\; \{v_k\}$
```

> **Worked example.** Compute the Householder QR of
>
> $$
> A = \begin{bmatrix} 3 & 0 \\ 4 & 5 \end{bmatrix}.
> $$
>
> The first column is $a = (3, 4)^\top$ with $\lVert a\rVert_2 = 5$. Choosing the
> sign opposite to $\alpha_1 = 3$ gives $\rho = -5$ and reflector vector $v = a -
> \rho e_1 = (8, 4)^\top$, with $v^\top v = 80$. The reflector is
>
> $$
> Q = I - \frac{2}{80}\, v v^\top = \begin{bmatrix} -0.6 & -0.8 \\ -0.8 & 0.6 \end{bmatrix},
> $$
>
> symmetric and orthogonal. Applying it,
>
> $$
> R = QA = \begin{bmatrix} -5 & -4 \\ 0 & 3 \end{bmatrix},
> $$
>
> upper triangular, with the first column reduced to $\rho e_1 = (-5, 0)^\top$ and
> its length $\lVert a\rVert_2 = 5 = \lvert\rho\rvert$ preserved. Since $Q^2 = I$,
> the factorization is $A = QR$.

The reflector is applied without ever forming $Q$: each application is a rank-one
update $Qx = x - (w^\top x)\,w$ with $w = v/\sqrt{\lvert\rho\rvert(\lvert\rho\rvert +
\lvert\alpha_1\rvert)}$, costing $O(m)$ flops. Storing the vectors $v_k$ instead
of the matrix keeps $Q$ implicit, and products $Qx$ or $Q^\top y$ cost only $O(mn)$.
A full Householder QR of $A$ costs

$$
\#\text{flop} \approx 2mn^2 - \frac{2}{3}n^3,
$$

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

## Choosing a least-squares method

| Method | Solve | Cost ($m \gg n$) | Conditioning | Use when |
| --- | --- | --- | --- | --- |
| Normal equations | $A^\top A x = A^\top b$ (Cholesky) | $mn^2 + \tfrac{1}{3}n^3$ | $\kappa_2(A)^2$ | $m \gg n$, well-conditioned $A$, large residual |
| Orthogonalization | $Rx = Q^\top b$ (Householder QR) | $2mn^2 - \tfrac{2}{3}n^3$ | $\kappa_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](/linear-algebra/numerical-linear-algebra/numerical-eigenvalues-and-svd).

[^bor-ne]: **Bornemann**, _Numerical Linear Algebra_, §16 — Normal Equation: the least-squares setup $b = Ap + e$, the normal equation $A^\top A x = A^\top b$ with $A^\top A$ s.p.d., its solution via Cholesky, and the condition relation $\kappa_2(A) = \sqrt{\kappa_2(A^\top A)}$.
[^bor-orth]: **Bornemann**, §17 — Orthogonalization: the stable reduction $Rx = Q^\top b$, backward stability from $\lVert Q\rVert_2\lVert R\rVert_2 = \lVert A\rVert_2$, the $Q$-free augmented-matrix variant, and the cost comparison with the normal equations.
[^bor-qr]: **Bornemann**, §9 — QR Decomposition (modified Gram–Schmidt, Givens rotations) and Appendix D — The Householder Method: the reflector $Q = I - 2vv^\top/(v^\top v)$, the cancellation-free sign choice $\rho = -\operatorname{sign}(\alpha_1)\lVert a\rVert_2$, and the $2mn^2 - \tfrac{2}{3}n^3$ cost.