New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Bug 1268 - Add LDLT with block diagonal
Summary: Add LDLT with block diagonal
Status: CONFIRMED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Cholesky (show other bugs)
Version: unspecified
Hardware: All All
: Normal Feature Request
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks:
 
Reported: 2016-07-30 16:00 UTC by w.h.greene
Modified: 2016-12-05 13:06 UTC (History)
4 users (show)



Attachments
simple test showing failure (1.03 KB, text/plain)
2016-07-30 16:00 UTC, w.h.greene
no flags Details
Detect faillure in LDLT (8.87 KB, patch)
2016-08-23 20:46 UTC, Gael Guennebaud
no flags Details | Diff

Description w.h.greene 2016-07-30 16:00:16 UTC
Created attachment 719 [details]
simple test showing failure

I have a simple symmetric matrix that I believe LDLT should solve. 
Instead it produces an incorrect solution. As the example code shows,
PartialPivLU produces a correct solution.

The incorrect solution is not reported; as can be seen info()
reports success.

Based on debugging, I can see that realAkk is flagged as zero during
the decomposition but this should not be the case. Note that
LU reports full-rank for this matrix.

In cases where a singularity is *really* encountered, I believe that
info() should report this rather than failing silently.
Comment 1 Christoph Hertzberg 2016-07-30 18:12:34 UTC
Ok, first of all, your matrix is indefinite -- you see this if you make an eigenvalue decomposition. Unfortunately, our ldlt decomposition often fails in that case.
Also realAkk==0 can happen for full-rank (but indefinite) matrices.
The easiest example is [0 1; 1 0]

I agree however, that silently returning Success in that case is not a good idea.
Perhaps, if realAkk==0 happens once, but we get abs(realAkk)>0 afterwards, declare the decomposition as failed (I have not thought through all cases ...)

Ideally, we'd support indefinite matrices as well, of course (would require allowing 2x2 blocks on the diagonal, i.e. an API change, though).
Comment 2 Gael Guennebaud 2016-08-23 17:24:10 UTC
This is rather a request for another LDLT variant.
Comment 3 Christoph Hertzberg 2016-08-23 17:29:09 UTC
I agree, but we definitely must not return Success, if the decomposition actually failed.
Comment 4 Gael Guennebaud 2016-08-23 20:46:28 UTC
Created attachment 722 [details]
Detect faillure in LDLT

right, I've been too fast sorry.

Here is an initial attempt. With this patch, LDLT returns NumericalIssue for a zero matrix, though the factorization is exact. To fix this we would need to explicitly check for a zero input in the case of a zero diagonal.
Comment 5 Gael Guennebaud 2016-08-23 21:17:41 UTC
ok, finally I've implemented a special path for a zero diagonal, and pushed the changes:

https://bitbucket.org/eigen/eigen/commits/c6f5c3c6ef20/
Summary:     Bug 1268: detect failure in LDLT and report them through info()
Comment 6 Allan Leal 2016-12-05 13:06:28 UTC
Out of curiosity, are there any plans on implementing LDLT solvers for indefinite matrices (i.e. when D has 1x1 and 2x2 blocks depending on the signs of the eigenvalues)?

Note You need to log in before you can comment on or make changes to this bug.