python - ¿Por qué la escritura lineal de lectura aleatoria no es más rápida que la escritura de lectura lineal aleatoria?

performance numpy (5)

El siguiente experimento corrobora que las escrituras aleatorias son más rápidas que las lecturas aleatorias. Para tamaños pequeños de los datos (cuando se ajusta completamente en cachés) el código de escritura aleatorio es más lento que el de lectura aleatoria (probablemente debido a ciertas peculiaridades de implementación numpy ), pero a medida que el tamaño de los datos aumenta la diferencia inicial de 1.7x en el tiempo de ejecución se elimina casi por completo (sin embargo, en caso de que numba haya una extraña reversión de esa tendencia al final).

$ cat test.py import numpy as np from timeit import timeit import numba def fwd(a,b,c): c = b[a] def inv(a,b,c): c[a] = b @numba.njit def fwd_numba(a,b,c): for i,j in enumerate(a): c[i] = b[j] @numba.njit def inv_numba(a,b,c): for i,j in enumerate(a): c[j] = b[i] for p in range(4, 8): N = 10**p n = 10**(9-p) a = np.random.permutation(N) b = np.random.random(N) c = np.empty_like(b) print(''---- N = %d ----'' % N) for f in ''fwd'', ''fwd_numba'', ''inv'', ''inv_numba'': print(f, timeit(f+''(a,b,c)'', number=n, globals=globals())) $ python test.py ---- N = 10000 ---- fwd 1.1199337750003906 fwd_numba 0.9052993479999714 inv 1.929507338001713 inv_numba 1.5510062070025015 ---- N = 100000 ---- fwd 1.8672701190007501 fwd_numba 1.5000483989970235 inv 2.509873716000584 inv_numba 2.0653326050014584 ---- N = 1000000 ---- fwd 7.639554155000951 fwd_numba 5.673054756000056 inv 7.685382894000213 inv_numba 5.439735023999674 ---- N = 10000000 ---- fwd 15.065879136000149 fwd_numba 12.68919651500255 inv 15.433822674000112 inv_numba 14.862108078999881

Actualmente estoy tratando de comprender mejor los problemas de rendimiento relacionados con la memoria / caché. Leí en algún lugar que la localidad de la memoria es más importante para la lectura que para la escritura, porque en el primer caso la CPU tiene que esperar los datos, mientras que en el último caso puede enviarlos y olvidarse de ellos.

Con eso en mente, hice la siguiente prueba rápida y sucia: escribí un script que crea una matriz de N flotantes aleatorios y una permutación, es decir, una matriz que contiene los números 0 a N-1 en orden aleatorio. Luego, repetidamente (1) lee la matriz de datos de forma lineal y la vuelve a escribir en una nueva matriz en el patrón de acceso aleatorio dado por la permutación o (2) lee la matriz de datos en el orden permutado y la escribe de forma lineal en una nueva matriz.

Para mi sorpresa (2) parecía consistentemente más rápido que (1). Hubo, sin embargo, problemas con mi guión.

El script está escrito en python / numpy. Al ser un lenguaje de alto nivel, no está claro cómo se implementa específicamente la lectura / escritura.
Probablemente no equilibré los dos casos correctamente.

Además, algunas de las respuestas / comentarios a continuación sugieren que mi expectativa original no es correcta y que, dependiendo de los detalles de la caché de la CPU, cualquier caso podría ser más rápido.

Mi pregunta es:

¿Cuál (si alguno) de los dos debería ser más rápido?
¿Cuáles son los conceptos de caché relvant aquí; ¿Cómo influyen en el resultado?

Se agradecería una explicación amigable para los principiantes. Cualquier código de soporte debe estar en C / cython / numpy / numba o python.

Opcionalmente:

Explique por qué las duraciones absolutas no son lineales en el tamaño del problema (consulte los tiempos a continuación).
Explique el comportamiento de mis experimentos con pitones claramente inadecuados.

Para referencia, mi plataforma es Linux-4.12.14-lp150.11-default-x86_64-with-glibc2.3.4 . La versión de Python es 3.6.5.

Aquí está el código que escribí:

