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

[ Thread Index | Date Index | More Archives ]

Hi Ben,

Thanks for your great test-case. It helped a lot. Thanks also for the
patch, but i didn't apply it, because the problems were deeper:

Normally, when we terminate without 'finalizing' the matrix, that's
because our termination criterion tells us that the remaining block of
the matrix is approximately zero anyway.

Unfortunately, there were 2 bugs in our code. Fixing either of them
was enough to make your test-case happy.

The first bug was that our second termination criterion didn't
actually say that the remaining block is zero. So it didn't tell us
that we could safely terminate. The fix was to quit using it as a
termination criterion, and only use it as a criterion to do the
division by Djj.

The second bug was that our cutoff was too big. It came from Higham's
formula that tells how to choose a cutoff for optimal rank
determination, but here in LDLt we're not even trying to be
rank-revealing. So a small cutoff (using machine precision) is better,
it gives better accuracy as it prevents triggering termination
criteria too early.

Both bugs fixed in the development branch.


2010/2/22 Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>:
> Hi,
> Thank you for creating eigen. It is great.
> I mentioned a bug about ldlt() with rank-deficient matrices and gave a
> test case here:
> I believe the problem is that when ldlt() decides to abort due to a
> diagonal entry of D being too close to zero, it does not change the
> remaining diagonal entries to be 0.0 . As a result, when the original
> matrix has a nullity of 2 or more, multiplying the factors together
> does not result in the original matrix.
> The attached patch fixes this issue for me in the test case by zeroing
> out the relevant diagonal entries. It applies to the development
> branch. I will reply to this message with a patch and a question about
> the stable branch.
> Thank you,
> Ben

Mail converted by MHonArc 2.6.19+