Re: [eigen] Getting Householder reflections right

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

BTW, JAMA used another scheme than LAPACK (or the GSL) for that. I
don't remember how they did it but I can have a look... perhaps there
is at least one good thing in their code.... Right now I can tell they
based their work on that book: Bowdler, Martin, Reinsch, and
Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra.

Also, I'm not considering the LAPACK compatibility as an argument
anymore. If there exists a better option, then let's go for it.

cheers,
Gael.

On Fri, May 8, 2009 at 9:22 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
> Hi,
>
> I'm wondering how to do Householder reflections right.
> Books (Golub&vanLoan) and libraries (LAPACK) seem to _always_ use the
> following packed storage for a Householder vector v:
> they normalize it so that v(0)=1, hence they don't need to store v(0)
> which is useful for storing decompositions in a packed way in a single
> matrix.
> However, they store separately the norm of v, in order to avoid having
> to recompute it everytime.
>
> I'm not saying this is always a bad idea, i'm saying that's not always
> a good idea.
>
> That's a good idea for the QR decomposition, because it allows to
> store the whole decomposition in a single matrix, plus a vector to
> store the norms of the v's.
>
> But right now i'm looking at bidiagonalization and i'm seriously
> wondering if this packed storage is a good idea.
> With the packed storage, i can store in a single matrix the bidiagonal
> matrix, and the essential parts of the householder vectors, but then i
> need to store in 2 separate vectors the norms of the householder
> vectors.
> While _without_ packed storage, i'd be able to store entirely in a
> single matrix (of the same size) the whole householder vectors, and
> then i'd store the bidiagonal matrix separately (obviously in a
> compact way).
> So here, I don't see any advantage in using packed Householder storage.
> On the other hand, NOT packing the householder vectors will give
> simpler code and slightly better speed.
>
> Maybe the only advantage is to make our bidiagonalization storage
> compatible with LAPACK. But if that's important, we can add conversion
> methods, and, having n^2 complexity, the conversion won't have a big
> overhead compared to the n^3 cost of the decomposition itself.
>
> Same remarks for these close relatives: tridiagonalization and
> hessenberg decomposition (which i'll rework along the way).
>
> Opinions?
>
> Cheers,
> Benoit
>
>
>