Computing log-determinant?

Question

Mathematica does-not have a function to compute the log-det of matrix? Naively computing Log[Det[M]] can be numerically unstable.

Answer 1

The determinant of a matrix may be outside the usual 64-bit floating point range, even if the log determinant isn't. Hence using Log[Det[...]] to compute the log determinant can be numerically problematic.

A simple method to directly compute the log-determinant of a symmetric positive definite matrix is:

logdet[m_] := 2 Total[Log[Diagonal[CholeskyDecomposition[m]]]];

(See the comments for dealing with more general matrices.) However, this approach is sensitive to slight asymmetries in the matrix, and near-singularity, so I use

   logdet2[m_] := Module[{mm}, 
      mm = (m + m\[Transpose])/2.;
      Do[mm[[i, i]] = mm[[i, i]] (1. + 10.^-8), {i, Length[m]}];
      2 Total[Log[Diagonal[CholeskyDecomposition[mm]]]] 
   ];

This is inelegant with the arbitrary small constant, but does the job for many applications.

A trick that I haven't tested well is to use Mathematica's internal algorithms using the multivariate normal log-likelihood:

logdet3[m_] := Module[{zero},
  zero = Table[0, Length[m]];
  -2 LogLikelihood[MultinormalDistribution[zero, m], {zero}] - Length[m] Log[2 \[Pi]]
  ]

It seems to be slower than the other methods.

To test these:

a = {{5, 2}, {2, 5}};
Log[Det[a]]
logdet[a] // Simplify
logdet2[a]
logdet3[a]

Using a = DiagonalMatrix[Table[10.^-10, {100}]];, Log[Det[a]] gives Indeterminate, while the other methods give the correct result -2302.59.