4.4 Rendimiento en R

“We should forget about small efficiencies, say about 97% of the time: premature optimization is the root of all evil. Yet we should not pass up our opportunities in that critical 3%. A good programmer will not be lulled into complacency by such reasoning, he will be wise to look carefully at the critical code; but only after that code has been identified.” -Donald Knuth

Diseña primero, luego optimiza. La optimización del código es un proceso iterativo:
1. Encuentra el cuello de botella más importante.
2. Intenta eliminarlo (no siempre se puede).
3. Repite hasta que tu código sea lo suficientemente rápido.

Diagnosticar

Una vez que tienes código que se puede leer y funciona, el perfilamiento (profiling) del código es un método sistemático que nos permite conocer cuanto tiempo se esta usando en diferentes partes del programa.

Comenzaremos con la función system.time (no es perfilamiento aún), esta calcula el tiempo en segundos que toma ejecutar una expresión (si hay un error, regresa el tiempo hasta que ocurre el error):

batting <- read.csv("data/batting.csv")
system.time(lm(R ~ AB + teamID, batting))
#>    user  system elapsed 
#>   2.996   0.048   3.043

• user time: Tiempo usado por el CPU(s) para evaluar esta expresión, tiempo que experimenta la computadora.
• elapsed time: tiempo en el reloj, tiempo que experimenta la persona.

Notemos que el tiempo de usuario (user) puede ser mayor al tiempo transcurrido (elapsed),

system.time(readLines("http://www.jhsph.edu"))
#>    user  system elapsed 
#>   0.024   0.000   1.161

o al revés:

library(parallel)
system.time(mclapply(2000:2007, 
  function(x){
    sub <- subset(batting, yearID == x)
    lm(R ~ AB + playerID, sub)
}, mc.cores = 7))
#>    user  system elapsed 
#>   9.540   0.384   7.218

Comparemos la velocidad de dplyr con funciones que se encuentran en R estándar y plyr.

# dplyr
dplyr_st <- system.time({
  batting %>%
  group_by(playerID) %>%
  summarise(total = sum(R, na.rm = TRUE), n = n()) %>%
  dplyr::arrange(desc(total))
})

# plyr
plyr_st <- system.time({
    batting %>% 
    plyr::ddply("playerID", plyr::summarise, total = sum(R, na.rm = TRUE), 
        n = length(R)) %>% 
    plyr::arrange(-total)
})

# estándar lento
est_l_st <- system.time({
  players <- unique(batting$playerID)
  n_players <- length(players)
  total <- rep(NA, n_players)
  n <- rep(NA, n_players)
  for(i in 1:n_players){
    sub_batting <- batting[batting$playerID == players[i], ]
    total[i] <- sum(sub_batting$R, na.rm = TRUE)
    n[i] <- nrow(sub_batting)
  }
  batting_2 <- data.frame(playerID = players, total = total, n = n)
  batting_2[order(batting_2$total, decreasing = TRUE), ]
})

# estándar rápido
est_r_st <- system.time({
  batting_2 <- aggregate(. ~ playerID, data = batting[, c("playerID", "R")], 
      sum)
  batting_ord <- batting_2[order(batting_2$R, decreasing = TRUE), ]
})

dplyr_st
#>    user  system elapsed 
#>   0.036   0.000   0.038
plyr_st
#>    user  system elapsed 
#>   4.360   0.012   4.374
est_l_st
#>    user  system elapsed 
#> 192.924   2.352 195.322
est_r_st
#>    user  system elapsed 
#>   0.128   0.008   0.136

La función system.time supone que sabes donde buscar, es decir, que expresiones debes evaluar, una función que puede ser más útil cuando uno desconoce cuál es la función que alenta un programa es profvis() del paquete con el mismo nombre.

library(profvis)
batting_recent <- filter(batting, yearID > 2006)
profvis({
  players <- unique(batting_recent$playerID)
  n_players <- length(players)
  total <- rep(NA, n_players)
  n <- rep(NA, n_players)
  for(i in 1:n_players){
    sub_batting <- batting_recent[batting_recent$playerID == players[i], ]
    total[i] <- sum(sub_batting$R, na.rm = TRUE)
    n[i] <- nrow(sub_batting)
  }
  batting_2 <- data.frame(playerID = players, total = total, n = n)
  batting_2[order(batting_2$total, decreasing = TRUE), ]
})

