计算 python 中每个点之间距离的最快方法
Fastest way to compute distance beetween each points in python
在我的项目中,我需要计算存储在数组中的每个点之间的欧几里得距离。
条目数组是一个二维 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 是初始条目数组的内存视图,但 flatten()
在函数调用之前由 numpy calcul1()
,R 是一维输出数组的内存视图,用于存储所有结果.
这是我的 Numpy 代码:
def calcul2(M):
return np.sqrt(((M[:,:,np.newaxis] - M[:,np.newaxis,:])**2).sum(axis=0))
此处 M 是初始条目数组,但 transpose()
在函数调用之前由 numpy 将坐标 (x,y,z) 作为行,将点作为列。
此外,这个 numpy 函数非常方便,因为它 returns 的数组组织得很好。它是一个 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(如并行计算)?
谢谢。
不确定你从哪里得到时间,但你可以使用 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 倍。
在我的项目中,我需要计算存储在数组中的每个点之间的欧几里得距离。 条目数组是一个二维 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 是初始条目数组的内存视图,但 flatten()
在函数调用之前由 numpy calcul1()
,R 是一维输出数组的内存视图,用于存储所有结果.
这是我的 Numpy 代码:
def calcul2(M):
return np.sqrt(((M[:,:,np.newaxis] - M[:,np.newaxis,:])**2).sum(axis=0))
此处 M 是初始条目数组,但 transpose()
在函数调用之前由 numpy 将坐标 (x,y,z) 作为行,将点作为列。
此外,这个 numpy 函数非常方便,因为它 returns 的数组组织得很好。它是一个 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(如并行计算)?
谢谢。
不确定你从哪里得到时间,但你可以使用 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 倍。