This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 1095 - log det for cholmod
Summary: log det for cholmod
Status: RESOLVED FIXED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Sparse (show other bugs)
Version: 3.3 (current stable)
Hardware: All All
: Normal Feature Request
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks:
 
Reported: 2015-10-21 08:55 UTC by Joshua Pritikin
Modified: 2019-12-04 15:05 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.
Comment 12 Nobody 2019-12-04 15:05:25 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/1095.

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