studio - r loop rows

¿Por qué “vectorizar” este simple bucle R da un resultado diferente? (4)

Tal vez una pregunta muy tonta.

Estoy tratando de "vectorizar" el siguiente bucle:

set.seed(0) x <- round(runif(10), 2) # [1] 0.90 0.27 0.37 0.57 0.91 0.20 0.90 0.94 0.66 0.63 sig <- sample.int(10) # [1] 1 2 9 5 3 4 8 6 7 10 for (i in seq_along(sig)) x[i] <- x[sig[i]] x # [1] 0.90 0.27 0.66 0.91 0.66 0.91 0.94 0.91 0.94 0.63

Creo que es simplemente x[sig] pero el resultado no coincide.

set.seed(0) x <- round(runif(10), 2) x[] <- x[sig] x # [1] 0.90 0.27 0.66 0.91 0.37 0.57 0.94 0.20 0.90 0.63

Que pasa

Observación

Obviamente, desde la salida vemos que el bucle for y x[sig] son diferentes. El significado de este último es claro: permutación, por lo tanto, muchas personas tienden a creer que el bucle simplemente está haciendo algo incorrecto. Pero nunca estés tan seguro; Puede ser algún proceso dinámico bien definido. El propósito de estas preguntas y respuestas no es juzgar qué es correcto, sino explicar por qué no son equivalentes. Esperemos que proporcione un caso de estudio sólido para comprender la "vectorización".

calentar

Como calentamiento, considere dos ejemplos más simples.

## example 1 x <- 1:11 for (i in 1:10) x[i] <- x[i + 1] x # [1] 2 3 4 5 6 7 8 9 10 11 11 x <- 1:11 x[1:10] <- x[2:11] x # [1] 2 3 4 5 6 7 8 9 10 11 11 ## example 2 x <- 1:11 for (i in 1:10) x[i + 1] <- x[i] x # [1] 1 1 1 1 1 1 1 1 1 1 1 x <- 1:11 x[2:11] <- x[1:10] x # [1] 1 1 2 3 4 5 6 7 8 9 10

"Vectorización" es exitosa en el primer ejemplo pero no en el segundo. ¿Por qué?

Aquí está el análisis prudente. La "vectorización" comienza con el desenrollado del bucle y luego ejecuta varias instrucciones en paralelo. Si un bucle se puede "vectorizar" depende de la dependencia de datos que lleve el bucle.

Desenrollando el bucle en el ejemplo 1 da

x[1] <- x[2] x[2] <- x[3] x[3] <- x[4] x[4] <- x[5] x[5] <- x[6] x[6] <- x[7] x[7] <- x[8] x[8] <- x[9] x[9] <- x[10] x[10] <- x[11]

La ejecución de estas instrucciones una por una y la ejecución dan simultáneamente un resultado idéntico. Así que este bucle puede ser "vectorizado".

El bucle en el ejemplo 2 es

x[2] <- x[1] x[3] <- x[2] x[4] <- x[3] x[5] <- x[4] x[6] <- x[5] x[7] <- x[6] x[8] <- x[7] x[9] <- x[8] x[10] <- x[9] x[11] <- x[10]

Desafortunadamente, ejecutar estas instrucciones una por una y ejecutarlas simultáneamente no daría un resultado idéntico. Por ejemplo, al ejecutarlos uno por uno, x[2] se modifica en la primera instrucción, luego este valor modificado se pasa a x[3] en la segunda instrucción. Entonces x[3] tendría el mismo valor que x[1] . Sin embargo, en ejecución paralela, x[3] es igual a x[2] . Como resultado, este bucle no puede ser "vectorizado".

En la teoría de la "vectorización",

El ejemplo 1 tiene una dependencia de "escritura después de leer" en los datos: x[i] se modifica después de que se lee;
El ejemplo 2 tiene una dependencia de "lectura después de escribir" en los datos: x[i] se lee después de su modificación.

Un bucle con la dependencia de datos de "escritura después de la lectura" puede ser "vectorizado", mientras que un bucle con la dependencia de datos de "lectura después de la lectura" no puede.

a fondo

