在 MATLAB 中优化/向量化马氏距离计算

Question

我有以下一段 Matlab 代码，它通过多次迭代计算向量和矩阵之间的 Mahalanobis distances。我正在尝试通过矢量化找到一种更快的方法来执行此操作，但没有成功。

S.data=0+(20-0).*rand(15000,3);
S.a=0+(20-0).*rand(2500,3);

S.resultat=ones(length(S.data),length(S.a))*nan;
S.b=ones(length(S.a),3,length(S.a))*nan;

for i=1:length(S.data)
    for j=1:length(S.a)
         S.a2=S.a;
         S.a2(j,:)=S.data(i,:);
         S.b(:,:,j)=S.a2;
           if j==length(S.a)
              for k=1:length(S.a);
                   S.resultat(i,k)=mahal(S.a(k,:),S.b(:,:,k));
              end
           end    
    end   
end

我现在已经修改了代码并避免了其中一个循环。但它仍然很长。如果有人有想法，我将不胜感激！

S.data=0+(20-0).*rand(15000,3);
S.a=0+(20-0).*rand(2500,3);

S.resultat=ones(length(S.data),length(S.a))*nan;
   for i=1:length(S.data)
       for j=1:length(S.a)
       S.a2=S.a;
       S.a2(j,:)=S.data(i,:);
       S.resultat(i,j)=mahal(S.a(j,:),S.a2);    
       end   
   end

Answer 1

简介及解决代码

您可以将使用 mahal 的最内层循环替换为位向量化的内容，因为它使用了一些mahal.

的循环缩短和破解版本中的预计算值（在 bsxfun 的帮助下）

基本上你有一个 2D 数组，我们称它为 A 以便于参考，还有一个 3D 数组，我们称它为 B。将输出存储到变量 out 中。因此，可以提取最里面的代码片段并基于假定的变量名称。

原始循环码

for k=1:size(A,1)
    out(k)=mahal(A(k,:),B(:,:,k));
end

所以，我所做的是侵入 mahal.m 并寻找当输入为 2D 和 3D 时可以矢量化的部分。现在，mahal 在其中使用 qr，无法对其进行矢量化。因此，我们最终得到了一个 hacked 代码。

黑客代码

%// Pre-calculate certain values that could be avoided than using into loop
meanB = mean(B,1); %// mean of B along dim-1
B_meanB = bsxfun(@minus,B,meanB); %// B minus mean values of B
A_B_meanB = A' - reshape(meanB,size(B,2),[]); %//'# A minus B_meanB 

%// QR calculations in a for-loop starts until the output is obtained
for k = 1:size(A,1)
    [~,R] = qr(B_meanB(:,:,k),0);
    out2(k) = sum((R'\A_B_meanB(:,k)).^2)*(size(A,1)-1);
end

现在，要将这个 hack 解决方案扩展到问题代码，可以引入更多的调整来预先计算更多正在使用那些嵌套循环的值。

最终解法代码

A = S.a; %// Get data from S
[rx,cx] = size(A); %// Get size parameters
Atr = A'; %//'# Pre-calculate transpose of A

%// Pre-calculate replicated B and the indices to be modified at each iteration
B_rep = repmat(S.a,1,1,rx);
B_idx = bsxfun(@plus,[(0:cx-1)*rx + 1]',[0:rx-1]*(rx*cx+1)); %//'

out = zeros(size(S.data,1),rx); %// initialize output array
for i=1:length(S.data)

    B = B_rep;
    B(B_idx) = repmat(S.data(i,:)',1,rx); %//'
    meanB = mean(B,1); %// mean of B along dim-1

    B_meanB = bsxfun(@minus,B,meanB); %// B minus mean values of B
    A_B_meanB = Atr - reshape(meanB,3,[]); %// A minus B_meanB
    for jj = 1:rx
        [~,R] = qr(B_meanB(:,:,jj),0);
        out(i,jj) = sum((R'\A_B_meanB(:,jj)).^2)*(rx-1); %//'
    end

end
S.resultat = out;

基准测试

这是基准测试代码，用于将建议的解决方案与问题中列出的代码进行比较 -

%// Random inputs
S.data=0+(20-0).*rand(1500,3); %(size 10x reduced for a quicker runtime test)
S.a=0+(20-0).*rand(250,3);

S.resultat=ones(length(S.data),length(S.a))*nan;
disp('----------------------------- With original code')
tic

S.b=ones(length(S.a),3,length(S.a))*nan;
for i=1:length(S.data)
    for j=1:length(S.a)
        S.a2=S.a;
        S.a2(j,:)=S.data(i,:);
        S.b(:,:,j)=S.a2;
        if j==length(S.a)
            for k=1:length(S.a);
                S.resultat(i,k)=mahal(S.a(k,:),S.b(:,:,k));
            end
        end
    end
end

toc, clear i j S.a2 k S.resultat

S.resultat=ones(length(S.data),length(S.a))*nan;
disp('----------------------------- With proposed solution code')
tic

[ ... Proposed solution code ...]

toc

运行时 -

----------------------------- With original code
Elapsed time is 17.734394 seconds.
----------------------------- With proposed solution code
Elapsed time is 6.602860 seconds.

因此，我们可能会通过建议的方法和一些调整来绕过 2.7x 加速！

在 MATLAB 中优化/向量化马氏距离计算

Optimize/ Vectorize Mahalanobis distance calculations in MATLAB

optimization

performance

matlab

matrix

vectorization

简介及解决代码

基准测试