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

Answer 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