Last active
May 13, 2023 12:49
-
-
Save redpony/fc8a0db6b20f7b1a3f23 to your computer and use it in GitHub Desktop.
Computing log(M.determinant()) in Eigen C++ is risky for large matrices since it may overflow or underflow. This gist uses LU (or, if applicable, Cholesky) decompositions to do the risky components in the log space.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// set use_cholesky if M is symmetric - it's faster and more stable | |
// for dep paring it won't be | |
template <typename MatrixType> | |
inline typename MatrixType::Scalar logdet(const MatrixType& M, bool use_cholesky = false) { | |
using namespace Eigen; | |
using std::log; | |
typedef typename MatrixType::Scalar Scalar; | |
Scalar ld = 0; | |
if (use_cholesky) { | |
LLT<Matrix<Scalar,Dynamic,Dynamic>> chol(M); | |
auto& U = chol.matrixL(); | |
for (unsigned i = 0; i < M.rows(); ++i) | |
ld += log(U(i,i)); | |
ld *= 2; | |
} else { | |
PartialPivLU<Matrix<Scalar,Dynamic,Dynamic>> lu(M); | |
auto& LU = lu.matrixLU(); | |
Scalar c = lu.permutationP().determinant(); // -1 or 1 | |
for (unsigned i = 0; i < LU.rows(); ++i) { | |
const auto& lii = LU(i,i); | |
if (lii < Scalar(0)) c *= -1; | |
ld += log(abs(lii)); | |
} | |
ld += log(c); | |
} | |
return ld; | |
} |
Line 11~14 can be rewritten as ld = 2 * chol.matrixL().toDenseMatrix().diagonal().array().log().sum()
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
@lood338 you can see the justification in the link above. The lower diagonal matrix L is easy to take the determinant (and log determinant of), but it is not the target matrix M. It is the square root of the target matrix (because M = LL^T).