8.3.

Householder QR factorization 231

8.3 Householder QR factorization

A fundamental problem to avoid in numerical codes is the situation where one starts with
large values and one ends up with small values with large relative errors in them. This
is known as catastrophic cancelation. The Gram-Schmidt algorithms can inherently fall
victim to this: column aj is successively reduced in length as components in the directions
of q0 , . . . , qj−1 are subtracted, leaving a small vector if aj was almost in the span of the first
j columns of A. Application of a unitary transformation to a matrix or vector inherently
preserves length. Thus, it would be beneficial if the QR factorization can be implementated
as the successive application of unitary transformations. The Householder QR factorization
accomplishes this.
The first fundamental insight is that the product of unitary matrices is itself unitary.
Thus, if, given A ∈ Cm×n one could find a sequence of unitary matrices H0 , . . . , Hn−1 such
that ! "
R
Hn−1 · · · , H0 A = ,
0
where R is upper triangular, then
! " ! " ! "! "
R R QL R
Hn−1 · · · , H0 A = H0 · · · Hn−1 =Q = = QL R,
# $% & 0 0 QR 0
Q
where QL equals the first n columns of A. Then A = QL R is the QR factorization of A. The
second fundamental insight is that the desired unitary transformations H0 , . . . , Hn−1 can be
computed and applied cheaply.

8.3.1 Householder transformations (reflectors)

In this section we discuss Householder transformations, also referred to as reflectors.
Definition 8.6 Let u ∈ Cn be a vector of unit length ("u"2 = 1). Then H = I − 2uuH is
said to be a reflector or Householder transformation.
We observe:
• Any vector z that is perpendicular to u is left unchanged:
(I − 2uuH )z = z − 2u(uH z) = z.

• Any vector x can be written as x = z + uH xu where z is perpendicular to u and uH xu

is the component of x in the direction of u. Then
(I − 2uuH )x = (I − 2uuH )(z + uH xu) = z + uH xu − 2u #$%&
uH z − 2uuH uH xu
0
= z + uH xu − 2uH x #$%&
uH u u = z − uH xu.
0
232 Chapter 8. QR factorization

x
H
x=z+u xu
H
u xu
v=x−y

z
H
u xu u
u

H H u
(I − 2 u u )x (I − 2 u u )x
y

Figure 8.8: ??

This can be interpreted as follows: The space perpendicular to u acts as a “mirror”: any
vector in that space (along the mirror) is not reflected, while any other vector has the
component that is orthogonal to the space (the component outside, orthogonal to, the mirror)
reversed in direction, as illustrated in Fig. 8.8. Notice that a reflection preserves the length
of the vector.
Exercise 8.7 Show that if H is a reflector, then

• HH = I (reflecting the reflection of a vector results in the original vector),

• H = H H , and

• H H H = I (a reflection is a unitary matrix and thus preserves the norm).

Next, let us ask the question of how to reflect a given x ∈ Cn into another vector y ∈ Cn ,
where "x"2 = "y"2 . In other words, how do we compute vector u so that (I − 2uuH )x = y.
From our discussion above, we need to find a vector u that is perpendicular to the space
with respect to which we will reflect. From Fig. 8.8(right) we notice that the vector from y
to x, v = x − y, is perpendicular to the desired space. Thus, u must equal a unit vector in
the direction v: u = v/"v"2 .
Remark 8.8 In subsequent discussion we will prefer to give Householder transformations
as I −uuH /τ , where τ = uH u/2 so that u needs no longer be a unit vector, just a direction.
The reason for this will become obvious later.

In the next subsection, we will wish to find a Householder transformation, H, for a vector,
x, such that Hx equals a vector with zeroes below the top element. More precisely, partition
! "
χ1
x= ,
x2
8.3. Householder QR factorization 233

where χ1 equals the first!element

