9
$\begingroup$

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)

$\endgroup$

migrated from math.stackexchange.com Feb 3 at 14:44

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

7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – geza Feb 5 at 11:25
2
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.