Re: [eigen] GivensQR

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

Wow, now I've started reading your code and it's quite subtle. I can
see where you're coming from, as sometimes one needs only Q, sometimes
one needs only solving, etc. I can see how on a single vector solve
it's much faster to apply the sequence of rotations as you do, than to
compute the matrix U and multiply by it. Very interesting, thank you.
But it'll take me some time to absorb into Eigen :) yes, if you can
write comments, it's probably useful ;) also, have you benchmarked a
single vector solve with your class against PartialLU? Because, since
you're optimizing so much for this kind of use, I wonder if you can
beat PartialLU; if not, then it's not worth optimizing so much.

Cheers,
Benoit

2009/8/18 Andrea Arteaga <yo.eres@xxxxxxxxx>:
>
>>
>> > Pros:
>> >
>> > The computation of Q is performed only if strictly needed; by default (by
>> > the constructor) no computation is performed. Q is computed only if the
>> > user requests it.
>>
>> that's something we'll have to uniformize across various
>> decompositions later on.
>> Actually I believe that the default should be "convenient but not
>> fast" i.e. compute Q, with an option to not compute Q. But we had yet
>> to discuss that on the list.
>
>
> I implemented the class paing attention to the perfomance, but in facts the result is quite
> "comfortable". The idea behind the method compute() is not to ask the user which matrices are to be
> computed, but only a tool to update the decomposition if the matrix has changed (since the original
> matrix is stored as alias). The decision behind the behaviour of solve() is simply: there is no need
> for the unitary matrix in order to solve the sistem.
> The only "problematic" issue is: if one needs both Q and R, he should request Q and then R and not
> viceversa. In the first case, Q and R are computed and Q is returned and then R is simply returned;
> in the second one R is computed and returned and then Q and R are both (re-)computed and Q is
> returned. In other words:
>
> Matrix<...> m;
> /* ... */
>
> GivensQR<Matrix<...> > qr1(m);
> qr1.matrixQ(); // Q and R are computed; Q is returned
> qr1.matrixR(); // R is returned (was computed before)
>
>
> GivensQR<Matrix<...> > qr2(m);
> qr2.matrixR(); // R is computed and returned; Q is not yet computed.
> qr2.matrixQ(); // Q and R are computed; Q is returned.
>
>
> The first case is faster than the second one, in spite of the same result. An alternative is:
>
> GivensQR<Matrix<...> > qr3(m);
> qr3.compute(); // Both matrices are computed and stored
> qr3.matrixR(); // R is returned; no computation
> qr3.matrixQ(); //Q is returned; no computation
>
> If the user doesn't care about performance, this issue is not a problem; if he does, the only think
> he must care is to call matrixQ() or compute() before matrixR() if he needs both Q and R.
>
>
>
> If you need help in order to wirk on my code, I can provide a fully-commented version.
>
> Just a note: I thought to store the R matrix in a row-major storage and Q in column-major and then I
> constructed the algorithm by updating R row-wise and Q column-wise. But I'm not sure about the
> performance and the stability, so I finally stored both in the default behaviour. What is your advice
> about that?
>
>
>