Aritmética de NumPy transmitida: ¿por qué es un método mucho más eficaz?

Esta pregunta es un seguimiento de mi respuesta de manera eficiente para calcular la matriz de Vandermonde .

Aquí está la configuración:

x = np.arange(5000) # an integer array N = 4 

Ahora, voy a calcular la matriz de Vandermonde de dos maneras diferentes:

 m1 = (x ** np.arange(N)[:, None]).T 

Y,

 m2 = x[:, None] ** np.arange(N) 

Prueba de cordura:

 np.array_equal(m1, m2) True 

Estos métodos son idénticos, pero su rendimiento no es:

 %timeit m1 = (x ** np.arange(N)[:, None]).T 42.7 µs ± 271 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) %timeit m2 = x[:, None] ** np.arange(N) 150 µs ± 995 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) 

Por lo tanto, el primer método, a pesar de requerir una transposición al final, aún es más de 3 veces más rápido que el segundo método.

La única diferencia es que en el primer caso, la matriz más pequeña se transmite, mientras que en el segundo caso, es la más grande .

Entonces, con una comprensión bastante decente de cómo funciona el entumecimiento, puedo adivinar que la respuesta involucraría al caché. El primer método es mucho más amigable con el caché que el segundo. Sin embargo, me gustaría una palabra oficial de alguien con más experiencia que yo.

¿Cuál podría ser la razón de este marcado contraste en los tiempos?

Aunque me temo que mi conclusión no será más fundamental que la suya (“probablemente almacenamiento en caché”), creo que puedo ayudar a centrar nuestra atención con un conjunto de pruebas más localizadas.

Considera tu problema de ejemplo:

 M,N = 5000,4 x1 = np.arange(M) y1 = np.arange(N)[:,None] x2 = np.arange(M)[:,None] y2 = np.arange(N) x1_bc,y1_bc = np.broadcast_arrays(x1,y1) x2_bc,y2_bc = np.broadcast_arrays(x2,y2) x1_cont,y1_cont,x2_cont,y2_cont = map(np.ascontiguousarray, [x1_bc,y1_bc,x2_bc,y2_bc]) 

Como puedes ver, definí un montón de matrices para comparar. x1 , y1 y x2 , y2 , respectivamente, corresponden a sus casos de prueba originales. ??_bc corresponde a las versiones de transmisión explícita de estos arreglos. Estos comparten los datos con los originales, pero tienen pasos 0 explícitos para obtener la forma adecuada. Finalmente, ??_cont son versiones contiguas de estos arreglos de transmisión, como si estuvieran construidos con np.tile .

Entonces, tanto x1_bc , y1_bc , x1_cont como y1_cont tienen forma (4, 5000) , pero mientras que los dos primeros tienen pasos de cero, los dos últimos son matrices contiguas. Para todos los efectos, tomar el poder de cualquiera de estos pares de arreglos correspondientes debería darnos el mismo resultado contiguo (como hpaulj señaló en un comentario, una transposición en sí misma es esencialmente gratuita, así que voy a ignorar esa transposición más externa en el seguimiento).

Aquí están los horarios correspondientes a su cheque original:

 In [143]: %timeit x1 ** y1 ...: %timeit x2 ** y2 ...: 52.2 µs ± 707 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) 96 µs ± 858 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) 

Aquí están los horarios para las matrices de transmisión explícita:

 In [144]: %timeit x1_bc ** y1_bc ...: %timeit x2_bc ** y2_bc ...: 54.1 µs ± 906 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) 99.1 µs ± 1.51 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each 

La misma cosa. Esto me dice que la discrepancia no se debe de alguna manera a la transición de las expresiones indexadas a las matrices de transmisión. Esto fue principalmente esperado, pero nunca está de más comprobar.

Finalmente, las matrices contiguas:

 In [146]: %timeit x1_cont ** y1_cont ...: %timeit x2_cont ** y2_cont ...: 38.9 µs ± 529 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) 45.6 µs ± 390 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) 

¡Una gran parte de la discrepancia desaparece!

Entonces, ¿por qué revisé esto? Existe una regla general de que puede trabajar con el almacenamiento en caché de la CPU si utiliza operaciones vectorizadas a lo largo de grandes dimensiones finales en Python. Para ser más específicos, para las matrices de fila mayor (“orden C”) las dimensiones finales son contiguas, mientras que para las matrices de columna mayor (“orden fortran”) las dimensiones iniciales son contiguas. Para dimensiones suficientemente grandes, arr.sum(axis=-1) debe ser más rápido que arr.sum(axis=0) para las matrices de números de fila mayor que dan o toman alguna impresión fina.

