R Language
Acelerar el código difícil de vectorizar

Aceleración de vectorización difícil para bucles con Rcpp

Considere el siguiente bucle for-vectorize for, que crea un vector de longitud len donde se especifica el primer elemento ( first ) y cada elemento x_i es igual a cos(x_{i-1} + 1) :

repeatedCosPlusOne <- function(first, len) {
  x <- numeric(len)
  x[1] <- first
  for (i in 2:len) {
    x[i] <- cos(x[i-1] + 1)
  }
  return(x)
}

Este código implica un bucle for con una operación rápida ( cos(x[i-1]+1) ), que a menudo se beneficia de la vectorización. Sin embargo, no es trivial vectorizar esta operación con la base R, ya que R no tiene una función de "coseno acumulativo de x + 1".

Un posible enfoque para acelerar esta función sería implementarla en C ++, utilizando el paquete Rcpp:

library(Rcpp)
cppFunction("NumericVector repeatedCosPlusOneRcpp(double first, int len) {
  NumericVector x(len);
  x[0] = first;
  for (int i=1; i < len; ++i) {
    x[i] = cos(x[i-1]+1);
  }
  return x;
}")

Esto a menudo proporciona aceleraciones significativas para cálculos grandes al mismo tiempo que produce los mismos resultados exactos:

all.equal(repeatedCosPlusOne(1, 1e6), repeatedCosPlusOneRcpp(1, 1e6))
# [1] TRUE
system.time(repeatedCosPlusOne(1, 1e6))
#    user  system elapsed 
#   1.274   0.015   1.310 
system.time(repeatedCosPlusOneRcpp(1, 1e6))
#    user  system elapsed 
#   0.028   0.001   0.030 

En este caso, el código Rcpp genera un vector de 1 millón de longitud en 0.03 segundos en lugar de 1.31 segundos con el enfoque de la base R.

Aceleración de vectores difíciles para la compilación de bytes por bucles

Siguiendo el ejemplo de Rcpp en esta entrada de documentación, considere la siguiente función difícil de vectorizar, que crea un vector de longitud len donde se especifica el primer elemento ( first ) y cada elemento x_i es igual a cos(x_{i-1} + 1) :

repeatedCosPlusOne <- function(first, len) {
  x <- numeric(len)
  x[1] <- first
  for (i in 2:len) {
    x[i] <- cos(x[i-1] + 1)
  }
  return(x)
}

Un enfoque simple para acelerar dicha función sin volver a escribir una sola línea de código es compilar el código mediante el paquete de compilación R:

library(compiler)
repeatedCosPlusOneCompiled <- cmpfun(repeatedCosPlusOne)

La función resultante a menudo será significativamente más rápida mientras sigue devolviendo los mismos resultados:

all.equal(repeatedCosPlusOne(1, 1e6), repeatedCosPlusOneCompiled(1, 1e6))
# [1] TRUE
system.time(repeatedCosPlusOne(1, 1e6))
#    user  system elapsed 
#   1.175   0.014   1.201 
system.time(repeatedCosPlusOneCompiled(1, 1e6))
#    user  system elapsed 
#   0.339   0.002   0.341 

En este caso, la compilación de bytes aceleró la operación de vectorización dura en un vector de 1 millón de longitud, de 1,20 segundos a 0,34 segundos.

Observación

La esencia de la repeatedCosPlusOne , como la aplicación acumulativa de una sola función, se puede expresar de forma más transparente con Reduce :

iterFunc <- function(init, n, func) {
  funcs <- replicate(n, func)
  Reduce(function(., f) f(.), funcs, init = init, accumulate = TRUE)
}
repeatedCosPlusOne_vec <- function(first, len) {
  iterFunc(first, len - 1, function(.) cos(. + 1))
}

repeatedCosPlusOne_vec puede considerarse como una "vectorización" de repeatedCosPlusOne . Sin embargo, se puede esperar que sea más lento por un factor de 2:

library(microbenchmark)
microbenchmark(
  repeatedCosPlusOne(1, 1e4),
  repeatedCosPlusOne_vec(1, 1e4)
)
#> Unit: milliseconds
#>                              expr       min        lq     mean   median       uq      max neval cld
#>      repeatedCosPlusOne(1, 10000)  8.349261  9.216724 10.22715 10.23095 11.10817 14.33763   100  a 
#>  repeatedCosPlusOne_vec(1, 10000) 14.406291 16.236153 17.55571 17.22295 18.59085 24.37059   100   b