Loopless Gaussian mixture model in Matlab

Question

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).

Divakar · Accepted Answer · 2015-12-08T11:34:35.590

1

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

edited Dec 08 '15 at 11:34

answered Dec 08 '15 at 11:19

Divakar

218,885
19
262
358

Thanks! The sum has to be done along the first dimension, though. The speed up is significant with higher `n` (up to 10x before my computer runs out of memory), but if I increase also `d` then it's the other way around. How to choose between the loop and `bsxfun` really seems to depend on the situation, I don't think there is a general rule. Nonetheless, thanks again! I don't think there is another way to kill the loop beside your suggestion. – Simon Dec 08 '15 at 15:21

score 0 · Answer 2 · answered Dec 07 '15 at 18:43

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

Thanks for the reply! Unfortunately, your solution is slower (at least 10x, even more increasing `n` and `m`) and also runs out of memory sooner than the loop (because of the preallocation by `zeros(...)`). — Simon, Dec 07 '15 at 22:00

Loopless Gaussian mixture model in Matlab

2 Answers2