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.

logdet.cc

// 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;
}

Hi, thanks for your implementations. I see that for symmetric matrices the formulation can be understood from http://blogs.sas.com/content/iml/2012/10/31/compute-the-log-determinant-of-a-matrix.html But, what is the formal formulation for the "else case"?

Thanks again.

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

Line 11~14 can be rewritten as ld = 2 * chol.matrixL().toDenseMatrix().diagonal().array().log().sum()