Skip to main content

Subsection 3.3.5 Forming \(Q \)

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

ALWAYS

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} } \\ ~~~=~~~~ \\ \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) } \quad {\rm where~ } \\ ~~~=~~~~ \\ \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{.}\)

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 ) \) ).
which, given \([ A, t ] = {\rm HQR\_unb\_var1}( A ) \) from 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.m. You will want to use Assignments/Week03/matlab/test_FormQ.m to check your implementation. Input is the \(m \times n \) matrix \(A \) and vector \(t \) that resulted from [ A, t ] = HQR( A ). Output is the matrix Q for the QR factorization. You may want to use Assignments/Week03/matlab/test_FormQ.m to check your implementation.

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

Modify the answer for Homework 3.3.4.2.

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=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*}
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