Tal vez muchas personas se han confundido por ahora. "Vectorización" es un "procesamiento paralelo"?

Sí. En la década de 1960, cuando la gente se preguntaba qué tipo de computadora de procesamiento paralelo se diseñaría para computación de alto rendimiento, Flynn clasificó las ideas de diseño en 4 tipos. La categoría "SIMD" (instrucción única, datos múltiples) se denomina "vectorización" y una computadora con capacidad "SIMD" se denomina "procesador de vectores" o "procesador de matrices".

En la década de 1960 no había muchos lenguajes de programación. La gente escribió el ensamblado (luego FORTRAN cuando se inventó un compilador) para programar los registros de la CPU directamente. Una computadora "SIMD" puede cargar datos múltiples en un registro vectorial con una sola instrucción y hacer la misma aritmética en esos datos al mismo tiempo. Así que el procesamiento de datos es de hecho paralelo. Considere nuestro ejemplo 1 otra vez. Supongamos que un registro vectorial puede contener dos elementos vectoriales, entonces el bucle se puede ejecutar con 5 iteraciones utilizando el procesamiento de vectores en lugar de 10 iteraciones como en el proceso escalar.

reg <- x[2:3] ## load vector register x[1:2] <- reg ## store vector register ------------- reg <- x[4:5] ## load vector register x[3:4] <- reg ## store vector register ------------- reg <- x[6:7] ## load vector register x[5:6] <- reg ## store vector register ------------- reg <- x[8:9] ## load vector register x[7:8] <- reg ## store vector register ------------- reg <- x[10:11] ## load vector register x[9:10] <- reg ## store vector register

Hoy en día hay muchos lenguajes de programación, como R. "Vectorización" ya no se refiere inequívocamente a "SIMD". R no es un lenguaje donde podamos programar registros de CPU. La "vectorización" en R es solo una analogía con "SIMD". En una sesión anterior de preguntas y respuestas: ¿El término "vectorización" significa cosas diferentes en contextos diferentes? He tratado de explicar esto. El siguiente mapa ilustra cómo se hace esta analogía:

single (assembly) instruction -> single R instruction CPU vector registers -> temporary vectors parallel processing in registers -> C/C++/FORTRAN loops with temporary vectors

Entonces, la R "vectorización" del bucle en el ejemplo 1 es algo así como

## the C-level loop is implemented by function "[" tmp <- x[2:11] ## load data into a temporary vector x[1:10] <- tmp ## fill temporary vector into x

La mayoría de las veces solo hacemos

x[1:10] <- x[2:10]

sin asignar explícitamente el vector temporal a una variable. El bloque de memoria temporal creado no es señalado por ninguna variable R y, por lo tanto, está sujeto a la recolección de basura.

una imagen completa

En lo anterior, la "vectorización" no se introduce con el ejemplo más simple. Muy a menudo, la "vectorización" se introduce con algo como

a[1] <- b[1] + c[1] a[2] <- b[2] + c[2] a[3] <- b[3] + c[3] a[4] <- b[4] + c[4]

donde a , b y c no tienen alias en la memoria, es decir, los bloques de memoria que almacenan los vectores a , c no se superponen. Este es un caso ideal, ya que ningún alias de memoria implica ninguna dependencia de datos.

Además de la "dependencia de datos", también existe la "dependencia de control", es decir, tratar con "if ... else ..." en la "vectorización". Sin embargo, por razones de tiempo y espacio no voy a dar más detalles sobre este tema.

Volviendo al ejemplo en la pregunta.

Ahora es el momento de investigar el bucle en la pregunta.

set.seed(0) x <- round(runif(10), 2) sig <- sample.int(10) # [1] 1 2 9 5 3 4 8 6 7 10 for (i in seq_along(sig)) x[i] <- x[sig[i]]

Desenrollando el bucle da

x[1] <- x[1] x[2] <- x[2] x[3] <- x[9] ## 3rd instruction x[4] <- x[5] x[5] <- x[3] ## 5th instruction x[6] <- x[4] x[7] <- x[8] x[8] <- x[6] x[9] <- x[7] x[10] <- x[10]

