[ Thread Index |
| More lists.tuxfamily.org/eigen Archives
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: Re: [eigen] GivensQR
- From: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
- Date: Tue, 18 Aug 2009 18:20:03 -0400
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:message-id:subject:from:to:content-type :content-transfer-encoding; bh=Z9bxJiThmF5MDHXbZMBdz1zNYEgkYU583JO49S+9fHg=; b=O1q+Qe4+2S6cGcs8QgmDfJ184nlm0ZY96sI8UDP98XJqidRMzeVdhuKqH6nS1FWFCt UaURBogY6Yead+uYVw4mgIp4VNXMK61d8wssL/oeMhrU6Wk/OYfc0DvJ5PYhDed+lBOF ckUTt0CWqZOZWGDZCHQ9V73Bvm7wEW15sF9jk=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type:content-transfer-encoding; b=bdp/lvsPklKjR31lxA+dmJVy1e+v50gkOxYQOaexwKW8pOj3JEOmtB/YSCgtw01t3E AtnHe6K09FAFnA/Ir3wM9YHUdCgrl5OhMPj86F5zauwL3BFrzg2wrny45rAnW3V/EI/1 y5EvUV9fOR3OP85VrUmjDaXNeAnpQIvEXO5Ig=
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.
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?