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?
Sure, patch welcome :)
Created attachment 636 [details]
How does this look?
Created attachment 637 [details]
How does this look?
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:
to support custom scalar types.
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.
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*.
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.
Created attachment 640 [details]
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))
Created attachment 644 [details]
add log det
Use fancy Eigen::Map instead of boring loop
thanks for the patch.
-- 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.