Il modo più veloce per calcolare la distanza tra ogni punto in python

Question

Nel mio progetto ho bisogno di calcolare distanza euclidea tra ogni punto memorizzato in una matrice. La matrice di immissione è una matrice numpia 2D con 3 colonne che sono le coordinate (x, y, z) e ogni riga definisce un nuovo punto.

Sono solito lavorare con 5000 - 6000 punti nei miei casi di test.

Il mio primo algoritmo utilizza Cython e il mio secondo numpy. Trovo che il mio algoritmo numpy sia più veloce di cython.

edit: con 6000 punti:

NumPy 1.76 s/Cython 4.36 s

Ecco il mio codice 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 è una vista della memoria della matrice di entrata iniziale, ma flatten() da numpy prima della chiamata della funzione calcul1(), R è una vista di memoria di una matrice di uscita 1D per memorizzare tutti i risultati.

Ecco mio codice Numpy:

def calcul2(M): 

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

Qui M è la matrice prima messa ma transpose() dal numpy prima della chiamata di funzione ha coordinate (x, y, z) come righe e punti come colonne.

Inoltre questa funzione numpy è abbastanza comoda perché la matrice restituita è ben organizzata. È una matrice di n con n il numero di punti e ogni punto ha una riga e una colonna. Così, per esempio la distanza AB è memorizzato l'indice intersezione di riga A e B. colonna

ecco come li (funzione Cython) chiamano:

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) 

sto facendo qualcosa di sbagliato qui? Per il mio progetto entrambi non sono abbastanza veloci.

1: C'è un modo per migliorare il mio codice cython al fine di battere la velocità di numpy?

2: C'è un modo per migliorare il mio codice di NumPy per il calcolo ancora più veloce?

3: O qualsiasi altra soluzione, ma deve essere un python/cython (come il calcolo parallelo)?

Grazie.

Answer 1

Non sai da dove hai trovato i tuoi tempi, ma è possibile utilizzare 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 

Timings:

%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 

È importante sottolineare che il suo utile anche per rendersi conto che l'output è simmetrica:

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

Un'alternativa è calcolare solo il triangolare superiore di questo array:

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

Modifica: Non sei sicuro di quale indice vuoi comprimere, sembra che tu stia facendo in entrambi i modi? Zippare l'altro indice per il confronto:

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

Ancora circa 10 volte più veloce dell'attuale implementazione NumPy.