Computing QR factorizations

It is possible to compute a thin QR factorization using the outer product formula (2.4.4), as we did with LU. However, to stably compute the factorization, a better strategy is to introduce zeros into the lower triangle, one column at a time, using orthogonal matrices. Thanks to Theorem 3.3.2, the product of orthogonal matrices will also be orthogonal.

3.4.1Householder reflections¶

We will use a particular type of orthogonal matrix.

In Exercise 3.4.2 you are asked to show that such a $\mathbf{P}$ is necessarily orthogonal. Note that for any vector $\mathbf{x}$ of appropriate dimension,

The reason $\mathbf{P}$ is called a reflector is sketched in Figure 3.4.1.

Given a vector $\mathbf{z}$, we can choose $\mathbf{v}$ so that $\mathbf{P}$ reflects $\mathbf{z}$ onto the $x_1$-axis—i.e., so that $\mathbf{P}\mathbf{z}$ is nonzero only in the first element. Because orthogonal matrices preserve the 2-norm, we must have

(Recall that $\mathbf{e}_k$ is the $k$th column of the identity matrix.) We choose the positive sign above for our discussion, but see Function 3.4.1 and Exercise 3.4.4 for important computational details. Let

If it turns out that $\mathbf{w}=\boldsymbol{0}$, then $\mathbf{z}$ is already in the target form and we can take $\mathbf{P}=\mathbf{I}$. Otherwise, we have the following.

3.4.2Factorization algorithm¶

The QR factorization is computed by using successive Householder reflections to introduce zeros in one column at a time. We first show the process for a small numerical example in Example 3.4.1.

You may be wondering what happened to $\mathbf{Q}$ in Example 3.4.1. Each Householder reflector is orthogonal but not full-size. We have to pad it out to represent algebraically the fact that a block of the first rows is left alone. Given a reflector $\mathbf{P}_k$ that is of square size $m-k+1$, we define

It is easy to show that $\mathbf{Q}_k$ is also orthogonal. Then the algorithm produces

But $\mathbf{Q}_n \mathbf{Q}_{n-1}\cdots \mathbf{Q}_1$ is orthogonal too, and we multiply on the left by its transpose to get $\mathbf{A}=\mathbf{Q}\mathbf{R}$, where $\mathbf{Q} = (\mathbf{Q}_n \mathbf{Q}_{n-1}\cdots \mathbf{Q}_1)^T$. We don’t even need to form these matrices explicitly. Writing

we can build $\mathbf{Q}^T$ iteratively by starting with the identity and doing the same row operations as on $\mathbf{A}$. That process uses much less memory than building the $\mathbf{Q}_k$ matrices explicitly.

The algorithm we have described is encapsulated in Function 3.4.1. There is one more refinement in it, however. As indicated by (3.4.2), the application of a reflector $\mathbf{P}$ to a vector does not require the formation of the matrix $\mathbf{P}$ explicitly.

3.4.3Q-less QR and least squares¶

In Example 3.3.1, it was seen that the $\mathbf{Q}$ output of Julia’s qr function is not a standard matrix. The reason is that Equation (3.3.8) shows that in order to solve the linear least-squares problem, all we need from $\mathbf{Q}$ is the computation of $\hat{\mathbf{Q}}^T\mathbf{b}$. Referring again to (3.4.10) and (3.4.2), the special structure of the reflectors is such that for this computation, we only need to apply code similar to lines 18 and 21 of Function 3.4.1 for each of the Householder vectors $\mathbf{v}$ that is constructed.

This observation leads to the idea of the Q-less QR factorization, in which the full or thin $\mathbf{Q}$ is never computed explicitly. This is the variant used by Julia’s qr. The returned value Q used within Function 3.3.2 is of a special type that allows Julia to perform Q'*b efficiently for the least-squares solution.

In Exercise 3.4.8 you are asked to derive the following result about the Q-less factorization.

The flop count quoted in Theorem 3.4.2 dominates the running time for least-squares solution via QR. Compared to the count from Theorem 3.2.3 for solution by the normal equations, the flops are essentially identical when $m=n$, but the QR solution is about twice the cost when $m\gg n$. The redeeming quality of the QR route is better stability, which we do not discuss here.

3.4.4Exercises¶