通过字节编译加速对循环进行难以矢量化
按照本文档条目中的 Rcpp 示例,考虑以下难以向量化的函数,该函数创建一个长度为 len
的向量,其中指定第一个元素(first
)并且每个元素 x_i
等于 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)
}
在不重写单行代码的情况下加速这种函数的一种简单方法是使用 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 百万的向量上的难以向量化的操作,从 1.20 秒到 0.34 秒。
备注
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 倍:
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