What is the reason that LAPACK uses $\tau$ in QR decomposition (instead of normalizing the reflection vector)?

Question

LAPACK's QR routine stores Q as Householder reflectors. It scales the reflection vector $v$ with $1/v_1$, so the first element of the result becomes $1$, so it doesn't have to be stored. And it stores a separate $\tau$ vector, which contains the needed scale factors. So a reflector matrix is like this: $$H=I-\tau v v^T,$$

where $v$ is not normalized. While, in textbooks, the reflector matrix is

$$H = I-2vv^T,$$

where $v$ is normalized.

Why does LAPACK scale $v$ with $1/v_1$, instead of normalizing it?

Needed storage is the same (instead of $\tau$, $v_1$ has to be stored), and afterwards, applying $H$ can be done faster, as there is no need to multiply with $\tau$ (multiplication with $2$ in the textbook version can be optimized, if instead of simple normalization, $v$ is scaled by $\sqrt 2/\|v\|$).

(The reason of my question is that I'm writing an QR and SVD routine, and I'd like to know the reason of this decision, whether I need to follow it or not)

Answer 1

It's the blocked variant of Householder-QR that is driving this design. If you look in Golub and Van Loan's book (Ch 5.2 or so) they talk about how k-iterations of the algorithm can be blocked together by accumulating the individual reflectors into a rank-k reflector of the form $\mathbf I + \mathbf W \mathbf Y^{\mathrm T}$, where both $\mathbf W$ and $\mathbf Y$ are "tall-skinny" matrices with size $n \times k$. This algorithm does more work but is faster in practice because it's rich in gemm() calls. Unfortunately it is wasteful in storage due to the need to represent $\mathbf W$ and $\mathbf Y$ independently.

In a later paper (cited below), Van Loan describes a more efficient "symmetrized" data structure, a block reflector of the form $\mathbf I + \mathbf Y \mathbf T \mathbf Y^{\mathrm T}$. Here $\mathbf Y$ is still $n \times k$, but the flop/storage requirement for forming $\mathbf W$ has been eliminated by introducing $\mathbf T$, a small $k \times k$ upper triangular matrix. Although the need to multiply by $\mathbf T$ introduces a small amount of extra work, it's typically a net gain because $k << n$.

Within LAPACK, the non-blocked algorithm is really just a limiting $k \rightarrow 1$ case of the block algorithm, all the way down to the choice of symbols (which leads us to $\tau$, a little $1\times1$ version of the $\mathbf T$ triangle).

Citation: Schreiber, Robert, and Charles Van Loan. "A storage-efficient WY representation for products of Householder transformations." SIAM Journal on Scientific and Statistical Computing 10.1 (1989): 53-57.

Answer 2

You don't have to store $\tau$, you can re-compute it from the rest of the vector. (You can recompute $v_1$ from the other entries also in the normalized version, but it's clearly an unstable computation because of those subtractions.)

Actually, you can reuse the lower-triangular part of $R$ to store $v_2,...v_n$, so that the factorization is computed fully in-place. Lapack cares a lot abut these in-place versions of algorithms.

Answer 3

My suggestion is based on the documentation for Intel MKL https://software.intel.com/en-us/mkl-developer-reference-c-geqrf. It looks like the values on and above the diagonal of the output store R so there is only lower triangle left for Q. The it seems natural to use additional storage for the scaling factors.