Re: [eigen] LDLt and LLt fixes

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

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 !

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