R Language
厳しいからベクトル化するコードの高速化

Rcppを使用したループのタフなベクトル化の高速化

最初の要素が指定され（ first ）、各要素x_iがcos(x_{i-1} + 1)等しい長さlenベクトルを作成する次のtough-to-vectorize forループを考えてみましょう。

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

このコードでは、高速化（ cos(x[i-1]+1) ）したforループが使用されます。しかし、Rに "x + 1の累積余弦"関数がないので、この演算を基底Rでベクトル化することは自明ではありません。

この機能を高速化する方法としては、Rcppパッケージを使用してC ++で実装する方法があります。

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;
}")

これはしばしば大規模な計算に大きなスピードアップをもたらし、まったく同じ結果をもたらします：

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 

この場合、Rcppコードは、基底Rのアプローチでは1.31秒ではなく、0.03秒で100万の長さのベクトルを生成します。

バイトコンパイルによるループの難易度の高いベクトル化の高速化

最初の要素が指定され（ first ）、各要素x_iがcos(x_{i-1} + 1)等しい長さlenベクトルを作成する次の難しいtoベクトル化関数を考えてみましょう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)
}

1行のコードを書き直さずにこのような関数を高速化する簡単な方法の1つは、Rコンパイルパッケージを使用してコードをバイトコンパイルすることです。

library(compiler)
repeatedCosPlusOneCompiled <- cmpfun(repeatedCosPlusOne)

結果の関数は、しばしば同じ結果を返しながら、かなり速くなります：

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 

この場合、バイトコンパイルは、1.20秒から0.34秒までの長さ1百万のベクトル上の厳しいto-vectorize操作をスピードアップしました。

リマーク

repeatedCosPlusOneの本質は、単一の関数の累積的なアプリケーションとして、 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は、 repeatedCosPlusOne 「ベクトル化」とみなすことができます。しかし、それは2分の1で遅くなることが予想される：

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