" of x and x2 is the rest of x. Then we will wish to find a
1
Householder vector u = so that
u2
' ! "! "H ( ! " ! "
1 1 1 χ1 $± "x"2
I− = .
τ u2 u2 x2 0
1
Here $± denotes a complex scalar on the complex
! unit circle
" Notice that this means that y
$
± "x"2
in the previous discussion equals the vector so that the direction of u is given
0
by ! "
χ1 − $± "x"2
v= .
x2
We now wish to normalize this vector so its first entry equals “1”:
! " ! "
v 1 χ1 − $ ± "x"2 1
u= = = .
"v"2 χ1 − $ ± "x"2 x2 x2 /ν1
where ν1 = χ1 − $
± "x"2 . (Note that if ν1 = 0 then u2 can be set to 0.)
Exercise 8.9 Verify that
' ! "! "H ( ! " ! "
1 1 1 χ1 ρ
I− =
τ u2 u2 x2 0

where τ = uH u/2 = (1 + uH 2 u2 )/2 and ρ = $ ± "x"2 .

2 2
Hint:
)z ρρ̄ = |rho| = "x" 2 since H preserves the norm. Also, "x"22 = |χ1 |2 + "x2 "22 and
z̄
= f racz|z|.
The choice $ ± is important when computer arithmetic is used and roundoff error is
a concern. In practice one chooses ν1 = χ1 + sign(χ1 )"x"2 to avoid creating catastrophic
cancelation so that $ ± = −sign(χ1 ) (−$ ± points in the opposite direction, in the complex
plane, as χ1 ). If χ1 is real valued, this simply means that sign(χ1 ) equals the sign of χ1 in
the usual sense. The reason is as follows (again restricting our discussion to real numbers):
If χ1 is positive and "x"2 is almost equal to χ1 , then χ1 − "x"2 is a small number and if
there is error in χ1 and/or "x"2 , this error becomes large relative to the result χ1 − "x"2 .
Regardless of whether χ1 is positive, negative, or complex valued, this can be avoided by
using ν1 = χ1 + sign(χ1 )"x"2 . It is not hard to see that sign(χ1 ) = χ1 /|χ1 |.
Let us introduce the notation
*! " + !! ""
ρ χ1
, τ := Housev
u2 x2
as the computation of the above mentioned vector u2 , and scalars ρ and τ , from vector x.
We will use the notation H(x) for the transformation I − τ1 uuH where u and τ are computed
by Housev(x).
1
For those who have problems thinking in terms of complex numbers, pretend this entire discussion deals
with real numbers and treat $± as ±.
234 Chapter 8. QR factorization

*! " + !! ""
ρ χ1
Algorithm: , τ = Housev
u2 x2
χ2 :=,"x
! 2 "2 ",
, χ1 ,
,
α := , , (= "x"2 )
χ2 ,2
ρ = −sign(χ1 )"x"2 ρ := −sign(χ1 )α
ν1 = χ1 + sign(χ1 )"x"2 ν1 := χ1 − ρ
u2 = x2 /ν1 u2 := x2 /ν1
χ2 = χ2 /|ν1 |(= "u2 "2 )
τ = (1 + uH
2 u2 )/2 τ = (1 + χ22 )/2

Figure 8.9: Computing the Householder transformation. Left: simple formulation. Right:
efficient computation. Note: I have not completely double-checked these formulas
for the complex case. They work for the real case.

8.3.2 Algorithms
Let A be an m × n with m ≥ n. We will now show how to compute A → QR, the QR
factorization, as a sequence of Householder transformations applied to A, which eventually
zeroes out all elements of that matrix below the diagonal.
In the first iteration, we partition
! "
α11 aT12
A→ .
a21 A22

Let *! " + ! "

ρ11 α11
, τ1 = Housev
u21 a21
be the Householder transform computed from the first column of A. Then applying this
Householder transform to A yields
! " ' ! "! "H ( ! "
α11 aT12 1 1 1 α11 aT12
:= I−
a21 A22 τ1 u 2 u2 a21 A22
! "
ρ11 aT12 − w12
T
= T ,
0 A22 − u21 w12

T
where w12 = (aT12 +uH
21 A22 )/τ1 . Computation of a full QR factorization of A will now proceed
with the updated matrix A22 .
Now let us assume that after k iterations of the algorithm matrix A contains
 
! " R00 r01 R02
RT L RT R
A→ = 0 α11 aT12  ,
0 ABR
0 a21 A22
8.3. Householder QR factorization 235

where RT L and R00 are k × k upper triangular matrices. Let

*! " + ! "
ρ11 α11
, τ1 = Housev .
u21 a21

and update
  
