LDLT is suppose to work for semi-definite symmetric matrices, but in some cases it returns a NAN result even when the system has a solution. Example code: #include <Eigen/Cholesky> #include <Eigen/Core> #include <iostream> using namespace Eigen; using namespace std; int main() { MatrixXd A(2,2); A << 1,1,1,1; VectorXd rhs(2); rhs << 0,0; VectorXd result = A.ldlt().solve(rhs); cout << result << endl; } Output: -nan -nan Expected Output: a -a, for any real number a.
I can confirm this bug. The problem is that the diagonal matrix D has zeros on the diagonal. In solve() we invert this matrix and this produces the NaN. Assuming that the documentation is right and that we do want to allow for the matrix being singular, one possibility would be to using the pseudo-inverse of D instead of the inverse (the pseudo-inverse of a nonzero scalar x is just 1/x but the pseudo-inverse of 0 is 0). This may kill vectorization in this step, but it's an O(N) operation so it shouldn't be a problem.
(In reply to comment #1) > Assuming that the documentation is right and that we do want to allow for the > matrix being singular, one possibility would be to using the pseudo-inverse of > D instead of the inverse I wouldn't expect the pseudo-inverse here, e.g. for rhs == Vector2d(1,0) the above example would give result == {1, -1}, which definitely is not a solution. Hm, ok I just re-read the documentation of it, and this seems to be the intended behavior -- question would be whether to change the implementation or the documentation ...
If the rhs is not in the column space of the matrix, like in your example, then clearly there is no reasonable possible output, so returning {1, -1} is not all that unreasonable. For now in my code I've replaced ldlt with fullPivLu, which seems to handle semidefinite matrices OK.
That seems to be the general behaviour of the solve() functions in Eigen: if there is a unique solution, try to compute it; if there are multiple solutions, try to compute one of them; if there are no solutions, return rubbish. Thus I implemented the solution with the pseudo-inverse in changesets 4d1457a1912e (dev branch) and c9e6c3df2893 (stable branch).
Using the biggest diagonal entries as threshold values to determine numerical "zeros" does not seem to be a good idea after all. Thus the following change: https://bitbucket.org/eigen/eigen/commits/f240bbe84f94/ Changeset: f240bbe84f94 User: ggael Date: 2014-07-02 23:04:46 Summary: LDLT is not rank-revealing, so we should not attempt to use the biggest diagonal elements as thresholds. The new unit tests with matrices having a large spectrum do not work without these changes. See also this thread where SimplicialLDLT (which does not perform pivoting) worked while the dense LDLT did not: http://forum.kde.org/viewtopic.php?f=74&t=121832#p314467 OK to backport to 3.2?
I saw cholesky_6 fail in 3 out of 5 cases: http://manao.inria.fr/CDash/testSummary.php?project=1&name=cholesky_6&date=2014-07-03 The fails all appear to be in the new unit test code. Maybe the spectra where too wide, or the required accuracy is too strict? Once the test passes reliably, I have no objections against backporting.
The problem was in LDLT itself: https://bitbucket.org/eigen/eigen/commits/ad2cfd97ef35/ Changeset: ad2cfd97ef35 User: ggael Date: 2014-07-08 10:04:27 Summary: Fix LDLT with semi-definite complex matrices: owing to round-off errors, the diagonal was not real. Also exploit the fact that the diagonal is real in the rest of LDLT A good reason to backport it!
Backport in 3.2: https://bitbucket.org/eigen/eigen/commits/619da854c7bd/ https://bitbucket.org/eigen/eigen/commits/d14370660187/
-- GitLab Migration Automatic Message -- This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity. You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/241.