寻找成对欧几里得距离（距离矩阵）的快速算法

Question

我知道Matlab有内置的pdist函数可以计算成对距离。但是，我的矩阵太大了，是60000乘以300，而且Matlab会用完内存。

这个问题是对Matlab欧几里得成对平方距离函数的后续提问。

有没有解决这种计算效率低下的方法？我试过手动编写成对距离计算的代码，通常需要一整天才能运行（有时需要6到7个小时）。

非常感谢任何帮助！

- user2521074

---

这永远不会快速。您需要计算大约20亿个结果，每个结果需要300次乘法和600次加减法。因此总共需要约2e12次操作。 - Oliver Charlesworth

话虽如此，如果代码经过充分优化，应该可以做得比6-7小时更好。 - Oliver Charlesworth

@OliCharlesworth - 唯一了解这个问题的方法是更多地了解正在使用的计算机。有多少内存？ - user85109

更好的投资时间和精力的方式可能是找到一台配备充足内存的64位机器，并在其上运行64位Matlab并运行您的代码（本地或远程）。您有没有拥有Mac和Matlab的朋友？在我的Retina MacBook Pro上，使用Matlab R2012b和OS X 10.8.4仅需要8分45秒（和14,399,760,000字节的内存）来运行d = pdist（rand（6e4,3e2））;。 - horchler

是的，了解限制是很好的，但指出它们并不一定解决问题或回答问题。一个人必须也要实用主义。用户有数据需要分析。也许下一次他们会有10倍的数据量。也可能不会。...所以，存在着限制。问题很大，但有些工具有点无力（如个人电脑）。你会如何解决这个问题？ - horchler

Answer 1

好吧，我忍不住玩了一下。我创建了一个名为pdistc的Matlab mex C文件，用于实现单精度和双精度的成对欧几里德距离计算。在我的机器上，使用Matlab R2012b和R2015a，对于大输入量（例如60,000乘以300），它比pdist（和底层的pdistmex辅助函数）快20-25％。

正如已经指出的那样，这个问题从根本上受到内存的限制，你正在请求很多内存。我的Mex C代码除了输出所需的内存外，使用了最少的内存。将其内存使用与pdist进行比较，它们看起来几乎相同。换句话说，pdist没有使用大量额外的内存。你的内存问题可能在调用pdist之前使用的内存中（你可以使用clear删除任何大数组吗？），或者仅仅因为你试图在微小的硬件上解决一个大型计算问题而产生的。因此，我的pdistc函数可能无法节省你的总内存，但你可以使用我内置的另一个功能。你可以计算你的整个成对距离向量的块。像这样：

m = 6e3;
n = 3e2;
X = rand(m,n);
sz = m*(m-1)/2;

for i = 1:m:sz-m
    D = pdistc(X', i, i+m); % mex C function, X is transposed relative to pdist
    ...                     % Process chunk of pairwise distances
end

这种方法速度较慢（大约慢10倍），而且我的C代码在这部分没有得到很好的优化，但它将允许更少的内存使用——假设您不需要一次性使用整个数组。请注意，您可以通过创建一个循环，在循环中直接传递X的子集，而不是全部传递到 pdist（或pdistc）中，以更高效地完成同样的事情。

如果你有一台64位的Intel Mac，你不需要编译，因为我已经包含了.mexmaci64二进制文件，但如果你使用其他机器，你需要找出如何为你的机器编译代码。我无法帮助你解决这个问题。你可能无法编译它，或者你需要自己编辑代码来解决兼容性问题。还有可能存在漏洞，导致代码崩溃Matlab。此外，请注意，相对于pdist，你可能会得到稍微不同的输出，两者之间的差异在机器epsilon范围内(eps)。 pdist可能会采取一些花哨的方法来避免大输入和其他数字问题的溢出，但请注意，我的代码不会这样做。

此外，我创建了一个简单的pure Matlab implementation。它比mex代码慢得多，但仍然比naïve实现或在pdist中找到的代码快。

所有文件在这里可以找到。ZIP压缩包包含了所有文件。它是BSD许可证。欢迎进行优化（我尝试使用C代码中的BLAS调用和OpenMP，但没有成功-也许一些指针魔法或GPU / OpenCL可以进一步加速）。我希望它对您或他人有所帮助。

- horchler

Answer 2

在我的系统上，以下方法是最快的（甚至比 @horchler 的 C 代码 pdistc 更快）：

function [ mD ] = CalcDistMtx ( mX )    
  vSsqX = sum(mX .^ 2);
  mD = sqrt(bsxfun(@plus, vSsqX.', vSsqX) - (2 * (mX.' * mX)));       
end

我认为你需要一份非常优秀的C代码才能打败它。 更新
自从MATLAB R2016b MATLAB支持隐式广播后，不再需要使用bsxfun()。

因此，代码可以这样编写：

function [ mD ] = CalcDistMtx ( mX )    
  vSsqX = sum(mX .^ 2, 1);
  mD = sqrt(vSsqX.'+ vSsqX - (2 * (mX.' * mX)));       
end

在我的计算距离矩阵项目中给出了一个概括。

附注：
使用MATLAB的pdist进行比较：squareform(pdist(mX.'))等价于CalcDistMtx(mX)。
也就是说，输入应该被转置。

- Royi

---

在当前函数中，mX 是一个 K 行 N 列的矩阵，其中 K 是变量的数量，N 是观察值的数量。通常的符号表示是其转置，所以要小心！ - PhABC

相关问题 - https://dev59.com/iazka4cB1Zd3GeqP_6eh。 - Royi

Answer 3

计算机并不是无限大或无限快的。人们认为它们有很多内存和快速的CPU，所以他们只是创造越来越大的问题，最终想知道为什么他们的问题运行缓慢。事实上，这不是计算效率低下，而只是CPU超载。

正如Oli在评论中指出的那样，即使假设您只计算距离矩阵的上半部分或下半部分，也有大约20亿个值需要计算。（6e4^2/2约为2e9。）这将需要大约16 GB的RAM来存储，假设内存中只创建了一个数组的副本。如果您的代码不够简洁，可能会轻易地将其加倍或加三倍。一旦进入虚拟内存，事情就会变得更加缓慢。

想要一个大问题快速运行是不够的。为了真正帮助您，我们需要知道可用的RAM有多少。这是虚拟内存问题吗？您是否在能处理所有所需RAM的CPU上使用64位MATLAB？

- user85109

---

注意措辞，有些词可能会误导人。实际上，我们一直在使用虚拟内存，而不是“进入”虚拟内存... ;) - Oliver Charlesworth

现代计算机总是使用虚拟内存（每个进程被赋予最大大小的连续虚拟地址空间，由操作系统/架构允许）。你需要注意的是页面置换/抖动。 - Amro

@OliCharlesworth等：你们可能会对我这个问题的新答案感兴趣。比Matlab的pdist稍微快一点。 - horchler