Log determinants frequently need to be computed as part of Bayesian inference in models with Gaussian priors or likelihoods. They can be quite troublesome to compute; done through the Cholesky decomposition, it’s an operation. It’s pretty much infeasible to compute the exact log det for arbitrary matrices > 10,000 x 10,000. Models with that many parameters frequently occur in non-parametric regression methods like Gaussian processes.

Zhang and Leithead (2007) offer an update to the power series expansion Monte Carlo method of Barry & Pace (1999). The old paper expanded log det(A) into a power series. All the terms in the expansion involve matrix traces; instead of computing the traces explicitly, they are approximated through a Monte Carlo method, which involves taking random matrix vector products. Zhang and Leithead expand upon the method by decreasing the number of iid Gaussian vectors needed to 1, not requiring any restrictions on the eigenvalues of the matrix A, and using extrapolation so the power series can be truncated earlier. They claim their method is , which is potentially much better than the Cholesky version for large N.

There’s Matlab code for the method in the appendix of the paper. Unfortunately, it’s mangled in the pdf and it’s missing a required function. For convenience, here is a cleaned up version:

function logdet=mcLogDet(C,pp,kk,beta)
% input arguments: covariance matrix C,
% (optional) the number of trace seeds selection pp,
% (optional) the number of power-series terms kk,
% (optional) the ensemble adjusting factor beta.
if nargin < 2
    pp = 20;
end
if nargin < 3
    kk = 50;
end
if nargin < 4
    beta = 1;
end

N = size(C,1);
Cn=norm(C,inf);
A=C/Cn;
if issparse(C)
    B=speye(N)-A;
else
    B=eye(N)-A;
end

clear C A;
traceBt=trace(B);
traceB2t=traceAxB(B,B); % a simple routine to get trace(B 2 ) with N 2 operations
deltaB122=inf;
for iii=1:pp % trace seeds selection for pp times
    xxt=randn(N,1);
    cct=B*xxt;
    xxBNxx=xxt'*cct;
    xxTxxt=xxt'*xxt;
    traceBe=N*xxBNxx/xxTxxt;
    deltaBt=traceBt-traceBe;
    cct=B*cct;
    xxBNxx=xxt'*cct;
    traceB2ee=N*xxBNxx/xxTxxt;
    deltaB2e=traceB2t-traceB2ee;
    deltaB122new=deltaBt+deltaB2e/2;
    if abs(deltaB122new)<=abs(deltaB122), % have found a better seed
        xx=xxt;
        cc=cct;
        xxTxx=xxTxxt; % save the workspace related to the new seed
        traceB1e=traceBe;
        traceB2e=traceB2ee;
        deltaB=deltaBt;
        deltaB2=deltaB2e;
        deltaB122=deltaB122new;
    end
end
clear xxt cct;
rhotr=deltaB2/deltaB;
traceB=traceB1e; traceB2=traceB2e;
if abs(rhotr)>=1, % disp(‘trace rho greater than or equal to 1?');
    compens2=deltaB/(1-rhotr/2); % a more robust yet loose trace-seed compensation asin Remark 3
elseif rhotr==0,
    compens2=0; % no need compensation as i = 0
else
    compens2=-deltaB*log(1-rhotr)/rhotr; % the standard trace-seed compensation as in Remark 3
end;
vv=traceB+traceB2/2;
preadd=traceB2; % save the computation of the first two power-series terms
compens1=0;
wthla1=0;
wthla2=0;
wthla3=0;
for jjj=3:kk, % start the computation from the third power-series term to the kth term
    cc=B*cc;
    xxBNxx=xx'*cc;
    newadd=N*xxBNxx/xxTxx;
    if preadd==0, % only for numerical-robustness consideration
        compens1=0; wthl=vv-compens1-compens2;
        logdet2=N*log(Cn); logdet=-wthl+logdet2;
        return
    end
    lambda=newadd/preadd;
    if lambda<1,
        sumlambdaik=0; lambdai=1; ndt5=0;
        for iii=1:jjj, % computing the second term of the truncation-error compensation as in Proposition 5
            lambdai=lambdai*lambda;
            ndt5=ndt5+lambdai/(jjj+iii);
            sumlambdaik=sumlambdaik+lambdai/iii;
        end
        if lambdai==0,
            compens1=-newadd*ndt5; % truncation-error compensation via (8)
        else
            compens1=newadd*(log(1-lambda)+sumlambdaik)/lambdai;
            % the final truncation-error compensation as in Proposition 5
        end
    end;
    preadd=newadd;
    newadd=newadd/jjj;
    vv=newadd+vv;
    wthl=vv-compens1-compens2; % for the finite sequence with the ensemble effect as in Proposition 6
    wthla1=wthla2; wthla2=wthla3; wthla3=wthl;
    % save the last three terms of such a sequence, i.e, k-2 , k-1 and k
end;
if lambda>=1,
% disp(‘final lambda greater than or equal to 1?');
    compens1=0; % only for numerical-robustness consideration
end
    aaa=wthla1;
    bbb=wthla2-wthla1;
    qqq=(wthla3-wthla2)/bbb;
if bbb==0,
    % disp(‘The last three terms of Gamma-sequence are equal: converge well');
else
    if abs(qqq)>= 1, % only for numerical-robustness consideration
        % disp(‘Gamma-sequence coefficient greater than or equal to 1?');
    else
        wthl=aaa+bbb*qqq/(1-qqq); % Geometric series based re-estimation as in Proposition 6
    end
end
logdet2=N*log(Cn);
logdet=-wthl+logdet2;
if logdet/max(abs([wthl,logdet2]))<1e-2
    wthl=wthl*1;
elseif bbb/aaa>1e-4 % use beta to adjust the final ensemble effect if necessary
    wthl=wthl*beta; % in general, beta:=1
end
logdet=-wthl+logdet2; % return the final logdet approximation
clear B; % the routine is thus complete
end

function r = traceAxB(A,B)
    r = sum(sum(A'.*B,1),2);
end

And some sample code that shows how to use it:

thestd = 0.2;
A = eye(1000);
A = A + randn(1000)*thestd;
A = A'*A;

L = chol(A);
thelogdet = 2*sum(log(diag(L)))
thelogdeta = mcLogDet(A,20,50,1)

In the paper, the authors used (10,30,1) as the parameters of the method by I prefer (20,50,1) for safety. Mind you, the method is not super precise, they say about 10%, but for model comparison it’s more than precise enough.