0

我有几个高斯分布，我想同时从它们中提取不同的值。因为这基本上是 GMM 所做的，所以我研究了 Matlab GMM 实现（gmrnd），我看到它对所有组件执行了一个简单的循环。

我想以更快的方式实现它，但问题是涉及到 3d 矩阵。一个简单的代码（带循环）将是

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

有没有可能加快速度？我不知道如何处理 3d 协方差矩阵（不使用cellfun，这要慢得多）。

4

2 回答 2

1

这里可以提出一些改进（希望是改进）。

PARTE #1您可以替换以下代码 -

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

单线bsxfun(@times,..)-

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

PARTE #2你可以像这样再次杀死代码的循环部分-bsxfun(@times,..)

samples = squeeze(sum(bsxfun(@times,U,permute(r,[1 3 2])),2)) + mu
于 2015-12-08T11:19:38.717 回答
0

我不完全相信这会更快，但它摆脱了循环。如果可以的话，看看基准测试结果会很有趣。我还认为这段代码相当丑陋，很难推断发生了什么，但我会让你在可读性和性能之间做出决定。

无论如何，我决定定义一个大n*d维高斯，其中每个d变量块彼此独立（如原件）。这允许将协方差定义为块对角矩阵，我使用blkdiag. 从那里开始，申请bsxfun消除循环的需要是一个问题。

使用相同的随机种子，我可以恢复与您的代码相同的样本：

%// 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
于 2015-12-07T18:43:58.617 回答