使用 Rcpp 为循环加速难以矢量化
考虑以下难以向量化的 for 循环,它创建一个长度为 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)
}
这段代码涉及一个快速操作的 for 循环(cos(x[i-1]+1)
),它通常受益于矢量化。然而,用基 R 对该操作进行矢量化并不是一件容易的事,因为 R 不具有“x + 1 的累积余弦”函数。
加速此功能的一种可能方法是使用 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 方法在 0.03 秒内生成长度为 100 万的向量,而不是 1.31 秒。