New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
Bug 631 - (Optionally) throw an exception when using an unsuccessful decomposition
Summary: (Optionally) throw an exception when using an unsuccessful decomposition
Status: DECISIONNEEDED
Alias: None
Product: Eigen
Classification: Unclassified
Component: General (show other bugs)
Version: unspecified
Hardware: All All
: Normal API Change
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on: 407
Blocks: 3.4
  Show dependency treegraph
 
Reported: 2013-07-17 18:39 UTC by Steeveness
Modified: 2016-01-28 22:14 UTC (History)
4 users (show)



Attachments

Description Steeveness 2013-07-17 18:39:44 UTC
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.
Comment 1 Gael Guennebaud 2013-07-17 23:12:37 UTC
Can you be more specific and provide an example. Thanks.
Comment 2 Steeveness 2013-07-18 10:12:14 UTC
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
Comment 3 Steeveness 2013-07-18 10:13:40 UTC
I forgot, this is the output...

0.0396232         0         0
0.0987804 0.0641749         0
  2.23354  0.152241         5
Comment 4 Christoph Hertzberg 2013-07-18 10:31:06 UTC
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?
Comment 5 Steeveness 2013-07-18 10:56:14 UTC
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.
Comment 6 Gael Guennebaud 2013-07-18 11:23:10 UTC
(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.
Comment 7 Jitse Niesen 2013-07-18 11:49:03 UTC
(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).
Comment 8 Christoph Hertzberg 2013-07-18 12:55:16 UTC
(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(...);
Comment 9 Jitse Niesen 2013-07-18 14:30:38 UTC
(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.
Comment 10 Christoph Hertzberg 2013-07-18 16:31:58 UTC
(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?
Comment 11 Christoph Hertzberg 2013-07-18 16:34:35 UTC
As this is more of a meta-question now, I changed the title.
Comment 12 Christoph Hertzberg 2014-06-15 04:28:30 UTC
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.
Comment 13 Gael Guennebaud 2014-06-16 14:13:12 UTC
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.
Comment 14 Christoph Hertzberg 2014-06-26 14:45:41 UTC
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.

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