Hay una dependencia de datos de "lectura después de la escritura" entre la tercera y la quinta instrucciones, por lo que el bucle no puede ser "vectorizado" (ver la Observación 1 ).

Bien, entonces, ¿qué hace x[] <- x[sig] ? Primero escribamos explícitamente el vector temporal:

tmp <- x[sig] x[] <- tmp

Como "[" se llama dos veces, en realidad hay dos bucles de nivel C detrás de este código "vectorizado":

tmp[1] <- x[1] tmp[2] <- x[2] tmp[3] <- x[9] tmp[4] <- x[5] tmp[5] <- x[3] tmp[6] <- x[4] tmp[7] <- x[8] tmp[8] <- x[6] tmp[9] <- x[7] tmp[10] <- x[10] x[1] <- tmp[1] x[2] <- tmp[2] x[3] <- tmp[3] x[4] <- tmp[4] x[5] <- tmp[5] x[6] <- tmp[6] x[7] <- tmp[7] x[8] <- tmp[8] x[9] <- tmp[9] x[10] <- tmp[10]

Entonces x[] <- x[sig] es equivalente a

for (i in 1:10) tmp[i] <- x[sig[i]] for (i in 1:10) x[i] <- tmp[i] rm(tmp); gc()

que no es en absoluto el bucle original dado en la pregunta.

Observación 1

Si la implementación del bucle en Rcpp se ve como una "vectorización", entonces déjelo. Pero no hay posibilidad de "vectorizar" aún más el bucle C / C ++ con "SIMD".

Observación 2

Este Q & A está motivado por este Q & A. OP presentó originalmente un bucle

for (i in 1:num) { for (j in 1:num) { mat[i, j] <- mat[i, mat[j, "rm"]] } }

Es tentador "vectorizarlo" como

mat[1:num, 1:num] <- mat[1:num, mat[1:num, "rm"]]

pero es potencialmente incorrecto. Más tarde OP cambió el bucle a

for (i in 1:num) { for (j in 1:num) { mat[i, j] <- mat[i, 1 + num + mat[j, "rm"]] } }

lo que elimina el problema del alias de memoria, porque las columnas que se reemplazarán son las primeras columnas num , mientras que las columnas que se buscarán están después de las primeras columnas num .

Observación 3

Recibí algunos comentarios sobre si el bucle en la pregunta está realizando una modificación "in situ" de x . Sí lo es. Podemos utilizar tracemem :

set.seed(0) x <- round(runif(10), 2) sig <- sample.int(10) tracemem(x) #[1] "<0x28f7340>" for (i in seq_along(sig)) x[i] <- x[sig[i]] tracemem(x) #[1] "<0x28f7340>"

Mi sesión R ha asignado un bloque de memoria señalado por la dirección <0x28f7340> para x y es posible que vea un valor diferente cuando ejecute el código. Sin embargo, la salida de tracemem no cambiará después del bucle, lo que significa que no se realiza ninguna copia de x . Así que el bucle está haciendo una modificación "in situ" sin usar memoria extra.

Sin embargo, el bucle no está haciendo permutación "in situ". La permutación "in situ" es una operación más complicada. No solo los elementos de x deben intercambiarse a lo largo del bucle, los elementos de sig también deben intercambiarse (y al final, sig sería 1:10 ).

Es interesante observar que aunque R "vectorización" es diferente de "SIMD" (como se explica muy bien OP), la misma lógica se aplica cuando se determina si un bucle es "vectorizable". Aquí hay una demostración usando ejemplos de la auto-respuesta de OP (con una pequeña modificación).

El ejemplo 1 con la dependencia de "escribir después de leer" es "vectorizable".

// "ex1.c" #include <stdlib.h> void ex1 (size_t n, size_t *x) { for (size_t i = 1; i < n; i++) x[i - 1] = x[i] + 1; } gcc -O2 -c -ftree-vectorize -fopt-info-vec ex1.c #ex1.c:3:3: note: loop vectorized

El ejemplo 2 con la dependencia de "lectura-escritura-posterior" no es "vectorizable".

