Calculating the log-determinant of a large sparse matrix

Like Brian Borchers says in his comment, the determinant can be determined from the diagonal entries of the Cholesky factor.

Theory:

Consider the Cholesky factorization $$R^TR = S^TMS,$$ where

• $R$ is an upper triangular matrix, and
• $S$ is a permutation matrix.

Then begin{align} text{det}(M) &= text{det}(SR^TRS^T) &= text{det}(S) text{det}(R^T) text{det}(R) text{det}(S^T) &= text{det}(R)^2. end{align}

Since $R$ is triangular, its eigenvalues are the diagonal entries, hence begin{align} text{log}~text{det}(M) &= 2 text{log}~text{det}(R) &= 2 sum_i text{log}(R_{ii}). end{align}

Code:

Matlab code to do this is, for example,

[R,~,S] = chol(M); 
log_det_chol = sum(2*log(diag(R)));

In python, the log determinant can be computed from the sparse cholesky factorization using the package scikits.sparse.cholmod. Paraphrasing from the linked page, the code to do this is:

from scikits.sparse.cholmod import cholesky
R = cholesky(M)
log_det_chol = R.logdet()

If you only want to consider the nonzero eigenvalues, you can just sum over the nonzero entries of the diagonal of $R$.

Test case:

I tested this against Matlab's build in det() function on the 2D finite difference Laplace operator with Dirichlet B.C.'s, and got exact agreement for matrices up to around 400-by-400, after which point using matlab's det() yielded Inf, whereas the cholesky method had no problem going up to million-by-million scale on my laptop. Here is the code, for reference:

n = 20; % M is n^2-by-n^2. Matlab's det() returns Inf for n=25+
M = gallery('poisson',n); 
[R,~,S] = chol(M); 
log_det_chol = sum(2*log(diag(R))), 
log_det = log(det(M))

Memory considerations:

Since this method is based on the Cholesky factorization, it inherits the speed and stability from that.

For large problems, the bottleneck using this method is probably going to be memory. In particular, using a Cholesky factorization with good pivoting (Matlab's is good), you can expect to need $O(s^2)$ memory, where $s$ is the number of degrees of freedom in the smallest "separator front" that partitions the degrees of freedom into distinct pieces.

For example, in a 2D problem with an n-by-n grid, the 1D separator front partitioning the domain into two pieces has $s = O(n)$ degrees of freedom. Thus amount of memory required is $O(n^2)$, the same order of magnitude as storing the right hand side.

A 3D n-by-n-by-n problem has a 2D separator front with $s = O(n^2)$ degrees of freedom, necessitating $O(n^4)$ memory.