Computing QR factorizations
Contents
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.
A Householder reflector is a matrix of the form
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,
The reason \(\mathbf{P}\) is called a reflector is sketched in Fig. 3.4.1.
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
(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
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.
The proofs of symmetry and orthogonality are left to Exercise 2. For the last fact, we use (3.4.2) to compute
Since \(\mathbf{e}_1^T\mathbf{z}=z_1\),
leading finally to
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.
A = rand(float(1:9),6,4)
m,n = size(A)
(6, 4)
z = A[:,1];
v = normalize(z - norm(z)*[1;zeros(m-1)])
P₁ = I - 2v*v' # reflector
6×6 Matrix{Float64}:
0.449719 0.449719 0.374766 0.599625 0.299813 0.0749532
0.449719 0.632466 -0.306279 -0.490046 -0.245023 -0.0612557
0.374766 -0.306279 0.744768 -0.408372 -0.204186 -0.0510464
0.599625 -0.490046 -0.408372 0.346606 -0.326697 -0.0816743
0.299813 -0.245023 -0.204186 -0.326697 0.836651 -0.0408372
0.0749532 -0.0612557 -0.0510464 -0.0816743 -0.0408372 0.989791
We check that this reflector introduces zeros as it should:
P₁*z
6-element Vector{Float64}:
13.341664064126332
1.942890293094024e-15
1.0269562977782698e-15
1.9984014443252818e-15
9.992007221626409e-16
2.220446049250313e-16
Now we replace \(\mathbf{A}\) by \(\mathbf{P}\mathbf{A}\).
A = P₁*A
6×4 Matrix{Float64}:
13.3417 10.7933 15.8151 15.4404
1.94289e-15 -1.00357 1.43032 2.10209
1.02696e-15 1.33036 0.358601 3.25174
1.9984e-15 -8.67143 -0.426239 -5.19721
9.99201e-16 2.66429 3.28688 3.4014
2.22045e-16 6.66607 5.07172 -0.149651
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.0881869 0.116902 -0.761985 0.234119 0.585768
0.116902 0.987441 0.081859 -0.0251511 -0.0629283
-0.761985 0.081859 0.466433 0.163938 0.410174
0.234119 -0.0251511 0.163938 0.94963 -0.126026
0.585768 -0.0629283 0.410174 -0.126026 0.684683
We now apply this reflector to rows 2 and below only.
A[2:m,:] = P₂*A[2:m,:]
A
6×4 Matrix{Float64}:
13.3417 10.7933 15.8151 15.4404
-1.21004e-15 11.3801 3.98095 4.86363
1.36567e-15 5.55112e-17 0.0845907 2.95508
-2.09384e-16 -2.22045e-15 1.35979 -3.26349
1.67754e-15 6.66134e-16 2.73812 2.80726
1.91926e-15 8.32667e-16 3.69873 -1.63618
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}:
13.3417 10.7933 15.8151 15.4404
0.0 11.3801 3.98095 4.86363
0.0 0.0 4.79939 -0.531909
0.0 0.0 0.0 5.44589
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
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.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.
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.
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#
⌨ 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}\]✍ Prove the unfinished items in Theorem 3.4.2, namely that a Householder reflector \(\mathbf{P}\) is symmetric and orthogonal.
✍ 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?
✍ 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.
✍ 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})\).
✍ 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}\).)
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.)
✍ Derive the result of Theorem 3.4.5 by analyzing Function 3.4.4 without lines 20–22.
✍ 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\).