为什么相对于简单的Rcpp实现，zoo::rollmean很慢？

Question

zoo::rollmean是一个有用的函数，它返回时间序列的滚动均值；对于长度为 n、窗口大小为 k 的向量 x，它返回向量 c(mean(x[1:k]), mean(x[2:(k+1)]), ..., mean(x[(n-k+1):n]))。

我注意到这个函数在我正在开发的一些代码中运行得比较慢，因此我使用 Rcpp 包和一个简单的 for 循环编写了自己的版本：

library(Rcpp)
cppFunction("NumericVector rmRcpp(NumericVector dat, const int window) {
  const int n = dat.size();
  NumericVector ret(n-window+1);
  double summed = 0.0;
  for (int i=0; i < window; ++i) {
    summed += dat[i];
  }
  ret[0] = summed / window;
  for (int i=window; i < n; ++i) {
    summed += dat[i] - dat[i-window];
    ret[i-window+1] = summed / window;
  }
  return ret;
}")

令我惊讶的是，这个函数版本比 zoo::rollmean 函数快得多：

# Time series with 1000 elements
set.seed(144)
y <- rnorm(1000)
x <- 1:1000
library(zoo)
zoo.dat <- zoo(y, x)

# Make sure our function works
all.equal(as.numeric(rollmean(zoo.dat, 3)), rmRcpp(y, 3))
# [1] TRUE

# Benchmark
library(microbenchmark)
microbenchmark(rollmean(zoo.dat, 3), rmRcpp(y, 3))
# Unit: microseconds
#                  expr     min       lq       mean    median        uq       max neval
#  rollmean(zoo.dat, 3) 685.494 904.7525 1776.88666 1229.2475 1744.0720 15724.321   100
#          rmRcpp(y, 3)   6.638  12.5865   46.41735   19.7245   27.4715  2418.709   100

即使对于更大的向量，速度提升仍然存在：

# Time series with 5 million elements
set.seed(144)
y <- rnorm(5000000)
x <- 1:5000000
library(zoo)
zoo.dat <- zoo(y, x)

# Make sure our function works
all.equal(as.numeric(rollmean(zoo.dat, 3)), rmRcpp(y, 3))
# [1] TRUE

# Benchmark
library(microbenchmark)
microbenchmark(rollmean(zoo.dat, 3), rmRcpp(y, 3), times=10)
# Unit: milliseconds
#                  expr        min         lq       mean     median         uq        max
#  rollmean(zoo.dat, 3) 2825.01622 3090.84353 3191.87945 3206.00357 3318.98129 3616.14047
#          rmRcpp(y, 3)   31.03014   39.13862   42.67216   41.55567   46.35191   53.01875

为什么一个简单的Rcpp实现比zoo::rollmean运行得快大约100倍？

- josliber

---

RcppRoll 软件包提供了 zoo::roll 函数更快的实现方式。 - ExperimenteR

Answer 1

感谢@DirkEddelbuettel指出问题，问题中的比较并不公平，因为我在将C++代码与纯R代码进行比较。下面是一个简单的基本R实现（没有来自zoo包的所有检查）；这与zoo::rollmean实现滚动均值的核心计算方式非常相似：

baseR.rollmean <- function(dat, window) {
  n <- length(dat)
  y <- dat[window:n] - dat[c(1, 1:(n-window))]
  y[1] <- sum(dat[1:window])
  return(cumsum(y) / window)
}

与zoo:rollmean相比，我们可以看出这个要快得多：

set.seed(144)
y <- rnorm(1000000)
x <- 1:1000000
library(zoo)
zoo.dat <- zoo(y, x)
all.equal(as.numeric(rollmean(zoo.dat, 3)), baseR.rollmean(y, 3), RcppRoll::roll_mean(y, 3), rmRcpp(y, 3))
# [1] TRUE
library(microbenchmark)
microbenchmark(rollmean(zoo.dat, 3), baseR.rollmean(y, 3), RcppRoll::roll_mean(y, 3), rmRcpp(y, 3), times=10)
# Unit: milliseconds
#                       expr        min         lq       mean     median         uq        max neval
#       rollmean(zoo.dat, 3) 507.124679 516.671897 646.813716 563.897005 593.861499 1220.08272    10
#       baseR.rollmean(y, 3)  46.156480  47.804786  53.923974  49.250144  55.061844   76.47908    10
#  RcppRoll::roll_mean(y, 3)   7.714032   8.513042   9.014886   8.693255   8.885514   11.32817    10
#               rmRcpp(y, 3)   7.729959   8.045270   8.924030   8.388931   8.996384   12.49042    10

