I am currently trying to use the llt decomposition of cholesky matrix. It seems to work perfectly, but there is a problem on the matrixL. For exemple when using a n*n matrix, the matrixL(n,n) is wrong. Did somebody noticed it? -- Thanks.
Can you be more specific and provide an example. Thanks.
Oh, excuse me, Here is a c++ code example unsigned int n = 3; Eigen::MatrixXd m(n,n); m << 0.00157, 0.003914, 0.0885, 0.003914, 0.013876, 0.2304, 0.0885, 0.2304, 5; Eigen::MatrixXd L(m.llt().matrixL()); std::cout << L << std::endl <<std::endl; The matrixL is almost good, but the last element - L(3,3) makes the product wrong. Meaning that L*L^T differs from the original matrix. -- Thanks
I forgot, this is the output... 0.0396232 0 0 0.0987804 0.0641749 0 2.23354 0.152241 5
Your matrix is not positive definite (it's almost singular). You need to check the info() of the decomposition. If llt() does not work, sometimes ldlt() does (even if your matrix is indefinite). Would it be a good idea, if we asserted that info()==Success whenever matrixL() or solve() is called for a decomposition? Or are there cases where a getting the components of a failed decomposition gives useful values? I.e. would we break working code?
Thank You for answering, Christoph Indeed, I did not check if the matrix is positive definite. I will now try to use the ldlt() method.
(In reply to comment #4) > Would it be a good idea, if we asserted that info()==Success whenever matrixL() > or solve() is called for a decomposition? > Or are there cases where a getting the components of a failed decomposition > gives useful values? I.e. would we break working code? We have assertions checking m_isInitialized, but currently m_isInitialized is set to true even if the factorization failed. That seems a bad idea, and I'm pretty sure that some other factorizations raises an assert in solve() if the factorization failed, so I'd be OK to set m_isInitialized to true only in case of numerical success.
(In reply to comment #6) > We have assertions checking m_isInitialized, but currently m_isInitialized is > set to true even if the factorization failed. That seems a bad idea, and I'm > pretty sure that some other factorizations raises an assert in solve() if the > factorization failed, so I'd be OK to set m_isInitialized to true only in case > of numerical success. I once put an assert in EigenSolver to prevent users from accessing eigenvalues if the computation did not converge, but I removed it again after Benoit argued against it, on the ground that we should not have asserts that depend on round-off error (see discussion at bug 354).
(In reply to comment #7) > [...], on the ground that we should not have asserts that depend on > round-off error (see discussion at bug 354). Yes, you are right. I forgot about that. Still, at the moment there is no comfortable way to check for successful decompositions, when using the decomposition functions of MatrixBase, i.e. if you don't want to store the decomposition, which you only use once. The following would compute the decomposition twice, and it is actually not completely guaranteed that both llt-calls lead to the same result: assert(M.llt().info()==Eigen::Success); M.llt().solve(...);
(In reply to comment #8) > Still, at the moment there is no comfortable way to check for successful > decompositions, when using the decomposition functions of MatrixBase, [...] We could add an output argument to MatrixBase::llt(), so that the user write something like: int info; /* not sure of type here */ M.llt(&info).solve(...); if (info != Eigen::Success) { /* handle error */ } In this case, the Eigen library should guarantee that you can call solve() on a decomposition that not computed successfully. I think the syntax is slightly awkward in that it seems more logical to check info before calling .solve(), but that is not possible as far as I can see without declaring an LLT decomposition object. As mentioned on IRC, exceptions are another possibility, leading to code like: try { M.llt().solve(); } catch (SomeExceptionType e) { ... } To me, and as I understand from the IRC discussion also to Christoph, this looks more like clean C++ code. It makes sure that the computation does not continue with wrong values, and it would probably have prevented the issue that started this report. On the other hand, we have never used exceptions in the Eigen library, and I think there are some users that like to compile without exceptions. This means that we should probably introduce another compile-time switch, EIGEN_USE_EXCEPTIONS, with all the maintainability problems that causes.
(In reply to comment #9) > but that is not possible as far as I can see without declaring an LLT > decomposition object. One alternative (not sure if that's really better): M.llt().assertSuccess().solve(...); Maybe alternatively with `throw` instead of assert. That would avoid declaring the decomposition object, but still is likely to be forgotten. > As mentioned on IRC, exceptions are another possibility, leading to code like: > > try { > M.llt().solve(); > } > catch (SomeExceptionType e) > { ... } Yes, exactly. And there is the additional advantage that you can catch the exception at any place you want (e.g. outside the current method), and still don't calculate with undefined values. (Big difference vs. checking info()) > [...] This means that we should probably introduce another compile-time > switch, EIGEN_USE_EXCEPTIONS, with all the maintainability problems that > causes. Yes, that's a problem indeed. So basically, the question is if avoiding undefined behavior or lots of info() checks by the user is worth that effort. The other problem is, are there cases where people don't care about wrong numeric results, e.g. if they need to check their result anyways afterwards? And should this be the default, in order to not break working code?
As this is more of a meta-question now, I changed the title.
I'd like to have a decision on that. I think throwing an exception is preferable to having silent undefined behavior when decompositions fail. In many cases it can make usage much easier than checking the info() return type. We actually already have an (internal) `EIGEN_EXCEPTIONS` which checks if exceptions are available. But I think that we could invest a user switch EIGEN_USE_EXCEPTIONS to let the user decide if exceptions shall be thrown.
I'm not big fan of c++ exceptions, but I've to admit that this is one of the very rare case in Eigen (not to say the only one) where they might the best answer. The remaining open question would then be what should be the default behavior? (in other words, shall we add a EIGEN_NO_EXCEPTIONS or a EIGEN_USE_EXCEPTIONS). I'd advocate for EIGEN_USE_EXCEPTIONS to be consistent with the "no-exceptions-policy" we had so far.
I'm ok with letting the current behavior be the default. Bug 601 would actually be another case where exceptions could be thrown: If something fails due to bad input and not due to programming errors, I'd always prefer exceptions over assertions.
-- 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/631.