R Language
Accélérer le code difficile à vectoriser

Accélérer la vectorisation des boucles avec Rcpp

Considérons la boucle difficile à vectoriser suivante, qui crée un vecteur de longueur len où le premier élément est spécifié (en first ) et chaque élément x_i est égal à 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)
}

Ce code implique une boucle for avec une opération rapide ( cos(x[i-1]+1) ), qui bénéficie souvent de la vectorisation. Cependant, il n'est pas trivial de vectoriser cette opération avec la base R, puisque R n'a pas de fonction "cosinus cumulatif de x + 1".

Une approche possible pour accélérer cette fonction serait de l'implémenter en C ++, en utilisant le package 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;
}")

Cela fournit souvent des accélérations significatives pour les gros calculs tout en donnant les mêmes résultats:

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 

Dans ce cas, le code Rcpp génère un vecteur de longueur 1 million en 0,03 seconde au lieu de 1,31 seconde avec l'approche base R.

Accélérer la vectorisation des boucles en compilant les octets

En suivant l'exemple de Rcpp dans cette entrée de documentation, considérez la fonction difficile à vectoriser suivante, qui crée un vecteur de longueur len où le premier élément est spécifié (en first ) et chaque élément x_i est égal à 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)
}

Une approche simple pour accélérer une telle fonction sans réécriture d'une seule ligne de code consiste à compiler l'octet en utilisant le package de compilation R:

library(compiler)
repeatedCosPlusOneCompiled <- cmpfun(repeatedCosPlusOne)

La fonction résultante sera souvent beaucoup plus rapide tout en renvoyant les mêmes résultats:

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 

Dans ce cas, la compilation des octets a accéléré l'opération difficile à vectoriser sur un vecteur de longueur 1 million de 1,20 seconde à 0,34 seconde.

Remarque

L'essence de repeatedCosPlusOne , en tant qu'application cumulative d'une fonction unique, peut être exprimée de manière plus transparente avec 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 peut être considéré comme une "vectorisation" de repeatedCosPlusOne . Cependant, on peut s'attendre à ce qu'il soit plus lent d'un facteur 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