¿Por qué Looping Beat Indexing aquí?

Hace unos años, alguien publicó en Active State Recipes con fines de comparación, tres funciones python / NumPy; cada uno de estos aceptó los mismos argumentos y devolvió el mismo resultado, una matriz de distancia .

Dos de estos fueron tomados de fonts publicadas; son ambos, o me parece que son, un código numérico idiomático. Los cálculos repetitivos necesarios para crear una matriz de distancia se basan en la elegante syntax de índice de numpy. Aquí hay uno de ellos:

from numpy.matlib import repmat, repeat def calcDistanceMatrixFastEuclidean(points): numPoints = len(points) distMat = sqrt(sum((repmat(points, numPoints, 1) - repeat(points, numPoints, axis=0))**2, axis=1)) return distMat.reshape((numPoints,numPoints)) 

El tercero creó la matriz de distancia utilizando un solo bucle (que, obviamente, es un montón de bucles, dado que una matriz de distancia de solo 1,000 puntos 2D, tiene un millón de entradas). A primera vista, esta función me parecía el código que solía escribir cuando estaba aprendiendo NumPy y escribiría código NumPy al escribir primero el código Python y luego traducirlo, línea por línea.

Varios meses después de la publicación de Active State, los resultados de las pruebas de rendimiento que comparan los tres se publicaron y discutieron en un hilo en la lista de correo NumPy.

La función con el bucle, de hecho, superó significativamente a los otros dos:

 from numpy import mat, zeros, newaxis def calcDistanceMatrixFastEuclidean2(nDimPoints): nDimPoints = array(nDimPoints) n,m = nDimPoints.shape delta = zeros((n,n),'d') for d in xrange(m): data = nDimPoints[:,d] delta += (data - data[:,newaxis])**2 return sqrt(delta) 

Un participante en el hilo (Keir Mierle) ofreció una razón por la que esto podría ser cierto:

La razón por la que sospecho que esto será más rápido es que tiene una mejor ubicación, completando por completo un cálculo en un conjunto de trabajo relativamente pequeño antes de pasar al siguiente. Los trazadores de líneas tienen que tirar de la matriz MxN potencialmente grande en el procesador repetidamente.

Según el propio relato de este póster, su comentario es solo una sospecha, y no parece que se haya discutido más.

¿Alguna otra idea sobre cómo dar cuenta de estos resultados?

En particular, ¿hay una regla útil, que se puede extraer de este ejemplo como una guía para escribir un código numpy?

Para aquellos que no están familiarizados con NumPy, o que no han visto el código, esta comparación no se basa en un caso de vanguardia: ciertamente no sería tan interesante para mí si lo fuera. En su lugar, esta comparación involucra una función que realiza una tarea común en el cálculo matricial (es decir, la creación de una matriz de resultados con dos antecedentes); además, cada función está a su vez comprendida entre las funciones numpy más comunes.

TL; DR El segundo código anterior es solo un bucle sobre el número de dimensiones de los puntos (3 veces a través del bucle for para puntos 3D), por lo que el bucle no es mucho. La verdadera aceleración en el segundo código anterior es que aprovecha mejor el poder de Numpy para evitar crear matrices adicionales al encontrar las diferencias entre los puntos. Esto reduce la memoria utilizada y el esfuerzo computacional.

Explicación más larga Creo que la función calcDistanceMatrixFastEuclidean2 está engañando con su bucle quizás. Solo se repite sobre el número de dimensiones de los puntos. Para los puntos 1D, el bucle solo se ejecuta una vez, para 2D, dos veces, y para 3D, tres veces. Esto realmente no es mucho bucle en absoluto.

Analicemos un poco el código para ver por qué el uno es más rápido que el otro. calcDistanceMatrixFastEuclidean llamaré fast1 y calcDistanceMatrixFastEuclidean2 será fast2 .

fast1 se basa en la forma de hacer las cosas de Matlab, como lo demuestra la función repmap . La función repmap crea una matriz en este caso que es solo los datos originales repetidos una y otra vez. Sin embargo, si miras el código de la función, es muy ineficiente. Utiliza muchas funciones de Numpy (3 reshape y 2 repeat ) para hacer esto. La función de repeat también se utiliza para crear una matriz que contiene los datos originales con cada elemento de datos repetido muchas veces. Si nuestros datos de entrada son [1,2,3] entonces estamos restando [1,2,3,1,2,3,1,2,3] de [1,1,1,2,2,2,3,3,3] . Numpy ha tenido que crear muchas matrices adicionales entre la ejecución del código C de Numpy, lo que podría haberse evitado.

