¿Por qué la transmisión de python en el siguiente ejemplo es más lenta que un simple bucle?

Tengo una gran variedad de vectores y calculo la norma de sus diferencias frente a la primera. Cuando se utiliza la transmisión de python, el cálculo es significativamente más lento que hacerlo a través de un simple bucle. ¿Por qué?

import numpy as np def norm_loop(M, v): n = M.shape[0] d = np.zeros(n) for i in range(n): d[i] = np.sum((M[i] - v)**2) return d def norm_bcast(M, v): n = M.shape[0] d = np.zeros(n) d = np.sum((M - v)**2, axis=1) return d M = np.random.random_sample((1000, 10000)) v = M[0] 

% timeit norm_loop (M, v) -> 25.9 ms

% timeit norm_bcast (M, v) -> 38.5 ms

Tengo Python 3.6.3 y Numpy 1.14.2

Para ejecutar el ejemplo en google colab: https://drive.google.com/file/d/1GKzpLGSqz9eScHYFAuT8wJt4UIZ3ZTru/view?usp=sharing

Acceso a la memoria.

En primer lugar, la versión de transmisión se puede simplificar para

 def norm_bcast(M, v): return np.sum((M - v)**2, axis=1) 

Esto todavía funciona un poco más lento que la versión en bucle. Ahora, la sabiduría convencional dice que el código vectorizado que usa la transmisión siempre debería ser más rápido, lo que en muchos casos no es cierto (aquí voy a poner descaradamente otra de mis respuestas). ¿Entonces que esta pasando?

Como he dicho, todo se reduce al acceso a la memoria.

En la versión de transmisión, cada elemento de M se resta de v. En el momento en que se procesa la última fila de M, los resultados del procesamiento de la primera fila se han desalojado de la memoria caché, por lo que para el segundo paso estas diferencias se cargan nuevamente en la memoria caché al cuadrado Finalmente, se cargan y procesan por tercera vez para la sum. Como M es bastante grande, se borran partes de la caché en cada paso para acomodar todos los datos.

En la versión en bucle, cada fila se procesa completamente en un paso más pequeño, lo que lleva a menos errores de caché y un código en general más rápido.

Por último, es posible evitar esto con algunas operaciones de matriz utilizando einsum . Esta función permite mezclar matrices y sums. Primero, señalaré que es una función que tiene una syntax poco intuitiva en comparación con el rest de números, y que las mejoras potenciales a menudo no valen el esfuerzo extra para entenderlas. La respuesta también puede ser ligeramente diferente debido a errores de redondeo. En este caso se puede escribir como

 def norm_einsum(M, v): tmp = Mv return np.einsum('ij,ij->i', tmp, tmp) 

Esto lo reduce a dos operaciones en toda la matriz: una resta y una llamada a einsum , que realiza la cuadratura y la sum. Esto le da una ligera mejora:

 %timeit norm_bcast(M, v) 30.1 ms ± 116 µs per loop (mean ± std. dev. of 7 runs, 10 loops each) %timeit norm_loop(M, v) 25.1 ms ± 37.3 µs per loop (mean ± std. dev. of 7 runs, 10 loops each) %timeit norm_einsum(M, v) 21.7 ms ± 65.3 µs per loop (mean ± std. dev. of 7 runs, 10 loops each) 

Exprimiendo el máximo rendimiento

En operaciones vectorizadas claramente tienes un mal comportamiento de caché. Pero el cálculo itsef también es lento debido a que no explota las instrucciones SIMD modernas (AVX2, FMA). Afortunadamente no es realmente complicado superar estos problemas.

Ejemplo

 import numpy as np import numba as nb @nb.njit(fastmath=True,parallel=True) def norm_loop_improved(M, v): n = M.shape[0] d = np.empty(n,dtype=M.dtype) #enables SIMD-vectorization #if the arrays are not aligned M=np.ascontiguousarray(M) v=np.ascontiguousarray(v) for i in nb.prange(n): dT=0. for j in range(v.shape[0]): dT+=(M[i,j]-v[j])*(M[i,j]-v[j]) d[i]=dT return d 

Actuación

 M = np.random.random_sample((1000, 1000)) norm_loop_improved: 0.11 ms**, 0.28ms norm_loop: 6.56 ms norm_einsum: 3.84 ms M = np.random.random_sample((10000, 10000)) norm_loop_improved:34 ms norm_loop: 223 ms norm_einsum: 379 ms 

** Tenga cuidado al medir el rendimiento

El primer resultado (0.11ms) viene de llamar a la función repetidamente con los mismos datos. Esto necesitaría 77 GB / s de lectura-grupo de RAM, que es mucho más de lo que mi DDR3 Dualchannel-RAM es capaz de hacer. Debido a que llamar a una función con los mismos parámetros de entrada sucesivamente no es realista en absoluto, tenemos que modificar la medición.

Para evitar este problema, debemos llamar a la misma función con datos diferentes al menos dos veces (8 MB de caché L3, 8 MB de datos) y luego dividir el resultado entre dos para borrar todos los cachés.

El rendimiento relativo de estos métodos también difiere en los tamaños de matriz (eche un vistazo a los resultados de einsum).