Matlab 中的无环高斯混合模型
Loopless Gaussian mixture model in Matlab
我有几个高斯分布,我想同时从所有高斯分布中提取不同的值。由于这基本上是 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
,速度要慢得多)。
我不完全相信这会更快,但它摆脱了循环。如果可以的话,看到基准测试结果会很有趣。我还认为这段代码相当难看,很难推断出发生了什么,但我会让你在可读性和性能之间做出决定。
无论如何,我决定定义一个大 n*d
维高斯分布,其中每个块 d
变量彼此独立(与原始版本一样)。这允许将协方差定义为块对角矩阵,为此我使用 blkdiag
. From there, it is a matter of applying 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
这里可以提出一些改进(希望是改进)。
PARTE #1 你可以替换下面这段代码-
I = repmat(triu(true(d,d)),[1,1,n]);
U(~I) = 0;
with bsxfun(@times,..)
一行 -
U = bsxfun(@times,triu(true(d,d)),U)
PARTE #2 你可以 kill 代码的循环部分再次使用 bsxfun(@times,..)
像这样 -
samples = squeeze(sum(bsxfun(@times,U,permute(r,[1 3 2])),2)) + mu
我有几个高斯分布,我想同时从所有高斯分布中提取不同的值。由于这基本上是 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
,速度要慢得多)。
我不完全相信这会更快,但它摆脱了循环。如果可以的话,看到基准测试结果会很有趣。我还认为这段代码相当难看,很难推断出发生了什么,但我会让你在可读性和性能之间做出决定。
无论如何,我决定定义一个大 n*d
维高斯分布,其中每个块 d
变量彼此独立(与原始版本一样)。这允许将协方差定义为块对角矩阵,为此我使用 blkdiag
. From there, it is a matter of applying 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
这里可以提出一些改进(希望是改进)。
PARTE #1 你可以替换下面这段代码-
I = repmat(triu(true(d,d)),[1,1,n]);
U(~I) = 0;
with bsxfun(@times,..)
一行 -
U = bsxfun(@times,triu(true(d,d)),U)
PARTE #2 你可以 kill 代码的循环部分再次使用 bsxfun(@times,..)
像这样 -
samples = squeeze(sum(bsxfun(@times,U,permute(r,[1 3 2])),2)) + mu