CODEWITHSUNDEEP
matlabmultidimensional-arraygaussiannormal-distributionmixture-model
Simon
Simon

Reputation: 5402

Loopless Gaussian mixture model in Matlab

I have several Gaussian distributions and I want to draw different values from all of them at the same time. Since this is basically what a GMM does, I have looked into Matlab GMM implementation (gmrnd) and I have seen that it performs a simple loop over all the components.

I would like to implement it in a faster way, but the problem is that 3d matrices are involved. A simple code (with loop) would be

n = 10; % number of Gaussians
d = 2; % dimension of each Gaussian
mu = rand(d,n); % init some means
U = rand(d,d,n); % init some covariances with their Cholesky decomposition (Cov = U'*U)
I = repmat(triu(true(d,d)),1,1,n);
U(~I) = 0;
r = randn(d,n); % random values for drawing samples

samples = zeros(d,n);
for i = 1 : n
    samples(:,i) = U(:,:,i)' * r(:,i) + mu(:,i);
end

Is it possible to speed it up? I do not know how to deal with the 3d covariances matrix (without using cellfun, which is much slower).

Upvotes: 0

Views: 169

Answers (2)

Divakar
Divakar

Reputation: 221504

Few improvements (hopefully are improvements) could be suggested here.

PARTE #1 You can replace the following piece of code -

I = repmat(triu(true(d,d)),[1,1,n]);
U(~I) = 0;

with bsxfun(@times,..) one-liner -

U = bsxfun(@times,triu(true(d,d)),U)

PARTE #2 You can kill the loopy portion of the code again with bsxfun(@times,..) like so -

samples = squeeze(sum(bsxfun(@times,U,permute(r,[1 3 2])),2)) + mu

Upvotes: 1

mikkola
mikkola

Reputation: 3476

I'm not fully convinced this is faster, but it gets rid of the loop. It would be interesting to see benchmarking results if you can do that. I also think this code makes is rather ugly and it's a bit hard to deduce what's going on, but I'll let you decide between readability and performance.

Anyway, I decided to define a big n*d dimensional Gaussian where each block d of variates are independent of each other (as in the original). This allows defining the covariance as a block diagonal matrix, for which I use blkdiag. From there, it is a matter of applying bsxfun to remove the need for looping.

Using the same random seed, I can recover the same samples as your code:

%// sampling with block diagonal covariance matrix
rng(1) %// set random seed
Ub = mat2cell(U, d, d, ones(n,1)); %// 1-by-1-by-10 cell of 2-by-2 matrices
C = blkdiag(Ub{:});
Ns = 1; %// number of samples
joint_samples = bsxfun(@plus, C'*randn(d*n, Ns), mu(:));
new_samples = reshape(joint_samples, [d n]); %// or [d n Ns] if Ns > 1

%//Compare to original
rng(1) %// set same seed for repeatability
r = randn(d,n); % random values for drawing samples
samples = zeros(d,n);
for i = 1 : n
    samples(:,i) = U(:,:,i)' * r(:,i) + mu(:,i);
end

isequal(samples, new_samples) %// true

Upvotes: 0

Related Questions