MATLAB 版本 7 中的 pdist2 等效项答案

【问题标题】：pdist2 equivalent in MATLAB version 7MATLAB 版本 7 中的 pdist2 等效项
【发布时间】：2011-10-08 12:18:06
【问题描述】：

我需要在 matlab 中计算 2 个矩阵之间的欧几里得距离。目前我正在使用 bsxfun 并计算距离如下（我附上了代码的 sn-p）：

for i=1:4754
test_data=fea_test(i,:);
d=sqrt(sum(bsxfun(@minus, test_data, fea_train).^2, 2));
end

fea_test 的大小是 4754x1024 和 fea_train 是 6800x1024 ，使用他的 for 循环导致 for 的执行需要大约 12 分钟，我认为这太高了。有没有办法更快地计算两个矩阵之间的欧式距离？

有人告诉我，通过删除不必要的 for 循环，我可以减少执行时间。我也知道 pdist2 可以帮助减少计算时间，但由于我使用的是 matlab 7 版，所以我没有 pdist2 函数。升级不是一种选择。

任何帮助。

问候，

佛系

【问题讨论】：

标签： matlab vectorization euclidean-distance bsxfun

【解决方案1】：

这里是计算欧几里得距离的矢量化实现，它比你所拥有的要快得多（甚至比我机器上的PDIST2 快得多）：

D = sqrt( bsxfun(@plus,sum(A.^2,2),sum(B.^2,2)') - 2*(A*B') );

基于以下事实：||u-v||^2 = ||u||^2 + ||v||^2 - 2*u.v

考虑以下两种方法之间的粗略比较：

A = rand(4754,1024);
B = rand(6800,1024);

tic
D = pdist2(A,B,'euclidean');
toc

tic
DD = sqrt( bsxfun(@plus,sum(A.^2,2),sum(B.^2,2)') - 2*(A*B') );
toc

在我运行 R2011b 的 WinXP 笔记本电脑上，我们可以看到 10 倍的时间改进：

Elapsed time is 70.939146 seconds.        %# PDIST2
Elapsed time is 7.879438 seconds.         %# vectorized solution

您应该知道，在最小精度下，它并没有完全给出与 PDIST2 相同的结果。通过比较结果，您会发现细微的差异（通常接近 eps浮点相对精度）：

>> max( abs(D(:)-DD(:)) )
ans =
  1.0658e-013

顺便说一句，我收集了大约 10 种不同的实现（有些只是彼此之间的微小变化），用于计算这种距离，并一直在对它们进行比较。与其他矢量化解决方案相比，您会惊讶于简单循环的速度有多快（感谢 JIT）...

【讨论】：

@Alex：谢谢。我看到您的解决方案基本相同，只是我避免使用 BSXFUN（而不是使用 REPMAT 等）在内存中创建临时矩阵
同意，bsxfun 是 MATLAB 中被低估和未知的函数之一。它不会像 repmat 那样浪费内存。更多用户需要学习如何在“矢量化”他们的代码时使用它。
快速说明一下，我比较了 PDIST2 和您使用 bsxfun 的方法，并且串行 bsxfun 的速度要快得多。但是，使用 DD 的 bsxfun 方法时，在 parfor 循环中的性能实际上变得更差了。
但是如果有 GPU 并且没有 for 循环，矢量化代码可以用来获得更好的性能......对吗？
@CharlieParker 我想，你必须进行基准测试才能看到.. 我的观点是 MATLAB 中的 for 循环在过去声名狼藉，但现在它们并没有那么糟糕。

【解决方案2】：

您可以通过重复fea_test 6800 次和fea_train 4754 次的行来完全向量化计算，如下所示：

rA = size(fea_test,1);
rB = size(fea_train,1);

[I,J]=ndgrid(1:rA,1:rB);

d = zeros(rA,rB);

d(:) = sqrt(sum(fea_test(J(:),:)-fea_train(I(:),:)).^2,2));

但是，这会导致大小为 6800x4754x1024 的中间数组（* 8 字节为双精度），这将占用约 250GB 的 RAM。因此，完全矢量化不起作用。

但是，您可以通过预分配来减少距离计算的时间，并且在必要之前不计算平方根：

rA = size(fea_test,1);
rB = size(fea_train,1);
d = zeros(rA,rB);

for i = 1:rA
    test_data=fea_test(i,:);
    d(i,:)=sum( (test_data(ones(nB,1),:) -  fea_train).^2, 2))';
end

d = sqrt(d);

【讨论】：

repmat 对性能没有好处。在这种情况下你最好内联：用test_data(ones(1,n8), :)替换repmat(test_data,nB,1)。
@Jones 感谢您的回复，但这并没有减少我的代码中的 for 循环数。我已经有一个for循环，在里面找到了这个循环，执行时间还是一样的。
@BhavyaPH：由于 RAM 不足，减少 for 循环数量的解决方案很可能无法在您的计算机上运行。另一种解决方案应该减少执行时间。您是否使用分析器进行比较？
由于 JIT，删除 for 循环并不总是最快的。请参阅stackoverflow.com/questions/7569368/…，以及关于 SO 的许多其他问题。

【解决方案3】：

试试这个矢量化版本，它应该非常有效。 编辑： 刚刚注意到我的回答与@Amro 的相似。

function K = calculateEuclideanDist(P,Q)
% Vectorized method to compute pairwise Euclidean distance
% Returns K(i,j) = sqrt((P(i,:) - Q(j,:))'*(P(i,:) - Q(j,:)))

[nP, d] = size(P);
[nQ, d] = size(Q);

pmag = sum(P .* P, 2);
qmag = sum(Q .* Q, 2);

K = sqrt(ones(nP,1)*qmag' + pmag*ones(1,nQ) - 2*P*Q');

end

【讨论】：