¿Por qué Looping Beat Indexing aquí?

Question

Hace unos años, alguien posted en Recetas de estado activo con fines de comparación, tres funciones de python/NumPy; cada uno de estos aceptó los mismos argumentos y arrojó el mismo resultado, una matriz de distancia .

Dos de estos fueron tomados de fuentes publicadas; ambos son, o me parecen ser, un código idiomático numpy. Los cálculos repetitivos necesarios para crear una matriz de distancia se basan en la sintaxis del índice elegante de numpy. He aquí 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)) 

La tercera creada la matriz de distancia utilizando un solo bucle (que, obviamente, una gran cantidad de looping, dado que una matriz de distancia de tan sólo 1.000 puntos 2D, tiene un millón de entradas). A primera vista, esta función me pareció como el código que solía escribir cuando estaba aprendiendo NumPy y escribía el código de NumPy escribiendo primero el código de Python y luego traduciéndolo, línea por línea.

Varios meses después de la publicación de estado activo, los resultados de las pruebas de rendimiento que comparaban los tres se publicaron y discutieron en un thread en la lista de correo de NumPy.

La función con el bucle de hecho significativamente superaron 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) 

Uno de los participantes en el hilo (Keir Mierle) ofreció una razón por qué esto podría ser verdad:

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

Por la propia cuenta de este cartel, su observación 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, ¿existe una regla útil sobre cuándo realizar un bucle y cuándo indexar que se puede extraer de este ejemplo como guía para escribir código numpy?

Para aquellos que no están familiarizados con NumPy, o que no han mirado el código, esta comparación no se basa en un caso marginal, sin duda no sería tan interesante para mí. En cambio, esta comparación implica una función que realiza una tarea común en el cálculo de la matriz (es decir, la creación de un conjunto de resultados dados dos antecedentes); además, cada función está compuesta a su vez por las más comunes numpy incorporadas.

Answer 1

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

más largo creo que la función calcDistanceMatrixFastEuclidean2 le está engañando con su bucle, tal vez. Solo se trata de recorrer el número de dimensiones de los puntos. Para los puntos 1D, el ciclo 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é uno es más rápido que el otro. calcDistanceMatrixFastEuclidean Llamaré fast1 y calcDistanceMatrixFastEuclidean2 será fast2.

fast1 se basa en la manera 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 la información original que se repite una y otra vez. Sin embargo, si observa el código de la función, es muy ineficiente. Utiliza muchas funciones de Numpy (3 reshape sy 2 repeat s) para hacer esto. La función repeat también se usa 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 ejecutar el código C de Numpy, que podrían haberse evitado.

fast2 utiliza más de la elevación pesada 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 acumulado de las diferencias cuadradas entre cada dimensión. Solo al final se completa la raíz cuadrada. Hasta ahora, esto puede no sonar tan eficiente como fast1, pero fast2 evita hacer las cosas repmat usando la indexación de Numpy. Veamos el caso 1D por simplicidad. fast2 hace 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 otros puntos sin tener que usar repmat y repeat y por lo tanto pasa por alto la creación de muchas matrices adicionales. Aquí es donde radica la verdadera diferencia de velocidad en mi opinión. fast1 crea un montón de extra entre matrices (y se crean costosamente computacionalmente) para encontrar las diferencias entre los puntos, mientras que fast2 aprovecha mejor la potencia de Numpy para evitar esto.

Por cierto, aquí es un poco más rápido versión 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.

Answer 2

dis para la 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ía que buscar más en el f Unciones que el primero está llamando para saber por qué tardan un poco. También hay una herramienta de perfilador de rendimiento con Python, cProfile.