Acelera la operación de bucle en R

Tengo un gran problema de rendimiento en R. Escribí una función que itera sobre un objeto data.frame . Simplemente agrega una nueva columna a un data.frame y acumula algo. (operación simple). El data.frame tiene aproximadamente 850K filas. Mi PC todavía funciona (aproximadamente 10 horas ahora) y no tengo idea del tiempo de ejecución.

 dayloop2 <- function(temp){ for (i in 1:nrow(temp)){ temp[i,10]  1) { if ((temp[i,6] == temp[i-1,6]) & (temp[i,3] == temp[i-1,3])) { temp[i,10] <- temp[i,9] + temp[i-1,10] } else { temp[i,10] <- temp[i,9] } } else { temp[i,10] <- temp[i,9] } } names(temp)[names(temp) == "V10"] <- "Kumm." return(temp) } 

¿Alguna idea de cómo acelerar esta operación?

El mayor problema y la raíz de la ineficacia es indexar data.frame, me refiero a todas estas líneas donde se usa temp[,] .
Intenta evitar esto tanto como sea posible. Tomé su función, cambio de indexación y aquí version_A

 dayloop2_A < - function(temp){ res <- numeric(nrow(temp)) for (i in 1:nrow(temp)){ res[i] <- i if (i > 1) { if ((temp[i,6] == temp[i-1,6]) & (temp[i,3] == temp[i-1,3])) { res[i] < - temp[i,9] + res[i-1] } else { res[i] <- temp[i,9] } } else { res[i] <- temp[i,9] } } temp$`Kumm.` <- res return(temp) } 

Como puedes ver, creo vectores res que recogen resultados. Al final lo agrego a data.frame y no necesito data.frame con los nombres. Entonces, ¿qué tan mejor es?

data.frame cada función para data.frame con nrow de 1,000 a 10,000 por 1,000 y system.time tiempo con system.time

 X < - as.data.frame(matrix(sample(1:10, n*9, TRUE), n, 9)) system.time(dayloop2(X)) 

El resultado es

actuación

Puedes ver que tu versión depende exponencialmente de nrow(X) . La versión modificada tiene relación lineal, y el modelo simple de lm predice que para 850,000 filas, el cálculo lleva 6 minutos y 10 segundos.

Poder de la vectorización

Como Shane y Calimo afirman en sus respuestas, la vectorización es la clave para un mejor rendimiento. Desde tu código, puedes moverte fuera del bucle:

  • acondicionamiento
  • inicialización de los resultados (que son temp[i,9] )

Esto lleva a este código

 dayloop2_B < - function(temp){ cond <- c(FALSE, (temp[-nrow(temp),6] == temp[-1,6]) & (temp[-nrow(temp),3] == temp[-1,3])) res <- temp[,9] for (i in 1:nrow(temp)) { if (cond[i]) res[i] <- temp[i,9] + res[i-1] } temp$`Kumm.` <- res return(temp) } 

Compare el resultado para estas funciones, esta vez para nrow de 10,000 a 100,000 por 10,000.

actuación

Sintonizando el afinado

Otro cambio es cambiar en una temp[i,9] indexación de bucle temp[i,9] a res[i] (que son exactamente iguales en la iteración de bucle i-th). data.frame vez es la diferencia entre indexar un vector e indexar un data.frame .
Lo segundo: cuando miras el ciclo puedes ver que no hay necesidad de recorrer todo i , sino solo para aquellos que se ajustan a la condición.
Así que, aquí vamos

 dayloop2_D < - function(temp){ cond <- c(FALSE, (temp[-nrow(temp),6] == temp[-1,6]) & (temp[-nrow(temp),3] == temp[-1,3])) res <- temp[,9] for (i in (1:nrow(temp))[cond]) { res[i] <- res[i] + res[i-1] } temp$`Kumm.` <- res return(temp) } 

El rendimiento que obtiene depende en gran medida de una estructura de datos. Precisamente, en porcentaje de valores TRUE en la condición. Para mis datos simulados, se necesita un tiempo de cálculo de 850,000 filas por debajo del segundo.

actuación

