Why is zoo::rollmean slow compared to a simple Rcpp implementation?

Question

zoo::rollmean is a helpful function that returns the rolling mean of a time series; for vector x of length n and window size k it returns the vector c(mean(x[1:k]), mean(x[2:(k+1)]), ..., mean(x[(n-k+1):n])).

I noticed that it seemed to be running slowly for some code I was developing, so I wrote my own version using the Rcpp package and a simple for loop:

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;
}")

To my surprise, this version of the function is much faster than the zoo::rollmean function:

# 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

The speedup holds even for much larger vectors:

# 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

Why does a simple Rcpp implementation run ~100x faster than zoo::rollmean?

Answer 1

Poking around in zoo it seem that the rollmean.* methods are all in implemented in R.

Whereas you implemented one in C++. The packaged R code probably also does a few more checks etc pp so maybe it is not so surprising that you beat it?