Re: [eigen] LDLt and LLt fixes

[ Thread Index | Date Index | More Archives ]

On Mon, Mar 30, 2009 at 7:22 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

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



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

Mail converted by MHonArc 2.6.19+