为什么与简单的 Rcpp 实现相比，zoo::rollmean 慢？

zoo::rollmean是一个有用的函数，它返回时间序列的滚动平均值；对于矢量x长度n和窗口大小k它返回向量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实现运行速度比以下快约 100 倍zoo::rollmean?

感谢@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 object:

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功能：

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 对象的内部结构，为其滚动平均操作增加了很多开销；与基本 R 实现相比，这会造成 10 倍的减速，与 Rcpp 实现相比，会造成 100 倍的减速。