Skip to main content

Subsection 3.3.6 Applying \(Q^H \)

In a future chapter, we will see that the QR factorization is used to solve the linear least-squares problem. To do so, we need to be able to compute \(\hat y = Q^H y \) where \(Q^H = H_{n-1} \cdots H_0 \text{.}\)

Let us start by computing \(H_0 y \text{:}\)

\begin{equation*} \begin{array}{l} \left( I - \frac{1}{\tau_1} \FlaTwoByOneSingleLine { 1 } { u_{2} } \FlaTwoByOneSingleLine { 1 } { u_{2} }^H \right) \FlaTwoByOneSingleLine { \psi_1 } { y_2 } \\ ~~~=~~~~ \\ \FlaTwoByOneSingleLine { \psi_1 } { y_2 } - \FlaTwoByOneSingleLine { 1 } { u_{2} } \begin{array}[t]{c} \underbrace{ \FlaTwoByOneSingleLine { 1 } { u_{2} }^H \FlaTwoByOneSingleLine { \psi_1 } { y_2 } / \tau_1 } \\ \omega_1 \end{array} \\ ~~~=~~~~ \\ \FlaTwoByOneSingleLine { \psi_1 } { y_2 } - \omega_1 \FlaTwoByOneSingleLine { 1 } { u_{2} } \\ ~~~=~~~~ \\ \FlaTwoByOneSingleLine { \psi_1 - \omega_1 } { y_2 - \omega_1 u_2 }. \end{array} \end{equation*}

More generally, let us compute \(H_k y \text{:}\)

\begin{equation*} \left( I - \frac{1}{\tau_1} \FlaThreeByOneB { 0 } { 1 } { u_{2} } \FlaThreeByOneB { 0 } { 1 } { u_{2} }^H \right) \FlaThreeByOneB { y_0 } { \psi_1 } { y_2 } = \FlaThreeByOneB { y_0 } { \psi_1 - \omega_1} { y_2 - \omega_1 u_2 }, \end{equation*}

where \(\omega_1 = ( \psi_1 + u_2^H y_2 ) / \tau_1 \text{.}\) This motivates the algorithm in Figure 3.3.6.1 for computing \(y := H_{n-1} \cdots H_0 y \) given the output matrix \(A \) and vector \(t \) from routine \({\rm HQR\_unb\_var1} \text{.}\)

\begin{equation*} \newcommand{\routinename}{[ y ] = {\rm Apply\_QH}( A, t, y )} \newcommand{\guard}{n( A_{BR} ) \lt 0} \newcommand{\partitionings}{ A \rightarrow \FlaTwoByTwo{ A_{TL} } {A_{TR}} { A_{BL} } { A_{BR} } , t \rightarrow \FlaTwoByOne{ t_{T} } { t_{B} } , y \rightarrow \FlaTwoByOne{ y_{T} } { y_{B} } } \newcommand{\partitionsizes}{ A_{TL} {\rm ~is~} 0 \times 0 {\rm ~and~} t_T, y_T {\rm ~have~} 0 {\rm ~elements}} \newcommand{\repartitionings}{ \FlaTwoByTwo{ A_{TL} } {A_{TR}} { A_{BL} } { A_{BR} } \rightarrow \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} } \rightarrow \FlaThreeByOneB{ t_{0} } { \tau_{1} } { t_{2} } , \FlaTwoByOne{ y_{T} } { y_{B} } \rightarrow \FlaThreeByOneB{ y_{0} } { \psi_{1} } { y_{2} } } \newcommand{\moveboundaries}{ \FlaTwoByTwo{ A_{TL} } {A_{TR}} { A_{BL} } { A_{BR} } \leftarrow \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} } \leftarrow \FlaThreeByOneT{ t_{0} } { \tau_{1} } { t_{2} } , \FlaTwoByOne{ y_{T} } { y_{B} } \leftarrow \FlaThreeByOneT{ y_{0} } { \psi_{1} } { y_{2} } } \newcommand{\update}{ \begin{array}{l} {\rm Update~} \left( \begin{array}{c} \psi_1 \\ y_{2} \end{array} \right) := \left( I - \frac{1}{\tau_1} \left( \begin{array}{c} 1 \\ u_{21} \end{array} \right) \left( \begin{array}{c c } 1 \amp u_{21}^H \end{array} \right) \right) \left( \begin{array}{c} \psi_1 \\ y_{2} \end{array} \right) \\ ~~~~~~~~~~{\rm via~the~steps} \\ ~~~~~~~~~~~~~~~ \omega_{1} := ( \psi_{1} + a_{21}^H y_{2} )/\tau_1 \\ ~~~~~~~~~~~~~ \FlaTwoByOneSingleLine { \psi_{1} } { y_{2} } := \FlaTwoByOneSingleLine { \psi_{1} - { \omega_{1} } } { y_{2} - \omega_1 { u_{2} } } \end{array} } \FlaAlgorithm \end{equation*}
Figure 3.3.6.1. Algorithm for computing \(y := Q^H y (= H_{n-1} \cdots H_0 y) \) given the output from the algorithm \({\rm HQR\_unb\_var1}\text{.}\)
Homework 3.3.6.1.

What is the approximate cost of algorithm in Figure 3.3.6.1 if \(Q \) (stored as Householder vectors in \(A \)) is \(m \times n \text{.}\)

Solution

The cost of this algorithm can be analyzed as follows: When \(y_T \) is of length \(k \text{,}\) the bulk of the computation is in a dot product with vectors of length \(m-k-1 \) (to compute \(\omega_1 \)) and an axpy operation with vectors of length \(m-k-1 \) to subsequently update \(\psi_1 \) and \(y_2 \text{.}\) Thus, the cost is approximately given by

\begin{equation*} \sum_{k=0}^{n-1} 4 (m-k-1) = 4 \sum_{k=0}^{n-1} m - 4 \sum_{k=0}^{n-1} (k-1) \approx 4 m n - 2 n^2. \end{equation*}

Notice that this is much cheaper than forming \(Q\) and then multiplying \(Q^H y \text{.}\)