import numpy as np from timeit import timeit def setup(): global a, b, c a = np.random.permutation(N) b = np.random.random(N) c = np.empty_like(b) def fwd(): c = b[a] def inv(): c[a] = b N = 10_000 setup() timeit(fwd, number=100_000) # 1.4942631321027875 timeit(inv, number=100_000) # 2.531870319042355 N = 100_000 setup() timeit(fwd, number=10_000) # 2.4054739447310567 timeit(inv, number=10_000) # 3.2365565397776663 N = 1_000_000 setup() timeit(fwd, number=1_000) # 11.131387163884938 timeit(inv, number=1_000) # 14.19817715883255

Como lo señalaron @Trilarion y @Yann Vernier, mis fragmentos no están bien equilibrados, así que los reemplacé con

def fwd(): c[d] = b[a] b[d] = c[a] def inv(): c[a] = b[d] b[a] = c[d]

donde d = np.arange(N) (Yo mezclo todas las formas para, con suerte, reducir a través de los efectos de almacenamiento en caché de prueba). También reemplacé el timeit con repeat y reduje el número de repeticiones por un factor de 10.

Entonces me pongo

[0.6757169323973358, 0.6705542299896479, 0.6702114241197705] #fwd [0.8183442652225494, 0.8382121799513698, 0.8173762648366392] #inv [1.0969422250054777, 1.0725746559910476, 1.0892365919426084] #fwd [1.0284497970715165, 1.025063106790185, 1.0247828317806125] #inv [3.073981977067888, 3.077839042060077, 3.072118630632758] #fwd [3.2967213969677687, 3.2996009718626738, 3.2817375687882304] #inv

Así que todavía parece haber una diferencia, pero es mucho más sutil y ahora puede ir en cualquier dirección dependiendo del tamaño del problema.

Este es un problema complejo estrechamente relacionado con las características arquitectónicas de los procesadores modernos y su intuición de que las lecturas aleatorias son más lentas que las escrituras aleatorias porque la CPU tiene que esperar a que no se verifiquen los datos de lectura (la mayoría de las veces). Hay varias razones para eso que detallaré.

Los procesadores modernos son muy eficientes para ocultar la latencia de lectura
mientras que las escrituras de memoria son más caras que las de memoria
especialmente en un entorno multinúcleo

Razón # 1 Los procesadores modernos son eficientes para ocultar la latencia de lectura.

El superscalar moderno puede ejecutar varias instrucciones simultáneamente y cambiar el orden de ejecución de la instrucción (ejecución fuera de orden ). Si bien la primera razón de estas características es aumentar el rendimiento de la instrucción, una de las consecuencias más interesantes es la capacidad de los procesadores para ocultar la latencia de las escrituras de memoria (o de operadores complejos, sucursales, etc.).

Para explicar eso, consideremos un código simple que copia la matriz en otro.

for i in a: c[i] = b[i]

Uno compilado, el código ejecutado por el procesador será así.

#1. (iteration 1) c[0] = b[0] 1a. read memory at b[0] and store result in register c0 1b. write register c0 at memory address c[0] #2. (iteration 2) c[1] = b[1] 2a. read memory at b[1] and store result in register c1 2b. write register c1 at memory address c[1] #1. (iteration 2) c[2] = b[2] 3a. read memory at b[2] and store result in register c2 3b. write register c2 at memory address c[2] # etc

(esto está terriblemente simplificado en exceso y el código real es más complejo y tiene que lidiar con la administración de bucles, el cálculo de direcciones, etc., pero este modelo simplista es actualmente suficiente).

Como se dijo en la pregunta, para las lecturas, el procesador tiene que esperar los datos reales. De hecho, 1b necesita los datos obtenidos por 1a y no puede ejecutarse mientras no se complete 1a. Dicha restricción se denomina dependencia y podemos decir que 1b depende de 1a. Las dependencias es una noción importante en los procesadores modernos. Las dependencias expresan el algoritmo (por ejemplo, escribo b a c) y deben respetarse absolutamente. Pero, si no hay dependencia entre las instrucciones, los procesadores intentarán ejecutar otras instrucciones pendientes para mantener el canal operativo siempre activo. Esto puede llevar a la ejecución fuera de orden, siempre y cuando se respeten las dependencias (similar a la regla del tipo de si mismo).

Para el código considerado, no hay dependencia entre las instrucciones de alto nivel 2. y 1. (o entre las instrucciones de asm 2a y 2b y las instrucciones anteriores). En realidad, el resultado final incluso sería idéntico si 2. se ejecuta antes de 1., y el procesador intentará ejecutar 2a y 2b, antes de completar 1a y 1b. Todavía hay una dependencia entre 2a y 2b, pero se pueden emitir ambas. Y de manera similar para 3a. y 3b., y así sucesivamente. Este es un medio poderoso para ocultar la latencia de la memoria . Si, por alguna razón, 2., 3. y 4. pueden terminar antes de que 1. cargue sus datos, es posible que incluso no note ninguna desaceleración.

