The theory behind computing a single Householder reflector to zero out part of a column of a matrix is pretty well described in Matrix Computations by Golub and Van Loan. However, the blocked algorithm in Lapack's xDLARFB routines do not match the description in Matrix Computations. The book describes computation of a non symmetric reflector of the form $I-WY^T$, where as Lapack computes something of the form $I-VTV^T$. Where is the algorithm as implemented in Lapack described? Or alternatively, can anyone provide a high level picture of what is going on?
$\begingroup$
$\endgroup$
1
-
2$\begingroup$ It's just a matrix version of the Householder reflector I - tau v v^H, where v is a unit vector. An inductive argument allows you to see that the scalar tau yields the triangular matrix T when applying the product of several reflectors. $\endgroup$– Jack PoulsonNov 15, 2012 at 7:03
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
I suppose I should give an actual answer instead of just a comment. The reason for the upper-triangular structure of the matrix $T$ in the accumulated Householder transformation $I - V T V^H$ can be seen from considering the product of two Householder transformations, say $I - \tau_1 v_1 v_1^H$ and $I - \tau_2 v_2 v_2^H$ (if $v_1$ and $v_2$ are unit vectors, then $\tau_1=\tau_2=2$, but it is often useful to scale the vectors such that their first entry is one and then to store the remaining entries directly within the matrix being modified). We see that
$$ (I - \tau_1 v_1 v_1^H) (I - \tau_2 v_2 v_2^H) = I - \tau_1 v_1 v_1^H - \tau_2 v_2 v_2^H + \tau_1 \tau_2 (v_1^H v_2) v_1 v_2^H, $$
which, as an exercise for the reader, can be rearranged into the form
$$ I - \begin{pmatrix} v_1 & v_2 \end{pmatrix} \begin{pmatrix} \tau_1 & -\tau_1 \tau_2 v_1^H v_2 \\ 0 & \tau_2 \end{pmatrix} \begin{pmatrix} v_1 & v_2 \end{pmatrix}^H. $$ This insight in fact generalizes to an arbitrary product of Householder reflections, and I recommend reading Accumulating Householder transformations, revisited, as well as the paper mentioned by Schreiber.
$\begingroup$
$\endgroup$
Jack's right, AFAIK, that what LAPACK uses is the "compact $WY$ form:" $I - \tau V T V^H$ where $T$ is triangular. It's equal to a product of elementary reflectors. And a product of reflectors is not in general a reflector.
($R$ is a reflector if $R$ is its own inverse.
$R_1 R_2 R_1 R_2$ is not $I$, even if $R_1 R_1$ and $R_2 R_2$ are.) But even though the compact $WY$ form
is not a reflector, it gets the usual jobs done, efficiently.
The theory of genuine (orthogonal, symmetric) block reflectors of the form $I - \tau V V^H$ (which can be used to make a matrix block upper triangular, for example) is contained in a paper by Parlett and me, published in SINUM 25:1 (1988), p. 189-205.