// "ex2.c" #include <stdlib.h> void ex2 (size_t n, size_t *x) { for (size_t i = 1; i < n; i++) x[i] = x[i - 1] + 1; } gcc -O2 -c -ftree-vectorize -fopt-info-vec-missed ex2.c #ex2.c:3:3: note: not vectorized, possible dependence between data-refs #ex2.c:3:3: note: bad data dependence

Use la palabra clave de restrict C99 para sugerir que el compilador no tenga alias de bloques de memoria entre tres arreglos.

// "ex3.c" #include <stdlib.h> void ex3 (size_t n, size_t * restrict a, size_t * restrict b, size_t * restrict c) { for (size_t i = 0; i < n; i++) a[i] = b[i] + c[i]; } gcc -O2 -c -ftree-vectorize -fopt-info-vec ex3.c #ex3.c:3:3: note: loop vectorized

Esto no tiene nada que ver con el alias de bloque de memoria (un término que nunca he encontrado antes). Tome un ejemplo de permutación particular y recorra las asignaciones que se producirían independientemente de la implementación en el nivel de lenguaje C o de ensamblador (o lo que sea); Es intrínseco a cómo se comportaría cualquier bucle for secuencial versus cómo se produciría cualquier permutación "verdadera" (lo que se obtiene con x[sig] ):

sample(10) [1] 3 7 1 5 6 9 10 8 4 2 value at 1 goes to 3, and now there are two of those values value at 2 goes to 7, and now there are two of those values value at 3 (which was at 1) now goes back to 1 but the values remain unchanged

... puede continuar, pero esto ilustra cómo esto usualmente no será una permutación "verdadera" y muy rara vez resultaría en una redistribución completa de los valores. Supongo que solo una permutación completamente ordenada (de la cual creo que solo hay una, es decir, 10:1 ) podría dar como resultado un nuevo conjunto de x que fuera único.

replicate( 100, {x <- round(runif(10), 2); sig <- sample.int(10); for (i in seq_along(sig)){ x[i] <- x[sig[i]]}; sum(duplicated(x)) } ) #[1] 4 4 4 5 5 5 4 5 6 5 5 5 4 5 5 6 3 4 2 5 4 4 4 4 3 5 3 5 4 5 5 5 5 5 5 5 4 5 5 5 5 4 #[43] 5 3 4 6 6 6 3 4 5 3 5 4 6 4 5 5 6 4 4 4 5 3 4 3 4 4 3 6 4 7 6 5 6 6 5 4 7 5 6 3 6 4 #[85] 8 4 5 5 4 5 5 5 4 5 5 4 4 5 4 5

Comencé a preguntarme cuál sería la distribución de los recuentos de duplicaciones en una serie grande. Parece bastante simétrico:

table( replicate( 1000000, {x <- round(runif(10), 5); sig <- sample.int(10); for (i in seq_along(sig)){ x[i] <- x[sig[i]]}; sum(duplicated(x)) } ) ) 0 1 2 3 4 5 6 7 8 1 269 13113 126104 360416 360827 125707 13269 294

Hay una explicación más simple. Con su bucle, está sobrescribiendo un elemento de x en cada paso, reemplazando su valor anterior por uno de los otros elementos de x . Así que obtienes lo que pediste. Esencialmente, es una forma complicada de muestrear con reemplazo ( sample(x, replace=TRUE) ): si necesita tal complicación, depende de lo que quiera lograr.

Con su código vectorizado, solo está solicitando una cierta permutación de x (sin reemplazo), y eso es lo que obtiene. El código vectorizado no está haciendo lo mismo que su bucle. Si desea obtener el mismo resultado con un bucle, primero deberá hacer una copia de x :

set.seed(0) x <- x2 <- round(runif(10), 2) # [1] 0.90 0.27 0.37 0.57 0.91 0.20 0.90 0.94 0.66 0.63 sig <- sample.int(10) # [1] 1 2 9 5 3 4 8 6 7 10 for (i in seq_along(sig)) x2[i] <- x[sig[i]] identical(x2, x[sig]) #TRUE

No hay peligro de alias aquí: x y x2 refieren inicialmente a la misma ubicación de memoria, pero cambiará tan pronto como cambie el primer elemento de x2 .