profvis() utiliza a su vez la función Rprof() de R base, este es un perfilador de muestreo que registra cambios en la pila de funciones, funciona tomando muestras a intervalos regulares y tabula cuánto tiempo se lleva en cada función.

Estrategias para mejorar desempeño

Algunas estrategias para mejorar desempeño:

1. Utilizar apropiadamente funciones de R, o funciones de paquetes que muchas veces están mejor escritas de lo que nosotros podríamos hacer.
2. Hacer lo menos posible.
3. Usar funciones vectorizadas en R (casi siempre). No hacer crecer objetos (es preferible definir su tamaño antes de operar en ellos).
4. Paralelizar.
5. La más simple y muchas veces la más barata: conseguir una máquina más grande (por ejemplo Amazon web services).

A continuación revisamos y ejemplificamos los puntos anteriores, los ejemplos de código se tomaron del taller EfficientR, impartido por Martin Morgan.

Utilizar apropiadamente funciones de R

Si el cuello de botella es la función de un paquete vale la pena buscar alternativas, CRAN task views es un buen lugar para buscar.

Hacer lo menos posible

Utiliza funciones más específicas, por ejemplo:
* rowSums(), colSums(), rowMeans() y colMeans() son más rápidas que las invocaciones equivalentes de apply().

• Si quieres checar si un vector contiene un valor any(x == 10) es más veloz que 10 %in% x, esto es porque examinar igualdad es más sencillo que examinar inclusión en un conjunto.
  Este conocimiento requiere que conozcas alternativas, para ello debes construir tu vocabulario, puedes comenzar por lo básico e ir incrementando conforme lees código.
  Otro caso es cuando las funciones son más rápidas cunado les das más información del problema, por ejemplo:

• read.csv(), especificar las clases de las columnas con colClasses.
  

• factor() especifica los niveles con el argumento levels.

Usar funciones vectorizadas en R

Es común escuchar que en R vectorizar es conveniente, el enfoque vectorizado va más allá que evitar ciclos for:

• Pensar en objetos, en lugar de enfocarse en las componentes de un vector, se piensa únicamente en el vector completo.

• Los ciclos en las funciones vectorizadas de R están escritos en C, lo que los hace más veloces.

Las funciones vectorizadas programadas en R pueden mejorar la interfaz de una función pero no necesariamente mejorar el desempeño. Usar vectorización para desempeño implica encontrar funciones de R implementadas en C.

Al igual que en el punto anterior, vectorizar requiere encontrar las funciones apropiadas, algunos ejemplos incluyen: _rowSums(), colSums(), rowMeans() y colMeans().

Ejemplo: iteración (for, lapply(), sapply(), vapply(), mapply(), apply(), …) en un vector de n elementos llama a R base n veces

compute_pi0 <- function(m) {
    s = 0
    sign = 1
    for (n in 0:m) {
        s = s + sign / (2 * n + 1)
        sign = -sign
    }
    4 * s
}

compute_pi1 <- function(m) {
    even <- seq(0, m, by = 2)
    odd <- seq(1, m, by = 2)
    s <- sum(1 / (2 * even + 1)) - sum(1 / (2 * odd + 1))
    4 * s
}
m <- 1e6

Utilizamos el paquete microbenchmark para medir tiempos varias veces.

library(microbenchmark)
m <- 1e4
result <- microbenchmark(
    compute_pi0(m),
    compute_pi0(m * 10),
    compute_pi0(m * 100),
    compute_pi1(m),
    compute_pi1(m * 10),
    compute_pi1(m * 100),
    compute_pi1(m * 1000),
    times = 20
)
result
#> Unit: microseconds
#>                   expr        min          lq        mean     median
#>         compute_pi0(m)    684.570    701.5120    743.2516    707.013
#>    compute_pi0(m * 10)   6904.621   6949.0440   7094.2657   6974.760
#>   compute_pi0(m * 100)  69539.332  70290.6035  70850.1057  70593.476
#>         compute_pi1(m)    246.899    264.0855    319.2403    297.886
#>    compute_pi1(m * 10)   2257.538   2310.1740   2873.6184   2329.766
#>   compute_pi1(m * 100)  22464.707  23488.9330  26011.5717  23742.731
#>  compute_pi1(m * 1000) 293987.507 298651.9915 357810.6187 333719.561
#>           uq        max neval
#>     716.5065   1174.749    20
#>    7215.8510   7647.342    20
#>   71265.4850  74877.953    20
#>     377.5405    400.109    20
#>    2411.1540  12930.811    20
#>   30343.4570  33618.310    20
#>  422141.7710 459192.225    20

