Python中最快的成对距离度量方法

Question

Python中最快的成对距离度量方法

21

我有一个一维数字数组，想要计算所有成对欧几里得距离。我有一种方法（感谢SO）可以使用广播来完成这个任务，但它效率低下，因为它会将每个距离计算两次。而且它不具有良好的扩展性。

这是一个示例，对于一个包含1000个数字的数组它能给我我想要的结果。

import numpy as np
import random
r = np.array([random.randrange(1, 1000) for _ in range(0, 1000)])
dists = np.abs(r - r[:, None])

在scipy/numpy/scikit-learn中，哪种实现方式是最快的？我需要处理的数据规模可能会超过10k个值。

注意：矩阵是对称的，因此我猜测至少可以通过处理对称性来加快速度，但我不知道具体方法。

- roblanf

6

有一个相应的功能：scipy.spatial.distance.pdist。我不确定这是否是最快的选项，因为它需要检查多维数据、非欧几里德范数和其他东西，但它是内置的。 - user2357112

你需要这个有多快？它永远不会比O(n^2)更好，因为你必须填充n^2个输出条目。你现有的解决方案是O(n^2)，似乎没有太多优化的空间。 - user2357112

赞同 @user2357112 和 @askewchan 的观点，但要确保你的 numpy 使用了 BLAS 或 MKL 编译，直接从 sourceforge 下载的可能没有使用这些编译选项。 - CT Zhu

3

我不认为它有...... 如果你跟随源代码，在最后调用的是这个函数。不仅没有花哨的优化，而且对于一维向量，它是通过求平方并取平方根来计算绝对值的。对于他特定的用例，可能比OP的代码更糟糕。 - Jaime

1

如果我没记错的话，scipy 总是使用 BLAS 编译的，这不像 numpy 那样是可选的。 - askewchan

显示剩余8条评论

3个回答

7

这里是一个Cython实现，在我的电脑上，该示例的速度提高了3倍以上。对于更大的数组，应该重新评估这个计时，因为BLAS例程可能比这个相当天真的代码更好地扩展。

我知道你要求一些在scipy/numpy/scikit-learn内部的东西，但也许这会为你打开新的可能性：

文件my_cython.pyx：

import numpy as np
cimport numpy as np
import cython

cdef extern from "math.h":
    double abs(double t)

@cython.wraparound(False)
@cython.boundscheck(False)
def pairwise_distance(np.ndarray[np.double_t, ndim=1] r):
    cdef int i, j, c, size
    cdef np.ndarray[np.double_t, ndim=1] ans
    size = sum(range(1, r.shape[0]+1))
    ans = np.empty(size, dtype=r.dtype)
    c = -1
    for i in range(r.shape[0]):
        for j in range(i, r.shape[0]):
            c += 1
            ans[c] = abs(r[i] - r[j])
    return ans

答案是一个包含所有非重复评估的一维数组。

要导入到Python中:

import numpy as np
import random

import pyximport; pyximport.install()
from my_cython import pairwise_distance

r = np.array([random.randrange(1, 1000) for _ in range(0, 1000)], dtype=float)

def solOP(r):
    return np.abs(r - r[:, None])

使用IPython进行时间控制：

In [2]: timeit solOP(r)
100 loops, best of 3: 7.38 ms per loop

In [3]: timeit pairwise_distance(r)
1000 loops, best of 3: 1.77 ms per loop

- Saullo G. P. Castro

1

你肯定是想用 fabs -- abs 是 int 变量。 - Fred Foo

7

相较于np.abs(r - r[:, None])，使用一半的内存，但速度慢6倍：

triu = np.triu_indices(r.shape[0],1)
dists2 = abs(r[triu[1]]-r[triu[0]])

- cyborg

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- roblanf · Accepted Answer

其他答案都没有完全回答这个问题，其中一个是用Cython编写的，另一个则比较慢。但它们都提供了非常有用的提示。按照这些提示进行后续研究表明 scipy.spatial.distance.pdist 是前进的方向。

以下是一些代码：

import numpy as np
import random
import sklearn.metrics.pairwise
import scipy.spatial.distance

r = np.array([random.randrange(1, 1000) for _ in range(0, 1000)])
c = r[:, None]

def option1(r):
    dists = np.abs(r - r[:, None])

def option2(r):
    dists = scipy.spatial.distance.pdist(r, 'cityblock')

def option3(r):
    dists = sklearn.metrics.pairwise.manhattan_distances(r)

使用IPython进行时间管理：

In [36]: timeit option1(r)
100 loops, best of 3: 5.31 ms per loop

In [37]: timeit option2(c)
1000 loops, best of 3: 1.84 ms per loop

In [38]: timeit option3(c)
100 loops, best of 3: 11.5 ms per loop

我没有尝试Cython实现（对于这个项目我无法使用它），但是与另一个已经尝试过的答案进行比较，看起来scipy.spatial.distance.pdist大约比Cython实现慢三分之一（通过在np.abs解决方案上进行基准测试考虑了不同机器）。