问求成对欧氏距离(距离矩阵)的快速算法

EN

好吧，我忍不住要玩了。我创建了一个名为pdistc的Matlab，它实现了单精度和双精度的成对欧几里德距离。在我使用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代码的这一部分没有得到很好的优化，但它将允许更少的内存使用-假设你不需要一次使用整个数组。请注意，使用pdist (或pdistc)可以更有效地完成相同的工作，方法是创建一个循环，直接传入X的子集，而不是全部。

如果您使用的是64位Intel Mac，则不需要编译，因为我已经包含了.mexmaci64二进制文件，否则您将需要弄清楚如何为您的机器编译代码。我帮不了你。您可能无法对其进行编译，或者可能存在兼容性问题，您需要通过自己编辑代码来解决这些问题。也有可能存在错误，代码会导致Matlab崩溃。此外，请注意，由于机器ε(eps)的范围不同，您可能会获得与pdist相关的略微不同的输出。pdist可能会做一些花哨的事情，也可能不会做一些花哨的事情来避免大输入溢出和其他数字问题，但要知道我的代码不会。

另外，我创建了一个简单的pure Matlab implementation。它比mex代码慢得多，但仍然比天真的实现或pdist中的代码快。

所有的文件都是can be found here。ZIP存档包括所有文件。这是BSD授权的。请随意优化(我在C代码中尝试了BLAS调用和OpenMP，但都无济于事--也许一些指针魔法或GPU/OpenCL可以进一步提高速度)。我希望它能对你或其他人有所帮助。

在我的系统上，下面是最快的(甚至比@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 supports implicit broadcasting没有使用bsxfun()。

因此，代码可以写成：

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

在我的Calculate Distance Matrix project中给出了一个概括。

另请注意：

使用MATLAB的pdist进行比较：squareform(pdist(mX.'))等同于CalcDistMtx(mX)。

也就是说，输入应该是转置的。

计算机不是无限大，也不是无限快。人们认为他们有大量的内存，一个快速的CPU，所以他们只会制造越来越大的问题，然后最终想知道为什么他们的问题运行缓慢。事实是，这并不是计算效率低下。它只是一个过载的CPU。

正如Oli在评论中指出的那样，即使假设您只计算距离矩阵的上半部分或下半部分，也有类似2e9的值需要计算。(6e4^2/2大约是2e9。)这将需要大约16 in的RAM来存储，假设只在内存中创建了数组的一个副本。如果你的代码是草率的，你可以很容易地将其增加一倍或三倍。一旦进入虚拟内存，事情就会变得慢得多。

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