Evitar copias

Otro aspecto importante es que generalmente conviene asignar objetos en lugar de hacerlos crecer (es más eficiente asignar toda la memoria necesaria antes del cálculo que asignarla sucesivamente). Esto es porque cuando se usan instrucciones para crear un objeto más grande (e.g. append(), cbind(), c(), rbind()) R debe primero asignar espacio a un nuevo objeto y luego copiar al nuevo lugar. Para leer más sobre esto Burns (2015) es una buena referencia.

Ejemplo: crecer un vector puede causar que R copie de manera repetida el vector chico en el nuevo vector, aumentando el tiempo de ejecución.

Solución: crear vector de tamaño final y llenarlo con valores. Las funciones como lapply() y map hacen esto de manera automática y son más sencillas que los ciclos for.

memory_copy1 <- function(n) {
    result <- numeric()
    for (i in seq_len(n))
        result <- c(result, 1/i)
    result
}

memory_copy2 <- function(n) {
    result <- numeric()
    for (i in seq_len(n))
        result[i] <- 1 / i
    result
}

pre_allocate1 <- function(n) {
    result <- numeric(n)
    for (i in seq_len(n))
        result[i] <- 1 / i
    result
}

pre_allocate2 <- function(n) {
    vapply(seq_len(n), function(i) 1 / i, numeric(1))
}

vectorized <- function(n) {
    1 / seq_len(n)
}

n <- 10000
microbenchmark(
    memory_copy1(n),
    memory_copy2(n),
    pre_allocate1(n),
    pre_allocate2(n),
    vectorized(n),
    times = 10, unit = "relative"
)
#> Unit: relative
#>              expr        min         lq       mean     median         uq
#>   memory_copy1(n) 4461.24790 4371.43198 518.426615 4008.84336 3268.67323
#>   memory_copy2(n)   77.17223   74.87649   9.981552   66.78001   57.63065
#>  pre_allocate1(n)   16.98927   15.93934   3.068832   13.94142   11.41144
#>  pre_allocate2(n)  166.12215  155.26760  24.136407  142.41561  125.37996
#>     vectorized(n)    1.00000    1.00000   1.000000    1.00000    1.00000
#>        max neval
#>  65.295191    10
#>   2.360671    10
#>   1.689472    10
#>   7.606423    10
#>   1.000000    10

Un caso común donde se hacen copias sin necesidad es al trabajar con data.frames.

Ejemplo: actualizar un data.frame copia el data.frame completo.

Solución: operar en vectores y actualiza el data.frame al final.

n <- 1e4
df <- data.frame(Index = 1:n, A = seq(10, by = 1, length.out = n))

f1 <- function(df) {
    ## constants
    cost1 <- 3
    cost2 <- 0.05
    cost3 <- 50

    ## update data.frame -- copies entire data frame each time!
    df$S[1] <- cost1
    for (j in 2:(n))
        df$S[j] <- df$S[j - 1] - cost3 + df$S[j - 1] * cost2 / 12

    ## return result
    df
}
.f2helper <- function(cost1, cost2, cost3, n) {
    ## create the result vector separately
    cost2 <- cost2 / 12   # 'hoist' common operations
    result <- numeric(n)
    result[1] <- cost1
    for (j in 2:(n))
        result[j] <- (1 + cost2) * result[j - 1] - cost3
    result
}

f2 <- function(df) {
    cost1 <- 3
    cost2 <- 0.05
    cost3 <- 50

    ## update the data.frame once
    df$S <- .f2helper(cost1, cost2, cost3, n)
    df
}

microbenchmark(
    f1(df),
    f2(df),
    times = 5, unit = "relative"
)
#> Unit: relative
#>    expr      min       lq     mean   median       uq      max neval
#>  f1(df) 355.3847 352.6145 129.6952 352.3906 104.6378 50.54361     5
#>  f2(df)   1.0000   1.0000   1.0000   1.0000   1.0000  1.00000     5

Paralelizar

Paralelizar usa varios cores para trabajar de manera simultánea en varias secciones de un problema, no reduce el tiempo computacional pero incrementa el tiempo del usuario pues aprovecha los recursos. Como referencia está [Parallel Computing for Data Science] de Norm Matloff.

Referencias

Burns, P. 2015. The R Inferno. The American Statistician. Vol. 4. 69.