how to sample a Gaussian given covariance or precision with symamd

gistfile1.matlab

% samples from N(mu, C) in case type=='cov', and from N(mu, C^(-1)) in case
% type=='prec'. z is optional and must be drawn from N(0,I), inds is a
% permutations vector that is applied to C before sampling (doesnt effect
% the output, just efficiency, smyamd is used as default)
function sample=gaussian_sampler_symamd(mu, C, type, z, inds)
if nargin<3
    type='prec';
end

if nargin <4
    z=randn(size(C,1),1);
end

if nargin < 5
    inds=symamd(C);
end

% inverse permutation
inds_inv(inds)=1:length(inds);

% permute input gaussian
z=z(inds);

% cholesky C_L'*C_L==C(inds,inds)
C_L=chol(C(inds,inds));

if strcmp(type, 'cov')
    % draw a sample from N(0, C(inds,inds))
    x=C_L'*z; 
else
    % draw a sample from N(0, C^(-1)(inds,inds))
    x=C_L\z;
end

% undo permutation, this is a sample from N(0, C^(-1)) or N(0, C), then add mean
sample=x(inds_inv)+mu;