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

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)

migrated from math.stackexchange.comFeb 3 at 14:44

This question came from our site for people studying math at any level and professionals in related fields.

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.

• Thanks for the answer! I don't see, that $\tau$ is just a $1 \times 1$-sized $\mathbf T$. In the cited paper, in Algorithm 5, $\mathbf Y$ is $v$, and $\mathbf T$ is -2. So it ends up as the textbook version, not as the LAPACK version. Do I miss something? – geza Feb 5 at 11:25

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.

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.