为了深入研究在使用基本R时我们为什么看到了10倍的加速，我使用了Hadley的lineprof工具，并在需要时从zoo包源代码中提取源代码：

lineprof(rollmean.zoo(zoo.dat, 3))
#     time  alloc release dups ref                  src
# 1  0.001  0.954       0   26 #27 rollmean.zoo/unclass
# 2  0.001  0.954       0    0 #28 rollmean.zoo/:      
# 3  0.002  0.954       0    1 #28 rollmean.zoo        
# 4  0.001  1.431       0    0 #28 rollmean.zoo/seq_len
# 5  0.001  0.000       0    0 #28 rollmean.zoo/c      
# 6  0.006  2.386       0    1 #28 rollmean.zoo        
# 7  0.002  0.954       0    2 #31 rollmean.zoo/cumsum 
# 8  0.001  0.000       0    0 #31 rollmean.zoo//      
# 9  0.005  1.912       0    1 #33 rollmean.zoo        
# 10 0.013  2.898       0   14 #33 rollmean.zoo/[<-    
# 11 0.299 28.941       0  127 #34 rollmean.zoo/na.fill

很明显，几乎所有时间都花在了na.fill函数上，实际上这个函数是在滚动均值已经计算完成之后才被调用的。

lineprof(na.fill.zoo(zoo.dat, fill=NA, 2:999999))
#     time  alloc release dups ref                  src
# 1  0.004  1.913       0   39 #26 na.fill.zoo/seq     
# 2  0.002  1.921       0    9 #33 na.fill.zoo/coredata
# 3  0.002  1.921       0    6 #37 na.fill.zoo/[<-     
# 4  0.001  0.955       0   10 #46 na.fill.zoo         
# 5  0.008  3.838       0   19 #46 na.fill.zoo/[<-     
# 6  0.003  0.959       0    2 #52 na.fill.zoo         
# 7  0.006  0.972       0   21 #52 na.fill.zoo/[<-     
# 8  0.001  0.486       0    0 #57 na.fill.zoo/seq_len 
# 9  0.005  0.959       0    6 #66 na.fill.zoo         
# 10 0.124 11.573       0   34 #66 na.fill.zoo/[ 

大部分时间都花在对zoo对象进行子集操作上：

lineprof("[.zoo"(zoo.dat, 2:999999))
#    time  alloc release dups          ref            src
# 1 0.004  0.004       0    0 character(0)               
# 2 0.002  1.922       0    4           #4 [.zoo/coredata
# 3 0.038 11.082       0   29          #19 [.zoo/zoo     
# 4 0.004  0.000       0    1          #28 [.zoo 

几乎所有时间都用于使用zoo函数构建新的zoo对象：

lineprof(zoo(y[2:999999], 2:999999))
#    time alloc release dups                ref        src
# 1 0.021 4.395       0    8 c("zoo", "unique") zoo/unique
# 2 0.012 0.477       0    8  c("zoo", "ORDER") zoo/ORDER 
# 3 0.001 0.477       0    1              "zoo" zoo       
# 4 0.001 0.954       0    0      c("zoo", ":") zoo/:     
# 5 0.015 3.341       0    5              "zoo" zoo      

设置新的zoo对象所需的各种操作（例如确定唯一时间点并对它们进行排序）。

总之，zoo包似乎通过构建新的zoo对象而不是使用当前zoo对象的内部来增加了滚动均值操作的开销；这导致与基本R实现相比慢了10倍，与Rcpp实现相比慢了100倍。

- josliber

Answer 2

在研究 zoo 的时候，似乎发现所有的 rollmean.* 方法都是用 R 实现的。

而你使用了 C++ 实现一个方法。打包的 R 代码可能也做了一些额外的检查等等，所以你超过它也许并不那么令人惊讶？

- Dirk Eddelbuettel