Si quieres puedes ir más lejos, veo al menos dos cosas que se pueden hacer:

  • escribe un código C para hacer cumsum condicional
  • si sabes que en tus datos la secuencia máxima no es grande, entonces puedes cambiar el ciclo a vectorizado mientras, algo así como

     while (any(cond)) { indx < - c(FALSE, cond[-1] & !cond[-n]) res[indx] <- res[indx] + res[which(indx)-1] cond[indx] <- FALSE } 

El código utilizado para simulaciones y figuras está disponible en GitHub .

Estrategias generales para acelerar el código R

Primero, averigua dónde está realmente la parte lenta. No es necesario optimizar el código que no se ejecuta lentamente. Para pequeñas cantidades de código, simplemente pensarlo puede funcionar. Si eso falla, RProf y herramientas de perfilado similares pueden ser útiles.

Una vez que descubra el cuello de botella, piense en algoritmos más eficientes para hacer lo que quiere. Los cálculos solo deben ejecutarse una vez, de ser posible, así que:

  • Almacene los resultados y acceda a ellos en lugar de recalcular repetidamente
  • Realice cálculos que no dependen del ciclo de bucles
  • Evite los cálculos que no son necesarios (p. Ej. , No usar expresiones regulares con búsquedas fijas )

El uso de funciones más eficientes puede producir ganancias de velocidad moderadas o grandes. Por ejemplo, paste0 produce una pequeña ganancia de eficiencia pero .colSums() y sus parientes producen ganancias algo más pronunciadas. mean es particularmente lenta .

Entonces puedes evitar algunos problemas particularmente comunes :

  • cbind lo retrasará muy rápido.
  • Inicialice sus estructuras de datos, luego instálelas, en lugar de expandirlas cada vez .
  • Incluso con la preasignación, puede cambiar a un enfoque de paso por referencia en lugar de un enfoque de paso por valor, pero puede no valer la pena la molestia.
  • Eche un vistazo al R Inferno para evitar más trampas.

Intente una mejor vectorización , que a menudo puede ayudar pero no siempre. A este respecto, los comandos inherentemente vectorizados como ifelse , diff y similares proporcionarán más mejoras que la familia de comandos apply (que proporcionan poco o ningún aumento de velocidad sobre un bucle bien escrito).

También puede intentar proporcionar más información a las funciones R. Por ejemplo, use vapply lugar de sapply , y especifique colClasses al leer en datos basados ​​en texto . Las ganancias de velocidad serán variables según la cantidad de conjeturas que elimines.

A continuación, considere los paquetes optimizados : el paquete data.table puede producir ganancias de velocidad masivas donde su uso es posible, en la manipulación de datos y en la lectura de grandes cantidades de datos ( fread ).

A continuación, intente obtener ganancias de velocidad a través de medios más eficientes para llamar a R :

  • Comstack tu guión R O utilice los paquetes Ra y jit en concierto para la comstackción justo a tiempo (Dirk tiene un ejemplo en esta presentación ).
  • Asegúrese de estar utilizando un BLAS optimizado. Estos proporcionan ganancias de velocidad generales. Honestamente, es una pena que R no use automáticamente la biblioteca más eficiente en la instalación. Es de esperar que Revolution R contribuya con el trabajo que han realizado aquí a la comunidad en general.
  • Radford Neal ha hecho un montón de optimizaciones, algunas de las cuales fueron adoptadas en R Core, y muchas otras que fueron incluidas en pqR .

Y, por último, si todo lo anterior aún no lo hace tan rápido como lo necesita, puede necesitar moverse a un lenguaje más rápido para el fragmento de código lento . La combinación de Rcpp y en inline aquí hace que reemplazar solo la parte más lenta del algoritmo con código C ++ sea particularmente fácil. Aquí, por ejemplo, es mi primer bash de hacerlo , y se lleva incluso las soluciones R altamente optimizadas.

