Subsection 3.3.6 Applying \(Q^H \)
¶ VIDEO 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 \(\widehat 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} ) \gt 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 \\
a_{21}
\end{array} \right)
\left(
\begin{array}{c c }
1 \amp a_{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 { a_{21} } }
\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{.}\)
