Pythonで各点間の距離を計算する最速の方法

Question

私のプロジェクトでは、配列に格納されている各点のユークリッド距離を計算する必要があります。 エントリ配列は、座標（x、y、z）である3列の2D numpy配列であり、各行は新しい点を定義します。

私は通常、テストケースで5000〜6000ポイントで作業しています。

私の最初のアルゴリズムはCythonと私の2番目のnumpyを使っています。私の手の込んだアルゴリズムは、cythonより高速です。

編集：

numpyの1.76秒/ 4.36秒

cythonここに私の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()によって6000ポイントで関数calcul1()の呼び出しの前にnumpy、Rはすべての結果を格納する1D出力配列のメモリビューです。ここ

は私numpyのコードは次のとおり

def calcul2(M): 

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

ここでMは、列としての行と点として座標（x、y、z）を有するように、関数呼び出しの前にnumpyのにより初期エントリアレイが、transpose()あります。

また、このnumpy関数は、それが返す配列がうまく整理されているので、非常に便利です。これはn個のn個の配列でn個の点を持ち、各点には行と列があります。

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) 

私はここで間違って何かをやっている？だから、例えばABは、A列の交差点インデックスに格納された距離と列B

は、ここで私は（cython機能）それらを呼び出す方法ですか私のプロジェクトでは、両方が十分に速いわけではありません。

1：numpyの速度を上げるために私のcythonコードを改善する方法はありますか？

2：私のnumpyコードをさらに高速に計算する方法はありますか？

3：または他の解決策ですが、Python/Cython（並列コンピューティングのような）でなければなりませんか？

ありがとうございます。

Answer 1

あなたのタイミングを取得しているが、あなたはscipy.spatial.distance使用できる場所わからない：

：あなたの出力が対称であることを理解することが重要なの

%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 

、そのも便利：

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 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倍速いです。