New user self-registration is currently disabled. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Bug 1095 - log det for cholmod
log det for cholmod
Status: RESOLVED FIXED
Product: Eigen
Classification: Unclassified
Component: Sparse
3.3 (current stable)
All All
: Normal Feature Request
Assigned To: Nobody
:
Depends on:
Blocks:
  Show dependency treegraph
 
Reported: 2015-10-21 08:55 UTC by Joshua Pritikin
Modified: 2016-01-22 15:06 UTC (History)
2 users (show)



Attachments
proposed change (6.17 KB, patch)
2015-12-27 20:24 UTC, Joshua Pritikin
no flags Details | Diff
proposed change (5.73 KB, patch)
2015-12-27 20:26 UTC, Joshua Pritikin
no flags Details | Diff
implement complex and address feedback (3.61 KB, patch)
2015-12-29 14:19 UTC, Joshua Pritikin
no flags Details | Diff
unremove tests (4.01 KB, patch)
2015-12-29 14:24 UTC, Joshua Pritikin
no flags Details | Diff
add log det (4.48 KB, patch)
2016-01-13 16:58 UTC, Joshua Pritikin
no flags Details | Diff

Description Joshua Pritikin 2015-10-21 08:55:36 UTC
Eigen::LDLT has an API vectorD from which it is easy to obtain the log determinant of the Cholesky factor. Eigen::CholmodDecomposition offers no similar API.

I found an implementation in python, https://github.com/njsmith/scikits-sparse/blob/master/scikits/sparse/cholmod.pyx (search for "def D(self)").
I plan translate this into C++. Once I get it working, I would be happy to contribute it. Would this be a welcome addition to Eigen?
Comment 1 Gael Guennebaud 2015-10-22 08:43:45 UTC
Sure, patch welcome :)
Comment 2 Joshua Pritikin 2015-12-27 20:24:31 UTC
Created attachment 636 [details]
proposed change

How does this look?
Comment 3 Joshua Pritikin 2015-12-27 20:26:32 UTC
Created attachment 637 [details]
proposed change

How does this look?
Comment 4 Gael Guennebaud 2015-12-28 10:41:59 UTC
Thanks for the patch, some comments:

- The documentation is a bit misleading, and I'd suggest something like:

/** \returns the determinant of the underlying matrix from the current factorization */

- I guess that int* should be replaced by StorageIndex* to support 64 bits integers.

- complex scalar type can easily be supported by observing that the diagonal entries are real. Use numext::real(val) to extract the real part (works on reals too).


- the nested for loop can probably be replaced by something like:

logDet += Array<Scalar,Dynamic,Dynamic>::Map(x+px[sx], nrows, ncols).diagonal().real().log().sum();


- Calls like std::bla(x) should be replaced by:

using std::bla;
...
bla(x)

to support custom scalar types.
Comment 5 Joshua Pritikin 2015-12-28 15:21:24 UTC
The int pointers point into the cholmod_factor type. This is a C structure provided by cholmod, not Eigen. I suspect it only implements integers.

I will address your other suggestions.
Comment 6 Gael Guennebaud 2015-12-29 11:27:38 UTC
I think that they can be long int if cholmod_sparse::itype == CHOLMOD_LONG, which is why those pointer members are void* and not int*.
Comment 7 Joshua Pritikin 2015-12-29 14:19:18 UTC
Created attachment 639 [details]
implement complex and address feedback

> the nested for loop can probably be replaced by something like:

If you look closely, we don't add up the diagonal of a dense matrix but entries in the first row. Maybe it is possible to use a Map'd Eigen type with TriangularView, but I didn't research it.
Comment 8 Joshua Pritikin 2015-12-29 14:24:06 UTC
Created attachment 640 [details]
unremove tests
Comment 9 Gael Guennebaud 2015-12-29 14:43:22 UTC
right, I'have been confused by the +1, then it's even simpler:

Map<const Matrix<Scalar,1,Dynamic>, 0, InnerStride<> >(x+px[sx], ncols, InnerStride<>(nrows+1))
Comment 10 Joshua Pritikin 2016-01-13 16:58:50 UTC
Created attachment 644 [details]
add log det

Use fancy Eigen::Map instead of boring loop
Comment 11 Gael Guennebaud 2016-01-22 15:06:09 UTC
thanks for the patch.

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