Computing QR factorizations

3.4. Computing QR factorizations#

It is possible to compute a thin QR factorization using the outer product formula (2.4.2), 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.6, the product of orthogonal matrices will also be orthogonal.

Householder reflections#

We will use a particular type of orthogonal matrix.

Definition 3.4.1 :  Householder reflector

A Householder reflector is a matrix of the form

(3.4.1)#\[ \mathbf{P} = \mathbf{I} - 2 \mathbf{v} \mathbf{v}^T,\]

where \(\mathbf{v}\) is any unit vector (in the 2-norm).

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

(3.4.2)#\[ \mathbf{P}\mathbf{x} = \mathbf{x} - 2 \mathbf{v} (\mathbf{v}^T\mathbf{x}).\]

The reason \(\mathbf{P}\) is called a reflector is sketched in Fig. 3.4.1.

../_images/hhreflect.svg

Fig. 3.4.1 A Householder reflector. Because \(\mathbf{v}\) is a unit vector, \(\mathbf{v}^T\mathbf{x}\) is the component of \(\mathbf{x}\) in the direction of \(\mathbf{v}\). Hence subtracting \((\mathbf{v}^T\mathbf{x})\mathbf{v}\) projects \(\mathbf{x}\) into a hyperplane orthogonal to \(\mathbf{v}\). By subtracting off twice as much, we get the reflection of \(\mathbf{x}\) through the hyperplane instead.#

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

(3.4.3)#\[\begin{split}\mathbf{P}\mathbf{z} = \begin{bmatrix} \pm \| \mathbf{z} \|\\0 \\ \vdots \\ 0 \end{bmatrix} = \pm \| \mathbf{z} \| \mathbf{e}_1.\end{split}\]

