Skip to content

Instantly share code, notes, and snippets.

@redpony
Last active May 13, 2023 12:49
Show Gist options
  • Save redpony/fc8a0db6b20f7b1a3f23 to your computer and use it in GitHub Desktop.
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.
// 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;
}
@trb5me
Copy link

trb5me commented Sep 5, 2017

@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).

@Alaya-in-Matrix
Copy link

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