Home > marsbar > @mardo_5 > private > pr_spm_logdet.m

pr_spm_logdet

PURPOSE ^

returns the log of the determinant of positive semi-definite matrix C

SYNOPSIS ^

function [H] = pr_spm_logdet(C)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Wed 11-May-2022 16:26:09 by m2html © 2003-2019