Re: [eigen] Getting Householder reflections right

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

LOL.
First I googled for "block householder", found this paper.
Then Rohit IMs me and points me to the same paper.
Then you. Small world!

Anyway, I was about to draw the same conclusion: The construction of
the block householder reflector relies on ....(drumroll).... SVD !

So yes, let's implement a plain SVD first.

Cheers,
Benoit

2009/5/8 Keir Mierle <mierle@xxxxxxxxx>:
> I found this:
>
> http://www.jstage.jst.go.jp/article/ipsjdc/2/0/298/_pdf
>
> Looks complicated. I suggest implementing the householder
> bidiagonalization in a straightforward way first, where code clarity
> and modularity goes before speed. Then we can look at blocking
> versions.
>
> Keir
>
> On Fri, May 8, 2009 at 12:33 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
>> thanks for the reminder, it's actually the right time to look at that
>> as it may impact the design.
>> Benoit
>>
>> 2009/5/8 Rohit Garg <rpg.314@xxxxxxxxx>:
>>> Possibly OT
>>>
>>> May be you should look at block wise House holder transformations. I'd
>>> infact recommend that all algorithms is-as-far-as-possible be done in
>>> a blockwise manner.
>>>
>>> On Sat, May 9, 2009 at 12:52 AM, 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
>>>>
>>>>
>>>>
>>>
>>>
>>>
>>> --
>>> Rohit Garg
>>>
>>> http://rpg-314.blogspot.com/
>>>
>>> Senior Undergraduate
>>> Department of Physics
>>> Indian Institute of Technology
>>> Bombay
>>>
>>>
>>>
>>
>>
>>
>
>
>