Subsection 11.2.3 Reduction to bidiagonal form
¶ VIDEO
Homework 11.2.3.1 .
Let \(B \in \mathbb R^{m \times m}\) be a bidiagonal matrix:
\begin{equation*}
B =
\left( \begin{array}{c c c c c c c c}
\beta_{0,0} \amp \beta_{0,1} \amp 0 \amp \cdots \amp 0
\amp 0\\
0 \amp \beta_{1,1} \amp \beta_{1,2} \amp \cdots \amp 0
\amp 0\\
\vdots \amp \ddots \amp \ddots \amp \ddots \amp \vdots
\amp \vdots \\
0 \amp 0 \amp 0 \amp \cdots \amp \beta_{m-2,m-2}
\amp \beta_{m-2,m-1}
\\
0 \amp 0 \amp 0 \amp \cdots \amp 0 \amp
\beta_{m-1,m-1}
\end{array}
\right).
\end{equation*}
Show that \(T = B^T B \) is a tridiagonal symmetric matrix.
Solution
\begin{equation*}
\begin{array}{l}
B^T B \\
~~~ = ~~~~ \\
\left( \begin{array}{c c c c }
\beta_{0,0} \amp \beta_{0,1} \amp 0 \amp \cdots \\
0 \amp \beta_{1,1} \amp \beta_{1,2} \amp \cdots \\
0 \amp 0 \amp \beta_{2,2} \amp \cdots \\
\vdots \amp \ddots \amp \ddots \amp \ddots \\
\end{array}
\right)^T
\left( \begin{array}{c c c c }
\beta_{0,0} \amp \beta_{0,1} \amp 0 \amp \cdots \\
0 \amp \beta_{1,1} \amp \beta_{1,2} \amp \cdots \\
0 \amp 0 \amp \beta_{2,2} \amp \cdots \\
\vdots \amp \ddots \amp \ddots \amp \ddots \\
\end{array}
\right)
~~~ = ~~~~ \\
\left( \begin{array}{c c c c }
\beta_{0,0} \amp 0\amp 0 \amp \cdots \\
\beta_{0,1} \amp \beta_{1,1} \amp 0 \amp \cdots \\
0 \amp \beta_{1,2} \amp \beta_{2,2} \amp \cdots \\
\vdots \amp \ddots \amp \ddots \amp \ddots \\
\end{array}
\right)
\left( \begin{array}{c c c c }
\beta_{0,0} \amp \beta_{0,1} \amp 0 \amp \cdots \\
0 \amp \beta_{1,1} \amp \beta_{1,2} \amp \cdots \\
0 \amp 0 \amp \beta_{2,2} \amp \cdots \\
\vdots \amp \ddots \amp \ddots \amp \ddots \\
\end{array}
\right) \\
~~~ = ~~~~ \\
\left( \begin{array}{c c c c }
\beta_{0,0}^2 \amp \beta_{0,1} \beta_{0,0} \amp 0 \amp \cdots \\
\beta_{0,1} \beta_{0,0}\amp \beta_{0,1}^2 + \beta_{1,1}^2 \amp \beta_{1,2} \beta_{1,1} \amp \cdots \\
0 \amp \beta_{1,2} \beta_{1,1}\amp \beta_{1,2}^2 + \beta_{2,2}^2 \amp \cdots \\
\vdots \amp \ddots \amp \ddots \amp \ddots \\
\end{array}
\right)
\end{array} \\
\end{equation*}
Given that we can preprocess our problem by computing its QR factorization, we focus now on the case where \(A \in
\mathbb C^{m \times m} \text{.}\) The next step is to reduce this matrix to bidiagonal form by multiplying the matrix from the left and right by two sequences of unitary matrices.
Once again, we employ Householder transformations. In Subsubsection 3.3.3.3 , we introduced the function
\begin{equation*}
\left[
\FlaTwoByOneSingleLine
{\rho}
{ u_2},
\tau \right] :=
{\rm Housev} \left(
\FlaTwoByOneSingleLine
{ \chi_1 }
{ x_2 } \right),
\end{equation*}
implemented by
function [ rho, ...
u2, tau ] = Housev( chi1, ...
x2 ),
to compute the vector \(u = \FlaTwoByOneSingleLine{1}{u_2}\) that reflects \(x \) into \(\complexone \| x \|_2 e_0 \) so that
\begin{equation*}
\left( I - \frac{1}{\tau}
\FlaTwoByOneSingleLine{1}{u_2}
\FlaTwoByOneSingleLine{1}{u_2}^H \right)
\FlaTwoByOneSingleLine
{ \chi_1 }
{ x_2 }
=
\complexone
\begin{array}[t]{c}
\underbrace{\| x \|_2}\\
\rho
\end{array} e_0.
\end{equation*}
In Subsection 10.3.1 , we introduced a variation on this function:
\begin{equation*}
\left[
u,
\tau \right] :=
{\rm Housev1} \left( x \right)
\end{equation*}
implemented by the function
function [ u, tau ] = Housev1( x ).
They differ only in how the input and output are passed to and from the function. We also introduce the notation \(H(u,\tau) \) for the transformation \(I - \frac{1}{\tau} u u^H \text{.}\)
We now describe an algorithm for reducing a square matrix to bidiagonal form:
Partition \(A \rightarrow
\FlaTwoByTwoSingleLine
{ \alpha_{11} }{ a_{12}^T }
{ a_{21} }{ A_{22} }\text{.}\)
Update \([
\left( \begin{array}{c}
\alpha_{11} \\
a_{21}
\end{array}
\right) , \tau_1 ] \becomes {\rm Housev}(
\left( \begin{array}{c}
\alpha_{11} \\
a_{21}
\end{array}
\right) ) \text{.}\) This overwrites \(\alpha_{11} \) with \(\complexone \left\|\left( \begin{array}{c}
\alpha_{11} \\
a_{21}
\end{array}
\right) \right\|_2 \) and \(a_{21} \) with \(u_{21} \text{.}\) Implicitly, \(a_{21} \) in the updated matrix equals the zero vector.
Update
\begin{equation*}
\left( \begin{array}{c}
a_{12}^T \\
A_{22}
\end{array}
\right)
:=
H( \left( \begin{array}{c}
1 \\
u_{21}
\end{array}
\right), \tau_1 )
\left( \begin{array}{c}
a_{12}^T \\
A_{22}
\end{array}
\right) .
\end{equation*}
This introduces zeroes below the first entry in the first column, as illustrated by
\begin{equation*}
\left( \begin{array}{| c c c c c} \hline
\times \amp \times \amp \times \amp \times \amp \times
\\
\times \amp \times \amp \times \amp \times \amp \times \\
\times \amp \times \amp \times \amp \times \amp \times \\
\times \amp \times \amp \times \amp \times \amp \times \\
\times \amp \times \amp \times \amp \times \amp \times
\\
\end{array} \right)~~
\longrightarrow
~~
\left( \begin{array}{| c c c c c} \hline
\times \amp \times \amp \times \amp \times \amp \times
\\
0 \amp \times \amp \times \amp \times \amp \times \\
0 \amp \times \amp \times \amp \times \amp \times \\
0 \amp \times \amp \times \amp \times \amp \times \\
0 \amp \times \amp \times \amp \times \amp \times
\\
\end{array} \right)
\end{equation*}
The Householder vector that introduced the zeroes is stored over those zeroes.
Next, we introduce zeroes in the first row of this updated matrix.
The matrix is still partitioned as \(A \rightarrow
\FlaTwoByTwoSingleLine
{ \alpha_{11} }{ a_{12}^T }
{ 0 }{ A_{22} }\text{,}\) where the zeroes have been overwritten with \(u_{21} \text{.}\)
We compute \([
u_{12}
, \rho_1 ] \becomes {\rm Housev1}(
(a_{12}^T )^T ) \text{.}\) The first element of \(u_{12} \) now holds \(\complexone{ \| ( a_{12}^T )^T \|_2 } \) and the rest of the elements define the Householder transformation that introduces zeroes in \((
a_{12}^T )^T \) below the first element. We store \(u_{12}^T \) in \(a_{12}^T \text{.}\)
After setting the first entry of \(u_{12} \) explicitly to one, we update \(A_{22} := A_{22} H( u_{12}, \rho_1 ) \text{.}\)
This introduces zeroes to the right of the first entry of \(a_{12}^T \text{,}\) as illustrated by
\begin{equation*}
\left( \begin{array}{| c c c c c} \hline
\times \amp \times \amp \times \amp \times \amp \times
\\
0 \amp \times \amp \times \amp \times \amp \times \\
0 \amp \times \amp \times \amp \times \amp \times \\
0 \amp \times \amp \times \amp \times \amp \times \\
0 \amp \times \amp \times \amp \times \amp \times
\\
\end{array} \right)~~
\longrightarrow
~~
\left( \begin{array}{| c c c c c} \hline
\times \amp \times \amp 0 \amp 0 \amp 0
\\
0 \amp \times \amp \times \amp \times \amp \times \\
0 \amp \times \amp \times \amp \times \amp \times \\
0 \amp \times \amp \times \amp \times \amp \times \\
0 \amp \times \amp \times \amp \times \amp \times
\\
\end{array} \right)
\end{equation*}
The Householder vector that introduced the zeroes is stored over those zeroes.
The algorithm continues this with the updated \(A_{22}\) as illustrated in Figure 11.2.3.1 .
\begin{equation*}
\begin{array}{| c c c c c} \hline
\times \amp { \times} \amp { \times} \amp { \times} \amp { \times} \\
\times \amp \times \amp { \times} \amp { \times} \amp { \times} \\
\times \amp \times \amp \times \amp { \times} \amp { \times} \\
\times \amp \times \amp \times \amp \times \amp { \times} \\
\times \amp \times \amp \times \amp \times \amp \times
\end{array}
\end{equation*}
\(\begin{array}{c}
\\
\\
\\
\longrightarrow
\end{array}\)
\begin{equation*}
\begin{array}{| c c c c c} \hline
\times \amp { \times} \amp { 0} \amp { 0} \amp { 0} \\
0 \amp \times \amp { \times} \amp { \times} \amp { \times} \\
0 \amp \times \amp \times \amp { \times} \amp { \times} \\
0 \amp \times \amp \times \amp \times \amp { \times} \\
0 \amp \times \amp \times \amp \times \amp \times
\end{array}
\end{equation*}
\(\begin{array}{c}
\\
\\
\\
\longrightarrow
\end{array}\)
\begin{equation*}
\begin{array}{c | c c c c}
\times \amp { \times} \amp { 0} \amp { 0} \amp { 0} \\ \hline
0 \amp \times \amp { \times} \amp { 0} \amp { 0} \\
0 \amp 0 \amp \times \amp { \times} \amp { \times} \\
0 \amp 0 \amp \times \amp \times \amp { \times} \\
0 \amp 0 \amp \times \amp \times \amp \times
\end{array}
\end{equation*}
\(\begin{array}{c}
\\
\\
\\
\longrightarrow
\end{array}\)
\begin{equation*}
\begin{array}{c c | c c c}
\times \amp { \times} \amp { 0} \amp
{ 0} \amp { 0} \\
0 \amp \times \amp { \times} \amp { 0} \amp { 0} \\ \hline
0 \amp 0 \amp \times \amp { \times} \amp { 0} \\
0 \amp 0 \amp 0 \amp \times \amp { \times} \\
0 \amp 0 \amp 0 \amp \times \amp \times
\end{array}
\end{equation*}
\(\begin{array}{c}
\\
\\
\\
\longrightarrow
\end{array}\)
\begin{equation*}
\begin{array}{c c c | c c}
\times \amp { \times} \amp
{ 0} \amp { 0} \amp { 0} \\
0 \amp \times \amp { \times} \amp { 0} \amp { 0} \\
0 \amp 0 \amp \times \amp { \times} \amp { 0} \\ \hline
0 \amp 0 \amp 0 \amp \times \amp { \times} \\
0 \amp 0 \amp 0 \amp 0 \amp \times
\end{array}
\end{equation*}
Figure 11.2.3.1. An illustration of the reduction of a square matrix to bidiagonal form. The \(\times \)s denote nonzero elements in the matrix.
Ponder This 11.2.3.2 .
Fill in the details for the above described algorithm that reduces a square matrix to bidiagonal form. In particular:
For the update
\begin{equation*}
\left( \begin{array}{c}
a_{12}^T \\
A_{22}
\end{array}
\right)
:=
H( \left( \begin{array}{c}
1 \\
u_{21}
\end{array}
\right), \tau_1 )
\left( \begin{array}{c}
a_{12}^T \\
A_{22}
\end{array}
\right) ,
\end{equation*}
describe how all the different parts of
\begin{equation*}
\left( \begin{array}{c}
a_{12}^T \\
A_{22}
\end{array}
\right)
\end{equation*}
are updated. (Hint: look at the QR factorization algorithm in Subsection 3.3.4 .)
For the update \(A_{22} := A_{22} H( u_{12}, \rho_1 ) \text{,}\) describe explicitly how \(A_{22} \) is updated. (Hint: look at Homework 10.3.1.1 .)
Next, state the algorithm by completing the skeleton in Figure 11.2.3.2 .
Finally, analyze the approximate cost of the algorithm, when started with a \(m \times m \) matrix.
\begin{equation*}
\newcommand{\routinename}{
\left[ A, t, r \right] := \mbox{BiRed-unb}( A )
}
\newcommand{\guard}{
m( A_{TL} ) \lt m( A )
}
\newcommand{\partitionings}{
A \rightarrow
\FlaTwoByTwo{A_{TL}}{A_{TR}}
{A_{BL}}{A_{BR}},
t \rightarrow
\FlaTwoByOne{ t_{T} }{t_B},
r \rightarrow
\FlaTwoByOne{ r_{T} }{r_B}
}
\newcommand{\partitionsizes}{
A_{TL} {\rm ~is~} 0 \times 0,
t_T, r_T \mbox{ have } 0 \mbox{ 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{r_T}{r_B}
\rightarrow
\FlaThreeByOneB{r_0}{\rho_1}{r_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{r_T}{r_B}
\leftarrow
\FlaThreeByOneT{r_0}{\rho_1}{r_2}
}
\newcommand{\update}{
\begin{array}{l}
~~~\\
\phantom{[ a_{21}, \tau_1 ] := {\rm Housev1}( a_{21} )} \\
\mbox{ Update } \phantom{A_{22} := H( u_{21} ) A_{22} H( u_{21}
)} \mbox{ via the steps } \\
~~~~~~
\left\{
\phantom{\begin{array}{ll}
y_{21} := A_{22} u_{21} \amp
\mbox{(Hermitian matrix-vector multiply!)}\\
\beta := u_{21}^H y_{21} / 2 \\
w_{21} := ( y_{21} - \beta u_{21} / \tau_1 ) / \tau_1 \\
A_{22} := A_{22} - \mbox{tril}(u_{21} w_{21}^H +
w_{21} u_{21}^H)
\amp
\mbox{(Hermitian~rank-2~update)}
\end{array}}
\right. \\
~~~\\
\phantom{[ a_{21}, \tau_1 ] := {\rm Housev1}( a_{21} )} \\
\mbox{ Update } \phantom{A_{22} := H( u_{21} ) A_{22} H( u_{21}
)} \mbox{ via the steps } \\
~~~~~~
\left\{
\phantom{\begin{array}{ll}
y_{21} := A_{22} u_{21} \amp
\mbox{(Hermitian matrix-vector multiply!)}\\
\beta := u_{21}^H y_{21} / 2 \\
w_{21} := ( y_{21} - \beta u_{21} / \tau_1 ) / \tau_1 \\
A_{22} := A_{22} - \mbox{tril}(u_{21} w_{21}^H +
w_{21} u_{21}^H)
\amp
\mbox{(Hermitian~rank-2~update)}
\end{array}}
\right.
\end{array}
}
\FlaAlgorithm
\end{equation*}
Figure 11.2.3.2. Algorithm skeleton for reduction of a square matrix to bidiagonal form.
Ponder This 11.2.3.3 .
Once you have derived the algorithm in Ponder This 11.2.3.2 , implement it.
You may want to start by executing git pull to update your directory Assignments.
In directory Assignments/Week11/matlab/, you will find the following files:
Housev.m and Housev1.m: Implementations of the function Housev and Housev1.
BiRed.m: A code skeleton for a function that reduces a square matrix to bidiagonal form.
[ B, t, r ] = BiRed( A, t, r )
returns the diagonal and first superdiagonal of the bidiagonal matrix in B, stores the Householder vectors below the subdiagonal and above the first superdiagonal, and returns the scalars \(\tau \) and \(\rho \) in vectors t and r .BiFromB.m: A function that extracts the bidiagonal matrix from matrix \(B \text{,}\) which also has the Householder vector information in it.
Bbi = BiFromB( B )
test_BiRed.m: A script that tests
BiRed.
These resources give you the tools to implement and test the reduction to bidiagonal form.
