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

【问题标题】：Why is zoo::rollmean slow compared to a simple Rcpp implementation?与简单的 Rcpp 实现相比，为什么 zoo::rollmean 慢？
【发布时间】：2015-07-17 09:25:21
【问题描述】：

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 倍？

【问题讨论】：

RcppRoll 包提供更快的zoo::rolls 实现。

标签： r zoo moving-average

【解决方案1】：

在zoo 中四处寻找，似乎rollmean.* 方法都在R 中实现。

而您在 C++ 中实现了一个。打包的 R 代码可能还会进行更多检查等 pp 所以也许你击败它并不奇怪？

【讨论】：

【解决方案2】：

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

为了深入研究为什么我们在使用 base 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 函数构造一个新的动物园对象上：

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 倍的速度下降。

【讨论】：