¿La familia “* apply” no está realmente vectorizada?

Así que estamos acostumbrados a decir a cada nuevo usuario de R que ” apply no está vectorizado, echa un vistazo al Patrick Burns R Inferno Circle 4 ” que dice (cito):

Un reflection común es usar una función en la familia de aplicaciones. Esto no es vectorización, es un bucle escondido . La función aplicar tiene un bucle for en su definición. La función de aplicar entierra el ciclo, pero los tiempos de ejecución tienden a ser aproximadamente iguales a un ciclo for explícito.

De hecho, una mirada rápida al código fuente de la apply revela el ciclo:

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

Ok hasta ahora, pero una mirada a lapply o vapply realidad revela una 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)) ## } ##  ##  

Entonces, aparentemente, no hay R for bucle oculto allí, sino que están llamando función escrita interna C.

Una mirada rápida en el agujero del conejo revela más o menos la misma imagen

Además, tomemos la función colMeans por ejemplo, que nunca fue acusada de no ser 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  length(dn) - 1L) # stop("invalid 'dims'") # n <- prod(dn[1L:dims]) # dn <- dn[-(1L:dims)] # z  1L) { # dim(z) <- dn # dimnames(z) <- dimnames(x)[-(1L:dims)] # } # else names(z) <- dimnames(x)[[dims + 1]] # z # } #  #  

