计算多元高斯变化均值的有效方法 - Matlab

Question

是否有一种有效的方法来计算返回矩阵 p 的多元高斯（如下所示），即利用某种矢量化？我知道矩阵 p 是对称的，但对于大小为 40000x3 的矩阵来说，这仍然需要相当长的时间。

Matlab 代码示例：

DataMatrix = [3 1 4; 1 2 3; 1 5 7; 3 4 7; 5 5 1; 2 3 1; 4 4 4];

[rows, cols ] = size(DataMatrix);

I = eye(cols);
p = zeros(rows);

for k = 1:rows

    p(k,:) = mvnpdf(DataMatrix(:,:),DataMatrix(k,:),I);

end

Answer 1

第一阶段：破解源代码

我们正在迭代地执行mvnpdf(DataMatrix(:,:),DataMatrix(k,:),I)

语法是：mvnpdf(X,Mu,Sigma)。

因此，与我们输入的对应关系变为：

X = DataMatrix(:,:);
Mu = DataMatrix(k,:);
Sigma = I

对于与我们的情况相关的大小，源代码mvnpdf.m 缩减为-

%// Store size parameters of X
[n,d] = size(X);

%// Get vector mean, and use it to center data
X0 = bsxfun(@minus,X,Mu);

%// Make sure Sigma is a valid covariance matrix
[R,err] = cholcov(Sigma,0);

%// Create array of standardized data, and compute log(sqrt(det(Sigma)))
xRinv = X0 / R;
logSqrtDetSigma = sum(log(diag(R)));

%// Finally get the quadratic form and thus, the final output
quadform = sum(xRinv.^2, 2);
p_out = exp(-0.5*quadform - logSqrtDetSigma - d*log(2*pi)/2)

现在，如果 Sigma 始终是一个单位矩阵，我们也会将 R 作为单位矩阵。因此，X0 / R 与X0 相同，保存为xRinv。所以，本质上是quadform = sum(X0.^2, 2);

因此，原始代码 -

for k = 1:rows
    p(k,:) = mvnpdf(DataMatrix(:,:),DataMatrix(k,:),I);
end

减少到 -

[n,d] = size(DataMatrix);
[R,err] = cholcov(I,0);
p_out = zeros(rows);
K = sum(log(diag(R))) + d*log(2*pi)/2;
for k = 1:rows  
    X0 = bsxfun(@minus,DataMatrix,DataMatrix(k,:));     
    quadform = sum(X0.^2, 2);
    p_out(k,:) = exp(-0.5*quadform - K);
end

现在，如果输入矩阵的大小为40000x3，您可能想在这里停下来。但在系统资源允许的情况下，您可以将所有内容矢量化，如下所述。

第 2 阶段：向量化所有内容

现在我们看到了实际发生的情况并且计算看起来可以并行化，是时候与他的好朋友 permute 一起在 3D 中使用 bsxfun 来实现矢量化解决方案，就像这样 -

%// Get size params and R
[n,d] = size(DataMatrix);
[R,err] = cholcov(I,0);

%// Calculate constants : "logSqrtDetSigma" and  "d*log(2*pi)/2`"
K1 = sum(log(diag(R)));
K2 = d*log(2*pi)/2;

%// Major thing happening here as we calclate "X0" for all iterations
%// in one go with permute and bsxfun
diffs = bsxfun(@minus,DataMatrix,permute(DataMatrix,[3 2 1]));

%// "Sigma" is an identity matrix, so it plays no in "/R" at "xRinv = X0 / R".
%// Perform elementwise squaring and summing rows to get vectorized "quadform"
quadform1 = squeeze(sum(diffs.^2,2))

%// Finally use "quadform1" and get vectorized output as a 2D array
p_out = exp(-0.5*quadform1 - K1 - K2)