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

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

2010/6/28 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. :)

Indeed his email was a bit unclear :) trust me, a separate class is
how we do everywhere else:
  FullPivLU vs PartialPivLU
  ColPivHouseholderQR vs FullPivHouseholderQR vs HouseholderQR
etc... not to mention that we chose to let LDLT and LLT be separate classes.

Having a separate class is also a very good idea because, depending on
whether we are pivoting or not, some methods are wanted or unwanted
(e.g. the accessor to get the permutation). Guarding that with asserts
isn't as good as guarding that at compile time...

Benoit


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