在Python中计算每个点之间距离的最快方法

Question

在Python中计算每个点之间距离的最快方法

pythonnumpyoptimizationcython

6

在我的项目中，我需要计算存储在数组中的每个点之间的欧几里得距离。输入数组是一个2D numpy数组，有3列，分别是坐标（x，y，z），每行定义一个新的点。

我通常在测试用例中使用5000-6000个点。

我的第一个算法使用Cython，第二个使用numpy。我发现numpy算法比cython更快。

编辑：使用6000个点：

numpy 1.76秒 / cython 4.36秒

这是我的cython代码：

cimport cython
from libc.math cimport sqrt
@cython.boundscheck(False)
@cython.wraparound(False)
cdef void calcul1(double[::1] M,double[::1] R):

  cdef int i=0
  cdef int max = M.shape[0]
  cdef int x,y
  cdef int start = 1

  for x in range(0,max,3):
     for y in range(start,max,3):

        R[i]= sqrt((M[y] - M[x])**2 + (M[y+1] - M[x+1])**2 + (M[y+2] - M[x+2])**2)
        i+=1  

     start += 1

M是初始入口数组的内存视图，但在调用函数calcul1()之前被numpy中的flatten()扁平化处理，R是一个1D输出数组的内存视图，用于存储所有结果。

这是我的Numpy代码：

def calcul2(M):

     return np.sqrt(((M[:,:,np.newaxis] - M[:,np.newaxis,:])**2).sum(axis=0))

这里的M是初始输入数组，但在函数调用之前通过numpy中的transpose()进行转置，使坐标(x,y,z)成为行，点成为列。

此外，这个numpy函数非常方便，因为它返回的数组组织得很好。它是一个n×n的数组，其中n是点的数量，每个点都有一行和一列。因此，例如距离AB存储在行A和列B的交叉索引处。

以下是我如何调用它们（cython函数）：

cpdef test():

  cdef double[::1] Mf 
  cdef double[::1] out = np.empty(17998000,dtype=np.float64) # (6000² - 6000) / 2

  M = np.arange(6000*3,dtype=np.float64).reshape(6000,3) # Example array with 6000 points
  Mf = M.flatten() #because my cython algorithm need a 1D array
  Mt = M.transpose() # because my numpy algorithm need coordinates as rows

  calcul2(Mt)

  calcul1(Mf,out)

我在这里做错了什么吗？对于我的项目，两者都不够快。

1：有没有办法改进我的Cython代码，以打败numpy的速度？

2：有没有办法改进我的numpy代码，使其计算速度更快？

3：或者其他的解决方案，但必须是Python / Cython（如并行计算）？

谢谢。

- UserAt

1

如果您不需要距离，只关心差异/排名，那么可以摆脱sqrt，这应该是计算中最慢的部分。也许您还可以使用更快的sqrt，它不太精确，或者使用其他度量标准（例如出租车距离）。 - sascha

2

有了5000到6000个点，您的矩阵将有大约3000万个条目。计算30m次平方根肯定会很慢。您真的需要完整的密集矩阵吗？在计算矩阵后，您要做什么？ - Sven Marnach

NumPy 比 Cython 快多少？ - sebacastroh

根据Sascha的建议，在计算L2范数时，使用x*x可能比使用x**2更快。 - AndyG

3

你看过 https://dev59.com/RILba4cB1Zd3GeqPc1Ku 吗？这不是同一个问题吗？ - sebacastroh

显示剩余2条评论

1个回答

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

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

- Daniel · Accepted Answer

不确定您从何处获取时间信息，但您可以使用 scipy.spatial.distance：

M = np.arange(6000*3, dtype=np.float64).reshape(6000,3)
np_result = calcul2(M)
sp_result = sd.cdist(M.T, M.T) #Scipy usage
np.allclose(np_result, sp_result)
>>> True

时序:

%timeit calcul2(M)
1000 loops, best of 3: 313 µs per loop

%timeit sd.cdist(M.T, M.T)
10000 loops, best of 3: 86.4 µs per loop

重要的是要意识到，输出是对称的：

np.allclose(sp_result, sp_result.T)
>>> True

另一种方法是仅计算该数组的上三角：

%timeit sd.pdist(M.T)
10000 loops, best of 3: 39.1 µs per loop

编辑：不确定您想要压缩哪个索引，看起来您可能会以两种方式进行？为了比较，压缩另一个索引：

%timeit sd.pdist(M)
10 loops, best of 3: 135 ms per loop

仍然比您当前的NumPy实现快10倍左右。