Este paralelismo de nivel de instrucción es administrado por un conjunto de "colas" en el procesador.

una cola de instrucciones pendientes en las estaciones de reserva RS (tipo 128 μinstructions in pentiums recientes). Tan pronto como los recursos requeridos por la instrucción estén disponibles (por ejemplo, el valor del registro c1 para la instrucción 1b), la instrucción puede ejecutarse.
una cola de accesos de memoria pendientes en el búfer de orden de memoria MOB antes de la caché L1. Esto es necesario para tratar con los alias de la memoria y para asegurar la secuencia en las escrituras o cargas de la memoria en la misma dirección (tipo 64 cargas, 32 tiendas)
una cola para imponer la secuencia al escribir resultados de resultados en registros (reordenar el búfer o ROB de 168 entradas) por razones similares.
y algunas otras colas en la búsqueda de instrucciones, para la generación de μops, escritura y pérdida de búferes en el caché, etc.

En un punto de ejecución del programa anterior habrá muchas instrucciones de tiendas pendientes en RS, varias cargas en MOB e instrucciones a la espera de retirarse en el ROB.

Tan pronto como un dato esté disponible (por ejemplo, una lectura termina) dependiendo de las instrucciones se pueden ejecutar y eso libera las posiciones en las colas. Pero si no se produce una terminación, y una de estas colas está llena, la unidad funcional asociada con esta cola se detiene (esto también puede ocurrir en el caso de una instrucción si al procesador le faltan nombres de registro). Las paradas son lo que crea pérdida de rendimiento y, para evitarlo, el llenado de la cola debe estar limitado.

Esto explica la diferencia entre los accesos de memoria lineales y aleatorios.
En un acceso lineal, 1 / el número de fallos será menor debido a la mejor ubicación espacial y porque los cachés pueden obtener accesos previos con un patrón regular para reducirlos aún más y 2 / cada vez que finalice una lectura, se tratará de una línea de caché completa y Puede liberar varias instrucciones de carga pendientes limitando el llenado de las colas de instrucciones. De esta manera, el procesador está permanentemente ocupado y la latencia de la memoria está oculta.
Para un acceso aleatorio, la cantidad de fallas será mayor y solo se podrá atender una carga cuando lleguen los datos. Por lo tanto, las colas de instrucciones se saturarán rápidamente, el procesador se detiene y la latencia de la memoria ya no se puede ocultar ejecutando otras instrucciones.

La arquitectura del procesador debe estar equilibrada en términos de rendimiento para evitar la saturación de la cola y los bloqueos. De hecho, en general, hay decenas de instrucciones en alguna etapa de ejecución en un procesador y el rendimiento global (es decir, la capacidad para atender solicitudes de instrucción por parte de la memoria (o unidades funcionales)) es el factor principal que determinará el rendimiento. El hecho de que algunas de estas instrucciones pendientes estén esperando un valor de memoria tiene un efecto menor ...

... excepto si tiene largas cadenas de dependencia.

Hay una dependencia cuando una instrucción tiene que esperar a que se complete una anterior. Usar el resultado de una lectura es una dependencia. Y las dependencias pueden ser un problema cuando se involucran en una cadena de dependencia.

Por ejemplo, considere el código for i in range(1,100000): s += a[i] . Todas las lecturas de memoria son independientes, pero hay una cadena de dependencia para la acumulación en s . No puede ocurrir ninguna adición hasta que la anterior haya terminado. Estas dependencias harán que las estaciones de reserva se llenen rápidamente y crearán paradas en la tubería.

Pero las lecturas rara vez están involucradas en las cadenas de dependencia. Todavía es posible imaginar un código patológico en el que todas las lecturas dependan del anterior (por ejemplo, for i in range(1,100000): s = a[s] ), pero no son comunes en el código real. Y el problema proviene de la cadena de dependencia, no del hecho de que es una lectura; la situación sería similar (e incluso probablemente peor) con el código dependiente del límite de cómputo como for i in range(1,100000): x = 1.0/x+1.0 .

