Bugzilla – Bug 1268
Add LDLT with block diagonal
Last modified: 2016-12-05 13:06:28 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()
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.
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).
This is rather a request for another LDLT variant.
I agree, but we definitely must not return Success, if the decomposition actually failed.
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.
ok, finally I've implemented a special path for a zero diagonal, and pushed the changes:
Summary: Bug 1268: detect failure in LDLT and report them through info()
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)?