¿La familia "*apply" realmente no está vectorizada?

Por eso estamos acostumbrados a decirle a cada nuevo usuario de R que " applyno está vectorizado, consulte el Patrick Burns R Inferno Circle 4 " que dice (cito):

Un reflejo común es utilizar una función de la familia apply. Esto no es vectorización, es ocultación de bucles . La función de aplicación tiene un bucle for en su definición. La función lapply oculta el bucle, pero los tiempos de ejecución tienden a ser aproximadamente iguales a los de un bucle for explícito.

De hecho, un vistazo rápido al applycódigo fuente revela el bucle:

grep("for", capture.output(getAnywhere("apply")), value = TRUE)
## [1] "        for (i in 1L:d2) {"  "    else for (i in 1L:d2) {"

Bien hasta ahora, pero un vistazo lapplyrevela vapplyuna imagen completamente diferente:

lapply
## function (X, FUN, ...) 
## {
##     FUN <- match.fun(FUN)
##     if (!is.vector(X) || is.object(X)) 
##        X <- as.list(X)
##     .Internal(lapply(X, FUN))
## }
## <bytecode: 0x000000000284b618>
## <environment: namespace:base>

Entonces, aparentemente no hay ningún forbucle R escondido allí, sino que están llamando a una función interna escrita en C.

Una mirada rápida a la madriguera del conejo revela prácticamente la misma imagen.

Además, tomemos colMeanspor ejemplo la función, a la que nunca se le acusó de no estar vectorizada.

colMeans
# function (x, na.rm = FALSE, dims = 1L) 
# {
#   if (is.data.frame(x)) 
#     x <- as.matrix(x)
#   if (!is.array(x) || length(dn <- dim(x)) < 2L) 
#     stop("'x' must be an array of at least two dimensions")
#   if (dims < 1L || dims > length(dn) - 1L) 
#     stop("invalid 'dims'")
#   n <- prod(dn[1L:dims])
#   dn <- dn[-(1L:dims)]
#   z <- if (is.complex(x)) 
#     .Internal(colMeans(Re(x), n, prod(dn), na.rm)) + (0+1i) * 
#     .Internal(colMeans(Im(x), n, prod(dn), na.rm))
#   else .Internal(colMeans(x, n, prod(dn), na.rm))
#   if (length(dn) > 1L) {
#     dim(z) <- dn
#     dimnames(z) <- dimnames(x)[-(1L:dims)]
#   }
#   else names(z) <- dimnames(x)[[dims + 1]]
#   z
# }
# <bytecode: 0x0000000008f89d20>
#   <environment: namespace:base>

