Calcolo della distanza efficiente tra N punti e un riferimento in numpy/scipy

Question

Ho appena iniziato a usare scipy/numpy. Ho una matrice 100000 * 3, ogni riga è una coordinata e un punto centrale 1 * 3. Voglio calcolare la distanza per ogni riga dell'array al centro e memorizzarla in un altro array. Qual è il modo più efficiente per farlo?

Answer 1

vorrei dare un'occhiata a scipy.spatial.distance.cdist:

http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html

import numpy as np 
import scipy 

a = np.random.normal(size=(10,3)) 
b = np.random.normal(size=(1,3)) 

dist = scipy.spatial.distance.cdist(a,b) # pick the appropriate distance metric 

dist per la metrica lontana predefinita è equivalente a:

np.sqrt(np.sum((a-b)**2,axis=1)) 

anche se cdist è molto più efficiente per grandi array (sulla mia macchina per il tuo problema di dimensioni, cdist è più veloce di un fattore di ~ 35x).

Answer 2

Potrebbe essere necessario specificare in modo più dettagliato la funzione di distanza a cui sei interessato, ma qui è un'implementazione molto semplice (ed efficiente) di Squared Euclidean Distance basata su inner product (che ovviamente può essere generalizzata, modo semplice, per altri tipi di distanza misure):

In []: P, c= randn(5, 3), randn(1, 3) 
In []: dot(((P- c)** 2), ones(3)) 
Out[]: array([ 8.80512, 4.61693, 2.6002, 3.3293, 12.41800]) 

Dove P sono i vostri punti e c è il centro.

Answer 3

È inoltre possibile utilizzare lo sviluppo della norma (simile a identità notevoli). Questo è probabilmente il modo più efficace per calcolare la distanza di una matrice di punti.

Ecco uno snippet di codice che ho utilizzato originariamente per un'implementazione k-Nearest-Neighbors, in Octave, ma è possibile adattarlo facilmente a numpy poiché utilizza solo moltiplicazioni di matrice (l'equivalente è numpy.dot()):

% Computing the euclidian distance between each known point (Xapp) and unknown points (Xtest) 
% Note: we use the development of the norm just like a remarkable identity: 
% ||x1 - x2||^2 = ||x1||^2 + ||x2||^2 - 2*<x1,x2> 
[napp, d] = size(Xapp); 
[ntest, d] = size(Xtest); 

A = sum(Xapp.^2, 2); 
A = repmat(A, 1, ntest); 

B = sum(Xtest.^2, 2); 
B = repmat(B', napp, 1); 

C = Xapp*Xtest'; 

dist = A+B-2.*C; 

Answer 4

Vorrei utilizzare l'implementazione sklearn della distanza euclidea. Il vantaggio è l'uso dell'espressione più efficiente utilizzando moltiplicazione di matrici:

dist(x, y) = sqrt(dot(x, x) - 2 * dot(x, y) + dot(y, y) 

Un semplice script sarebbe simile a questa:

import numpy as np 

x = np.random.rand(1000, 3) 
y = np.random.rand(1000, 3) 

dist = np.sqrt(np.dot(x, x)) - (dot(x, y) + dot(x, y)) + dot(y, y) 

Il vantaggio di questo approccio è stato ben descritto nella documentazione sklearn : http://scikit-learn.org/stable/modules/generated/sklearn.metrics.pairwise.euclidean_distances.html#sklearn.metrics.pairwise.euclidean_distances

Sto utilizzando questo approccio per il crunch di grandi datamatrices (10000, 10000) con alcune modifiche minori come l'utilizzo della funzione np.einsum.

Answer 5

Questo potrebbe non rispondere alla domanda direttamente, ma se si è dopo tutte le permutazioni di coppie di particelle, ho trovato la seguente soluzione per essere più veloce della funzione pdist in alcuni casi.

import numpy as np 

L = 100  # simulation box dimension 
N = 100  # Number of particles 
dim = 2   # Dimensions 

# Generate random positions of particles 
r = (np.random.random(size=(N,dim))-0.5)*L 

# uti is a list of two (1-D) numpy arrays 
# containing the indices of the upper triangular matrix 
uti = np.triu_indices(100,k=1)  # k=1 eliminates diagonal indices 

# uti[0] is i, and uti[1] is j from the previous example 
dr = r[uti[0]] - r[uti[1]]   # computes differences between particle positions 
D = np.sqrt(np.sum(dr*dr, axis=1)) # computes distances; D is a 4950 x 1 np array 

Vedi this per uno sguardo più approfondito su questo tema, nel mio blog.