Householder QR Factorization

Definition Let $v \in \mathbb{R}^m$ be nonzero. The Householder reflector determined by $v$ is the $m \times m$ matrix

The matrix $P = vv^{\mathsf{T}}/(v^{\mathsf{T}}v)$ is the orthogonal projector onto the line spanned by $v$. So $H = I - 2P$ reverses the component of a vector along $v$ and keeps the component perpendicular to $v$. In particular $Hv = -v$, and $Hw = w$ whenever $v^{\mathsf{T}}w = 0$. In $\mathbb{R}^3$ this is the mirror reflection in the plane perpendicular to $v$, and the same formula works in every dimension. The matrix $H$ is symmetric and orthogonal, and it is its own inverse. Indeed $P^{\mathsf{T}} = P$, and $P^2 = P$ because the scalar $v^{\mathsf{T}}v$ in the middle cancels one factor of the denominator,

Hence $H^{\mathsf{T}} = H$ and

Consequently $\lVert Hz \rVert^2 = z^{\mathsf{T}}H^{\mathsf{T}}Hz = z^{\mathsf{T}}z$ for every $z$, inner products are preserved as well, and $H^2 = I$. Multiplying $v$ by a nonzero number does not change $H$. Over $\mathbb{C}$ one uses $I - 2vv^{*}/(v^{*}v)$, which is Hermitian and unitary, with the phase of the first entry in place of the sign chosen below. This lesson stays with real matrices.

Definition Let $a \in \mathbb{R}^m$ with $m \ge 2$, and suppose its tail $(a_2, \dots, a_m)$ is not zero. Let $\alpha = \lVert a \rVert$, let $s = 1$ if $a_1 \ge 0$ and $s = -1$ if $a_1 < 0$, so that the sign of $0$ is taken to be $1$, and set

Then $v \ne 0$, and the reflector $H$ of $v$ satisfies $Ha = -s\alpha\,e_1$. To check this, note that $sa_1 = \lvert a_1 \rvert$, so

Both are positive because $\alpha > 0$, and therefore

The other sign, $v = a - s\alpha\,e_1$, also gives a reflector that sends $a$ to $s\alpha\,e_1$ in exact arithmetic. Its first entry, however, has absolute value $\alpha - \lvert a_1 \rvert$. When $a$ is close to the $e_1$ axis this is a difference of two nearly equal numbers, and such a difference can lose most of its significant digits. The stable choice adds the two nonnegative numbers $\lvert a_1 \rvert$ and $\alpha$, so this step has no cancellation. That explains the design choice. It is not by itself a proof that the whole algorithm is stable, which needs an error analysis of every operation, including the computation of $\alpha$.

When the tail is already zero, including the case $a = 0$, there is nothing to eliminate, and the algorithm uses $H = I$. A reflector is never formed from the zero vector, whose denominator $v^{\mathsf{T}}v$ would be $0$. A zero leading entry with a nonzero tail is not such a case. There the formula applies with $s = 1$.

Definition Let $A \in \mathbb{R}^{m \times n}$ with $m \ge n$. For $k = 1, \dots, n$, look at column $k$ of the current matrix in rows $k$ through $m$. If its entries below row $k$ are all zero, take $\tilde{H}_k = I$. Otherwise take the stable reflector $\tilde{H}_k$ of size $m - k + 1$ that sends this piece of the column to a multiple of its first unit vector. Embed it as

which acts only on rows $k$ through $m$. When $m = n$ the last step is the identity, since a single entry has no tail. The zeros already created in columns $1, \dots, k - 1$ survive, because in those columns rows $k$ through $m$ are zero and $\tilde{H}_k$ sends zero to zero. After the last step

where $R \in \mathbb{R}^{n \times n}$ is upper triangular and the zero block has $m - n$ rows. Each $H_k$ is symmetric and orthogonal, so $Q = H_1 H_2 \cdots H_n$ is orthogonal with $Q^{\mathsf{T}} = H_n \cdots H_2 H_1$, and

