如何在没有循环（矢量化）的情况下构建距离矩阵？答案

【问题标题】：How to Build a Distance Matrix without a Loop (Vectorization)?如何在没有循环（矢量化）的情况下构建距离矩阵？
【发布时间】：2013-10-14 12:15:19
【问题描述】：

我有很多点，我想建立距离矩阵，即每个点与所有其他点的距离，但我不想使用 from 循环，因为时间太长...... 构建这个矩阵有更好的方法吗？这是我的循环：对于大小为 10000x3 的 setl，此方法需要花费我很多时间 :(

 for i=1:size(setl,1)
     for j=1:size(setl,1)        
         dist = sqrt((xl(i)-xl(j))^2+(yl(i)-yl(j))^2+...
             (zl(i)-zl(j))^2);
         distanceMatrix(i,j) = dist;
     end
 end

【问题讨论】：

你有统计工具箱吗？如果是，那么：mathworks.com/help/stats/pdist.html
您只需要计算矩阵的一半，因为两点之间的距离是对称的（例如 d(x,y)=d(y,x) ）。否则会计算两次。
不，如果你使用pdist，距离只计算一次。然后可以使用squareform 来构建对称距离矩阵。见my answer here。您也可以输入edit squareform 来查看使用的代码（没有for 循环）。
如果你关心速度，pdist（一个原生 C 函数）和squareform 是唯一的选择，除非你想try compiling mex code 更快一点。

标签： matlab matrix vectorization distance linear-algebra

【解决方案1】：

使用一些线性代数怎么样？两点的距离可以通过它们的位置向量的内积来计算，

D(x, y) = ∥y – x∥ = √ ( x^T x + y^T y – 2 x^T y ),

所有点对的内积都可以通过简单的矩阵运算得到。

x = [xl(:)'; yl(:)'; zl(:)'];
IP = x' * x;
d = sqrt(bsxfun(@plus, diag(IP), diag(IP)') - 2 * IP);

对于 10000 点，我得到以下计时结果：

ahmad 的循环 + shoelzer 的预分配：7.8 秒
Dan 的矢量化索引：5.3 秒
Mohsen 的 bsxfun：1.5 秒
我的解决方案：1.3 秒

【讨论】：

在我的机器上计时，即 Intel Core i7-3770 CPU @ 3.40GHz（四核 + 超线程），在 Debian Linux 下运行的 64 位 Matlab。
我赢了。 ;-) 另外，它真的很优雅。
刚刚注意到您的帖子。 +1 我一直最喜欢这种方法，因为它速度快，而且是对问题的巧妙表述。其实在a different answer中，我对比了几个方法，其中一个和这个是等价的。这是最好的！ ;) 该答案中有一些额外的数学描述，但我得出了相同的结论。另一方面，另一个问题是一团糟。
@chappjc：谢谢！我将根据您的答案提出一个解释公式。可惜这里没有 LaTeX 数学......
对不起，我迟到的评论。非常有用的帖子，但 xl、yl 和 zl 变量是什么？谢谢。

【解决方案2】：

您可以使用bsxfun，这通常是一种更快的解决方案：

s = [xl(:) yl(:) zl(:)];
d = sqrt(sum(bsxfun(@minus, permute(s, [1 3 2]), permute(s, [3 1 2])).^2,3));

【讨论】：

【解决方案3】：

您可以像这样完全矢量化：

n = numel(xl);
[X, Y] = meshgrid(1:n,1:n);
Ix = X(:)
Iy = Y(:)
reshape(sqrt((xl(Ix)-xl(Iy)).^2+(yl(Ix)-yl(Iy)).^2+(zl(Ix)-zl(Iy)).^2), n, n);

如果您查看 Ix 和 Iy（尝试像 3x3 数据集一样尝试），它们使每个矩阵的线性索引的每种组合都成为可能。现在您可以一次完成每个减法！

但是，将 shoelzer 和 Jost 的建议混合起来会给您带来几乎相同的性能提升：

n = 50;

xl = rand(n,1);
yl = rand(n,1);
zl = rand(n,1);

tic
for t = 1:100
    distanceMatrix = zeros(n); %// Preallocation
    for i=1:n
       for j=min(i+1,n):n %// Taking advantge of symmetry
           distanceMatrix(i,j) = sqrt((xl(i)-xl(j))^2+(yl(i)-yl(j))^2+(zl(i)-zl(j))^2);
       end
    end
    d1 = distanceMatrix + distanceMatrix';           %'
end
toc


%// Vectorized solution that creates linear indices using meshgrid
tic
for t = 1:100
    [X, Y] = meshgrid(1:n,1:n);
    Ix = X(:);
    Iy = Y(:);
    d2 = reshape(sqrt((xl(Ix)-xl(Iy)).^2+(yl(Ix)-yl(Iy)).^2+(zl(Ix)-zl(Iy)).^2), n, n);
end
toc

返回：

Elapsed time is 0.023332 seconds.
Elapsed time is 0.024454 seconds.

但如果我将 n 更改为 500 那么我会得到

Elapsed time is 1.227956 seconds.
Elapsed time is 2.030925 seconds.

这只是表明您应该始终在 Matlab 中对解决方案进行基准测试，然后再将循环记为慢！在这种情况下，根据解决方案的规模，循环可能会明显更快。

【讨论】：

非常好的答案，但我质疑你的计时结果。当我在我的机器上运行这段代码时，两种技术都给出了相似的时间，大约 0.02 秒。
@shoelzer 我在线运行它，可能没有 JIT。我会检查真正的 Matlab 并更新
@shoelzer 我在 Matlab 中运行它，你是对的，它们给出了非常相似的结果！已更新。
感谢您的更新！我没想到结果会这样改变。 “始终作为基准”确实是个好建议。
coords = [x,y,z]; d2 = sqrt(sum( bsxfun(@minus,permute(coords,[1,3,2]),permute(coords,[3,1,2])).^2, 3)); 可能更有效。创建两个数组而不是permute创建一个数组更好。

【解决方案4】：

一定要preallocatedistanceMatrix。您的循环将运行得更快，并且可能不需要矢量化。即使你这样做了，也可能不会进一步提高速度。

【讨论】：

我不知道有preallocate的东西（!!!）但是我不能使用这个项目，我对它的使用很感兴趣......我相信我可以使用这个属性给出我的答案...有没有人告诉我：我怎样才能预先分配 distanceMatrix？
@ahmadhosseini distanceMatrix = zeros(size(setl,1),size(setl,1)) 在两个 for 循环之前预分配矩阵。见mathworks.co.uk/support/solutions/en/data/1-18150

【解决方案5】：

MATLAB support Implicit Broadcasting 的最新版本（自 R2016b 起）（另请参阅 bsxfun()）。

因此距离矩阵最快的方法是：

function [ mDistMat ] = CalcDistanceMatrix( mA, mB )

mDistMat = sum(mA .^ 2).' - (2 * mA.' * mB) + sum(mB .^ 2);


end

点在集合的列中的位置。
在你的情况下mA = mB。

看看我的Calculate Distance Matrix Project。

【讨论】：

相关问题 - stackoverflow.com/questions/51181392.