Re: [eigen] LDLt and LLt fixes

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

Thanks, Jitse and Keir, for your insights! I didn't know about
2x2-block-Cholesky but that seems to make sense.
Benoit

2009/3/31 Keir Mierle <mierle@xxxxxxxxx>:
> On Mon, Mar 30, 2009 at 7:22 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
> wrote:
>>
>> Bill,
>>
>> Thanks for the testcase. First, I just want to clarify that our LDLt
>> does diagonal pivoting which is all the pivoting that a Cholesky dec
>> can do.
>>
>> In our unit tests (which succeed) we have random positive and negative
>> semidef matrices, so I'm confident that our code handles them well.
>> Your matrix is:
>>
>> 0   -2
>> -2  1
>>
>> which has trace 1 and determinant -4, hence is indefinite (it has one
>> positive and one negative eigenvalue). So it is not surprising that
>> our LDLt fails on it: we only claim to support semidefinite matrices
>> (pos or neg).
>>
>> What is more surprising is how other libraries can do, I'd like to
>> know, but I wont be able to investigate this as I have a ton other
>> things to do in Eigen first !
>
> Note that the "Cholesky" of an indefinite matrix results in a block-diagonal
> D with 1x1 and 2x2  blocks. This is known as the block Cholesky. If we want
> to support that, it should be in an entirely different routine because it
> will ruin the performance of normal LDLt due to all the special casing.
>
> Also, in the case of a semidefinite matrix, a pivoting cholesky
> implementation will produce a PLDLPt decomposition such that PL has a
> triangular part that is zeros at the bottom.
>
> See page 2 of http://eprints.ma.man.ac.uk/1193/01/high90c.pdf
>
> Keir
>
>
>>
>>
>> Benoit
>>
>> 2009/3/30 Bill Greene <w.h.greene@xxxxxxxxx>:
>> > Benoit,
>> >
>> >    I think I accidentally sent you a test case that actually *requires*
>> > pivoting because the (0,0) term in the
>> > matrix is zero. So please ignore it if you like.
>> >
>> > Bill
>> >
>> > On Mon, Mar 30, 2009 at 9:12 PM, Bill Greene <w.h.greene@xxxxxxxxx>
>> > wrote:
>> >>
>> >> Benoit,
>> >>
>> >>    Wow, thanks for the quick work!
>> >>
>> >>    I just downloaded the latest code from svn and tried a quick test.
>> >> Surprisingly, it didn't work. After some
>> >> debugging and investigation, I think this may be an unusual case. My
>> >> other
>> >> test cases work correctly--
>> >> exactly as I expect. Here is the code for my failed case
>> >>
>> >> void bathe2DOF()
>> >> {
>> >>   const int n = 2;
>> >>   Eigen::Matrix2d K, M;
>> >>   K << 6., -2., -2., 4.;
>> >>   M << 2., 0., 0., 1.;
>> >>   const double omega2 = 3.0;
>> >>   Eigen::Matrix2d A = K - omega2*M;
>> >>   Eigen::LDLT<Eigen::Matrix2d> factA(A);
>> >>   std::cout << "Is invertible: " << factA.isInvertible() << std::endl;
>> >>   Eigen::Vector2d b, x;
>> >>   b << 0.0, 10.0;
>> >>   factA.solve(b, &x);
>> >>   prtMat(x, "x");
>> >> }
>> >>
>> >> The correct value for x, by the way, is [-5.0, 0].
>> >>
>> >> The lapack routines dsytrf and dsytrs compute this solution. I promised
>> >> to
>> >> provide more details of the lapack
>> >> solution-- and I intend to-- but I'm somewhat puzzled by the factored
>> >> matrix returned by dsytrf. Even though
>> >> the solution from dsytrs is correct, the factored matrix is not what I
>> >> expected-- or perhaps its just that the
>> >> documentation is so cryptic I don't know what I'm seeing. At any rate,
>> >> I'll provide more lapack details as
>> >> soon as I understand them.
>> >>
>> >> Bill
>> >>
>> >> On Mon, Mar 30, 2009 at 5:55 PM, Benoit Jacob
>> >> <jacob.benoit.1@xxxxxxxxx>
>> >> wrote:
>> >>>
>> >>> Hi,
>> >>>
>> >>> the thread with Bill inspired me... so now (trunk revision 947097),
>> >>> LDLt supports the negative semidefinite case.
>> >>>
>> >>> There's a couple of things to discuss.
>> >>>
>> >>> First, the floating-point comparisons in both LLt and LDLt had issues.
>> >>> LLt did
>> >>>
>> >>>   if(x < eps)
>> >>>
>> >>> and LDLt did
>> >>>
>> >>>   if(Djj <= 0)
>> >>>
>> >>> and replacing these by proper fuzzy comparisons (experimenting with a
>> >>> few different ones) resulted in a big improvement in the accuracy of
>> >>> some methods, like LDLt::rank() which used to be untested and had
>> >>> typically 3% imprecision, and now is exact in 95% of cases (However
>> >>> for an exact rank computation, it's still not nearly as reliable as
>> >>> the full-pivoting LU which can be explained by the different between
>> >>> bi-directional pivoting and the diagonal pivoting that's all one can
>> >>> do with Cholesky).
>> >>>
>> >>> The other issue is that LLt and LDLt have claims in their dox that
>> >>> they use only a triangular half of the matrix and don't need the other
>> >>> half to be initialized, well that's not currently the case, one can
>> >>> check that by uncommenting a line I added to the cholesky.cpp
>> >>> unit-test. Help fixing this welcome...
>> >>>
>> >>> Cheers,
>> >>> Benoit
>> >>>
>> >>>
>> >>
>> >
>> >
>>
>>
>
>