Numpy Broadcast para realizar la distancia euclídea vectorizada

Tengo matrices que son 2 x 4 y 3 x 4. Quiero encontrar la distancia euclidiana a través de las filas y obtener una matriz de 2 x 3 al final. Aquí está el código con one for loop que calcula la distancia euclidiana para cada vector de fila en contra de todos los vectores de fila b. ¿Cómo hago lo mismo sin usar para bucles?

import numpy as np a = np.array([[1,1,1,1],[2,2,2,2]]) b = np.array([[1,2,3,4],[1,1,1,1],[1,2,1,9]]) dists = np.zeros((2, 3)) for i in range(2): dists[i] = np.sqrt(np.sum(np.square(a[i] - b), axis=1)) 

Simplemente use np.newaxis en el lugar correcto:

  np.sqrt((np.square(a[:,np.newaxis]-b).sum(axis=2))) 

Aquí están las variables de entrada originales:

 A = np.array([[1,1,1,1],[2,2,2,2]]) B = np.array([[1,2,3,4],[1,1,1,1],[1,2,1,9]]) A # array([[1, 1, 1, 1], # [2, 2, 2, 2]]) B # array([[1, 2, 3, 4], # [1, 1, 1, 1], # [1, 2, 1, 9]]) 

A es una matriz de 2×4. B es una matriz de 3×4.

Queremos calcular la operación de la matriz de distancia euclidiana en una operación completamente vectorizada, donde dist[i,j] contiene la distancia entre la instancia ith en A y la instancia jth en B. Así que dist es 2×3 en este ejemplo.

La distancia

introduzca la descripción de la imagen aquí

ostensiblemente podría escribirse con numpy como

 dist = np.sqrt(np.sum(np.square(AB))) # DOES NOT WORK # Traceback (most recent call last): # File "", line 1, in  # ValueError: operands could not be broadcast together with shapes (2,4) (3,4) 

Sin embargo, como se muestra arriba, el problema es que la operación de sustracción AB involucra tamaños de matriz incompatibles, específicamente los 2 y 3 en la primera dimensión.

 A has dimensions 2 x 4 B has dimensions 3 x 4 

Para realizar la resta de elementos, tenemos que rellenar A o B para satisfacer las reglas de transmisión de números. Elegiré rellenar A con una dimensión adicional para que se convierta en 2 x 1 x 4, lo que permite que las dimensiones de los arreglos se alineen para la transmisión. Para obtener más información sobre la difusión numpy, consulte el tutorial en el manual de scipy y el ejemplo final en este tutorial .

Puede realizar el relleno con el valor np.newaxis o con el comando np.reshape . Os muestro a continuación:

 # First approach is to add the extra dimension to A with np.newaxis A[:,np.newaxis,:] has dimensions 2 x 1 x 4 B has dimensions 3 x 4 # Second approach is to reshape A with np.reshape np.reshape(A, (2,1,4)) has dimensions 2 x 1 x 4 B has dimensions 3 x 4 

Como puede ver, el uso de cualquiera de los dos enfoques permitirá que las dimensiones se alineen. np.newaxis el primer enfoque con np.newaxis . Así que ahora, esto funcionará para crear AB, que es una matriz de 2x3x4:

 diff = A[:,np.newaxis,:] - B # Alternative approach: # diff = np.reshape(A, (2,1,4)) - B diff.shape # (2, 3, 4) 

Ahora podemos poner esa expresión de diferencia en la instrucción de la ecuación dist para obtener el resultado final:

 dist = np.sqrt(np.sum(np.square(A[:,np.newaxis,:] - B), axis=2)) dist # array([[ 3.74165739, 0. , 8.06225775], # [ 2.44948974, 2. , 7.14142843]]) 

Tenga en cuenta que la sum está sobre el axis=2 , lo que significa tomar la sum sobre el tercer eje de la matriz 2x3x4 (donde la identificación del eje comienza con 0).

