Acelerar la operación de bucle en R

193

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

dayloop2 <- function(temp){
    for (i in 1:nrow(temp)){    
        temp[i,10] <- i
        if (i > 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?

Kay
fuente

Respuestas:

433

El mayor problema y la raíz de la ineficacia es indexar data.frame, me refiero a todas estas líneas donde las usa temp[,].
Intenta evitar esto tanto como sea posible. Tomé su función, cambio de indexación y aquí versión_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 puede ver, creo un vector resque reúne resultados. Al final lo agrego data.framey no necesito meterme con los nombres. Entonces, ¿qué tan mejor es?

Corro para cada función data.framecon nrowentre 1.000 y 10.000 antes de 1000 y medir el tiempo consystem.time

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

El resultado es

actuación

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

Poder de vectorización

Como Shane y Calimo afirman en sus respuestas, la vectorización es la clave para un mejor rendimiento. Desde su código puede moverse 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 de estas funciones, esta vez nrowde 10,000 a 100,000 por 10,000.

actuación

Afinando lo sintonizado

Otro ajuste es cambiar en una indexación de bucle temp[i,9]a res[i](que son exactamente iguales en la iteración del i-ésimo bucle). Nuevamente es la diferencia entre indexar un vector e indexar a data.frame.
Segunda cosa: cuando observa el bucle, puede ver que no hay necesidad de recorrer todo i, sino solo los 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 TRUEvalores en la condición. Para mis datos simulados, toma tiempo de cálculo para 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:

  • escribir un Ccódigo para hacer cumsum condicional
  • si sabe que en su secuencia máxima de datos no es grande, puede cambiar el ciclo a vectorizado, 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 .

Marek
fuente
2
Como no puedo encontrar una manera de preguntarle a Marek en privado, ¿cómo se generaron esos gráficos?
carbontwelve
@carbontwelve ¿Estás preguntando por datos o parcelas? Las parcelas se hicieron con paquete de celosía. Si tengo tiempo, pongo el código en algún lugar de la web y le doy aviso.
Marek
@carbontwelve Ooops, me equivoqué :) Estas son tramas estándar (de la base R).
Marek
@ Gregor Desafortunadamente no. Es acumulativo, por lo que no puede vectorizarlo. Un simple ejemplo: res = c(1,2,3,4)y condes todo TRUE, entonces resultado final debe ser: 1, 3(causa 1+2), 6(causa segunda es ahora 3, y tercero es 3también), 10( 6+4). Haciendo simple suma que tienes 1, 3, 5, 7.
Marek
Ah, debería haberlo pensado más detenidamente. Gracias por mostrarme el error.
Gregor Thomas
132

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 creación de perfiles 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 si es posible, por lo tanto:

El uso de funciones más eficientes puede producir ganancias de velocidad moderadas o grandes. Por ejemplo, paste0produce una pequeña ganancia de eficiencia pero .colSums()sus parientes producen ganancias algo más pronunciadas. meanEs particularmente lento .

Entonces puede evitar algunos problemas particularmente comunes :

  • cbind te ralentizará muy rápido.
  • Inicialice sus estructuras de datos, luego complételas, 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 que no valga la pena la molestia.
  • Echa un vistazo al R Inferno para ver más trampas que debes evitar.

Intente una mejor vectorización , que a menudo puede ayudar, pero no siempre. En este sentido, los comandos inherentemente vectorizados como ifelse, diffy similares proporcionarán más mejoras que la applyfamilia de comandos (que proporcionan poco o ningún aumento de velocidad en un bucle bien escrito).

También puede tratar de proporcionar más información a las funciones R . Por ejemplo, use en vapplylugar desapply , y especifique colClassesal leer datos basados ​​en texto . Las ganancias de velocidad serán variables dependiendo de cuantas conjeturas elimines.

A continuación, considere paquetes optimizados : el data.tablepaquete puede producir ganancias masivas de velocidad donde es posible su uso, 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 :

  • Compila tu script R. O utilice los paquetes Ray jiten concierto para la compilación justo a tiempo (Dirk tiene un ejemplo en esta presentación ).
  • Asegúrese de utilizar un BLAS optimizado. Estos proporcionan ganancias de velocidad en todos los ámbitos. Honestamente, es una pena que R no use automáticamente la biblioteca más eficiente en la instalación. Esperemos que Revolution R contribuya el trabajo que han hecho aquí a la comunidad en general.
  • Radford Neal ha realizado varias optimizaciones, algunas de las cuales se adoptaron en R Core, y muchas otras se bifurcaron en pqR .

Y, por último, si todo lo anterior aún no le da la velocidad que necesita, es posible que deba pasar a un idioma más rápido para el fragmento de código lento . La combinación de Rcppy inlineaquí 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 intento de hacerlo , y deslumbra incluso las soluciones R altamente optimizadas.

Si aún tiene problemas después de todo esto, solo necesita más potencia informática. Examine la paralelización ( http://cran.r-project.org/web/views/HighPerformanceComputing.html ) o incluso soluciones basadas en GPU ( gpu-tools).

Enlaces a otras orientaciones

Ari B. Friedman
fuente
35

Si está utilizando forbucles, lo más probable es que esté codificando R como si fuera C o Java u otra cosa. El código R que se vectoriza correctamente es extremadamente rápido.

Tome por ejemplo estos dos simples bits de código para generar una lista de 10,000 enteros en secuencia:

El primer ejemplo de código es cómo se codificaría un bucle utilizando un paradigma de codificación 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 por la simple acción de preasignar 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 usando la operación del vector base R usando el operador de colon, :esta operación es prácticamente instantánea:

system.time(a <- 1:1e5)

   user  system elapsed 
      0       0       0 
Andrie
fuente
+1 aunque consideraría su segundo ejemplo como poco convincente ya a[i]que no cambia. Pero system.time({a <- NULL; for(i in 1:1e5){a[i] <- 2*i} }); system.time({a <- 1:1e5; for(i in 1:1e5){a[i] <- 2*i} }); system.time({a <- NULL; a <- 2*(1:1e5)})tiene un resultado similar.
Henry
@Henry, comentario justo, pero como señalas, los resultados son los mismos. He modificado el ejemplo para inicializar un to rep(1, 1e5): los tiempos son idénticos.
Andrie
17

Esto podría hacerse mucho más rápido omitiendo los bucles utilizando índices o ifelse()sentencias 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."
Shane
fuente
Gracias por la respuesta. Trato de entender tus declaraciones. La línea 4: "temp [idx1,10] <- temp [idx1,9] + temp [which (idx1) -1,10]" causó un error porque la longitud del objeto más largo no es un múltiplo de la longitud del objeto más corto "temp [idx1,9] = num [1: 11496]" y "temp [which (idx1) -1,10] = int [1: 11494]" por lo que faltan 2 filas.
Kay
Si proporciona una muestra de datos (use dput () con algunas filas), lo arreglaré por usted. Debido a which () - 1 bit, los índices son desiguales. Pero debería ver cómo funciona desde aquí: no hay necesidad de realizar ningún bucle o aplicación; solo use funciones vectorizadas.
Shane
1
¡Guauu! ¡Acabo de cambiar un bloque de función if..else anidado y mapply, a una función ifelse anidada y obtuve una aceleración de 200x!
James
Su consejo general es correcto, pero en el código que omitió el hecho, ese ivalor depende de i-1-th, por lo que no se pueden configurar de la manera en que lo hace (usando which()-1).
Marek
8

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

Con frecuencia 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 versiones:

   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 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 esto en class(d)='list'lugar de d=as.list(d)yclass(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)
Chris
fuente
1
Probablemente sea gracias a la sobrecarga de [<-.data.frame, que de alguna manera se llama cuando lo hace d$foo[i] = marky puede terminar haciendo una nueva copia del vector de posiblemente el data.frame completo en cada <-modificación. Haría una pregunta interesante sobre SO.
Frank
2
@Frank It (i) debe asegurarse de que el objeto modificado sigue siendo un data.frame válido y (ii) afaik realiza al menos una copia, posiblemente más de una. Se sabe que la subasignación de marco de datos es lenta y, si observa el código fuente largo, no es realmente sorprendente.
Roland
@Frank, @Roland: ¿La df$var[i]notación pasa por la misma [<-.data.framefunción? Me di cuenta de que es bastante largo. Si no, ¿qué función usa?
Chris
@Chris creo que d$foo[i]=markse traduce aproximadamente d <- `$<-`(d, 'foo', `[<-`(d$foo, i, mark)), pero con algún uso de variables temporales.
Tim Goodman
7

Como Ari mencionó al final de su respuesta, los paquetes Rcppy inlinehacen que sea increíblemente fácil hacer las cosas rápido. Como ejemplo, intente este inlinecódigo (advertencia: no probado):

body <- 'Rcpp::NumericMatrix nm(temp);
         int nrtemp = Rccp::as<int>(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 #includelas cosas, donde simplemente pasa un parámetro

inc <- '#include <header.h>

a cxxfunction, como include=inc . Lo realmente bueno de esto es que hace todo el enlace y la compilación por ti, por lo que la creación de prototipos es realmente rápida.

Descargo de responsabilidad: no estoy totalmente seguro de que la clase de tmp deba ser numérica y no matriz numérica u otra cosa. Pero estoy casi seguro.

Editar: si aún necesita más velocidad después de esto, OpenMP es un recurso de paralelización bueno para C++. No he intentado usarlo inline, pero debería funcionar. La idea sería, en el caso de los nnúcleos, hacer que la iteración de bucle kse lleve a cabo k % n. Una introducción adecuada se encuentra en Matloff es el arte de la programación R , disponible aquí , en el capítulo 16, Recurrir a C .

jclancy
fuente
3

Las respuestas aquí son geniales. Un aspecto menor no cubierto es que la pregunta dice " Mi PC todavía funciona (aproximadamente 10 horas ahora) y no tengo idea sobre el tiempo de ejecución ". Siempre pongo el siguiente código en bucles cuando desarrollo 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 bucle es bastante rápida pero el número de bucles es grande, considere imprimir de vez en cuando, ya que imprimir en la consola tiene una sobrecarga. 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)
}
novato
fuente
Una opción similar, imprimir la fracción i / n. Siempre tengo algo así, cat(sprintf("\nNow running... %40s, %s/%s \n", nm[i], i, n))ya que generalmente estoy recorriendo cosas con nombre (con nombres en nm).
Frank
2

En R, a menudo puede acelerar el procesamiento del bucle utilizando las applyfunciones familiares (en su caso, probablemente lo sería replicate). Echa un vistazo al plyrpaquete que proporciona barras de progreso.

Otra opción es evitar los bucles por completo 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.

Calimo
fuente
14
Estás en el punto en que las funciones vectoriales serán más rápidas que los bucles o apply (), pero no es cierto que apply () sea más rápido que los bucles. En muchos casos, apply () es simplemente abstraer el bucle del usuario, pero aún continúa. Vea esta pregunta anterior: stackoverflow.com/questions/2275896/…
JD Long
0

Procesar con data.tablees 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, es útil.

Bulat
fuente
2
¿Por qué repite la sugerencia de usar data.table? Ya se ha hecho varias veces en las respuestas anteriores.
IRTFM