向量之间的向量化距离计算答案

【问题标题】：Vectorizing distance calculation between vectors向量之间的向量化距离计算
【发布时间】：2014-06-21 20:45:07
【问题描述】：

我有一个 3 X 1000（以及后来的 3 X 10 000）矩阵cord，其中包含我的像素的三维坐标。

我的目的是计算所有像素之间的距离，我用一个 for 循环（见下文）来计算，但我很快就必须为巨大的矩阵计算这个，我想知道我是否可以对代码进行矢量化让它更快...？

dist = zeros(size(cord,2),size(cord,2)); 
for i = 1:size(cord,2)
    for j = 1:size(cord,2)
        dist(i,j) = norm(cord(:,i)-cord(:,j));
        dist(j,i) = dist(i,j);
    end 
end

【问题讨论】：

标签： image matlab for-loop vectorization

【解决方案1】：

pdist 正是这样做的。需要squareform 才能得到方阵对称矩阵形式的结果：

dist = squareform(pdist(cord.'));

【讨论】：

【解决方案2】：

方法 1（bsxfun 的矢量化方法）-

squeeze(sqrt(sum(bsxfun(@minus,cord,permute(cord,[1 3 2])).^2)))

不确定这是否会更快。

方法 2 -

受此 very smart approach 的启发，以及海报的所有功劳。此处发布的代码只是针对您的情况稍微定制的，希望在运行时稍微好一些。在这里——

A = cord'; %//'
numA = size(cord,2);
helpA = ones(numA,9);
helpB = ones(numA,9);
for idx = 1:3
    sqA_idx = A(:,idx).^2;
    helpA(:,3*idx-1:3*idx) = [-2*A(:,idx), sqA_idx ];
    helpB(:,3*idx-2:3*idx-1) = [sqA_idx , A(:,idx)];
end
dist1 = sqrt(helpA * helpB'); %// desired output

【讨论】：

【解决方案3】：

从您的代码中，您已经认识到 dist 矩阵是对称的

dist(i,j) = norm(cord(:,i)-cord(:,j));
dist(j,i) = dist(i,j);

您可以更改内部循环以解决此问题，并将所需计算次数减少大约一半

for j = i:size(cord,2)

此外，我们可以在每次迭代中避免dist(j,i) = dist(i,j);，并在最后通过提取dist 的上三角部分并将其转置添加到dist 矩阵中来解决对称性

dist = zeros(size(cord,2),size(cord,2)); 
for i = 1:size(cord,2)
    for j = i:size(cord,2)
        dist(i,j) = norm(cord(:,i)-cord(:,j));
    end 
end
dist = dist + triu(dist)';

上面的加法很好，因为主对角线全为零。

虽然它仍然表现不佳，所以我们应该利用矢量化。我们可以针对内部循环执行以下操作

dist = zeros(size(cord,2),size(cord,2)); 
for i = 1:size(cord,2)
    dist(i,i+1:end) = sum((repmat(cord(:,i),1,size(cord,2)-i)-cord(:,i+1:end)).^2);
end
dist = dist + triu(dist)';
dist = sqrt(dist);

对于cord 中的每个元素，我们需要计算它与跟随它的所有其他元素的距离。我们用repmat 复制元素，这样我们就可以从后面的每个元素中减去它，而无需循环。将差值平方和求和并分配给dist 矩阵。我们处理对称性，然后取矩阵的平方根，完成norm运算。

使用tic 和toc，使用随机绳索 (cord = rand(3,num);) 的原始距离计算耗时约 93 秒。这个版本花了~2.8。

【讨论】：