Por lo tanto, excepto en algunas situaciones, el tiempo de cálculo está más relacionado con el rendimiento que con la dependencia de lectura, gracias al hecho de que la superscalar o la ejecución de una orden ocultan la latencia. Y para lo que se refiere al rendimiento, las escrituras son peores que las leídas.

Razón # 2: las escrituras en memoria (especialmente las aleatorias) son más caras que las lecturas en memoria

Esto está relacionado con la forma en que se comportan los caches . La memoria caché es una memoria rápida que almacena una parte de la memoria (llamada línea ) por el procesador. Las líneas de caché tienen actualmente 64 bytes y permiten explotar la localidad espacial de las referencias de memoria: una vez que se almacena una línea, todos los datos en la línea están disponibles de inmediato. El aspecto importante aquí es que todas las transferencias entre el caché y la memoria son líneas .

Cuando un procesador realiza una lectura en un dato, el caché comprueba si la línea a la que pertenece el dato está en el caché. Si no, la línea se recupera de la memoria, se almacena en la memoria caché y los datos deseados se envían de vuelta al procesador.

Cuando un procesador escribe datos en la memoria, el caché también verifica la presencia de la línea. Si la línea no está presente, el caché no puede enviar sus datos a la memoria (porque todas las transferencias se basan en la línea) y realiza los siguientes pasos:

caché recupera la línea de la memoria y la escribe en la línea de caché.
los datos se escriben en el caché y la línea completa se marca como modificada (sucia)
cuando se suprime una línea de la memoria caché, comprueba el indicador modificado, y si la línea se ha modificado, la vuelve a escribir en la memoria (caché de escritura)

Por lo tanto, cada escritura de memoria debe ir precedida por una lectura de memoria para obtener la línea en el caché. Esto agrega una operación adicional, pero no es muy costoso para escrituras lineales. Habrá una falta de memoria caché y una lectura de memoria para la primera palabra escrita, pero las escrituras sucesivas solo se referirán a la memoria caché y serán aciertos.

Pero la situación es muy diferente para las escrituras al azar. Si el número de errores es importante, cada error de caché implica una lectura seguida de solo un pequeño número de escrituras antes de que la línea sea expulsada del caché, lo que aumenta significativamente el costo de escritura. Si una línea se expulsa después de una sola escritura, incluso podemos considerar que una escritura es el doble del costo temporal de una lectura.

Es importante tener en cuenta que el aumento del número de accesos a la memoria (tanto de lectura como de escritura) tiende a saturar la ruta de acceso a la memoria y ralentiza globalmente todas las transferencias entre el procesador y la memoria.

En cualquier caso, las escrituras son siempre más caras que las leídas. Y los multicores aumentan este aspecto.

Razón # 3: Escrituras aleatorias crean fallos de caché en multinúcleos

No estoy seguro de que esto se aplique realmente a la situación de la pregunta. Si bien las rutinas de BLAS numpy son multiproceso, no creo que la copia de matriz básica sea. Pero está estrechamente relacionado y es otra razón por la que las escrituras son más caras.

El problema con los multicores es garantizar la coherencia adecuada de la memoria caché de tal manera que los datos compartidos por varios procesadores se actualicen correctamente en la memoria caché de cada núcleo. Esto se hace por medio de un protocolo como MESI que actualiza una línea de caché antes de escribirla, e invalida otras copias de caché (leer para el propietario).

Si bien ninguno de los datos se comparte realmente entre los núcleos de la pregunta (o una versión paralela), tenga en cuenta que el protocolo se aplica a las líneas de caché . Cada vez que se modifica una línea de caché, se copia del caché que contiene la copia más reciente, se actualiza localmente y todas las demás copias se invalidan. Incluso si los núcleos están accediendo a diferentes partes de la línea de caché. Dicha situación se denomina intercambio falso y es un tema importante para la programación multinúcleo.

Con respecto al problema de las escrituras aleatorias, las líneas de caché tienen 64 bytes y pueden contener 8 int64, y si la computadora tiene 8 núcleos, cada núcleo procesará en el promedio de 2 valores. Por lo tanto, hay un intercambio falso importante que ralentizará las escrituras.

Hicimos algunas evaluaciones de desempeño. Se realizó en C para incluir una evaluación del impacto de la paralelización. Comparamos 5 funciones que procesan arreglos int64 de tamaño N.

