Processing math: 100%
Skip to main content
πŸ”—

Subsection 3.3.4 Householder QR factorization algorithm

fit width

πŸ”—
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. The process is illustrated in Figure 3.3.4.1.

πŸ”—
Original matrix[(ρ11u21),Ο„1]=    Housev(Ξ±11a21)(Ξ±11aT12a21A22):=(ρ11aT12βˆ’wT120A22βˆ’u21wT12)β€˜β€˜Move forwardβ€³Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—Γ—0Γ—Γ—Γ—0Γ—Γ—Γ—0Γ—Γ—Γ—0Γ—Γ—Γ—Γ—Γ—Γ—Γ—0Γ—Γ—Γ—0Γ—Γ—Γ—0Γ—Γ—Γ—0Γ—Γ—Γ—Γ—Γ—Γ—Γ—0Γ—Γ—Γ—0Γ—Γ—Γ—0Γ—Γ—Γ—0Γ—Γ—Γ—Γ—Γ—Γ—Γ—0Γ—Γ—Γ—00Γ—Γ—00Γ—Γ—00Γ—Γ—Γ—Γ—Γ—Γ—0Γ—Γ—Γ—00Γ—Γ—00Γ—Γ—00Γ—Γ—Γ—Γ—Γ—Γ—0Γ—Γ—Γ—00Γ—Γ—00Γ—Γ—00Γ—Γ—Γ—Γ—Γ—Γ—0Γ—Γ—Γ—00Γ—Γ—000Γ—000Γ—Γ—Γ—Γ—Γ—0Γ—Γ—Γ—00Γ—Γ—000Γ—000Γ—Γ—Γ—Γ—Γ—0Γ—Γ—Γ—00Γ—Γ—000Γ—000Γ—Γ—Γ—Γ—Γ—0Γ—Γ—Γ—00Γ—Γ—000Γ—0000Γ—Γ—Γ—Γ—0Γ—Γ—Γ—00Γ—Γ—000Γ—0000
Figure 3.3.4.1. Illustration of Householder QR factorization.

πŸ”—
In the first iteration, we partition

A→(α11aT12a21A22).

πŸ”—
Let

[(ρ11u21),Ο„1]=Housev((Ξ±11a21))

πŸ”—
be the Householder transform computed from the first column of A. Then applying this Householder transform to A yields

(Ξ±11aT12a21A22):=(Iβˆ’1Ο„1(1u21)(1u21)H)(Ξ±11aT12a21A22)=(ρ11aT12βˆ’wT120A22βˆ’u21wT12),

πŸ”—
where wT12=(aT12+uH21A22)/Ο„1. Computation of a full QR factorization of A will now proceed with the updated matrix A22.

πŸ”—

fit width

πŸ”—

πŸ”—
Homework 3.3.4.1.

Show that

(I00Iβˆ’1Ο„1(1u21)(1u21)H)=(Iβˆ’1Ο„1(01u21)(01u21)H).
Solution
\begin{equation*} \begin{array}{rcl} \FlaTwoByTwo { I } { 0 } { 0 } { I - \frac{1}{\tau_1} \FlaTwoByOneSingleLine { 1 } { u_{21} } \FlaTwoByOneSingleLine { 1 } { u_{21} }^H } \amp=\amp I - \FlaTwoByTwo {0 } { 0 } { 0 } { \frac{1}{\tau_1} \FlaTwoByOneSingleLine { 1 } { u_{21} } \FlaTwoByOneSingleLine { 1 } { u_{21} }^H } \\ \amp=\amp I - \frac{1}{\tau_1} \FlaTwoByTwo {0 } { 0 } { 0 } { \FlaTwoByOneSingleLine { 1 } { u_{21} } \FlaTwoByOneSingleLine { 1 } { u_{21} }^H } \\ \amp=\amp I - \frac{1}{\tau_1} \FlaThreeByThreeBR {0 } { 0 }{0} { 0 }{1}{u_{21}^H} { 0}{u_{21}}{ u_{21} u_{21}^H} \\ \amp=\amp \left( I - \frac{1}{\tau_1} \FlaThreeByOneB { 0 } { 1 } { u_{21} } \FlaThreeByOneB { 0 } { 1 } { u_{21} }^H \right). \end{array} \end{equation*}

πŸ”—
More generally, let us assume that after k iterations of the algorithm matrix A contains

A→(RTLRTR0ABR)=(R00r01R020α11aT120a21A22),

πŸ”—
where RTL and R00 are kΓ—k upper triangular matrices. Let

[(ρ11u21),Ο„1]=Housev((Ξ±11a21)).

πŸ”—
and update

