如何构建没有循环的距离矩阵(矢量化)

Question

我有很多点，我想建立距离矩阵，即每个点与所有其他点的距离，但我想不要使用 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

Answer 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 点，我得到以下计时结果：

• 艾哈迈德循环 + 舒尔泽的预分配：7.8 秒
• Dan 的矢量化索引：5.3 秒
• 莫森的bsxfun：1.5秒
• 我的解决方案：1.3秒

Answer 2

您可以使用

bsxfun这通常是一个更快的解决方案：

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

Answer 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 中对解决方案进行基准测试！在这种情况下，根据解决方案的规模，循环可能会明显更快。

Answer 4

请务必预先分配distanceMatrix 。您的循环将运行得更快，并且可能不需要矢量化。即使你这样做，也可能不会有任何进一步的速度增加。

Answer 5

MATLAB 的最新版本（自 R2016b 起）支持隐式广播（另见 bsxfun() 上注明）。

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

function [ mDistMat ] = CalcDistanceMatrix( mA, mB )
mDistMat = sum(mA .^ 2).' - (2 * mA.' * mB) + sum(mB .^ 2);

end

点沿集合的列的位置。
在您的情况下mA = mB.

看看我的计算距离矩阵项目。