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

[Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>]

Hi Benoit,

On Mon, Jun 28, 2010 at 8:51 AM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
> 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

That is all fine with me; I don't have any strong reason to prefer a
runtime option. But I interpreted Gael's email on June 16 as asking
for a runtime option and then later he asked about the possibly of a
separate class. So, it seems we now have a consensus on the separate
class at least. :)

Ben

> 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
>>
>
>
>