This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 241 - LDLT gives NAN result for linear solve using semi-definite matrix
LDLT gives NAN result for linear solve using semi-definite matrix
 Status: RESOLVED FIXED None Eigen Unclassified Cholesky (show other bugs) 3.0 x86 - 64-bit Linux --- Unknown Nobody 3.1 Show dependency tree / graph

 Reported: 2011-03-31 21:53 UTC by Etienne Vouga 2019-12-04 10:36 UTC (History) 4 users (show) chtz frank.heckel gael.guennebaud jitseniesen

Attachments

 Etienne Vouga 2011-03-31 21:53:17 UTC ```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 #include #include 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.``` Jitse Niesen 2011-04-01 11:35:35 UTC ```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.``` Christoph Hertzberg 2011-04-28 18:56:26 UTC ```(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 ...``` Etienne Vouga 2011-04-28 19:09:34 UTC ```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.``` Jitse Niesen 2011-09-11 07:36:57 UTC ```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).``` Gael Guennebaud 2014-07-02 23:10:51 UTC ```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?``` Christoph Hertzberg 2014-07-03 10:54:46 UTC ```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.``` Gael Guennebaud 2014-07-08 10:14:12 UTC ```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!``` Gael Guennebaud 2014-07-08 10:22:54 UTC ```Backport in 3.2: https://bitbucket.org/eigen/eigen/commits/619da854c7bd/ https://bitbucket.org/eigen/eigen/commits/d14370660187/``` Nobody 2019-12-04 10:36:45 UTC ```-- 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.```

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