向量化依赖嵌套循环答案

【问题标题】：Vectorizing Dependent Nested Loops向量化依赖嵌套循环
【发布时间】：2014-08-28 13:34:00
【问题描述】：

我正在寻找一种有效的方法来使用 MATLAB 计算大小为 (nodeCount x 2) 的坐标矩阵中所有点之间的成对距离。我不想计算两次成对距离（例如，节点 1-2 之间和节点 2-1 之间）。我构建了一个外部“for”循环，该循环通过每个节点递增，内部循环仅评估索引号较高的节点。结果是由节点间隔距离填充的上三角矩阵。我想向量化这些计算，或者至少提高这个操作的效率。任何帮助，将不胜感激。

gap = 10;

for s = 1:(nodeCount); 
  for ss = s+1:(nodeCount);
    if abs(nodeCoord(s,1)-nodeCoord(ss,1)) < gap;
      sep(s,ss) = sqrt((nodeCoord(s,1)-nodeCoord(ss,1))^2+(nodeCoord(s,2)-nodeCoord(ss,2))^2);
    end
  end
end

【问题讨论】：

很想看看这里发布的解决方案在效率方面如何为您服务！
大家好，感谢伟大的 cmets！我已经修改了我的原始代码，以显示一个我用来加快应用程序计算速度的技巧。我最感兴趣的是计算节点对之间的距离，这些节点对之间的距离由变量“gap”定义。我的代码中有一个额外的屏幕，可以过滤掉比“x”方向间隙更远的节点坐标，因此我进一步减少了计算时间和内存需求。我正在处理 10 的数千个节点，而这段代码是瓶颈。
您是否使用零预分配sep？我问是因为如果条件 - if abs(nodeCoord(s,1)-nodeCoord(ss,1)) < gap 不满足，那么 sep(s,ss) 必须是什么？
你是对的。矩阵“sep”预先分配了零。
@MattM：介绍gap 你从根本上改变了问题。高效的实现基于 kD 树。这里讨论的所有解决方案基本上都具有 O(n^2) 运行时间，使用 kd-trees 和恒定数量的 neibours，您可以达到 O(nlogn)

标签： matlab vectorization

【解决方案1】：

循环并不真正依赖于那个角度。我猜你想找到到所有其他坐标的距离试试这个：

xCoord = [1;2;3;4;5];
yCoord = [1;2;3;4;5]:
xSquare = bsxfun(@(x,y) power((x-y),2),xCoord,xCoord.');
ySquare = bsxfun(@(x,y) power((x-y),2),yCoord,yCoord.');
dist = sqrt(xSquare+ySquare);

【讨论】：

感谢您的回复。这段代码似乎比我构建的代码慢，而且内存不足（我在笔记本电脑上运行，大约有 50,000 个节点对）。

【解决方案2】：

xCoord = [1;2;3;4;5];
yCoord = [1;2;3;4;5];

dist = sqrt(pdist2(xCoord,yCoord,'euclidean'));

可以使用函数pdist2

【讨论】：

谢谢！ pdist2 实际上比我的代码慢，大概是因为每对之间的距离计算了两次。

【解决方案3】：

与其尝试使用不需要下三角元素的事实，因为它们在输出中为零，我认为您最好使用基于快速矩阵乘法的技术，如very smart solution 中讨论的那样得到欧几里得距离的完整矩阵。要得到想要的上三角矩阵输出，可以用triu包裹输出。

接下来的代码是它的一个稍微修改的版本，我们正在计算来自nodeCoord的同一对坐标之间的距离。

代码

numA = size(nodeCoord,1);
helpA = ones(numA,6);
helpB = ones(numA,6);
for idx = 1:2
    sqA_idx = nodeCoord(:,idx).^2;
    helpA(:,3*idx-1:3*idx) = [-2*nodeCoord(:,idx), sqA_idx ];
    helpB(:,3*idx-2:3*idx-1) = [sqA_idx , nodeCoord(:,idx)];
end
sep = triu(sqrt(helpA(:,1:3) * helpB(:,1:3)')<gap).* sqrt(helpA * helpB');

【讨论】：

【解决方案4】：

pdist(nodeCoord) 以一种快速的方式执行此操作，但以向量的形式返回数据。将其映射回矩阵的成本与计算距离大致相同：

sep3=zeros(nodeCount,nodeCount);
sep3(tril(true(nodeCount),-1))=pdist(nodeCoord);
sep3=sep3+sep3.';

如果您对下三角矩阵感到满意，可以省略最后一行。

【讨论】：