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):

data("Batting", package = "Lahman") 
system.time(lm(R ~ AB + teamID, Batting))
#>    user  system elapsed 
#>   3.782   0.132   3.914

• 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 menor al tiempo transcurrido (elapsed),

system.time(readLines("http://www.jhsph.edu"))
#>    user  system elapsed 
#>   0.017   0.004   0.680

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 
#>  12.734   0.749   9.747

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.113   0.000   0.113
plyr_st
#>    user  system elapsed 
#>   4.081   0.015   4.097
est_l_st
#>    user  system elapsed 
#>  61.097   0.156  61.255
est_r_st
#>    user  system elapsed 
#>    0.39    0.00    0.39

La función system.time supone que sabes donde buscar, es decir, que sabes 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)    784.123    792.4700    808.6164    799.2345
#>    compute_pi0(m * 10)   7870.525   7916.8640   8321.5265   7947.4775
#>   compute_pi0(m * 100)  78685.802  79211.9975  80101.5360  79777.1910
#>         compute_pi1(m)    159.604    251.3475    267.6655    292.2355
#>    compute_pi1(m * 10)   1238.680   1313.3290   1763.8384   1388.0875
#>   compute_pi1(m * 100)  12729.967  12886.8455  23304.1597  18982.4395
#>  compute_pi1(m * 1000) 253927.749 358522.2820 357314.8858 364344.4460
#>           uq        max neval
#>     824.0710    871.498    20
#>    8155.9330  14220.185    20
#>   80910.3130  83693.327    20
#>     304.7665    321.351    20
#>    1442.7955   9186.193    20
#>   22289.3580 128275.896    20
#>  373671.3445 498513.471    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) 5537.66878 4543.10043 668.948791 4050.45337 3627.59050
#>   memory_copy2(n)   91.23861   76.10998  11.923358   70.27946   67.66902
#>  pre_allocate1(n)   19.86093   16.45651   3.738255   14.68634   13.53066
#>  pre_allocate2(n)  194.01336  163.16694  23.303504  148.40098  136.98439
#>     vectorized(n)    1.00000    1.00000   1.000000    1.00000    1.00000
#>         max neval
#>  136.454492    10
#>    2.920912    10
#>    2.038541    10
#>    4.136380    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) 246.1295 243.2241 81.37993 243.5072 61.90639 32.68433     5
#>  f2(df)   1.0000   1.0000  1.00000   1.0000  1.00000  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.