¿Eh? También simplemente llama .Internal(colMeans(...y también lo podemos encontrar en la madriguera del conejo . Entonces, ¿en qué se diferencia esto de .Internal(lapply(..?

En realidad, una prueba comparativa rápida revela que sapplysu rendimiento no es peor colMeansy mucho mejor que un forbucle para un gran conjunto de datos.

m <- as.data.frame(matrix(1:1e7, ncol = 1e5))
system.time(colMeans(m))
# user  system elapsed 
# 1.69    0.03    1.73 
system.time(sapply(m, mean))
# user  system elapsed 
# 1.50    0.03    1.60 
system.time(apply(m, 2, mean))
# user  system elapsed 
# 3.84    0.03    3.90 
system.time(for(i in 1:ncol(m)) mean(m[, i]))
# user  system elapsed 
# 13.78    0.01   13.93 

En otras palabras, ¿es correcto decir que lapplyy vapply en realidad están vectorizados (en comparación con applycuál es un forbucle que también llama a lapply) y qué quiso decir realmente Patrick Burns?

Jan 01 '70 08:01 David Arenburg

Para mí, la vectorización se trata principalmente de hacer que el código sea más fácil de escribir y de entender.

El objetivo de una función vectorizada es eliminar la contabilidad asociada con un bucle for. Por ejemplo, en lugar de:

means <- numeric(length(mtcars))
for (i in seq_along(mtcars)) {
  means[i] <- mean(mtcars[[i]])
}
sds <- numeric(length(mtcars))
for (i in seq_along(mtcars)) {
  sds[i] <- sd(mtcars[[i]])
}

Puedes escribir:

means <- vapply(mtcars, mean, numeric(1))
sds   <- vapply(mtcars, sd, numeric(1))

Eso hace que sea más fácil ver qué es igual (los datos de entrada) y qué es diferente (la función que estás aplicando).

Una ventaja secundaria de la vectorización es que el bucle for a menudo se escribe en C, en lugar de R. Esto tiene importantes beneficios de rendimiento, pero no creo que sea la propiedad clave de la vectorización. La vectorización se trata fundamentalmente de salvar tu cerebro, no de salvar el trabajo de la computadora.

Mar 11 '2015 12:03 hadley

Estoy de acuerdo con la opinión de Patrick Burns de que se trata más bien de ocultación de bucles y no de vectorización de código . Este es el por qué:

Considere este Cfragmento de código:

for (int i=0; i<n; i++)
  c[i] = a[i] + b[i]

Lo que nos gustaría hacer está bastante claro. Pero en realidad no lo es cómo se realiza la tarea o cómo se podría realizar. Un bucle for por defecto es una construcción en serie. No informa si se pueden hacer cosas en paralelo ni cómo.

La forma más obvia es que el código se ejecuta de forma secuencial . Cargue a[i]y b[i]continúe con los registros, agréguelos, almacene el resultado en c[i]y haga esto para cada uno i.

Sin embargo, los procesadores modernos tienen un conjunto de instrucciones vectoriales o SIMD que son capaces de operar en un vector de datos durante la misma instrucción cuando se realiza la misma operación (por ejemplo, sumar dos vectores como se muestra arriba). Dependiendo del procesador/arquitectura, podría ser posible sumar, digamos, cuatro números desde ay bbajo la misma instrucción, en lugar de uno a la vez.

Nos gustaría explotar los datos múltiples de instrucción única y realizar paralelismo a nivel de datos , es decir, cargar 4 cosas a la vez, agregar 4 cosas a la vez, almacenar 4 cosas a la vez, por ejemplo. Y esta es la vectorización de código .

Tenga en cuenta que esto es diferente de la paralelización de código, donde se realizan múltiples cálculos simultáneamente.

Sería fantástico si el compilador identificara dichos bloques de código y los vectorizara automáticamente , lo cual es una tarea difícil. La vectorización automática de códigos es un tema de investigación desafiante en Ciencias de la Computación. Pero con el tiempo, los compiladores han mejorado. Puede consultar las capacidades de vectorización automática GNU-gcc aquí . Lo mismo ocurre LLVM-clang aquí . Y también puede encontrar algunos puntos de referencia en el último enlace comparados con gccy ICC(compilador Intel C++).

gcc(Estoy en v4.9), por ejemplo, no vectoriza el código automáticamente en el -O2nivel de optimización. Entonces, si tuviéramos que ejecutar el código que se muestra arriba, se ejecutaría secuencialmente. Este es el momento para sumar dos vectores enteros de longitud 500 millones.

Necesitamos agregar la bandera -ftree-vectorizeo cambiar la optimización a nivel -O3. (Tenga en cuenta que también -O3realiza otras optimizaciones adicionales ). La bandera -fopt-info-veces útil ya que informa cuando un bucle se vectorizó con éxito).

# compiling with -O2, -ftree-vectorize and  -fopt-info-vec
# test.c:32:5: note: loop vectorized
# test.c:32:5: note: loop versioned for vectorization because of possible aliasing
# test.c:32:5: note: loop peeled for vectorization to enhance alignment    

Esto nos dice que la función está vectorizada. Estos son los tiempos que comparan las versiones vectorizadas y no vectorizadas en vectores enteros de longitud 500 millones:

x = sample(100L, 500e6L, TRUE)
y = sample(100L, 500e6L, TRUE)
z = vector("integer", 500e6L) # result vector

# non-vectorised, -O2
system.time(.Call("Csum", x, y, z))
#    user  system elapsed 
#   1.830   0.009   1.852

# vectorised using flags shown above at -O2
system.time(.Call("Csum", x, y, z))
#    user  system elapsed 
#   0.361   0.001   0.362

# both results are checked for identicalness, returns TRUE

Esta parte se puede omitir de forma segura sin perder continuidad.

Los compiladores no siempre tendrán suficiente información para vectorizar. Podríamos usar la especificación OpenMP para programación paralela , que también proporciona una directiva de compilador simd para indicar a los compiladores que vectoricen el código. Es esencial asegurarse de que no haya superposiciones de memoria, condiciones de carrera, etc. al vectorizar el código manualmente; de ​​lo contrario, se obtendrán resultados incorrectos.

#pragma omp simd
for (i=0; i<n; i++) 
  c[i] = a[i] + b[i]

Al hacer esto, le pedimos específicamente al compilador que lo vectorice pase lo que pase. Necesitaremos activar las extensiones OpenMP usando el indicador de tiempo de compilación -fopenmp. Haciendo eso:

# timing with -O2 + OpenMP with simd
x = sample(100L, 500e6L, TRUE)
y = sample(100L, 500e6L, TRUE)
z = vector("integer", 500e6L) # result vector
system.time(.Call("Cvecsum", x, y, z))
#    user  system elapsed 
#   0.360   0.001   0.360

which is great! This was tested with gcc v6.2.0 and llvm clang v3.9.0 (both installed via homebrew, MacOS 10.12.3), both of which support OpenMP 4.0.

In this sense, even though Wikipedia page on Array Programming mentions that languages that operate on entire arrays usually call that as vectorised operations, it really is loop hiding IMO (unless it is actually vectorised).

In case of R, even rowSums() or colSums() code in C don't exploit code vectorisation IIUC; it is just a loop in C. Same goes for lapply(). In case of apply(), it's in R. All of these are therefore loop hiding.

In short, wrapping an R function by:

just writing a for-loop in C != vectorising your code.
just writing a for-loop in R != vectorising your code.

Intel Math Kernel Library (MKL) for example implements vectorised forms of functions.

HTH

References:

1. Talk by James Reinders, Intel (this answer is mostly an attempt to summarise this excellent talk)

Mar 13 '2015 00:03 Arun