Bugzilla – Bug 1095
log det for cholmod
Last modified: 2016-01-22 15:06:09 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?
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.