Cálculo de distancia eficiente entre N puntos y una referencia en numpy/scipy

Question

Acabo de empezar a usar scipy/numpy. Tengo una matriz de 100000 * 3, cada fila es una coordenada y un punto central de 1 * 3. Quiero calcular la distancia para cada fila en el conjunto al centro y almacenarlos en otro conjunto. ¿Cuál es la forma más eficiente de hacerlo?

Answer 1

Me gustaría echar un vistazo 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 para la métrica de distancia predeterminado es equivalente a:

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

aunque cdist es mucho más eficiente para grandes conjuntos (En mi máquina para su problema de tamaño, cdist es más rápido por un factor de ~ 35x).

Answer 2

Puede que necesite especificar de una manera más detallada la función de distancia que le interesa, pero aquí hay una implementación muy simple (y eficiente) de Squared Euclidean Distance basada en inner product (que obviamente puede ser generalizada, de manera directa, a otro tipo de distancia): medidas

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]) 

Dónde P son sus puntos y c es el centro.

Answer 3

También puede usar el desarrollo de la norma (similar a las identidades notables). Esta es probablemente la forma más eficiente de calcular la distancia de una matriz de puntos.

Aquí hay un fragmento de código que utilicé originalmente para una implementación de k-Nearest-Neighbors, en Octave, pero puedes adaptarlo fácilmente a numpy ya que solo usa multiplicaciones de matrices (el equivalente es 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

Utilizaría la implementación de sklearn de la distancia euclidiana. La ventaja es el uso de la expresión más eficiente mediante el uso de la multiplicación de matrices:

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

Un simple script se vería así:

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) 

La ventaja de este enfoque ha sido muy bien descrito en la documentación sklearn : http://scikit-learn.org/stable/modules/generated/sklearn.metrics.pairwise.euclidean_distances.html#sklearn.metrics.pairwise.euclidean_distances

Estoy usando este enfoque para procesar las grandes matrices de datos (10000, 10000) con algunas modificaciones menores como usar la función np.einsum.

Answer 5

Esto podría no responder su pregunta directamente, pero si lo es después de todas las permutaciones de pares de partículas, en algunos casos he encontrado que la siguiente solución es más rápida que la función pdist.

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 

Ver this para una mirada más a fondo sobre este asunto, en mi blog.