fast2 usa más del trabajo pesado de Numpy sin crear tantas matrices entre las llamadas de Numpy. fast2 recorre cada dimensión de los puntos, realiza la resta y mantiene un total fast2 de las diferencias cuadradas entre cada dimensión. Sólo al final se hace la raíz cuadrada. Hasta ahora, esto puede no sonar tan eficiente como fast1 , pero fast2 evita hacer las cosas de repmat mediante el uso de la indexación de Numpy. Veamos el caso 1D para la simplicidad. fast2 una matriz 1D de los datos y la resta de una matriz 2D (N x 1) de los datos. Esto crea la matriz de diferencia entre cada punto y todos los demás puntos sin tener que usar repmat y repeat y, por lo tanto, evita crear muchos arreglos adicionales. Aquí es donde la diferencia de velocidad real radica en mi opinión. fast1 crea una gran cantidad de extras entre las matrices (y se crean de manera costosa computacional) para encontrar las diferencias entre los puntos, mientras que fast2 aprovecha mejor el poder de Numpy para evitarlos.

Por cierto, aquí hay una versión un poco más rápida de fast2 :

 def calcDistanceMatrixFastEuclidean3(nDimPoints): nDimPoints = array(nDimPoints) n,m = nDimPoints.shape data = nDimPoints[:,0] delta = (data - data[:,newaxis])**2 for d in xrange(1,m): data = nDimPoints[:,d] delta += (data - data[:,newaxis])**2 return sqrt(delta) 

La diferencia es que ya no estamos creando delta como una matriz de ceros.

dis por diversión

dis.dis(calcDistanceMatrixFastEuclidean)

  2 0 LOAD_GLOBAL 0 (len) 3 LOAD_FAST 0 (points) 6 CALL_FUNCTION 1 9 STORE_FAST 1 (numPoints) 3 12 LOAD_GLOBAL 1 (sqrt) 15 LOAD_GLOBAL 2 (sum) 18 LOAD_GLOBAL 3 (repmat) 21 LOAD_FAST 0 (points) 24 LOAD_FAST 1 (numPoints) 27 LOAD_CONST 1 (1) 30 CALL_FUNCTION 3 4 33 LOAD_GLOBAL 4 (repeat) 36 LOAD_FAST 0 (points) 39 LOAD_FAST 1 (numPoints) 42 LOAD_CONST 2 ('axis') 45 LOAD_CONST 3 (0) 48 CALL_FUNCTION 258 51 BINARY_SUBTRACT 52 LOAD_CONST 4 (2) 55 BINARY_POWER 56 LOAD_CONST 2 ('axis') 59 LOAD_CONST 1 (1) 62 CALL_FUNCTION 257 65 CALL_FUNCTION 1 68 STORE_FAST 2 (distMat) 5 71 LOAD_FAST 2 (distMat) 74 LOAD_ATTR 5 (reshape) 77 LOAD_FAST 1 (numPoints) 80 LOAD_FAST 1 (numPoints) 83 BUILD_TUPLE 2 86 CALL_FUNCTION 1 89 RETURN_VALUE 

dis.dis(calcDistanceMatrixFastEuclidean2)

  2 0 LOAD_GLOBAL 0 (array) 3 LOAD_FAST 0 (nDimPoints) 6 CALL_FUNCTION 1 9 STORE_FAST 0 (nDimPoints) 3 12 LOAD_FAST 0 (nDimPoints) 15 LOAD_ATTR 1 (shape) 18 UNPACK_SEQUENCE 2 21 STORE_FAST 1 (n) 24 STORE_FAST 2 (m) 4 27 LOAD_GLOBAL 2 (zeros) 30 LOAD_FAST 1 (n) 33 LOAD_FAST 1 (n) 36 BUILD_TUPLE 2 39 LOAD_CONST 1 ('d') 42 CALL_FUNCTION 2 45 STORE_FAST 3 (delta) 5 48 SETUP_LOOP 76 (to 127) 51 LOAD_GLOBAL 3 (xrange) 54 LOAD_FAST 2 (m) 57 CALL_FUNCTION 1 60 GET_ITER >> 61 FOR_ITER 62 (to 126) 64 STORE_FAST 4 (d) 6 67 LOAD_FAST 0 (nDimPoints) 70 LOAD_CONST 0 (None) 73 LOAD_CONST 0 (None) 76 BUILD_SLICE 2 79 LOAD_FAST 4 (d) 82 BUILD_TUPLE 2 85 BINARY_SUBSCR 86 STORE_FAST 5 (data) 7 89 LOAD_FAST 3 (delta) 92 LOAD_FAST 5 (data) 95 LOAD_FAST 5 (data) 98 LOAD_CONST 0 (None) 101 LOAD_CONST 0 (None) 104 BUILD_SLICE 2 107 LOAD_GLOBAL 4 (newaxis) 110 BUILD_TUPLE 2 113 BINARY_SUBSCR 114 BINARY_SUBTRACT 115 LOAD_CONST 2 (2) 118 BINARY_POWER 119 INPLACE_ADD 120 STORE_FAST 3 (delta) 123 JUMP_ABSOLUTE 61 >> 126 POP_BLOCK 8 >> 127 LOAD_GLOBAL 5 (sqrt) 130 LOAD_FAST 3 (delta) 133 CALL_FUNCTION 1 136 RETURN_VALUE 

No soy un experto en dis , pero parece que tendrías que mirar más las funciones que llama primero para saber por qué tardan un tiempo. También hay una herramienta de perfilador de rendimiento con Python, cProfile .