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

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

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.

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

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

yes I'd much prefer the other way, such code duplication is prone to bugs.

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

You should also make the solve function skips the transpositions when
no pivoting has been computed.

cheers,
gael