# Computing log-determinant?

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

• If speed isn't an issue, you could just do a SingularValueDecomposition, take the log of the singular values and take the Total Commented Feb 13, 2020 at 11:48
• Yes, speed (and memory) are an issue. Having to compute and store all the singular vectors is very wasteful.
– a06e
Commented Feb 13, 2020 at 12:01
• SingularValueList is about 2-3 times faster, and I would like to know too if there are significantly better ways to compute the log-determinant. Commented Feb 13, 2020 at 12:08
• LUDecomposition is much faster, but potentially less stable. So, something like Total[Log[Abs[Diagonal[First[LUDecomposition[m]]]]]] should work (this removes the sign of the determinant, though) Commented Feb 13, 2020 at 13:33
• For general $\mathbf M$, I don't expect that you can do better than $O(n^3)$ effort for $\log\det\mathbf M$, even if it's SPD. But if the matrix has some sort of exploitable structure, perhaps... Commented May 24, 2020 at 6:44

## 1 Answer

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.

• These seem OK for symmetric positive definite matrices. For general square matrices one might use logDet[mat_] := With[{evals = Eigenvalues[mat]}, Total[Log[evals]]] Commented Apr 10 at 18:23
• @DanielLichtblau some of the eigenvalues may be negative, and so applying the function Re[#] + I*Mod[Im[#], 2π, -π] & to the result may prettify it a bit. Commented Apr 10 at 18:31
• Yes, bringing the imaginary part into the correct range should also be done. Commented Apr 10 at 23:43
• Good points @DanielLichtblau and Roman, I qualified the answer to say it is for symmetric p.d. matrices. Commented Apr 11 at 10:11
• @DanielLichtblau (or others), the Do loop in logdet2 above seems to mess up the storing of the matrix in some way, causing CholeskyDecomposition to slow down enormously (for large matrices). I tried applying Developer`ToPackedArray but to no avail. Any ideas what is happening? Commented Apr 11 at 10:57