Solo una copia de b a c ( c[i] = b[i] ) (implementado por el compilador con memcpy() )
Copie con un índice lineal c[i] = b[d[i]] donde d[i]==i ( read_linear )
Copie con un índice aleatorio c[i] = b[a[i]] donde a es una permutación aleatoria de 0..N-1 ( read_random es equivalente a fwd en la pregunta original)
Escribe c[d[i]] = b[i] lineal c[d[i]] = b[i] donde d[i]==i ( write_linear )
Escriba al azar c[a[i]] = b[i] con a permutación aleatoria de 0..N-1 ( write_random es equivalente a inv en la pregunta)

El código se ha compilado con gcc -O3 -funroll-loops -march=native -malign-double en un procesador skylake. Las actuaciones se miden con _rdtsc() y se dan en ciclos por iteración. La función se ejecuta varias veces (1000-20000 dependiendo del tamaño de la matriz), se realizan 10 experimentos y se mantiene el menor tiempo.

El tamaño de los arreglos varía de 4000 a 1200000. Todo el código se ha medido con una versión secuencial y una versión paralela con openmp.

Aquí hay una gráfica de los resultados. Las funciones son con diferentes colores, con la versión secuencial en líneas gruesas y la paralela con las delgadas.

La copia directa es (obviamente) la más rápida y es implementada por gcc con el memcpy() altamente optimizado. Es un medio para obtener una estimación del rendimiento de datos con memoria. Va desde 0,8 ciclos por iteración (CPI) para matrices pequeñas hasta 2,0 CPI para las grandes.

Las interpretaciones lineales de lectura son aproximadamente dos veces más largas que memcpy, pero hay 2 lecturas y una escritura, frente a 1 lectura y una escritura para la copia directa. Más el índice añade alguna dependencia. El valor mínimo es 1.56 CPI y el valor máximo 3.8 CPI. Escribir lineal es ligeramente más largo (5-10%).

Las lecturas y escrituras con un índice aleatorio son el propósito de la pregunta original y merecen un comentario más largo. Aquí están los resultados.

size 4000 6000 9000 13496 20240 30360 45536 68304 102456 153680 230520 345776 518664 777992 1166984 rd-rand 1.86821 2.52813 2.90533 3.50055 4.69627 5.10521 5.07396 5.57629 6.13607 7.02747 7.80836 10.9471 15.2258 18.5524 21.3811 wr-rand 7.07295 7.21101 7.92307 7.40394 8.92114 9.55323 9.14714 8.94196 8.94335 9.37448 9.60265 11.7665 15.8043 19.1617 22.6785

valores pequeños (<10k): el caché L1 es 32k y puede contener una matriz de 4k de uint64. Tenga en cuenta que, debido a la aleatoriedad del índice, después de ~ 1/8 de iteraciones, la memoria caché L1 se llenará completamente con los valores de la matriz de índice aleatorio (ya que las líneas de la memoria caché son de 64 bytes y pueden contener 8 elementos de matriz). Al acceder a las otras matrices lineales, generaremos rápidamente muchos errores de L1 y tendremos que utilizar el caché de L2. El acceso al caché L1 es de 5 ciclos, pero está segmentado y puede servir un par de valores por ciclo. El acceso a L2 es más largo y requiere 12 ciclos. La cantidad de fallos es similar para las lecturas y escrituras aleatorias, pero vemos que pagamos completamente el doble acceso requerido para las escrituras cuando el tamaño de la matriz es pequeño.
valores medios (10k-100k): el caché L2 es 256k y puede contener una matriz de 32k int64. Después de eso, tenemos que ir a la caché L3 (12Mo). A medida que aumenta el tamaño, aumenta el número de fallos en L1 y L2 y el tiempo de cálculo correspondiente. Ambos algoritmos tienen un número similar de fallas, principalmente debido a lecturas o escrituras aleatorias (otros accesos son lineales y pueden ser captados de manera muy eficiente por los cachés). Recuperamos el factor dos entre lecturas aleatorias y escrituras ya anotadas en la respuesta de BM. Se puede explicar en parte por el doble costo de las escrituras.
Valores grandes (> 100k): la diferencia entre métodos se reduce progresivamente. Para estos tamaños, una gran parte de la información se almacena en el caché L3. El tamaño de L3 es suficiente para mantener una matriz completa de 1.5M y es menos probable que las líneas sean expulsadas. Por lo tanto, para las escrituras, después de la lectura inicial, se puede hacer un mayor número de escrituras sin expulsión de línea, y se reduce el costo relativo de las escrituras frente a las leídas. Para estos tamaños grandes, también hay muchos otros factores que deben ser considerados. Por ejemplo, las memorias caché solo pueden servir un número limitado de errores (tipo 16) y cuando el número de errores es grande, este puede ser el factor limitante.