This shows that every real $m \times n$ matrix with $m \ge n$ has such a factorization, with no assumption on its rank. If $A$ has full column rank $n$, then $R$ has the same rank, so it is invertible and its diagonal entries are nonzero. They may be negative. Changing the sign of column $j$ of $Q$ together with row $j$ of $R_{\text{full}}$ gives another factorization, so a nonnegative diagonal can always be arranged if a convention requires it. In practice $Q$ is often stored as the list of vectors that define the reflectors rather than multiplied out.

Definition Let $R \in \mathbb{R}^{n \times n}$ be upper triangular with nonzero diagonal entries $r_{11}, \dots, r_{nn}$, and let $c \in \mathbb{R}^n$. The system $Rx = c$ has exactly one solution, found by back substitution from the last equation upward,

for $i = n - 1, \dots, 1$ in turn. It takes about $n^2$ arithmetic operations and never forms $R^{-1}$. A lower triangular system is solved the same way from the top down, by forward substitution. If some $r_{ii} = 0$, the matrix $R$ is singular and the substitution breaks down at that step. For a QR factorization this happens exactly when $A$ does not have full column rank, and a rank-revealing method such as QR with column pivoting or the SVD is then needed. To solve a square invertible system $Ax = b$ with $A = QR$, apply $Q^{\mathsf{T}}$ and solve $Rx = Q^{\mathsf{T}}b$. Least squares uses the same step on the first $n$ coordinates of $Q^{\mathsf{T}}b$.

Definition Let $A \in \mathbb{R}^{m \times n}$ with $m \ge n$. A full QR factorization is $A = QR_{\text{full}}$ with $Q \in \mathbb{R}^{m \times m}$ orthogonal and $R_{\text{full}} \in \mathbb{R}^{m \times n}$ upper triangular in its top $n \times n$ block $R$ and zero in its last $m - n$ rows. Split $Q = [\,Q_1 \ \ Q_2\,]$ into its first $n$ columns $Q_1 \in \mathbb{R}^{m \times n}$ and its last $m - n$ columns $Q_2$. The zero rows of $R_{\text{full}}$ multiply $Q_2$, so the reduced, or thin, QR factorization is

When $m > n$ the matrix $Q_1$ is not orthogonal, since $Q_1Q_1^{\mathsf{T}} \ne I_m$. If $A$ has full column rank, then $R$ is invertible, the columns of $Q_1$ form an orthonormal basis of the column space of $A$, the matrix $Q_1Q_1^{\mathsf{T}}$ is the orthogonal projector onto that column space, and the columns of $Q_2$ form an orthonormal basis of its orthogonal complement. In that case the reduced factorization with positive diagonal entries in $R$ is unique. When $m > n$ the block $Q_2$ is not unique, since any orthonormal basis of the complement will do. When $m = n$ there is no $Q_2$, and the full and reduced factorizations coincide.

Definition Let $A \in \mathbb{R}^{m \times n}$ with $m \ge n$ have full column rank and a full QR factorization $A = QR_{\text{full}}$ with $Q = [\,Q_1 \ \ Q_2\,]$ as above, and let $b \in \mathbb{R}^m$. Because $Q^{\mathsf{T}}$ preserves lengths and $Q^{\mathsf{T}}A = R_{\text{full}}$, for every $x \in \mathbb{R}^n$

The two terms come from the first $n$ and the last $m - n$ coordinates, which are perpendicular groups. The second term does not depend on $x$. The first is zero exactly for $\hat{x} = R^{-1}Q_1^{\mathsf{T}}b$, computed by a triangular solve, so $\hat{x}$ is the unique least-squares solution, and the minimal residual norm is

The residual itself is $r = b - A\hat{x} = b - Q_1Q_1^{\mathsf{T}}b = Q_2Q_2^{\mathsf{T}}b$, the component of $b$ perpendicular to the column space, so $A^{\mathsf{T}}r = 0$. It is zero exactly when $b$ lies in the column space of $A$. Watch for this when only the reduced factorization is kept. The triangular system $R\hat{x} = Q_1^{\mathsf{T}}b$ is then solved exactly, which says nothing about whether $b$ lies in the column space. The omitted component still exists and can be computed as $b - Q_1(Q_1^{\mathsf{T}}b)$.