returns the log of the determinant of positive semi-definite matrix C FORMAT [H] = pr_spm_logdet(C) H = log(det(C)) spm_logdet is a computationally efficient operator that can deal with sparse matrices _______________________________________________________________________ Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
0001 function [H] = pr_spm_logdet(C) 0002 % returns the log of the determinant of positive semi-definite matrix C 0003 % FORMAT [H] = pr_spm_logdet(C) 0004 % H = log(det(C)) 0005 % 0006 % spm_logdet is a computationally efficient operator that can deal with 0007 % sparse matrices 0008 %_______________________________________________________________________ 0009 % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience 0010 0011 % Karl Friston 0012 % $Id: spm_logdet.m 309 2005-11-24 16:24:04Z karl $ 0013 0014 % assume diagonal form 0015 %----------------------------------------------------------------------- 0016 TOL = 1e-8; % c.f. n*max(s)*eps 0017 n = length(C); 0018 s = diag(C); 0019 i = find(s > TOL & s < 1/TOL); 0020 C = C(i,i); 0021 H = sum(log(diag(C))); 0022 0023 % invoke det if non-diagonal 0024 %----------------------------------------------------------------------- 0025 w = warning; 0026 warning off 0027 [i j] = find(C); 0028 if any(i ~= j) 0029 n = length(C); 0030 a = exp(H/n); 0031 H = H + log(det(C/a)); 0032 end 0033 warning(w) 0034 0035 % invoke svd is rank deficient 0036 %----------------------------------------------------------------------- 0037 if imag(H) | isinf(H) 0038 s = svd(full(C)); 0039 H = sum(log(s(s > TOL & s < 1/TOL))); 0040 end