Lo que sucede aquí es que existe una gran diferencia entre las dos dimensiones (tamaño 4 y 5000, respectivamente), pero la gran asimetría de rendimiento entre los dos casos transpuestos solo ocurre en el caso de la transmisión. Mi impresión es que la transmisión utiliza pasos de 0 para construir vistas del tamaño apropiado. Estos 0 pasos implican que, en el caso más rápido, el acceso a la memoria se ve así para la matriz x larga:

 [mem0,mem1,mem2,...,mem4999, mem0,mem1,mem2,...,mem4999, ...] # and so on 

donde mem* solo denota un valor float64 de x sentado en algún lugar de la RAM. Compare esto con el caso más lento donde trabajamos con la forma (5000,4) :

 [mem0,mem0,mem0,mem0, mem1,mem1,mem1,mem1, mem2,mem2,mem2,mem2, ...] 

Mi noción ingenua es que trabajar con el primero permite que la CPU almacene en caché porciones más grandes de los valores individuales de x a la vez, por lo que el rendimiento es excelente. En este último caso, los pasos en 0 hacen que la CPU salte en la misma dirección de memoria de x cuatro veces, haciendo esto 5000 veces. Me parece razonable creer que esta configuración funciona contra el almacenamiento en caché, lo que lleva a un mal rendimiento general. Esto también estaría de acuerdo con el hecho de que los casos contiguos no muestran este impacto de rendimiento: allí la CPU tiene que trabajar con todos los 5000 * 4 valores float64 únicos, y el almacenamiento en caché puede no verse impedido por estas lecturas extrañas.

Yo también traté de mirar broadcast_arrays :

 In [121]: X,Y = np.broadcast_arrays(np.arange(4)[:,None], np.arange(1000)) In [122]: timeit X+Y 10.1 µs ± 31.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each) In [123]: X,Y = np.broadcast_arrays(np.arange(1000)[:,None], np.arange(4)) In [124]: timeit X+Y 26.1 µs ± 30.6 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) In [125]: X.shape, X.strides Out[125]: ((1000, 4), (4, 0)) In [126]: Y.shape, Y.strides Out[126]: ((1000, 4), (0, 4)) 

np.ascontiguousarray convierte las dimensiones de 0 pasos en unas completas

 In [132]: Y1 = np.ascontiguousarray(Y) In [134]: Y1.strides Out[134]: (16, 4) In [135]: X1 = np.ascontiguousarray(X) In [136]: X1.shape Out[136]: (1000, 4) 

Operar con los arreglos completos es más rápido:

 In [137]: timeit X1+Y1 4.66 µs ± 161 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each) 

Por lo tanto, hay algún tipo de penalización de tiempo para usar las matrices de 0 pasos, incluso si no se expanden explícitamente las matrices primero. Y el costo está vinculado a las formas, y posiblemente a qué dimensión se expande.

No estoy convencido de que el almacenamiento en caché sea realmente el factor más influyente aquí.

Tampoco soy un experto en informática, por lo que podría estar equivocado, pero déjame explicarte un par de observaciones. Para simplificar, estoy usando la llamada de @ hpaulj de que ‘+’ muestra esencialmente el mismo efecto que ‘**’.

Mi hipótesis de trabajo es que es la sobrecarga de los bucles externos lo que creo que son sustancialmente más caros que los bucles internos más contiguos vectorizables.

Así que primero minimicemos la cantidad de datos que se repiten, por lo que es poco probable que el almacenamiento en caché tenga mucho impacto:

 >>> from timeit import repeat >>> import numpy as np >>> >>> def mock_data(k, N, M): ... x = list(np.random.randint(0, 10000, (k, N, M))) ... y = list(np.random.randint(0, 10000, (k, M))) ... z = list(np.random.randint(0, 10000, (k, N, 1))) ... return x, y, z ... >>> k, N, M = 500, 5000, 4 >>> >>> repeat('x.pop() + y.pop()', setup='x, y, z = mock_data(k, M, N)', globals=globals(), number=k) [0.017986663966439664, 0.018148145987652242, 0.018077059998176992] >>> repeat('x.pop() + y.pop()', setup='x, y, z = mock_data(k, N, M)', globals=globals(), number=k) [0.026680009090341628, 0.026304758968763053, 0.02680662798229605] 

Aquí, ambos escenarios tienen datos contiguos, el mismo número de adiciones, pero la versión con 5000 iteraciones externas es sustancialmente más lenta. Cuando recuperamos el almacenamiento en caché aunque a través de las pruebas, la diferencia se mantiene aproximadamente igual, pero la proporción se vuelve aún más pronunciada:

 >>> repeat('x[0] + y[0]', setup='x, y, z = mock_data(k, M, N)', globals=globals(), number=k) [0.011324503924697638, 0.011121788993477821, 0.01106808998156339] >>> repeat('x[0] + y[0]', setup='x, y, z = mock_data(k, N, M)', globals=globals(), number=k) [0.020170683041214943, 0.0202067659702152, 0.020624138065613806] 

