在 Matlab 中有效地计算两两平方欧几里得距离答案

【问题标题】：Efficiently compute pairwise squared Euclidean distance in Matlab在 Matlab 中有效地计算两两平方欧几里得距离
【发布时间】：2014-07-17 16:18:45
【问题描述】：

给定两组d维点。如何在 Matlab 中最有效地计算 成对平方欧式距离矩阵？

符号： 第一组由(numA,d)-矩阵A 给出，第二组由(numB,d)-矩阵B 给出。生成的距离矩阵的格式应为(numA,numB)。

示例点：

d = 4;            % dimension
numA = 100;       % number of set 1 points
numB = 200;       % number of set 2 points
A = rand(numA,d); % set 1 given as matrix A
B = rand(numB,d); % set 2 given as matrix B

【问题讨论】：

你看过pdist2函数吗？ mathworks.com/help/stats/pdist2.html
@rayryeng 是的，请在我的回答中查看我的评估部分:)

标签： performance matlab matrix distance euclidean-distance

【解决方案1】：

这里通常给出的答案是基于bsxfun（参见例如[1]）。我提出的方法基于矩阵乘法，结果证明比我能找到的任何可比算法都要快：

helpA = zeros(numA,3*d);
helpB = zeros(numB,3*d);
for idx = 1:d
    helpA(:,3*idx-2:3*idx) = [ones(numA,1), -2*A(:,idx), A(:,idx).^2 ];
    helpB(:,3*idx-2:3*idx) = [B(:,idx).^2 ,    B(:,idx), ones(numB,1)];
end
distMat = helpA * helpB';

请注意： 对于常量d，可以通过硬编码实现替换for-loop，例如

helpA(:,3*idx-2:3*idx) = [ones(numA,1), -2*A(:,1), A(:,1).^2, ... % d == 2
                          ones(numA,1), -2*A(:,2), A(:,2).^2 ];   % etc.

评估：

%% create some points
d = 2; % dimension
numA = 20000;
numB = 20000;
A = rand(numA,d);
B = rand(numB,d);

%% pairwise distance matrix
% proposed method:
tic;
helpA = zeros(numA,3*d);
helpB = zeros(numB,3*d);
for idx = 1:d
    helpA(:,3*idx-2:3*idx) = [ones(numA,1), -2*A(:,idx), A(:,idx).^2 ];
    helpB(:,3*idx-2:3*idx) = [B(:,idx).^2 ,    B(:,idx), ones(numB,1)];
end
distMat = helpA * helpB';
toc;

% compare to pdist2:
tic;
pdist2(A,B).^2;
toc;

% compare to [1]:
tic;
bsxfun(@plus,dot(A,A,2),dot(B,B,2)')-2*(A*B');
toc;

% Another method: added 07/2014
% compare to ndgrid method (cf. Dan's comment)
tic;
[idxA,idxB] = ndgrid(1:numA,1:numB);
distMat = zeros(numA,numB);
distMat(:) = sum((A(idxA,:) - B(idxB,:)).^2,2);
toc;

结果：

Elapsed time is 1.796201 seconds.
Elapsed time is 5.653246 seconds.
Elapsed time is 3.551636 seconds.
Elapsed time is 22.461185 seconds.

如需更详细的评估 w.r.t.数据点的维度和数量遵循下面的讨论 (@cmets)。事实证明，在不同的设置中应该首选不同的算法。在非时间紧迫的情况下，只需使用pdist2 版本。

进一步发展： 可以考虑用基于相同原理的任何其他度量来替换平方欧几里得：

help = zeros(numA,numB,d);
for idx = 1:d
    help(:,:,idx) = [ones(numA,1), A(:,idx)     ] * ...
                    [B(:,idx)'   ; -ones(1,numB)];
end
distMat = sum(ANYFUNCTION(help),3);

尽管如此，这非常耗时。将较小的d 3 维矩阵help 替换为d 2 维矩阵可能很有用。特别是对于d = 1，它提供了一种通过简单的矩阵乘法计算成对差异的方法：

pairDiffs = [ones(numA,1), A ] * [B'; -ones(1,numB)];

你还有什么想法吗？

【讨论】：

真的很有趣！+1 在另一个故事中：在我的机器上大约从 d>30 开始，bsxfun 将再次获胜，因为内存开销较低。
@knedlsepp 感谢您花时间将所有这些放在一起！好吧，我再次对这两个矢量化版本进行了基准测试，这里提出了基于循环的版本，我没有看到很多差异，至少对于小到体面的 dims 来说不是。
@Divakar：就像在我的机器上一样：如果我们想要平方距离，你的Vec1 版本对于较低维度来说是最快的，直到它被bsxfun 击败。如果我们想要实际的sqrt-距离pdist2 更快，直到它最终也被bsxfun 击败。在做了所有这些比较之后：我想，尽管很高兴知道我们可以从所有这些中榨取最后一点速度，但我还是觉得简单地使用 pdist2 是不费吹灰之力的，如果你有的话安装了统计工具箱，因为它很灵活但仍然非常快速。
@knedlsepp 非常感谢 - 这是一个非常有趣的评估！我想补充一点，log10 中的时间尺度有点误导，因为计算时间的相关性并不存在于对数尺度上（例如，因子 2 对节省时间非常有趣，但在 log10 上看起来什么都没有——规模）。我的事实是：为时间紧迫的实现测试不同的算法是值得的（尤其是对于大量点而言）。例如。对于大量 2d 数据点，使用我的实现很有用。我真的很喜欢我们的算法集合！ :)
这是一个非常有趣的建议和比较。似乎 pdist2 版本缺乏效率主要是由于元素平方，而 Matlab 现在提供了 'squaredeuclidean' 选项来直接获得它。有了这个，所提出的方法和 pdist2 似乎非常接近（也许 pdist2 在某些情况下更快）。该选项可能比发布的答案更新。

【解决方案2】：

对于平方欧几里得距离，也可以使用以下公式

||a-b||^2 = ||a||^2 + ||b||^2 - 2<a,b>

其中<a,b> 是a 和b 之间的点积

nA = sum( A.^2, 2 ); %// norm of A's elements
nB = sum( B.^2, 2 ); %// norm of B's elements
distMat = bsxfun( @plus, nA, nB' ) - 2 * A * B' ;

最近，我一直 told 表示，从 R2016b 开始，这种计算平方欧几里得距离的方法比公认的方法更快。

【讨论】：