Si aún te quedan problemas después de todo esto, solo necesitas más poder de cómputo. Consulte la paralelización ( http://cran.r-project.org/web/views/HighPerformanceComputing.html ) o incluso las soluciones basadas en GPU ( gpu-tools ).

Enlaces a otra guía

Si está utilizando bucles for , lo más probable es que esté codificando R como si fuera C o Java o algo más. El código R que está correctamente vectorizado es extremadamente rápido.

Tomemos por ejemplo estos dos bits simples de código para generar una lista de 10.000 enteros en secuencia:

El primer ejemplo de código es cómo se codifica un bucle utilizando un paradigma de encoding tradicional. Tarda 28 segundos en completarse

 system.time({ a < - NULL for(i in 1:1e5)a[i] <- i }) user system elapsed 28.36 0.07 28.61 

Puede obtener una mejora de casi 100 veces mediante la simple acción de preasignar la memoria:

 system.time({ a < - rep(1, 1e5) for(i in 1:1e5)a[i] <- i }) user system elapsed 0.30 0.00 0.29 

Pero utilizando la operación de vector base R utilizando el operador de dos puntos : esta operación es prácticamente instantánea:

 system.time(a < - 1:1e5) user system elapsed 0 0 0 

Esto podría hacerse mucho más rápido omitiendo los bucles usando índices o declaraciones ifelse() anidadas.

 idx < - 1:nrow(temp) temp[,10] <- idx idx1 <- c(FALSE, (temp[-nrow(temp),6] == temp[-1,6]) & (temp[-nrow(temp),3] == temp[-1,3])) temp[idx1,10] <- temp[idx1,9] + temp[which(idx1)-1,10] temp[!idx1,10] <- temp[!idx1,9] temp[1,10] <- temp[1,9] names(temp)[names(temp) == "V10"] <- "Kumm." 

Como Ari mencionó al final de su respuesta, los paquetes Rcpp y en inline hacen que sea increíblemente fácil hacer las cosas rápido. Como ejemplo, pruebe este código en inline (advertencia: no probado):

 body < - 'Rcpp::NumericMatrix nm(temp); int nrtemp = Rccp::as(nrt); for (int i = 0; i < nrtemp; ++i) { temp(i, 9) = i if (i > 1) { if ((temp(i, 5) == temp(i - 1, 5) && temp(i, 2) == temp(i - 1, 2) { temp(i, 9) = temp(i, 8) + temp(i - 1, 9) } else { temp(i, 9) = temp(i, 8) } } else { temp(i, 9) = temp(i, 8) } return Rcpp::wrap(nm); ' settings < - getPlugin("Rcpp") # settings$env$PKG_CXXFLAGS <- paste("-I", getwd(), sep="") if you want to inc files in wd dayloop <- cxxfunction(signature(nrt="numeric", temp="numeric"), body-body, plugin="Rcpp", settings=settings, cppargs="-I/usr/include") dayloop2 <- function(temp) { # extract a numeric matrix from temp, put it in tmp nc <- ncol(temp) nm <- dayloop(nc, temp) names(temp)[names(temp) == "V10"] <- "Kumm." return(temp) } 

Hay un procedimiento similar para #include cosas, donde solo pasas un parámetro

 inc < - '#include  

a cxxfunction, como include=inc . Lo bueno de esto es que hace todo el enlace y la comstackción para usted, por lo que el prototipado es realmente rápido.

Descargo de responsabilidad: no estoy totalmente seguro de que la clase de tmp debe ser numérica y no matriz numérica o algo más. Pero estoy casi seguro.

Editar: si aún necesita más velocidad después de esto, OpenMP es una instalación de paralelización adecuada para C++ . No he intentado usarlo en inline , pero debería funcionar. La idea sería, en el caso de n núcleos, que la iteración de ciclo k se lleve a cabo por k % n . Una introducción adecuada se encuentra en The Art of R Programming de Matloff , disponible aquí , en el capítulo 16, Recurrir a C.

No me gusta volver a escribir el código … También, por supuesto, ifelse y lapply son mejores opciones, pero a veces es difícil hacerlo.

Frecuentemente uso data.frames como uno usaría listas como df$var[i]

Aquí hay un ejemplo inventado:

 nrow=function(x){ ##required as I use nrow at times. if(class(x)=='list') { length(x[[names(x)[1]]]) }else{ base::nrow(x) } } system.time({ d=data.frame(seq=1:10000,r=rnorm(10000)) d$foo=d$r d$seq=1:5 mark=NA for(i in 1:nrow(d)){ if(d$seq[i]==1) mark=d$r[i] d$foo[i]=mark } }) system.time({ d=data.frame(seq=1:10000,r=rnorm(10000)) d$foo=d$r d$seq=1:5 d=as.list(d) #become a list mark=NA for(i in 1:nrow(d)){ if(d$seq[i]==1) mark=d$r[i] d$foo[i]=mark } d=as.data.frame(d) #revert back to data.frame }) 

versión de data.frame:

  user system elapsed 0.53 0.00 0.53 

lista de la versión:

  user system elapsed 0.04 0.00 0.03 

17 veces más rápido para usar una lista de vectores que un data.frame.

¿Algún comentario sobre por qué internamente los data.frames son tan lentos a este respecto? Uno pensaría que operan como listas …

Para un código aún más rápido, haga esta class(d)='list' lugar de d=as.list(d) y class(d)='data.frame'

 system.time({ d=data.frame(seq=1:10000,r=rnorm(10000)) d$foo=d$r d$seq=1:5 class(d)='list' mark=NA for(i in 1:nrow(d)){ if(d$seq[i]==1) mark=d$r[i] d$foo[i]=mark } class(d)='data.frame' }) head(d) 

Las respuestas aquí son geniales. Un aspecto menor no cubierto es que la pregunta dice ” Mi PC todavía está funcionando (alrededor de 10 horas ahora) y no tengo idea del tiempo de ejecución “. Siempre introduzco el siguiente código en bucles al desarrollar para tener una idea de cómo los cambios parecen afectar la velocidad y también para monitorear cuánto tiempo tomará completar.

 dayloop2 < - function(temp){ for (i in 1:nrow(temp)){ cat(round(i/nrow(temp)*100,2),"% \r") # prints the percentage complete in realtime. # do stuff } return(blah) } 

Funciona con lapply también.

 dayloop2 < - function(temp){ temp <- lapply(1:nrow(temp), function(i) { cat(round(i/nrow(temp)*100,2),"% \r") #do stuff }) return(temp) } 

Si la función dentro del ciclo es bastante rápida pero la cantidad de ciclos es grande, entonces considere imprimir solo cada cierto tiempo, ya que la impresión en la consola tiene un costo adicional. p.ej

 dayloop2 < - function(temp){ for (i in 1:nrow(temp)){ if(i %% 100 == 0) cat(round(i/nrow(temp)*100,2),"% \r") # prints every 100 times through the loop # do stuff } return(temp) } 

En R, a menudo puede acelerar el procesamiento de bucles utilizando las funciones de apply familiar (en su caso, probablemente se replicate ). Eche un vistazo al paquete plyr que proporciona barras de progreso.

Otra opción es evitar por completo los bucles y reemplazarlos con aritmética vectorizada. No estoy seguro exactamente de lo que está haciendo, pero probablemente pueda aplicar su función a todas las filas a la vez:

 temp[1:nrow(temp), 10] < - temp[1:nrow(temp), 9] + temp[0:(nrow(temp)-1), 10] 

Esto será mucho más rápido, y luego puede filtrar las filas con su condición:

 cond.i < - (temp[i, 6] == temp[i-1, 6]) & (temp[i, 3] == temp[i-1, 3]) temp[cond.i, 10] <- temp[cond.i, 9] 

La aritmética vectorizada requiere más tiempo y pensar en el problema, pero a veces puede ahorrar varios órdenes de magnitud en el tiempo de ejecución.

Procesar con data.table es una opción viable:

 n < - 1000000 df <- as.data.frame(matrix(sample(1:10, n*9, TRUE), n, 9)) colnames(df) <- paste("col", 1:9, sep = "") library(data.table) dayloop2.dt <- function(df) { dt <- data.table(df) dt[, Kumm. := { res <- .I; ifelse (res > 1, ifelse ((col6 == shift(col6, fill = 0)) & (col3 == shift(col3, fill = 0)) , res < - col9 + shift(res) , # else res <- col9 ) , # else res <- col9 ) } ,] res <- data.frame(dt) return (res) } res <- dayloop2.dt(df) m <- microbenchmark(dayloop2.dt(df), times = 100) #Unit: milliseconds # expr min lq mean median uq max neval #dayloop2.dt(df) 436.4467 441.02076 578.7126 503.9874 575.9534 966.1042 10 

Si ignora las posibles ganancias del filtrado de condiciones, es muy rápido. Obviamente, si puede hacer el cálculo en el subconjunto de datos, ayuda.