4

我有很多点，我想建立距离矩阵，即每个点与所有其他点的距离，但我不想使用循环，因为太花时间了......是建立这个矩阵的更好方法吗？这是我的循环：对于大小为 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
4

5 回答 5

19

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

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 秒
于 2013-10-14T13:43:23.660 回答
3

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

s = [xl(:) yl(:) zl(:)];
d = sqrt(sum(bsxfun(@minus, permute(s, [1 3 2]), permute(s, [3 1 2])).^2,3));
于 2013-10-14T13:50:40.580 回答
2

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

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);

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

于 2013-10-14T12:56:42.383 回答
1

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

于 2013-10-14T12:30:41.337 回答
0

MATLAB的最新版本（自 R2016b 起）支持隐式广播（另请参阅 参考资料bsxfun()）。

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

function [ mDistMat ] = CalcDistanceMatrix( mA, mB )

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


end

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

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

于 2018-04-20T16:56:29.183 回答