pro - graficos en r

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

Por lo tanto, estamos acostumbrados a decirle a cada nuevo usuario de R que "la apply no está vectorizada, vea el Círculo Infernal 4 de Patrick Burns R " que dice (cito):

Un reflejo común es usar una función en la familia de postulantes. Esto no es vectorización, está ocultando bucles . La función de aplicación tiene un bucle for en su definición. La función lapply entierra el ciclo, pero los tiempos de ejecución tienden a ser aproximadamente iguales a un ciclo for explícito.

De hecho, un vistazo rápido al código fuente de 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)) ## } ## <bytecode: 0x000000000284b618> ## <environment: namespace:base>

Entonces, aparentemente no hay una R for bucle que se esconde allí, sino que están llamando a la función escrita interna de C.

Una mirada rápida en la rabbit del rabbit revela casi 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 < 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 solo llama .Internal(colMeans(... que también podemos encontrar en la madriguera del conejo . Entonces, ¿cómo es esto diferente de .Internal(lapply(.. ?

En realidad, un punto de referencia rápido revela que sapply funciona no peor que colMeans y mucho mejor que un bucle for 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 lapply y vapply están realmente vectorizados (en comparación con apply que es un bucle for que también llama lapply ) y qué quería decir realmente Patrick Burns?

En primer lugar, en su ejemplo, realiza pruebas en un "data.frame" que no es justo para colMeans , 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

Regulando la parte principal de la pregunta, la principal diferencia entre lapply / mapply / etc y R-loops directos es donde se realiza el bucle. Como señala Roland, los bucles C y R deben evaluar una función R en cada iteración, que es la más costosa. Las funciones C realmente rápidas son aquellas que hacen todo en C, así que, supongo, ¿esto debería ser de lo que se trata la "vectorización"?

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

( EDITAR 11 de mayo de 16 : creo que el ejemplo de encontrar la "media" no es una buena configuración para las diferencias entre evaluar una función R de forma iterativa y un código compilado, (1) debido a la particularidad del algoritmo medio 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 "list" s con length(x) >> lengths(x) . Entonces, el ejemplo "mean" se mueve hasta el final y reemplazado por otro.)

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

En un archivo tmp.c :

#include <R.h> #define USE_RINTERNALS #include <Rinternals.h> #include <Rdefines.h> /* 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 medio):

#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

Entonces, para resumir las excelentes respuestas / comentarios en una respuesta general y proporcionar algunos antecedentes: R tiene 4 tipos de bucles ( desde el orden no vectorizado hasta el vectorizado )

R for loop que llama repetidamente a las funciones R en cada iteración ( no vectorizado )
C loop que llama repetidamente a las funciones R en cada iteración ( no vectorizado )
C loop que llama a la función R solo una vez ( algo vectorizado )
Un bucle C simple que no llama a ninguna función R y usa sus propias funciones compiladas ( Vectorizado )

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

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

/ * .Internal (lapply (X, FUN)) * /

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

Eso significa que el código C de lapply s acepta una función no evaluada de R y luego la evalúa dentro del propio código C. 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 el colMeans .Internal Llamada .Internal que no tiene un argumento .Internal

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

colMeans , a diferencia de lapply sabe exactamente qué función necesita usar, por lo que 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 lapply

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); }

Para resumir, lapply no está vectorizado , aunque tiene dos posibles ventajas sobre el bucle R simple.

El acceso y la asignación en un bucle parece ser más rápido en C (es decir, en la lapply de una función) Aunque la diferencia parece grande, nosotros, aún, permanecemos en el nivel de microsegundos 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
Como mencionó @Roland, ejecuta un bucle C compilado en lugar de un bucle R interpretado

Aunque al vectorizar su código, hay algunas cosas que debe tener en cuenta.

Si su conjunto de datos (llamémoslo df ) es de clase data.frame , algunas funciones vectorizadas (como colMeans , colSums , rowSums , etc.) rowSums tendrán que convertirlo en una matriz, simplemente porque así es como fueron diseñados . Esto significa que para un gran df esto puede crear una gran sobrecarga. Si bien lapply no tendrá que hacer esto, ya que extrae los vectores reales de df (ya que 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) .
Otra cosa para recordar es que R tiene una gran variedad de diferentes tipos de funciones, como .Primitive y genérico ( S3 , S4 ). Consulte here para obtener información adicional. La función genérica tiene que hacer un envío de método que a veces es una operación costosa. Por ejemplo, la mean es la función genérica S3 mientras que la sum es Primitive . Por lo tanto, algunas veces lapply(df, sum) podría ser muy eficiente en comparación con colSums por las razones enumeradas anteriormente

Estoy de acuerdo con la opinión de Patrick Burns de que es más bien un bucle oculto 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<n; i++) c[i] = a[i] + b[i]

Lo que nos gustaría hacer es bastante claro. Pero cómo se realiza la tarea o cómo se podría realizar no es realmente. Un bucle for 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 manera secuencial . Cargue a[i] 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 al realizar la misma operación (por ejemplo, agregando dos vectores como se muestra arriba). Dependiendo del procesador / arquitectura, podría ser posible agregar, por ejemplo, cuatro números de a y 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 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 genial si el compilador identifica dichos bloques de código y los vectoriza automáticamente , lo cual es una tarea difícil. La vectorización automática de código es un tema de investigación desafiante en informática. Pero con el tiempo, los compiladores han mejorado en eso. Puede verificar las capacidades de vectorización automática de GNU-gcc here . Del mismo modo para LLVM-clang here . Y también puede encontrar algunos puntos de referencia en el último enlace en comparación con gcc e ICC (compilador Intel C ++).

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

Necesitamos agregar la bandera -ftree-vectorize o cambiar la optimización al nivel -O3 . (Tenga en cuenta que -O3 realiza otras optimizaciones adicionales ). La bandera -fopt-info-vec es útil, ya que informa cuando un ciclo 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 no vectorizadas y vectorizadas en vectores enteros de 500 millones de longitud:

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 información suficiente para vectorizar. Podríamos usar la especificación OpenMP para programación paralela , que también proporciona una directiva de compilación simd para instruir a los compiladores 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<n; i++) c[i] = a[i] + b[i]

Al hacer esto, le pedimos específicamente al compilador que lo vectorice sin importar qué. Tendremos que activar las extensiones de OpenMP utilizando 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

¡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), los cuales son compatibles con OpenMP 4.0.

En este sentido, a pesar de que la página de Wikipedia sobre la programación de matrices menciona que los lenguajes que operan en matrices enteras generalmente lo llaman operaciones vectorizadas , en realidad se trata de ocultar IMO (a menos que en realidad esté vectorizado).

En el caso de R, incluso el rowSums() o colSums() en C no explota la vectorización de código IIUC; es solo un bucle en C. Lo mismo ocurre con lapply() . En el caso de apply() , está en R. Todos estos son, por lo tanto, ocultación de bucles .

En resumen, envolviendo una función R:

simplemente escribiendo un bucle for en C ! = vectorizando su código.
simplemente escribiendo un bucle for en R ! = vectorizando su código.

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

HTH

Referencias

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

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 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 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 bucle for 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 su cerebro, no de salvar el trabajo de la computadora.