(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.4 and Exercise 4 for important computational details. Let

(3.4.4)#\[ \mathbf{w} = \| \mathbf{z} \| \mathbf{e}_1-\mathbf{z}, \quad \mathbf{v} = \frac{\mathbf{w}}{\|\mathbf{w}\|}.\]

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.

Theorem 3.4.2 :  Householder reflector

Let \(\mathbf{v}\) be defined by (3.4.4) and let \(\mathbf{P}\) be given by (3.4.1). Then \(\mathbf{P}\) is symmetric and orthogonal, and \(\mathbf{P}\mathbf{z}=\| \mathbf{z} \|\mathbf{e}_1\).

Proof

The proofs of symmetry and orthogonality are left to Exercise 2. For the last fact, we use (3.4.2) to compute

\[\mathbf{P}\mathbf{z} = \mathbf{z} - 2 \frac{\mathbf{w}^T \mathbf{z}}{\mathbf{w}^T\mathbf{w}} \mathbf{w}.\]

Since \(\mathbf{e}_1^T\mathbf{z}=z_1\),

\[\begin{split}\begin{split} \mathbf{w}^T\mathbf{w} &= \| \mathbf{z} \|^2 - 2 \| \mathbf{z} \| z_1 + \mathbf{z}^T\mathbf{z} = 2\| \mathbf{z} \|(\| \mathbf{z} \|-z_1),\\ \mathbf{w}^T\mathbf{z} &= \| \mathbf{z} \|z_1 - \mathbf{z}^T\mathbf{z} = -\| \mathbf{z} \|\bigl(\| \mathbf{z} \|-z_1\bigr), \end{split}\end{split}\]

leading finally to

\[\mathbf{P}\mathbf{z} = \mathbf{z} - 2\cdot \frac{-\| \mathbf{z} \| \bigl(\| \mathbf{z} \|-z_1\bigr)}{2\| \mathbf{z} \| \bigl(\| \mathbf{z} \|-z_1\bigr)} \mathbf{w} = \mathbf{z} + \mathbf{w} = \| \mathbf{z} \|\mathbf{e}_1.\]

Factorization 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 Demo 3.4.3.

Demo 3.4.3

We will use Householder reflections to produce a QR factorization of a random matrix.

The rand function can select randomly from within the interval \([0,1]\), or from a vector or range that you specify.

A = rand(float(1:9),6,4)
m,n = size(A)

(6, 4)

Our first step is to introduce zeros below the diagonal in column 1 by using (3.4.4) and (3.4.1).

I can stand for an identity matrix of any size, inferred from the context when needed.

z = A[:,1];
v = normalize(z - norm(z)*[1;zeros(m-1)])
P₁ = I - 2v*v'   # reflector

6×6 Matrix{Float64}:
 0.367194   0.550791   0.244796    0.428393    0.122398    0.550791
 0.550791   0.520594  -0.213069   -0.372871   -0.106535   -0.479406
 0.244796  -0.213069   0.905303   -0.16572    -0.0473487  -0.213069
 0.428393  -0.372871  -0.16572     0.709989   -0.0828602  -0.372871
 0.122398  -0.106535  -0.0473487  -0.0828602   0.976326   -0.106535
 0.550791  -0.479406  -0.213069   -0.372871   -0.106535    0.520594

We check that this reflector introduces zeros as it should:

P₁*z

6-element Vector{Float64}:
 16.340134638368195
 -1.942890293094024e-15
 -3.0531133177191805e-16
 -1.1102230246251565e-15
 -2.636779683484747e-16
 -2.4424906541753444e-15

Now we replace \(\mathbf{A}\) by \(\mathbf{P}\mathbf{A}\).

A = P₁*A

6×4 Matrix{Float64}:
 16.3401       13.219       9.48585    9.05745
 -1.94289e-15  -0.0241574  -6.38604   -3.40198
 -3.05311e-16  -0.566292   -0.282683   4.04356
 -1.11022e-15   0.758989    3.2553     5.57624
 -2.63678e-16   5.21685    -0.641341   4.02178
 -2.44249e-15  -3.02416    -0.386037  -2.40198

We are set to put zeros into column 2. We must not use row 1 in any way, lest it destroy the zeros we just introduced. So we leave it out of the next reflector.

z = A[2:m,2]
v = normalize(z - norm(z)*[1;zeros(m-2)])
P₂ = I - 2v*v'

5×5 Matrix{Float64}:
 -0.00395766  -0.0927744   0.124343    0.854666   -0.495441
 -0.0927744    0.991427    0.0114904   0.0789786  -0.0457831
  0.124343     0.0114904   0.9846     -0.105853    0.061362
  0.854666     0.0789786  -0.105853    0.272426    0.421768
 -0.495441    -0.0457831   0.061362    0.421768    0.755506

We now apply this reflector to rows 2 and below only.

A[2:m,:] = P₂*A[2:m,:]
A

6×4 Matrix{Float64}:
 16.3401       13.219         9.48585     9.05745
  8.82719e-16   6.10397       0.0994012   4.95901
 -4.42004e-17  -1.00293e-16   0.316628    4.81619
 -1.46018e-15   7.66754e-17   2.45206     4.5407
 -2.66911e-15   1.11068e-15  -6.16237    -3.09591
 -1.04809e-15  -5.61114e-16   2.81445     1.72407

We need to iterate the process for the last two columns.

for j in 3:n
    z = A[j:m,j]
    v = normalize(z - norm(z)*[1;zeros(m-j)])
    P = I - 2v*v'
    A[j:m,:] = P*A[j:m,:]
end

We have now reduced the original to an upper triangular matrix using four orthogonal Householder reflections:

R = triu(A)

6×4 Matrix{Float64}:
 16.3401  13.219    9.48585    9.05745
  0.0      6.10397  0.0994012  4.95901
  0.0      0.0      7.21172    5.07361
  0.0      0.0      0.0        5.53436
  0.0      0.0      0.0        0.0
  0.0      0.0      0.0        0.0

You may be wondering what happened to \(\mathbf{Q}\) in Demo 3.4.3. 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

\[\begin{split}\mathbf{Q}_k = \begin{bmatrix} \mathbf{I}_{k-1} & \boldsymbol{0} \\ \boldsymbol{0} & \mathbf{P}_k \end{bmatrix}.\end{split}\]

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

\[ \mathbf{Q}_n \mathbf{Q}_{n-1}\cdots \mathbf{Q}_1 \mathbf{A} = \mathbf{R}.\]

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

(3.4.5)#\[\mathbf{Q}^T = \mathbf{Q}_n \mathbf{Q}_{n-1}\cdots \mathbf{Q}_1 = \mathbf{Q}_n \Bigl( \mathbf{Q}_{n-1}\bigl(\cdots (\mathbf{Q}_1\mathbf{I})\cdots\bigr)\Bigr),\]

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.4. 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.

Function 3.4.4 :  qrfact

QR factorization by Householder reflections

 1"""
 2    qrfact(A)
 3
 4QR factorization by Householder reflections. Returns Q and R.
 5"""
 6function qrfact(A)
 7    m,n = size(A)
 8    Qt = diagm(ones(m))
 9    R = float(copy(A))
10    for k in 1:n
11        z = R[k:m,k]
12        w = [ -sign(z[1])*norm(z) - z[1]; -z[2:end] ]
13        nrmw = norm(w)
14        if nrmw < eps() continue; end    # skip this iteration
15        v = w / nrmw;
16        # Apply the reflection to each relevant column of A and Q
17        for j in 1:n
18            R[k:m,j] -= v*( 2*(v'*R[k:m,j]) )
19        end
20        for j in 1:m
21            Qt[k:m,j] -= v*( 2*(v'*Qt[k:m,j]) )
22        end
23    end
24    return Qt',triu(R)
25end

Q-less QR and least squares#

In Demo 3.3.9 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.4) 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.5) 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.4 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.11 is of a special type that allows Julia to perform Q'*b efficiently for the least-squares solution.

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

Theorem 3.4.5

Q-less QR factorization by Householder reflections takes \(\sim(2mn^2-\frac{2}{3}n^3)\) flops.

The flop count quoted in Theorem 3.4.5 dominates the running time for least-squares solution via QR. Compared to the count from Theorem 3.2.7 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.

Exercises#

  1. ⌨ Find a Householder reflector \(\mathbf{P}\) such that

    \[\begin{split}\mathbf{P} \begin{bmatrix} 2 \\ 9 \\ -6 \end{bmatrix} = \begin{bmatrix} 11\\0\\0 \end{bmatrix}.\end{split}\]

  2. ✍ Prove the unfinished items in Theorem 3.4.2, namely that a Householder reflector \(\mathbf{P}\) is symmetric and orthogonal.

  3. ✍ Let \(\mathbf{P}\) be a Householder reflector as in (3.4.1).

    (a) Find a vector \(\mathbf{u}\) such that \(\mathbf{P}\mathbf{u} = -\mathbf{u}\). (Fig. 3.4.1 may be of help.)

    (b) What algebraic condition is necessary and sufficient for a vector \(\mathbf{x}\) to satisfy \(\mathbf{P}\mathbf{x}=\mathbf{x}\)? In \(n\) dimensions, how many such linearly independent vectors are there?

  4. ✍ Under certain circumstances, computing the vector \(\mathbf{v}\) in (3.4.4) could lead to subtractive cancellation, which is why line 12 of Function 3.4.4 reads as it does. Devise an example that causes subtractive cancellation if (3.4.4) is used.

  5. ✍ Suppose QR factorization is used to compute the solution of a square linear system, \(\mathbf{A}\mathbf{x}=\mathbf{b}\), i.e., let \(m=n\).

    (a) Find an asymptotic flop count for this procedure, and compare it to the LU factorization algorithm.

    (b) Show that \(\kappa_2(\mathbf{A}) = \kappa_2(\mathbf{R})\).

  6. ✍ Prove that \(\kappa_2(\mathbf{A})=\kappa_2(\mathbf{R})\) when \(\mathbf{A}\) is not square. (Be careful! You can’t take an inverse of \(\mathbf{A}\) or \(\mathbf{R}\).)

  7. Another algorithmic technique for orthogonally introducing zeros into a matrix is the Givens rotation. Given a 2-vector \([\alpha,\, \beta]\), it defines an angle \(\theta\) such that

    \[\begin{split}\begin{bmatrix} \cos(\theta) & \sin(\theta) \\ -\sin(\theta) & \cos(\theta) \end{bmatrix} \begin{bmatrix} \alpha \\ \beta \end{bmatrix} = \begin{bmatrix} \sqrt{\alpha^2 + \beta^2} \\ 0 \end{bmatrix}.\end{split}\]

    (a) ✍ Given \(\alpha\) and \(\beta\), show how to compute \(\cos \theta\) and \(\sin \theta\) without evaluating any trig functions.

    (b) ⌨ Given the vector \(\mathbf{z}=[1\;2\;3\;4\;5]^T\), use Julia to find a sequence of Givens rotations that transforms \(\mathbf{z}\) into the vector \(\| \mathbf{z} \|\mathbf{e}_1\). (Hint: You can operate only on pairs of elements at a time, introducing a zero at the lower of the two positions.)

  8. ✍ Derive the result of Theorem 3.4.5 by analyzing Function 3.4.4 without lines 20–22.

  9. ✍ Suppose \(m=Kn\) for constant \(K \ge 1\) as both \(m\) and \(n\) go to infinity. Show that the flop counts from Theorem 3.4.5 and Theorem 3.2.7 have a ratio of 1 when \(K=1\) and approaches 2 as \(K\to \infty\).

Contents