Si sus arreglos son pequeños, entonces el comando anterior funcionará bien. Sin embargo, si tiene matrices grandes, es posible que tenga problemas de memoria. Tenga en cuenta que en el ejemplo anterior, numpy creó internamente una matriz de 2x3x4 para realizar la transmisión. Si generalizamos A para tener las dimensiones axz y B para tener las dimensiones bxz , entonces numpy creará internamente una matriz axbxz para la difusión.

Podemos evitar crear esta matriz intermedia haciendo alguna manipulación matemática. Debido a que está calculando la distancia euclidiana como una sum de diferencias al cuadrado, podemos aprovechar el hecho matemático de que la sum de las diferencias al cuadrado puede reescribirse.

introduzca la descripción de la imagen aquí

Tenga en cuenta que el término medio implica la sum sobre la multiplicación por elementos . Esta sum sobre multiplicaciones es mejor conocida como un producto de puntos. Debido a que A y B son cada una matriz, entonces esta operación es en realidad una multiplicación de matrices. Así podemos reescribir lo anterior como:

introduzca la descripción de la imagen aquí

Luego podemos escribir el siguiente código numpy:

 threeSums = np.sum(np.square(A)[:,np.newaxis,:], axis=2) - 2 * A.dot(BT) + np.sum(np.square(B), axis=1) dist = np.sqrt(threeSums) dist # array([[ 3.74165739, 0. , 8.06225775], # [ 2.44948974, 2. , 7.14142843]]) 

Tenga en cuenta que la respuesta anterior es exactamente la misma que la implementación anterior. Nuevamente, la ventaja aquí es que no necesitamos crear la matriz intermedia de 2x3x4 para la transmisión.

Para completar, verifiquemos que las dimensiones de cada sumndo en threeSums permitieron la transmisión.

 np.sum(np.square(A)[:,np.newaxis,:], axis=2) has dimensions 2 x 1 2 * A.dot(BT) has dimensions 2 x 3 np.sum(np.square(B), axis=1) has dimensions 1 x 3 

Entonces, como se esperaba, la matriz dist final tiene dimensiones 2×3.

Este uso del producto punto en lugar de la sum de la multiplicación de elementos también se discute en este tutorial .

Recientemente tuve el mismo problema al trabajar con el aprendizaje profundo (stanford cs231n, Asignación1), pero cuando lo usé

  np.sqrt((np.square(a[:,np.newaxis]-b).sum(axis=2))) 

Hubo un error

 MemoryError 

Eso significa que me quedé sin memoria (de hecho, eso produjo una matriz de 500 * 5000 * 1024 en el medio. ¡Es tan enorme!)

Para evitar ese error, podemos usar una fórmula para simplificar:

código:

 import numpy as np aSumSquare = np.sum(np.square(a),axis=1); bSumSquare = np.sum(np.square(b),axis=1); mul = np.dot(a,bT); dists = np.sqrt(aSumSquare[:,np.newaxis]+bSumSquare-2*mul) 

Esta funcionalidad ya está incluida en el módulo espacial de scipy y recomiendo su uso ya que será vectorizado y altamente optimizado bajo el capó. Pero, como lo demuestra la otra respuesta, hay formas de hacerlo usted mismo.

 import numpy as np a = np.array([[1,1,1,1],[2,2,2,2]]) b = np.array([[1,2,3,4],[1,1,1,1],[1,2,1,9]]) np.sqrt((np.square(a[:,np.newaxis]-b).sum(axis=2))) # array([[ 3.74165739, 0. , 8.06225775], # [ 2.44948974, 2. , 7.14142843]]) from scipy.spatial.distance import cdist cdist(a,b) # array([[ 3.74165739, 0. , 8.06225775], # [ 2.44948974, 2. , 7.14142843]]) 

El uso de numpy.linalg.norm también funciona bien con la difusión. La especificación de un valor entero para el axis utilizará una norma vectorial, que por defecto es la norma euclidiana.

 import numpy as np a = np.array([[1,1,1,1],[2,2,2,2]]) b = np.array([[1,2,3,4],[1,1,1,1],[1,2,1,9]]) np.linalg.norm(a[:, np.newaxis] - b, axis = 2) # array([[ 3.74165739, 0. , 8.06225775], # [ 2.44948974, 2. , 7.14142843]])