|Summary:||log det for cholmod|
|Product:||Eigen||Reporter:||Joshua Pritikin <jpritikin>|
|Severity:||Feature Request||CC:||chtz, gael.guennebaud|
|Version:||3.3 (current stable)|
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 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.