Volviendo al escenario original de “sum externa”, vemos que en el caso de dimensión larga no contigua estamos empeorando. Como no tenemos que leer más datos que en el escenario contiguo, esto no se puede explicar porque los datos no se almacenan en la memoria caché.

 >>> repeat('z.pop() + y.pop()', setup='x, y, z = mock_data(k, M, N)', globals=globals(), number=k) [0.013918839977122843, 0.01390116906259209, 0.013737019035033882] >>> repeat('z.pop() + y.pop()', setup='x, y, z = mock_data(k, N, M)', globals=globals(), number=k) [0.0335254140663892, 0.03351909795310348, 0.0335453050211072] 

Además, ambos se benefician del almacenamiento en caché de prueba:

 >>> repeat('z[0] + y[0]', setup='x, y, z = mock_data(k, M, N)', globals=globals(), number=k) [0.012061356916092336, 0.012182610924355686, 0.012071475037373602] >>> repeat('z[0] + y[0]', setup='x, y, z = mock_data(k, N, M)', globals=globals(), number=k) [0.03265167598146945, 0.03277428599540144, 0.03247103898320347] 

Desde el punto de vista de un cachista, esto no es concluyente en el mejor de los casos.

Así que echemos un vistazo a la fuente. Después de construir un NumPy actual a partir del tarball, encontrará en algún lugar del árbol casi 15000 líneas de código generado por computadora en un archivo llamado ‘loops.c’. Estos bucles son los bucles más internos de ufuncs, el bit más relevante para nuestra situación parece ser

 #define BINARY_LOOP\ char *ip1 = args[0], *ip2 = args[1], *op1 = args[2];\ npy_intp is1 = steps[0], is2 = steps[1], os1 = steps[2];\ npy_intp n = dimensions[0];\ npy_intp i;\ for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op1 += os1) /* * loop with contiguous specialization * op should be the code working on `tin in1`, `tin in2` and * storing the result in `tout * out` * combine with NPY_GCC_OPT_3 to allow autovectorization * should only be used where its worthwhile to avoid code bloat */ #define BASE_BINARY_LOOP(tin, tout, op) \ BINARY_LOOP { \ const tin in1 = *(tin *)ip1; \ const tin in2 = *(tin *)ip2; \ tout * out = (tout *)op1; \ op; \ } etc. 

En nuestro caso, la carga útil parece lo suficientemente sencilla, especialmente si interpreto correctamente el comentario sobre la especialización contigua y la autovectorización. Ahora, si hacemos solo 4 iteraciones, la relación entre los gastos generales y la carga útil empieza a parecer un poco preocupante y no termina aquí.

En el archivo ufunc_object.c encontramos el siguiente fragmento de código

 /* * If no trivial loop matched, an iterator is required to * resolve broadcasting, etc */ NPY_UF_DBG_PRINT("iterator loop\n"); if (iterator_loop(ufunc, op, dtypes, order, buffersize, arr_prep, arr_prep_args, innerloop, innerloopdata) < 0) { return -1; } return 0; 

el bucle real parece

  NPY_BEGIN_THREADS_NDITER(iter); /* Execute the loop */ do { NPY_UF_DBG_PRINT1("iterator loop count %d\n", (int)*count_ptr); innerloop(dataptr, count_ptr, stride, innerloopdata); } while (iternext(iter)); NPY_END_THREADS; 

innerloop es el bucle interno que vimos arriba. ¿Cuánta sobrecarga viene con el siguiente?

Para esto necesitamos pasar al archivo nditer_templ.c.src donde encontramos

 /*NUMPY_API * Compute the specialized iteration function for an iterator * * If errmsg is non-NULL, it should point to a variable which will * receive the error message, and no Python exception will be set. * This is so that the function can be called from code not holding * the GIL. */ NPY_NO_EXPORT NpyIter_IterNextFunc * NpyIter_GetIterNext(NpyIter *iter, char **errmsg) { etc. 

Esta función devuelve un puntero a una de las cosas que hace el preprocesamiento

 /* Specialized iternext (@const_itflags@,@tag_ndim@,@tag_nop@) */ static int npyiter_iternext_itflags@tag_itflags@_dims@tag_ndim@_iters@tag_nop@( NpyIter *iter) { etc. 

El análisis de esto está fuera de mi scope, pero en cualquier caso es un puntero de función que debe llamarse en cada iteración del bucle externo, y por lo que sé, los punteros de función no pueden estar en línea, por lo que se comparan con 4 iteraciones de un cuerpo de bucle trivial. Será sustantivo.

Probablemente debería hacer un perfil de esto, pero mis habilidades son insuficientes.