I 0 R00 r01 R02
 ' ! "! (
"H 
A :=  1 1 1   0 α11 aT12 
0 I − τ1
u2 u2 0 a21 A22
   H   
0 0 R00 r01 R02
 1 
= I −  1   1    0 α11 aT12 
τ1
u2 u2 0 a21 A22
 
R00 r01 R02
=  0 ρ11 T
a12 − w12T ,
T
0 0 A22 − u21 w12
T
where again w12 = (aT12 + uH
21 A22 )/τ1 .
Let    H 
0k 0k
 1 
Hk = I −  1   1  
τ1
u21 u21
be the Householder transform so computed during the (k + 1)st iteration. Then upon com-
pletion matrix A contains
! "
RT L
R= = Hn−1 · · · H1 H0 Â
0

where  denotes the original contents of A and RT L is an upper triangular matrix. Rear-
ranging this we find that
 = H0 H1 · · · Hn−1 R
which shows that if Q = H0 H1 · · · Hn−1 then  = QR.
Exercise 8.10 Show that
     H 
I 0 0 0
 ' ! "! "H (   1 
 1 1  = I −  1   1   .
0 I − τ11 τ1
u2 u2 u2 u2

Typically, the algorithm overwrites the original matrix A with the upper triangular ma-
trix, and at each step u21 is stored over the elements that become zero, thus overwriting a21 .
(It is for this reason that the first element of u was normalized to equal “1”.) In this case Q
236 Chapter 8. QR factorization

Algorithm: [A, t] = URt(A)