Una palabra en la versión omp paralela de lecturas y escrituras aleatorias. A excepción de los tamaños pequeños, donde tener la matriz de índice aleatorio repartida en varios cachés puede no ser una ventaja, son sistemáticamente ~ dos veces más rápidos. Para tamaños grandes, vemos claramente que la brecha entre lecturas y escrituras aleatorias aumenta debido a la falsa compartición.

Es casi imposible hacer predicciones cuantitativas con la complejidad de las arquitecturas informáticas actuales, incluso para códigos simples, e incluso las explicaciones cualitativas del comportamiento son difíciles y deben tener en cuenta muchos factores. Como se mencionó en otras respuestas, los aspectos de software relacionados con python también pueden tener un impacto. Pero, aunque puede suceder en algunas situaciones, la mayoría de las veces, no se puede considerar que las lecturas son más caras debido a la dependencia de los datos.

Su función fwd no está tocando la variable global c . No le dijiste global c (solo en la setup ), por lo que tiene su propia variable local, y usa STORE_FAST en cpython:

>>> import dis >>> def fwd(): ... c = b[a] ... >>> dis.dis(fwd) 2 0 LOAD_GLOBAL 0 (b) 3 LOAD_GLOBAL 1 (a) 6 BINARY_SUBSCR 7 STORE_FAST 0 (c) 10 LOAD_CONST 0 (None) 13 RETURN_VALUE

Ahora, vamos a intentarlo con un global:

>>> def fwd2(): ... global c ... c = b[a] ... >>> dis.dis(fwd2) 3 0 LOAD_GLOBAL 0 (b) 3 LOAD_GLOBAL 1 (a) 6 BINARY_SUBSCR 7 STORE_GLOBAL 2 (c) 10 LOAD_CONST 0 (None) 13 RETURN_VALUE

Aun así, puede diferir en el tiempo en comparación con la función inv que llama a setitem para un global.

De cualquier manera, si desea que escriba en c , necesita algo como c[:] = b[a] o c.fill(b[a]) . La asignación reemplaza la variable (nombre) con el objeto del lado derecho, por lo que la antigua c podría ser desasignada en lugar de la nueva b[a] , y ese tipo de mezcla de memoria puede ser costoso.

En cuanto al efecto, creo que quería medir, básicamente si las permutaciones directas o inversas son más costosas, eso sería altamente dependiente de la caché. La permutación directa (el almacenamiento en índices ordenados aleatoriamente a partir de una lectura lineal) en principio podría ser más rápida porque puede usar el enmascaramiento de escritura y nunca obtener la nueva matriz, asumiendo que el sistema de caché es lo suficientemente inteligente como para conservar las máscaras de bytes en el búfer de escritura. Al revés, existe un alto riesgo de colisiones de caché mientras se realiza la lectura aleatoria si la matriz es lo suficientemente grande.

Esa fue mi impresión inicial; Los resultados, como usted dice, son opuestos. Esto podría ser el resultado de una implementación de caché que no tiene un búfer de escritura grande o que no puede explotar escrituras pequeñas. De todos modos, si los accesos fuera del caché requieren la misma hora del bus de memoria, el acceso de lectura tendrá la posibilidad de cargar datos que no serán eliminados del caché antes de que sea necesario. Con un caché de múltiples vías, las líneas parcialmente escritas también tendrán la posibilidad de no ser elegidas para la expulsión; y solo las líneas de caché sucias requieren tiempo de bus de memoria para soltar. Un programa de nivel inferior escrito con otro conocimiento (por ejemplo, que la permutación es completa y no se superpone) podría mejorar el comportamiento utilizando sugerencias tales como escrituras SSE no temporales.

Sus dos fragmentos NumPy b[a] y c[a] = b parecen ser heurísticas razonables para medir velocidades de lectura / escritura aleatorias / lineales, como trataré de argumentar mirando el código NumPy subyacente en la primera sección a continuación.

Respecto a la pregunta de cuál debería ser más rápido, parece plausible que la escritura de lectura aleatoria barajada pueda ganar (como parecen mostrarse los puntos de referencia), pero la diferencia en la velocidad puede verse afectada por la "barajada" del índice barajado. , y uno o más de:

Las políticas de lectura / actualización de la memoria caché de la CPU ( reescritura frente a escritura directa , etc.).
Cómo la CPU elige (re) ordenar las instrucciones que necesita ejecutar (canalización).
La CPU reconoce los patrones de acceso a la memoria y los datos de búsqueda previa.
Lógica de desalojo de caché.

Incluso haciendo suposiciones sobre qué políticas existen, estos efectos son difíciles de modelar y razonar de manera analítica, por lo que no estoy seguro de que sea posible una respuesta general aplicable a todos los procesadores (aunque no soy un experto en hardware).

Sin embargo, en la segunda sección a continuación, trataré de razonar por qué la escritura aleatoria de lectura lineal es aparentemente más rápida, dadas algunas suposiciones.

Indexación de lujo "trivial"

El propósito de esta sección es revisar el código fuente de NumPy para determinar si hay explicaciones obvias para los tiempos, y también obtener una idea lo más clara posible de lo que sucede cuando se ejecuta A[B] o A[B] = C .

La rutina de iteración que sustenta la indexación de getitem para las operaciones getitem y setitem en esta pregunta es " trivial ":

B es una matriz de indexación única con un solo paso
A y B tienen el mismo orden de memoria (tanto C-contiguo como Fortran-contiguo)

Además, en nuestro caso, tanto A como B están alineados en Uint :

Código de copia estriada: aquí, en su lugar, se utiliza "alineación uint". Si el tamaño de elemento [N] de una matriz es igual a 1, 2, 4, 8 o 16 bytes y la matriz está alineada en uint, en su lugar [de usar buffering] numpy hará *(uintN*)dst) = *(uintN*)src) para el N. apropiado. De lo contrario, tendrá que hacer memcpy(dst, src, N) copias haciendo memcpy(dst, src, N) .

El punto aquí es que el uso de un búfer interno para evitar la alineación se evita. La copia subyacente implementada con *(uintN*)dst) = *(uintN*)src) es tan directa como "poner los X bytes de offset src en los X bytes en offset dst".

Los compiladores probablemente traducirán esto muy simplemente en instrucciones de mov (en x86 por ejemplo), o similar.

El código básico de bajo nivel que realiza la obtención y configuración de elementos se encuentra en las funciones mapiter_trivial_get y mapiter_trivial_set . Estas funciones se producen en lowlevel_strided_loops.c.src , donde las plantillas y las macros hacen que sea un poco difícil de leer (una ocasión para agradecer los idiomas de nivel superior).

Perseverando, podemos ver que hay poca diferencia entre getitem y setitem. Aquí hay una versión simplificada del bucle principal para la exposición. Las líneas de macro determinan si se estaban ejecutando getitem o setitem:

while (itersize--) { char * self_ptr; npy_intp indval = *((npy_intp*)ind_ptr); #if @isget@ if (check_and_adjust_index(&indval, fancy_dim, 0, _save) < 0 ) { return -1; } #else if (indval < 0) { indval += fancy_dim; } #endif self_ptr = base_ptr + indval * self_stride; /* offset into array being indexed */ #if @isget@ *(npy_uint64 *)result_ptr = *(npy_uint64 *)self_ptr; #else *(npy_uint64 *)self_ptr = *(npy_uint64 *)result_ptr; #endif ind_ptr += ind_stride; /* move to next item of index array */ result_ptr += result_stride; /* move to next item of result array */

Como podríamos esperar, esto simplemente equivale a cierta aritmética para obtener el desplazamiento correcto en las matrices, y luego copiar los bytes de una ubicación de memoria a otra.

Indices adicionales para el setitem

Una cosa que vale la pena mencionar es que para setitem, la validez de los índices (ya sea que estén todos dentro de la matriz de destino) se verifica antes de que comience la copia (a través de check_and_adjust_index ), que también reemplaza los índices negativos con los índices positivos correspondientes.

En el fragmento anterior se puede ver check_and_adjust_index llamado getitem en el bucle principal, mientras que para el setitem se produce una verificación más simple (posiblemente redundante) de índices negativos.

Esta comprobación preliminar adicional podría tener un impacto pequeño pero negativo en la velocidad de setitem ( A[B] = C ).

Falta de caché

Debido a que el código para ambos fragmentos de código es muy similar, la sospecha recae en la CPU y en cómo maneja el acceso a las matrices subyacentes de la memoria.

La CPU almacena en la memoria caché pequeños bloques de memoria (líneas de caché) a los que se ha accedido recientemente con la anticipación de que probablemente pronto necesitará acceder nuevamente a esa región de la memoria.

Para el contexto, las líneas de caché son generalmente de 64 bytes. La memoria caché de datos L1 (la más rápida) en la CPU de mi computadora portátil más antigua es de 32 KB (suficiente para contener alrededor de 500 valores int64 de la matriz, pero tenga en cuenta que la CPU hará otras cosas que requieren otra memoria mientras se ejecuta el fragmento NumPy):

$ cat /sys/devices/system/cpu/cpu0/cache/index0/coherency_line_size 64 $ cat /sys/devices/system/cpu/cpu0/cache/index0/size 32K

Como probablemente ya sepa, para la lectura / escritura de la memoria el almacenamiento en caché secuencialmente funciona bien porque se recuperan bloques de memoria de 64 bytes según sea necesario y se almacenan más cerca de la CPU. El acceso repetido a ese bloque de memoria es más rápido que obtener de la RAM (o un caché de nivel superior más lento). De hecho, la CPU puede incluso obtener la siguiente línea de caché antes de que el programa lo solicite.

Por otro lado, es probable que el acceso aleatorio a la memoria cause frecuentes fallas en la memoria caché. Aquí, la región de la memoria con la dirección requerida no se encuentra en la memoria caché rápida cerca de la CPU y, en su lugar, se debe acceder desde una memoria caché de nivel superior (más lenta) o desde la memoria real (mucho más lenta).

Entonces, ¿qué es más rápido para que la CPU maneje: errores frecuentes en la lectura de datos o errores en la escritura de datos?

Asumamos que la política de escritura de la CPU es de retroceso, lo que significa que una memoria modificada se escribe de nuevo en el caché. La memoria caché se marca como modificada (o "sucia"), y el cambio solo se volverá a escribir en la memoria principal una vez que la línea se haya desalojado de la memoria caché (la CPU aún puede leer desde una línea de memoria caché sucia).

Si estamos escribiendo puntos aleatorios en una gran matriz, la expectativa es que muchas de las líneas de caché en el caché de la CPU se ensucien. Será necesario escribir en la memoria principal ya que cada uno es desalojado, lo que puede ocurrir con frecuencia si el caché está lleno.

Sin embargo, esta escritura simultánea debería ocurrir con menos frecuencia cuando se escriben datos de forma secuencial y se leen de forma aleatoria, ya que esperamos que menos líneas de caché se ensucien y los datos se vuelvan a escribir en la memoria principal o que los cachés más lentos sean menos frecuentes.

Como se mencionó, este es un modelo simplificado y puede haber muchos otros factores que influyen en el rendimiento de la CPU. Alguien con más experiencia que yo podría mejorar este modelo.

Primero, una refutación de su intuición: fwd beats inv incluso sin mecanismo numpy.

Es el caso de esta versión numba :

import numba @numba.njit def fwd_numba(a,b,c): for i in range(N): c[a[i]]=b[i] @numba.njit def inv_numba(a,b,c): for i in range(N): c[i]=b[a[i]]

Tiempos para N = 10 000:

%timeit fwd() %timeit inv() %timeit fwd_numba(a,b,c) %timeit inv_numba(a,b,c) 62.6 µs ± 3.84 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each) 144 µs ± 2 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each) 16.6 µs ± 1.52 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each) 34.9 µs ± 1.57 µs per loop (mean ± std. dev. of 7 runs, 100000 loops each)

En segundo lugar, Numpy tiene que lidiar con temibles problemas de alineación y localidad (caché).

Es esencialmente una envoltura en procedimientos de bajo nivel de BLAS / ATLAS / MKL ajustados para eso. La indexación sofisticada es una buena herramienta de alto nivel pero herética para estos problemas; No hay una traducción directa de este concepto a bajo nivel.

En tercer lugar, Numpy Dev Docs: detalles de indexación de lujo. En particular:

A menos que haya una sola matriz de indexación durante la obtención del elemento, la validez de los índices se verifica de antemano. De lo contrario, se maneja en el propio bucle interno para la optimización.

Estamos en este caso aquí. Creo que esto puede explicar la diferencia, y por qué establecer es más lento que obtener.

También explica por qué el numba hecho a mano suele ser más rápido: no comprueba nada y se bloquea en un índice inconsistente.