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

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

Hi Ben,

Let's get this sorted :) So far I didn't participate much in the
discussion so I may be missing a lot of things.

When it was question of making pivoting/nopivoting an option of the
LDLT class, I thought we were talking about a template parameter. If
fact, in your patch, it is a runtime option. In all other
decompositions in Eigen, we have made this be separate classes, e.g.
HouseholderQR vs ColPivHouseholderQR. Frankly I would like LDLT to do
the same. Keep LDLT pivoting, and add a NoPivLDLT class. This will
actually make your NoPivLDLT be more optimized (the biggest advantage
will be smaller executable code). The only advantage of having that as
runtime param is to share code if a user is using both flavors of
LDLT. I believe that the most common case is that a program will use
at most one flavor of LDLT.

If you send a patch with this change, I promise I'll review it :)

Benoit

2010/6/25 Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>:
> Hi Gael and everyone else,
>
> On Wed, Jun 16, 2010 at 1:05 PM, Gael Guennebaud
> <gael.guennebaud@xxxxxxxxx> wrote:
>> On Sun, Jun 6, 2010 at 5:34 AM, 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 :( .
>>>
>>> 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.
>>
>> When I redesigned LDLT to support blocking and to work fully inplace,
>> I also thought it might be cool to have such a pivoting option,
>> however I was unsure whether this should be a runtime and or compile
>> time option. And if it is a runtime option, shall it be set once for
>> all by the ctor, or modifiable when calling compute, or via a specific
>> function... Finally, we can consider this option to be in pair with
>> the compute* options for the other decompositions, so let's do the
>> same here, i.e.:
>>
>> add a "int option" parameter to all LDLT ctor and to the compute
>> method defaulting to "Pivoting". To disable pivoting will use the
>> NoPivoting constant. That's better than meaningless booleans.
>
> Here is a patch that uses the options=Pivoting or options=NoPivoting
> syntax instead of pivoting=true or pivoting=false and avoids code
> duplication.
>
>>> 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.
>>
>> maybe this is fixed now ?
>
> I am not sure I understand you. As far as I know, it was not broken in
> the Pivoting case, but the decomposition is done in-place now. So, if
> we finish early in the NoPivoting case, then we have to put zeros on
> the diagonal to get vectorD() to extract the correct thing, right?
>
>>> 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?
>>
>> why not.
>
> I have not added the singular tests yet, but I can do that soon.
>
>> You should also make the solve function skips the transpositions when
>> no pivoting has been computed.
>
> I have not done this yet either. What did you decide about making a
> new class versus putting a flag in the class definition?
>
> Ben
>