A:=(I00(Iβˆ’1Ο„1(1u21)(1u21)H))(R00r01R020Ξ±11aT120a21A22)=(Iβˆ’1Ο„1(01u21)(01u21)H)(R00r01R020Ξ±11aT120a21A22)=(R00r01R020ρ11aT12βˆ’wT1200A22βˆ’u21wT12),

πŸ”—
where, again, wT12=(aT12+uH21A22)/Ο„1.

πŸ”—
Let

Hk=(Iβˆ’1Ο„1(0k1u21)(0k1u21)H)

πŸ”—
be the Householder transform so computed during the (k+1)st iteration. Then upon completion matrix A contains

R=(RTL0)=Hnβˆ’1β‹―H1H0Λ†A

πŸ”—
where Λ†A denotes the original contents of A and RTL is an upper triangular matrix. Rearranging this we find that

Λ†A=H0H1β‹―Hnβˆ’1R

πŸ”—
which shows that if Q=H0H1β‹―Hnβˆ’1 then Λ†A=QR.

πŸ”—
Typically, the algorithm overwrites the original matrix A with the upper triangular matrix, 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 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 Figure 3.3.4.2.

πŸ”—
[A,t]=HQR_unb_var1(A)Aβ†’(ATLATRABLABR) and tβ†’(tTtB)   ATL is 0Γ—0 and tT has 0 elementswhile n(ABR)>0   (ATLATRABLABR)β†’(A00a01A02aT10Ξ±11aT12A20a21A22) and (tTtB)β†’(t0Ο„1t2)   [(Ξ±11a21),Ο„1]:=[(ρ11u21),Ο„1]=Housev(Ξ±11a21)Update (aT12A22):=(Iβˆ’1Ο„1(1u21)(1uH21))(aT12A22)      via the steps                wT12:=(aT12+aH21A22)/Ο„1               (aT12A22):=(aT12βˆ’wT12A22βˆ’a21wT12)   (ATLATRABLABR)←(A00a01A02aT10Ξ±11aT12A20a21A22) and (tTtB)←(t0Ο„1t2)endwhile
Figure 3.3.4.2. Unblocked Householder transformation based QR factorization.

πŸ”—
In that figure,

[A,t]=HQR_unb_var1(A)

πŸ”—
denotes 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.

πŸ”—
Homework 3.3.4.2.

Given A∈RmΓ—n show that the cost of the algorithm in Figure 3.3.4.2 is given by

CHQR(m,n)β‰ˆ2mn2βˆ’23n3 flops.
Solution

The bulk of the computation is in

\begin{equation*} w_{12}^{T} = ( a_{12}^{T} + u_{21}^{H} A_{22} ) / \tau_1 \end{equation*}

and

\begin{equation*} A_{22} - u_{21} w_{12}^{T} \text{.} \end{equation*}

During the \(k \)th iteration (when \(R_{TL} \) is \(k \times k \)), this means a matrix-vector multiplication (\(u_{21}^H A_{22} \)) and rank-1 update with matrix \(A_{22} \) which is of size approximately \((m-k) \times (n-k) \) for a cost of \(4 (m-k)(n-k) \) flops. Thus the total cost is approximately

\begin{equation*} \begin{array}{l} \sum_{k=0}^{n-1} 4 (m-k)(n-k) \\ ~~~ = ~~~~\\ 4 \sum_{j=0}^{n-1} (m - n + j )j \\ ~~~ = ~~~~ \\ 4 ( m-n) \sum_{j=0}^{n-1} j + 4 \sum_{j=0}^{n-1} j^2\\ ~~~ = ~~~~ \\ 2 ( m-n) n ( n-1 ) + 4 \sum_{j=0}^{n-1} j^2 \\ ~~~ \approx ~~~~ \\ 2 ( m-n) n^2 + 4 \int_0^n x^2 dx \\ ~~~ = ~~~~ \\ 2 m n^2 - 2 n^3 + \frac{4}{3} n^3 \\ ~~~ = ~~~~ \\ 2 m n^2 - \frac{2}{3} n^3. \end{array} \end{equation*}
πŸ”—
Homework 3.3.4.3.

Implement the algorithm given in Figure 3.3.4.2 as 

function [ A_out, t ] = HQR( A  )

by completing the code in Assignments/Week03/matlab/HQR.m. Input is an mΓ—n matrix A. Output is the matrix A_out with the Householder vectors below its diagonal and R in its upper triangular part. You may want to use Assignments/Week03/matlab/test_HQR.m to check your implementation.

Solution

See Assignments/Week03/answers/HQR.m. Warning: it only checks if \(R \) is computed correctly.

login