! " ! "
{U \R}T L RT R tT
Partition A → and t →
UBL ABR tB
where {U \R}T L is 0 × 0 and tT has 0 elements
while n(ABR ) (= 0 do
Repartition    
! " {U \R}00 r01 R02 ! " t0
{U \R}T L RT R t T
→ uT10 α11 aT12  and →  τ1 
UBL ABR tB
U20 a21 A22 t2
where α11 and τ1 are scalars
*! " + ! "
ρ11 α11
, τ := Housev
u21 ! 1 " ! a21
! " "! "
aT12 1 3 4 aT12
Update := I − τ1
1
1 uH
21
A22 u21 A22
via the steps
T
• w12 := (aT12 + uH
21 A22 )/τ1
! T " ! "
a12 aT12 − w12
T
• := T
A22 A22 − u21 w12

Continue with    
! " {U \R}00 r01 R02 ! " t0
{U \R}T L RT R tT
← uT10 ρ11 T 
r12 and ←  τ1 
UBL ABR tB
U20 u21 A22 t2
endwhile

Figure 8.10: Unblocked Householder transformation based QR factorization.

is usually not explicitly formed as it can be stored as the separate Householder vectors below
the diagonal of the overwritten matrix. The algorithm that overwrites A in this manner is
given in Fig. 8.10.
We will let
[{U \R}, t] = URt(A)
denote the operation that computes the QR factorization of m × n matrix A, with m ≥ n,
via Householder transformations. It returns the Householder vectors and matrix R in the
first argument and the vector of scalars “τi ” that are computed as part of the Householder
transformations in t.
Theorem 8.11 Given A ∈ Cm×n the cost of the algorithm in Figure 8.10 is given by

2
CHQR (m, n) ≈ 2mn2 − n3 flops.
3

T
Proof: The bulk of the computation is in w12 = (aT12 + uH T
21 A22 )/τ1 and A22 − u21 w12 . During
the kth iteration (when RT L is k × k), this means a matrix-vector multiplication (uH 21 A22 )
8.3. Householder QR factorization 237

and rank-1 update with matrix A22 which is of size approximately (m − k) × (n − k) for a
cost of 4(m − k)(n − k) flops. Thus the total cost is approximately
n−1
5 n−1
5 n−1
5 n−1
5
4(m − k)(n − k) = 4 (m − n + j)j = 4(m − n) j+4 j2
k=0 j=0 j=0 j=0
n−1
5
= 2(m − n)n(n − 1) + 4 j2
j=0
6 n
4 2
≈ 2(m − n)n2 + 4 x2 dx = 2mn2 − 2n3 + n3 = 2mn2 − n3 .
0 3 3

8.3.3 Forming Q
Given A ∈ Cm×n , let [A, t] = URt(A) yield the matrix A with the Householder vectors stored
below the diagona, R stored on and above the diagonal, and the τi stored in vector t. We
now discuss how to form the first n columns of Q = H0 H1 · · · Hn−1 . Notice that to pick out
the first n columns we must form
! " ! " ! "
In×n In×n In×n
Q = H0 · · · Hn−1 = H0 · · · Hk−1 Hk · · · Hn−1 .
0 0 0
# $% &
Bk
where Bk is defined as indicated.
Lemma 8.12 Bk has the form
! " ! "
In×n Ik×k 0
Bk = Hk · · · Hn−1 = .
0 0 B̃k

Proof: The proof of this is by induction on k:

! "
In×n
• Base case: k = n. Then Bn = , which has the desired form.
0
• Inductive step: Assume the result is true for Bk . We show it is true for Bk−1 :
! " ! "
In×n Ik×k 0
Bk−1 = Hk−1 Hk · · · Hn−1 Hk−1 Bk = Hk−1 .
0 0 B̃k
  
I(k−1)×(k−1) 0 I(k−1)×(k−1) 0 0
! "
=  1 3 4   0 1 0 
0 I − τ1k 1 uHk
uk 0 0 B̃k
238 Chapter 8. QR factorization

 
I(k−1)×(k−1) 0
! ! " "! "
=  1 1 3 4 1 0 
0 I− τk
1 uH
k
uk 0 B̃k
 
I(k−1)×(k−1) 0
! " " !
=  1 0 1 3 4  where ykT = uHk B̃k /τk
0 − 1/τk tTk
0 B̃k uk
 
I(k−1)×(k−1) 0
! "
=  1 − 1/τk −ykT 
0
−uk /τk Bk − uk ykT
 
I(k−1)×(k−1) 0 0 ! "
I(k−1)×(k−1) 0
=  0 1 − 1/τk −ykT = .
T 0 B̃k−1
0 −uk /τk Bk − uk yk

• By the Principle of Mathematical Induction the result holds for B0 , . . . , Bn .

Theorem 8.13 Given [A, t] = URt(A) from Figure 8.10, the algorithm in Figure 8.11 over-
writes A with the first n = n(A) columns of Q as defined by the Householder transformations
stored below the diagonal of A and in the vector t.

Proof: The algorithm is justified by the proof of Lemma 8.12.

Theorem 8.14 Given A ∈ Cm×n the cost of the algorithm in Figure 8.11 is given by
2
CFormQ (m, n) ≈ 2mn2 − n3 flops.
3

Proof: Hence the proof for Theorem 8.11 can be easily modified to establish this result.

Exercise 8.15 If m = n then Q could be accumulated by the sequence

Q = (· · · ((IH0 )H1 ) · · · Hn−1 ).

Give a high-level reason why this would be (much) more expensive than the algorithm in
Figure 8.11.
8.4. Solving Linear Least-Squares Problems 239

Algorithm: [A] = FormQ(A, t)

! " ! "
{U \R}T L RT R tT
Partition A → and t →
UBL ABR tB
where {U \R}T L is n(A) × n(A) and tT has n(A) elements
while n(AT R ) (= 0 do
Repartition    
! " {U \R}00 r01 R02 ! " t0
{U \R}T L RT R t T
→ uT10 ρ11 T 
r12 and →  τ1 
UBL ABR tB
U20 u21 A22 t2
where ρ11 and τ1 are scalars
! " ! ! " "! "
α11 aT12 1 3 4 1 0
Update := I− 1
1 uH
a21 A22 τ1 u21 21
0 A22
via the steps
• α11 := 1 − 1/τ1
• aT12 := −(uH
21 A22 )/τ1

• A22 := A22 + u21 w12

T

• a21 := −u21 /τ1

Continue with    
! " {U \R}00 r01 R02 ! " t0
{U \R}T L RT R t T
← uT10 α11 aT12  and ←  τ1 
UBL ABR tB
U20 a21 A22 t2
endwhile

Figure 8.11: Algorithm for overwriting A with Q from the Householder transformations
stored as Householder vectors below the diagonal of A (as produced by [A, t] = URt(A) ).

8.4 Solving Linear Least-Squares Problems