¿Huh? También llama a .Internal(colMeans(... que también podemos encontrar en el agujero del conejo . Entonces, ¿cómo es esto diferente de .Internal(lapply(.. ?

En realidad, un benchmark rápido revela que sapply no es peor que colMeans y mucho mejor que un bucle for para un big data set

 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 lapply y vapply están realmente vectorizados (en comparación con apply que es un ciclo for que también llama lapply ) y qué fue lo que realmente quiso decir Patrick Burns?

En primer lugar, en su ejemplo realiza pruebas en un “data.frame” que no es justo para colMeans , se apply y "[.data.frame" ya que tienen una sobrecarga:

 system.time(as.matrix(m)) #called by `colMeans` and `apply` # user system elapsed # 1.03 0.00 1.05 system.time(for(i in 1:ncol(m)) m[, i]) #in the `for` loop # user system elapsed # 12.93 0.01 13.07 

En una matriz, la imagen es un poco diferente:

 mm = as.matrix(m) system.time(colMeans(mm)) # user system elapsed # 0.01 0.00 0.01 system.time(apply(mm, 2, mean)) # user system elapsed # 1.48 0.03 1.53 system.time(for(i in 1:ncol(mm)) mean(mm[, i])) # user system elapsed # 1.22 0.00 1.21 

Regading la parte principal de la pregunta, la principal diferencia entre lapply / mapply / etc y bucles R directos es donde se realiza el bucle. Como señala Roland, los bucles C y R necesitan evaluar una función R en cada iteración, que es la más costosa. Las funciones realmente rápidas de C son aquellas que hacen todo en C, así que, supongo, ¿de qué se tratará “vectorizado”?

Un ejemplo donde encontramos la media en cada uno de los elementos de una “lista”:

( EDITAR 11 de mayo de16 : creo que el ejemplo para encontrar la “media” no es una buena configuración para las diferencias entre evaluar una función R de forma iterativa y código comstackdo, (1) debido a la particularidad del algoritmo de la media de R en “numérico” s sobre una sum(x) / length(x) simple sum(x) / length(x) y (2) debería tener más sentido probar en “listas” con length(x) >> lengths(x) . Entonces, se mueve el ejemplo “malo” hasta el final y reemplazado por otro.)

Como un simple ejemplo, podríamos considerar el hallazgo del opuesto de cada length == 1 elemento de una “lista”:

En un archivo tmp.c :

 #include  #define USE_RINTERNALS #include  #include  /* call a C function inside another */ double oppC(double x) { return(ISNAN(x) ? NA_REAL : -x); } SEXP sapply_oppC(SEXP x) { SEXP ans = PROTECT(allocVector(REALSXP, LENGTH(x))); for(int i = 0; i < LENGTH(x); i++) REAL(ans)[i] = oppC(REAL(VECTOR_ELT(x, i))[0]); UNPROTECT(1); return(ans); } /* call an R function inside a C function; * will be used with 'f' as a closure and as a builtin */ SEXP sapply_oppR(SEXP x, SEXP f) { SEXP call = PROTECT(allocVector(LANGSXP, 2)); SETCAR(call, install(CHAR(STRING_ELT(f, 0)))); SEXP ans = PROTECT(allocVector(REALSXP, LENGTH(x))); for(int i = 0; i < LENGTH(x); i++) { SETCADR(call, VECTOR_ELT(x, i)); REAL(ans)[i] = REAL(eval(call, R_GlobalEnv))[0]; } UNPROTECT(2); return(ans); } 

Y en el lado R:

 system("R CMD SHLIB /home/~/tmp.c") dyn.load("/home/~/tmp.so") 

con datos:

 set.seed(007) myls = rep_len(as.list(c(NA, runif(3))), 1e7) #a closure wrapper of `-` oppR = function(x) -x for_oppR = compiler::cmpfun(function(x, f) { f = match.fun(f) ans = numeric(length(x)) for(i in seq_along(x)) ans[[i]] = f(x[[i]]) return(ans) }) 

Benchmarking:

 #call a C function iteratively system.time({ sapplyC = .Call("sapply_oppC", myls) }) # user system elapsed # 0.048 0.000 0.047 #evaluate an R closure iteratively system.time({ sapplyRC = .Call("sapply_oppR", myls, "oppR") }) # user system elapsed # 3.348 0.000 3.358 #evaluate an R builtin iteratively system.time({ sapplyRCprim = .Call("sapply_oppR", myls, "-") }) # user system elapsed # 0.652 0.000 0.653 #loop with a R closure system.time({ forR = for_oppR(myls, "oppR") }) # user system elapsed # 4.396 0.000 4.409 #loop with an R builtin system.time({ forRprim = for_oppR(myls, "-") }) # user system elapsed # 1.908 0.000 1.913 #for reference and testing system.time({ sapplyR = unlist(lapply(myls, oppR)) }) # user system elapsed # 7.080 0.068 7.170 system.time({ sapplyRprim = unlist(lapply(myls, `-`)) }) # user system elapsed # 3.524 0.064 3.598 all.equal(sapplyR, sapplyRprim) #[1] TRUE all.equal(sapplyR, sapplyC) #[1] TRUE all.equal(sapplyR, sapplyRC) #[1] TRUE all.equal(sapplyR, sapplyRCprim) #[1] TRUE all.equal(sapplyR, forR) #[1] TRUE all.equal(sapplyR, forRprim) #[1] TRUE 

(Sigue el ejemplo original de hallazgo malo):

 #all computations in C all_C = inline::cfunction(sig = c(R_ls = "list"), body = ' SEXP tmp, ans; PROTECT(ans = allocVector(REALSXP, LENGTH(R_ls))); double *ptmp, *pans = REAL(ans); for(int i = 0; i < LENGTH(R_ls); i++) { pans[i] = 0.0; PROTECT(tmp = coerceVector(VECTOR_ELT(R_ls, i), REALSXP)); ptmp = REAL(tmp); for(int j = 0; j < LENGTH(tmp); j++) pans[i] += ptmp[j]; pans[i] /= LENGTH(tmp); UNPROTECT(1); } UNPROTECT(1); return(ans); ') #a very simple `lapply(x, mean)` C_and_R = inline::cfunction(sig = c(R_ls = "list"), body = ' SEXP call, ans, ret; PROTECT(call = allocList(2)); SET_TYPEOF(call, LANGSXP); SETCAR(call, install("mean")); PROTECT(ans = allocVector(VECSXP, LENGTH(R_ls))); PROTECT(ret = allocVector(REALSXP, LENGTH(ans))); for(int i = 0; i < LENGTH(R_ls); i++) { SETCADR(call, VECTOR_ELT(R_ls, i)); SET_VECTOR_ELT(ans, i, eval(call, R_GlobalEnv)); } double *pret = REAL(ret); for(int i = 0; i < LENGTH(ans); i++) pret[i] = REAL(VECTOR_ELT(ans, i))[0]; UNPROTECT(3); return(ret); ') R_lapply = function(x) unlist(lapply(x, mean)) R_loop = function(x) { ans = numeric(length(x)) for(i in seq_along(x)) ans[i] = mean(x[[i]]) return(ans) } R_loopcmp = compiler::cmpfun(R_loop) set.seed(007); myls = replicate(1e4, runif(1e3), simplify = FALSE) all.equal(all_C(myls), C_and_R(myls)) #[1] TRUE all.equal(all_C(myls), R_lapply(myls)) #[1] TRUE all.equal(all_C(myls), R_loop(myls)) #[1] TRUE all.equal(all_C(myls), R_loopcmp(myls)) #[1] TRUE microbenchmark::microbenchmark(all_C(myls), C_and_R(myls), R_lapply(myls), R_loop(myls), R_loopcmp(myls), times = 15) #Unit: milliseconds # expr min lq median uq max neval # all_C(myls) 37.29183 38.19107 38.69359 39.58083 41.3861 15 # C_and_R(myls) 117.21457 123.22044 124.58148 130.85513 169.6822 15 # R_lapply(myls) 98.48009 103.80717 106.55519 109.54890 116.3150 15 # R_loop(myls) 122.40367 130.85061 132.61378 138.53664 178.5128 15 # R_loopcmp(myls) 105.63228 111.38340 112.16781 115.68909 128.1976 15 

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

El objective 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 lo mismo (los datos de entrada) y qué es diferente (la función que está aplicando).

Una ventaja secundaria de la vectorización es que el for-loop a menudo se escribe en C, en lugar de en 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.

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

Considere este fragmento de código C :

 for (int i=0; i 

Lo que nos gustaría hacer es bastante claro. Pero la forma en que se realiza la tarea o cómo se puede realizar no es realmente así. Un for-loop por defecto es una construcción en serie. No informa si o cómo se pueden hacer las cosas en paralelo.

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

Sin embargo, los procesadores modernos tienen un conjunto de instrucciones vectoriales o SIMD que es capaz de operar en un vector de datos durante la misma instrucción cuando se realiza la misma operación (por ejemplo, agregar dos vectores como se muestra arriba). Dependiendo del procesador / architecture, podría ser posible agregar, por ejemplo, cuatro números de b bajo 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 el paralelismo de nivel de datos , es decir, cargar 4 elementos a la vez, agregar 4 elementos al tiempo, almacenar 4 elementos a la vez, por ejemplo. Y esto es vectorización de código .

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

Sería genial si el comstackdor identifica tales bloques de código y los vectoriza 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 informática. Pero con el tiempo, los comstackdores han mejorado. Puede verificar las capacidades de vectorización automática de GNU-gcc aquí . Del mismo modo para LLVM-clang aquí . Y también puede encontrar algunos puntos de referencia en el último enlace en comparación con gcc e ICC (comstackdor Intel C ++).

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

O bien necesitamos agregar el indicador -ftree-vectorize o cambiar la optimización al nivel -O3 . (Tenga en cuenta que -O3 realiza otras optimizaciones adicionales también). La bandera -fopt-info-vec es ú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. Aquí están los tiempos que comparan versiones no vectorizadas y 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 sin perder continuidad.

Los comstackdores no siempre tendrán información suficiente para vectorizar. Podríamos usar la especificación OpenMP para progtwigción paralela , que también proporciona una directiva de comstackdor simd para instruir a los comstackdores a vectorizar 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 

Al hacer esto, le pedimos específicamente al comstackdor que lo vectorice pase lo que pase. Tendremos que activar las extensiones de OpenMP utilizando el indicador de tiempo de comstackció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 

¡Lo cual es genial! Esto se probó con gcc v6.2.0 y llvm clang v3.9.0 (ambos instalados a través de homebrew, MacOS 10.12.3), que admiten OpenMP 4.0.


En este sentido, aunque la página de Wikipedia en Array Programming menciona que los lenguajes que operan en arreglos completos usualmente lo llaman operaciones vectorizadas , realmente es un bucle que oculta IMO (a menos que realmente se vectorice).

En el caso de R, incluso los rowSums() o colSums() en C no explotan la vectorización de código IIUC; es solo un ciclo en C. Lo mismo ocurre con lapply() . En el caso de apply() , está en R. Por lo tanto, todos estos son bucles ocultos .

En resumen, envolviendo una función R por:

simplemente escribiendo un for-loop en C ! = vectorizando tu código.
simplemente escribiendo un for-loop en R ! = vectorizando tu código.

Intel Math Kernel Library (MKL), por ejemplo, implementa formas vectorizadas de funciones.

HTH


Referencias

  1. Charla de James Reinders, Intel (esta respuesta es principalmente un bash de resumir esta excelente charla)

Por lo tanto, para sumr las grandes respuestas / comentarios en una respuesta general y proporcionar algunos antecedentes: R tiene 4 tipos de bucles ( desde vector no vectorizado a vectorizado )

  1. R for bucle que llama repetidamente funciones R en cada iteración ( No vectorizado )
  2. C loop que llama repetidamente funciones R en cada iteración ( No vectorizado )
  3. C loop que llama a la función R solo una vez ( algo vectorizado )
  4. Un bucle C simple que no llama ninguna función R en absoluto y usa sus propias funciones comstackdas ( Vectorizado )

Entonces, la familia *apply es el segundo tipo. Excepto apply que es más del primer tipo

Puedes entender esto desde el comentario en su código fuente

/ * .Internal (aplicar (X, DIVERSIÓN)) * /

/ * Esto es un especial .Interno, por lo que tiene argumentos no evaluados. Es
llamado desde un contenedor de cierre, por lo que X y FUN son promesas. FUN debe estar sin evaluar para su uso en, por ejemplo, bquote. * /

Esto significa que el código C de la aplicación acepta una función no evaluada de R y luego la evalúa dentro del código C mismo. Esta es básicamente la diferencia entre lapply s .Internal Llamada .Internal

 .Internal(lapply(X, FUN)) 

Que tiene un argumento FUN que tiene una función R

Y los colMeans .Internal Llamada .Internal que no tiene un argumento FUN

 .Internal(colMeans(Re(x), n, prod(dn), na.rm)) 

colMeans , a diferencia de lapply sabe exactamente qué función necesita usar, por lo tanto, calcula la media internamente dentro del código C.

Puede ver claramente el proceso de evaluación de la función R en cada iteración dentro del código C de la aplicación

  for(R_xlen_t i = 0; i < n; i++) { if (realIndx) REAL(ind)[0] = (double)(i + 1); else INTEGER(ind)[0] = (int)(i + 1); tmp = eval(R_fcall, rho); // <----------------------------- here it is if (MAYBE_REFERENCED(tmp)) tmp = lazy_duplicate(tmp); SET_VECTOR_ELT(ans, i, tmp); } 

En resumen, lapply no está vectorizado , aunque tiene dos ventajas posibles sobre el R simple for ciclo

  1. El acceso y la asignación en un bucle parece ser más rápido en C (es decir, al lapply una función). Aunque la diferencia parece grande, aún nos mantenemos en el nivel de microsegundo y lo costoso es la valoración de una función R en cada iteración. Un simple ejemplo:

     ffR = function(x) { ans = vector("list", length(x)) for(i in seq_along(x)) ans[[i]] = x[[i]] ans } ffC = inline::cfunction(sig = c(R_x = "data.frame"), body = ' SEXP ans; PROTECT(ans = allocVector(VECSXP, LENGTH(R_x))); for(int i = 0; i < LENGTH(R_x); i++) SET_VECTOR_ELT(ans, i, VECTOR_ELT(R_x, i)); UNPROTECT(1); return(ans); ') set.seed(007) myls = replicate(1e3, runif(1e3), simplify = FALSE) mydf = as.data.frame(myls) all.equal(ffR(myls), ffC(myls)) #[1] TRUE all.equal(ffR(mydf), ffC(mydf)) #[1] TRUE microbenchmark::microbenchmark(ffR(myls), ffC(myls), ffR(mydf), ffC(mydf), times = 30) #Unit: microseconds # expr min lq median uq max neval # ffR(myls) 3933.764 3975.076 4073.540 5121.045 32956.580 30 # ffC(myls) 12.553 12.934 16.695 18.210 19.481 30 # ffR(mydf) 14799.340 15095.677 15661.889 16129.689 18439.908 30 # ffC(mydf) 12.599 13.068 15.835 18.402 20.509 30 
  2. Como menciona @Roland, ejecuta un bucle C comstackdo en lugar de un bucle R interpretado


Aunque al vectorizar tu código, hay algunas cosas que debes tener en cuenta.

  1. Si su conjunto de datos (llamémoslo df ) es de clase data.frame , algunas funciones vectorizadas (como colMeans , colSums , rowSums , etc.) deberán convertirlo primero a una matriz, simplemente porque así es como fueron diseñadas . Esto significa que para un gran df esto puede generar una gran sobrecarga. Si bien lapply no tendrá que hacer esto, ya que extrae los vectores reales de df (como data.frame es solo una lista de vectores) y, por lo tanto, si no tiene tantas columnas sino muchas filas, lapply(df, mean) a veces puede ser una mejor opción que colMeans(df) .
  2. Otra cosa para recordar es que R tiene una gran variedad de tipos de funciones diferentes, como .Primitive y genérico ( S3 , S4 ), consulte aquí para obtener información adicional. La función genérica tiene que hacer un despacho de método que a veces es una operación costosa. Por ejemplo, mean es la función genérica S3 mientras que sum es Primitive . Por lo tanto, algunas veces lapply(df, sum) podría ser muy eficiente comparado con colSums por los motivos enumerados anteriormente