Subsection 3.3.5 Forming \(Q \)
¶ VIDEO Given \(A \in \C^{m \times n} \text{,}\) let \([ A, t ] = {\rm HQR\_unb\_var1(A)} \) yield the matrix \(A \) with the Householder vectors stored below the diagonal, \(R \) stored on and above the diagonal, and the scalars \(\tau_i \text{,}\) \(0 \leq i \lt n \text{,}\) stored in vector \(t \text{.}\) We now discuss how to form the first \(n \) columns of \(Q = H_0 H_1 \cdots H_{n-1} \text{.}\) The computation is illustrated in Figure 3.3.5.1 .
\begin{equation*}
\begin{array}{ c | c | c }
{\rm Original~matrix}
\amp
\begin{array}{l}
\FlaTwoByTwo
{ \alpha_{11} }{ a_{12}^T }
{ a_{21} }{ A_{22} }
:= \\
\FlaTwoByTwo
{ \color{red}{1 - 1/\tau_1} }
{ - ( u_{21}^H A_{22} )/\tau_1 }
{ \color{blue} {- u_{21}/\tau_1} }
{ \color{green} {A_{22} + { u_{21} } { a_{12}^T }} }
\end{array}
\amp
{\rm ``Move~forward''}
\\ \hline
\begin{array}{ c c c c | }
1 \amp 0 \amp 0 \amp 0 \\
0 \amp 1 \amp 0 \amp 0 \\
0 \amp 0 \amp 1 \amp 0 \\
0 \amp 0 \amp 0 \amp 1 \\ \hline
0 \amp 0 \amp 0 \amp 0
\end{array}
\amp
\begin{array}{ c c c c |}
1 \amp 0 \amp 0 \amp 0 \\
0 \amp 1 \amp 0 \amp 0 \\
0 \amp 0 \amp 1 \amp 0 \\
0 \amp 0 \amp 0 \amp \color{red}\times \\ \hline
0 \amp 0 \amp 0 \amp \color{blue}\times
\end{array}
\amp
\begin{array}{c c c | c}
1 \amp 0 \amp 0 \amp 0 \\
0 \amp 1 \amp 0 \amp 0 \\
0 \amp 0 \amp 1 \amp 0 \\ \hline
0 \amp 0 \amp 0 \amp \times \\
0 \amp 0 \amp 0 \amp \times
\end{array}
\\ \hline
\amp
\begin{array}{c c c | c}
1 \amp 0 \amp 0 \amp 0 \\
0 \amp 1 \amp 0 \amp 0 \\
0 \amp 0 \amp \color{red} \times \amp \times \\ \hline
0 \amp 0 \amp \color{blue} \times \amp \color{green} \times \\
0 \amp 0 \amp \color{blue} \times \amp \color{green} \times
\end{array}
\amp
\begin{array}{c c | c c}
1 \amp 0 \amp 0 \amp 0 \\
0 \amp 1 \amp 0 \amp 0 \\ \hline
0 \amp 0 \amp \times \amp \times \\
0 \amp 0 \amp \times \amp \times \\
0 \amp 0 \amp \times \amp \times
\end{array}
\\ \hline
\amp
\begin{array}{c c | c c}
1 \amp 0 \amp 0 \amp 0 \\
0 \amp \color{red} \times \amp \times \amp \times \\ \hline
0 \amp \color{blue} \times \amp \color{green} \times \amp \color{green} \times \\
0 \amp \color{blue} \times \amp \color{green} \times \amp \color{green} \times \\
0 \amp \color{blue} \times \amp \color{green} \times \amp \color{green} \times
\end{array}
\amp
\begin{array}{c | c c c}
1 \amp 0 \amp 0 \amp 0 \\ \hline
0 \amp \times \amp \times \amp \times \\
0 \amp \times \amp \times \amp \times \\
0 \amp \times \amp \times \amp \times \\
0 \amp \times \amp \times \amp \times
\end{array}
\\ \hline
\amp
\begin{array}{ c | c c c}
\color{red} \times \amp \times \amp \times \amp \times \\ \hline
\color{blue} \times \amp \color{green} \times \amp \color{green} \times \amp \color{green} \times \\
\color{blue} \times \amp \color{green} \times \amp \color{green} \times \amp \color{green} \times \\
\color{blue} \times \amp \color{green} \times \amp \color{green} \times \amp \color{green} \times \\
\color{blue} \times \amp \color{green} \times \amp \color{green} \times \amp \color{green} \times
\end{array}
\amp
\begin{array}{| c c c c} \hline
\times \amp \times \amp \times \amp \times \\
\times \amp \times \amp \times \amp \times \\
\times \amp \times \amp \times \amp \times \\
\times \amp \times \amp \times \amp \times \\
\times \amp \times \amp \times \amp \times
\end{array}
\end{array}
\end{equation*}
Figure 3.3.5.1. Illustration of the computation of \(Q \text{.}\)Notice that to pick out the first \(n \) columns we must form
\begin{equation*}
Q
\left( \begin{array}{c}
I_{n \times n} \\ \hline
0
\end{array}
\right) \\
~~~=~~~~ \\
H_0 \cdots H_{n-1}
\left( \begin{array}{c}
I_{n \times n} \\ \hline
0
\end{array}
\right) \\
~~~ = ~~~~ \\
H_0 \cdots H_{k-1}
\begin{array}[t]{c}
\underbrace{
H_k \cdots H_{n-1}
\left( \begin{array}{c}
I_{n \times n} \\ \hline
0
\end{array}
\right)
}
\\
B_k
\end{array}
\end{equation*}
so that \(Q = B_0 \text{,}\) where \(B_k = H_k \cdots H_{n-1}
\left( \begin{array}{c}
I_{n \times n} \\ \hline
0
\end{array}
\right) \text{.}\)
Homework 3.3.5.1 .
ALWAYS/SOMETIMES/NEVER:
\begin{equation*}
B_k = H_k \cdots H_{n-1}
\left( \begin{array}{c}
I_{n \times n} \\ \hline
0
\end{array}
\right)
=
\FlaTwoByTwo{ I_{k \times k} }{ 0 }{ 0 }{ \widetilde{B}_k }.
\end{equation*}
for some \((m-k) \times (n-k) \) matrix \(\widetilde{B}_k \text{.}\)
Answer Solution
The proof of this is by induction on \(k \text{:}\)
Base case: \(k = n \text{.}\) Then \(B_n = \left( \begin{array}{c}
I_{n \times n} \\ \hline
0
\end{array}
\right) \text{,}\) which has the desired form.
Inductive step: Assume the result is true for \(B_k \text{.}\) We show it is true for \(B_{k-1} \text{:}\)
\begin{equation*}
\begin{array}{l}
B_{k-1} \\
~~~=~~~~ \\
H_{k-1} H_k \cdots H_{n-1} \left( \begin{array}{c}
I_{n \times n} \\ \hline
0
\end{array}
\right) \\
~~~=~~~~ \\
H_{k-1} B_k \\
~~~=~~~~ \\
H_{k-1} \FlaTwoByTwo{ I_{k \times k} }{ 0 }{ 0 }{ \widetilde{B}_k } \\
~~~=~~~~ \\
\FlaTwoByTwoSingleLine
{ I_{(k-1) \times (k-1)} }{ 0 }
{ 0 }{ I - \frac{1}{\tau_k}
\left( \begin{array}{c}
1 \\ \hline
u_k
\end{array}
\right)
\left( \begin{array}{c | c}
1 \amp u_k^H
\end{array}
\right)
}
\FlaThreeByThreeTL
{ I_{(k-1) \times (k-1)} }{ 0 }{ 0 }
{0 }{ 1 }{ 0 }
{0 }{ 0 }{ \widetilde{B}_k } \\
~~~=~~~~ \\
\FlaTwoByTwoSingleLine
{ I_{(k-1) \times (k-1)} }{ 0 }
{ 0 }{ \left( I - \frac{1}{\tau_k}
\left( \begin{array}{c}
1 \\ \hline
u_k
\end{array}
\right)
\left( \begin{array}{c | c}
1 \amp u_k^H
\end{array}
\right) \right)
\FlaTwoByTwo{1}{0}{0}{\widetilde{B}_k}
}
\\
~~~=~~~~ \lt \mbox{ choose }
y_k^T = u_k^H \widetilde{B}_k/\tau_k \gt
\\
\FlaTwoByTwoSingleLine
{ I_{(k-1) \times (k-1)} }{ 0 }
{ 0 }{
\FlaTwoByTwo{1}{0}{0}{\widetilde{B}_k}
-
\left( \begin{array}{c}
1 \\ \hline
u_k
\end{array}
\right)
\left( \begin{array}{c | c}
1/\tau_k \amp y_k^T
\end{array}\right)
}
\\
~~~=~~~~ \\
\FlaTwoByTwoSingleLine
{ I_{(k-1) \times (k-1)} }{ 0 }
{ 0 }{
\FlaTwoByTwo{1-1/\tau_k}{- y_k^T}{- u_k / \tau_k}{\widetilde B_k - u_k y_k^T}
} \\
~~~=~~~~ \\
\FlaThreeByThreeTL
{ I_{(k-1) \times (k-1)} }{ 0 }{0}
{ 0 }{1-1/\tau_k}{- y_k^T }
{ 0 }{- u_k / \tau_k}{\widetilde B_k - u_k y_k^T} \\
~~~=~~~~ \\
\FlaTwoByTwoSingleLine
{ I_{(k-1) \times (k-1)} }{ 0 }
{0}{ \widetilde{B}_{k-1} }.
\end{array}
\end{equation*}
By the Principle of Mathematical Induction the result holds for \(B_0, \ldots, B_{n} \text{.}\)
VIDEO
The last exercise justifies the algorithm in Figure 3.3.5.2 ,
\begin{equation*}
\newcommand{\routinename}{[ A ] = {\rm FormQ}( A, t )}
\newcommand{\guard}{n( A_{TL} ) > 0}
\newcommand{\partitionings}{
A \rightarrow
\FlaTwoByTwo{ A_{TL} } { A_{TR}}
{ A_{BL} } { A_{BR} },
t \rightarrow
\FlaTwoByOne{ t_{T} }
{ t_{B} }
}
\newcommand{\partitionsizes}{
A_{TL} {\rm ~is~} n(A) \times n(A) {\rm ~and~} t_T {\rm ~has~} n(A)
{\rm ~elements}
}
\newcommand{\repartitionings}{
\FlaTwoByTwo{ A_{TL} } {A_{TR}}
{ A_{BL} } { A_{BR} } \rightarrow
\FlaThreeByThreeTL{ A_{00} }{ a_{01} }{ A_{02} }
{ a_{10}^T }{ \alpha_{11} }{ a_{12}^T }
{ A_{20} }{ a_{21} }{ A_{22} },
\FlaTwoByOne{ t_{T} }
{ t_{B} } \rightarrow
\FlaThreeByOneT{ t_{0} }
{ \tau_{1} }
{ t_{2} }
}
\newcommand{\moveboundaries}{
\FlaTwoByTwo{ A_{TL} } {A_{TR}}
{ A_{BL} } { A_{BR} } \leftarrow
\FlaThreeByThreeBR{ A_{00} }{ a_{01} }{ A_{02} }
{ a_{10}^T }{ \alpha_{11} }{ a_{12}^T }
{ A_{20} }{ a_{21} }{ A_{22} },
\FlaTwoByOne{ t_{T} }
{ t_{B} } \leftarrow
\FlaThreeByOneB{ t_{0} }
{ \tau_{1} }
{ t_{2} }
}
\newcommand{\update}{
\begin{array}{l}
{\rm Update~}
\left(
\begin{array}{c | c}
\alpha_{11} \amp a_{12}^T \\ \hline
a_{21} \amp A_{22}
\end{array}
\right) := \\
~~~~~~~~~~~~~~~ \left(
I -
\frac{1}{\tau_1}
\left(
\begin{array}{c}
1 \\ \hline
u_{21}
\end{array} \right)
\left(
\begin{array}{c | c }
1 \amp u_{21}^H
\end{array}
\right)
\right)
\left(
\begin{array}{c | c}
1 \amp 0 \\ \hline
0 \amp A_{22}
\end{array}
\right) \\
~~~~{\rm via~the~steps} \\
~~~~~~~~~~{ \alpha_{11} }
:=
{ 1 - 1/\tau_1 } \\
~~~~~~~~~~a_{12}^T := - ( a_{21}^H A_{22} )/\tau_1 \\
~~~~~~~~~~{ A_{22} }
:=
{ A_{22} + { a_{21} } { a_{12}^T } }\\
~~~~~~~~~~{ a_{21} }
:=
{ - a_{21}/\tau_1 }
\end{array}
}
\FlaAlgorithm
\end{equation*}
Figure 3.3.5.2. Algorithm for overwriting \(A \) with \(Q \) from the Householder transformations stored as Householder vectors below the diagonal of \(A \) (as produced by \([ A, t ] = {\rm
HQR\_unb\_var1}( A, t ) \) ).Figure 3.3.4.2 , overwrites \(A \) with the first \(n = n(A) \) columns of \(Q \text{.}\)
Homework 3.3.5.2 .
Implement the algorithm in Figure 3.3.5.2 as
function [ A_out ] = FormQ( A, t )
by completing the code in Assignments/Week03/matlab/FormQ.mAssignments/Week03/matlab/test_FormQ.m [ A, t ] = HQR( A ). Output is the matrix Q for the QR factorization. You may want to use Assignments/Week03/matlab/test_FormQ.m
Homework 3.3.5.3 .
Given \(A \in \C^{m \times n} \text{,}\) show that the cost of the algorithm in Figure 3.3.5.2 is given by
\begin{equation*}
C_{\rm FormQ}( m, n ) \approx
2 m n^2 - \frac{2}{3} n^3
{\rm ~flops}.
\end{equation*}
Hint Solution
When computing the Householder QR factorization, the bulk of the cost is in the computations
\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*}
When forming \(Q \text{,}\) the cost is in computing
\begin{equation*}
a_{12}^T := - ( a_{21}^H A_{22} / \tau_1
\end{equation*}
and
\begin{equation*}
A_{22} := A_{22} + u_{21} w_{12}^{T} \text{.}
\end{equation*}
During the when \(A_{TL} \) is \(k \times k \)), these represent, essentially, identical costs: p the 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=n-1}^{0} 4 (m-k)(n-k) \\
~~~ = ~~~~ \lt
\mbox{ reverse the order of the summation } \gt \\
\sum_{k=0}^{n-1} 4 (m-k)(n-k) \\
~~~ = ~~~~\\
4 \sum_{j=1}^{n} (m - n + j )j \\
~~~ = ~~~~ \\
4 ( m-n) \sum_{j=1}^{n} j +
4 \sum_{j=1}^{n} j^2\\
~~~ = ~~~~ \\
2 ( m-n) n ( n+1 ) +
4 \sum_{j=1}^{n} 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*}
Ponder This 3.3.5.4 .
If \(m = n \) then \(Q \) could be accumulated by the sequence
\begin{equation*}
Q = ( \cdots ( ( I H_{0} ) H_1 ) \cdots H_{n-1} ).
\end{equation*}
Give a high-level reason why this would be (much) more expensive than the algorithm in Figure 3.3.5.2
