Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices

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

Hi Ben,

Sorry that nobody replied to you. As you can see, this is a busy
period, hopefully that excuses us.
Gael, since you rewrote the LDLT recently, do you want to handle this?
I don't want to give you homework, just trying to optimize the load
balancing according and memory locality, as this is hopefully still in
your L2 cache ;-)

Benoit

2010/6/15 Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>:
> Hi,
>
> On Sat, Jun 5, 2010 at 11:34 PM, Ben Goodrich <bgokgm@xxxxxxxxxxxxxx> wrote:
>> Hi,
>>
>> On Sun, Feb 28, 2010 at 1:44 AM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
>>> If I understand correctly your other email, you'd be happy with
>>> instead a non-pivoting LDLt. That sounds more reasonable than
>>> "unpivoting" a pivoting LDLt.
>>>
>>> So: OK to add a non-pivoting option to LDLt.
>>>
>>> Benoit
>>
>> Here is a unified patch to reimplement the non-pivoting option. Sorry
>> for the more than 3 month delay; had to write words in dissertation
>> rather than code :( .
>
> Has anyone tried my patch or has comments on the things below? I
> rebased the patch in the attachment to make sure it would still apply
> and to improve the documentation a little more. -- Ben
>
>> Some thoughts:
>>
>> 0) If you changed your mind and don't want this option anymore, that's
>> fine. I can apply the patch locally. That said, having a pivot option
>> makes it possible to be compatible with eigen2 and is probably useful
>> to someone besides just me. And it is easy to do.
>>
>> 1) Gael made some big changes to LDLT a couple of days ago that I
>> think I understand but might have misunderstood something when I
>> reimplemented the non-pivoting option on top of them.
>>
>> 2) For example, in the positive semi-definite but (numerically)
>> rank-deficient case, I had to put exact zeros on the diagonal of D to
>> get the reconstruction right. This contravenes Benoit's statement
>> earlier in this thread that Eigen does not "finish" decompositions of
>> rank-deficient matrices. But I didn't immediately see any other good
>> way to do it, in light of Gael's recent changes to LDLT.
>>
>> 3) I implemented the option like
>>
>> if(pivot)
>> {
>> ....
>> }
>> else
>> {
>> ... // mostly copy-and-paste from above
>> }
>>
>> which resulted in some code duplication but avoids the need to
>> repeatedly evaluate if(pivot) inside the loop. I guess that is
>> slightly faster and clearer, but if you would rather not have the code
>> duplication, I can do it the other way.
>>
>> 4) I made some changes to the documentation explaining the pivot
>> option, but just noticed (and haven't really read yet) today's long
>> documentation thread.
>>
>> 5) I copied-and-pasted a block inside /test/cholesky.cpp and exercised
>> the pivot=false option. It seems to work when you do ./check.sh
>> cholesky. I did some other tests locally with singular matrices, but
>> /test/cholesky.cpp does not seem to have any tests with singular
>> matrices, so maybe some should be added?
>>
>> Thanks,
>> Ben
>>
>