自适应移动平均线 - R 中的顶级性能

Question

我正在寻找 R 中滚动/滑动窗口函数方面的一些性能提升。这是一项非常常见的任务，可用于任何有序的观察数据集。我想分享一些我的发现，也许有人可以提供反馈以使其更快。
重要的一点是我关注align="right" 和自适应滚动窗口的情况，所以width 是一个向量（与我们的观察向量长度相同）。如果我们将width 作为标量，那么zoo 和TTR 包中已经有非常完善的功能，这些功能将很难被击败（4 年后：这比我容易预期）因为其中一些甚至使用 Fortran（但使用下面提到的wapply 仍然可以更快地使用用户定义的 FUN）。
RcppRoll 包由于其出色的性能而值得一提，但到目前为止还有没有任何功能可以回答这个问题。如果有人可以扩展它来回答这个问题，那就太好了。

假设我们有以下数据：

x = c(120,105,118,140,142,141,135,152,154,138,125,132,131,120)
plot(x, type="l")

并且我们想对具有可变滚动窗口width 的x 向量应用滚动函数。

set.seed(1)
width = sample(2:4,length(x),TRUE)

在这种特殊情况下，我们将滚动函数自适应于 sample 或 c(2,3,4)。
我们将应用mean 函数，预期结果：

r = f(x, width, FUN = mean)
print(r)
##  [1]       NA       NA 114.3333 120.7500 141.0000 135.2500 139.5000
##  [8] 142.6667 147.0000 146.0000 131.5000 128.5000 131.5000 127.6667
plot(x, type="l")
lines(r, col="red")

任何指标都可以用来产生width 参数作为自适应移动平均线的不同变体，或任何其他函数。

寻求最佳表现。

Answer 1

2018 年 12 月更新

自适应滚动功能的高效实现已在 data.table 最近 - ?froll 手册中的更多信息。此外，已经确定了使用基础 R 的有效替代解决方案（下面的fastama）。不幸的是，Kevin Ushey 的回答没有解决这个问题，因此它不包含在基准测试中。 由于比较微秒毫无意义，因此增加了基准规模。

set.seed(108)
x = rnorm(1e6)
width = rep(seq(from = 100, to = 500, by = 5), length.out=length(x))
microbenchmark(
  zoo=rollapplyr(x, width = width, FUN=mean, fill=NA),
  mapply=base_mapply(x, width=width, FUN=mean, na.rm=T),
  wmapply=wmapply(x, width=width, FUN=mean, na.rm=T),
  ama=ama(x, width, na.rm=T),
  fastama=fastama(x, width),
  frollmean=frollmean(x, width, na.rm=T, adaptive=TRUE),
  frollmean_exact=frollmean(x, width, na.rm=T, adaptive=TRUE, algo="exact"),
  times=1L
)
#Unit: milliseconds
#            expr          min           lq         mean       median           uq          max neval
#             zoo 32371.938248 32371.938248 32371.938248 32371.938248 32371.938248 32371.938248     1
#          mapply 13351.726032 13351.726032 13351.726032 13351.726032 13351.726032 13351.726032     1
#         wmapply 15114.774972 15114.774972 15114.774972 15114.774972 15114.774972 15114.774972     1
#             ama  9780.239091  9780.239091  9780.239091  9780.239091  9780.239091  9780.239091     1
#         fastama   351.618042   351.618042   351.618042   351.618042   351.618042   351.618042     1
#       frollmean     7.708054     7.708054     7.708054     7.708054     7.708054     7.708054     1
# frollmean_exact   194.115012   194.115012   194.115012   194.115012   194.115012   194.115012     1

ama = function(x, n, na.rm=FALSE, fill=NA, nf.rm=FALSE) {
  # more or less the same as previous forloopply
  stopifnot((nx<-length(x))==length(n))
  if (nf.rm) x[!is.finite(x)] = NA_real_
  ans = rep(NA_real_, nx)
  for (i in seq_along(x)) {
    ans[i] = if (i >= n[i])
      mean(x[(i-n[i]+1):i], na.rm=na.rm)
    else as.double(fill)
  }
  ans
}
fastama = function(x, n, na.rm, fill=NA) {
  if (!missing(na.rm)) stop("fast adaptive moving average implemented in R does not handle NAs, input having NAs will result in incorrect answer so not even try to compare to it")
  # fast implementation of adaptive moving average in R, in case of NAs incorrect answer
  stopifnot((nx<-length(x))==length(n))
  cs = cumsum(x)
  ans = rep(NA_real_, nx)
  for (i in seq_along(cs)) {
    ans[i] = if (i == n[i])
      cs[i]/n[i]
    else if (i > n[i])
      (cs[i]-cs[i-n[i]])/n[i]
    else as.double(fill)
  }
  ans
}

---

旧答案：

我选择了 4 个不需要 C++ 的可用解决方案，很容易找到或 google。

# 1. rollapply
library(zoo)
?rollapplyr
# 2. mapply
base_mapply <- function(x, width, FUN, ...){
  FUN <- match.fun(FUN)
  f <- function(i, width, data){
    if(i < width) return(NA_real_)
    return(FUN(data[(i-(width-1)):i], ...))
  }
  mapply(FUN = f, 
         seq_along(x), width,
         MoreArgs = list(data = x))
}
# 3. wmapply - modified version of wapply found: https://rmazing.wordpress.com/2013/04/23/wapply-a-faster-but-less-functional-rollapply-for-vector-setups/
wmapply <- function(x, width, FUN = NULL, ...){
  FUN <- match.fun(FUN)
  SEQ1 <- 1:length(x)
  SEQ1[SEQ1 <  width] <- NA_integer_
  SEQ2 <- lapply(SEQ1, function(i) if(!is.na(i)) (i - (width[i]-1)):i)
  OUT <- lapply(SEQ2, function(i) if(!is.null(i)) FUN(x[i], ...) else NA_real_)
  return(base:::simplify2array(OUT, higher = TRUE))
}
# 4. forloopply - simple loop solution
forloopply <- function(x, width, FUN = NULL, ...){
  FUN <- match.fun(FUN)
  OUT <- numeric()
  for(i in 1:length(x)) {
    if(i < width[i]) next
    OUT[i] <- FUN(x[(i-(width[i]-1)):i], ...)
  }
  return(OUT)
}

以下是prod 函数的时间安排。 mean 函数可能已经在 rollapplyr 内部进行了优化。所有结果均等。

library(microbenchmark)
# 1a. length(x) = 1000, window = 5-20
x <- runif(1000,0.5,1.5)
width <- rep(seq(from = 5, to = 20, by = 5), length(x)/4)
microbenchmark(
  rollapplyr(data = x, width = width, FUN = prod, fill = NA),
  base_mapply(x = x, width = width, FUN = prod, na.rm=T),
  wmapply(x = x, width = width, FUN = prod, na.rm=T),
  forloopply(x = x, width = width, FUN = prod, na.rm=T),
  times=100L
)
Unit: milliseconds
                                                       expr       min        lq    median       uq       max neval
 rollapplyr(data = x, width = width, FUN = prod, fill = NA) 59.690217 60.694364 61.979876 68.55698 153.60445   100
   base_mapply(x = x, width = width, FUN = prod, na.rm = T) 14.372537 14.694266 14.953234 16.00777  99.82199   100
       wmapply(x = x, width = width, FUN = prod, na.rm = T)  9.384938  9.755893  9.872079 10.09932  84.82886   100
    forloopply(x = x, width = width, FUN = prod, na.rm = T) 14.730428 15.062188 15.305059 15.76560 342.44173   100

# 1b. length(x) = 1000, window = 50-200
x <- runif(1000,0.5,1.5)
width <- rep(seq(from = 50, to = 200, by = 50), length(x)/4)
microbenchmark(
  rollapplyr(data = x, width = width, FUN = prod, fill = NA),
  base_mapply(x = x, width = width, FUN = prod, na.rm=T),
  wmapply(x = x, width = width, FUN = prod, na.rm=T),
  forloopply(x = x, width = width, FUN = prod, na.rm=T),
  times=100L
)
Unit: milliseconds
                                                       expr      min       lq   median       uq      max neval
 rollapplyr(data = x, width = width, FUN = prod, fill = NA) 71.99894 74.19434 75.44112 86.44893 281.6237   100
   base_mapply(x = x, width = width, FUN = prod, na.rm = T) 15.67158 16.10320 16.39249 17.20346 103.6211   100
       wmapply(x = x, width = width, FUN = prod, na.rm = T) 10.88882 11.54721 11.75229 12.19790 106.1170   100
    forloopply(x = x, width = width, FUN = prod, na.rm = T) 15.70704 16.06983 16.40393 17.14210 108.5005   100

# 2a. length(x) = 10000, window = 5-20
x <- runif(10000,0.5,1.5)
width <- rep(seq(from = 5, to = 20, by = 5), length(x)/4)
microbenchmark(
  rollapplyr(data = x, width = width, FUN = prod, fill = NA),
  base_mapply(x = x, width = width, FUN = prod, na.rm=T),
  wmapply(x = x, width = width, FUN = prod, na.rm=T),
  forloopply(x = x, width = width, FUN = prod, na.rm=T),
  times=100L
)
Unit: milliseconds
                                                       expr       min       lq   median       uq       max neval
 rollapplyr(data = x, width = width, FUN = prod, fill = NA) 753.87882 781.8789 809.7680 872.8405 1116.7021   100
   base_mapply(x = x, width = width, FUN = prod, na.rm = T) 148.54919 159.9986 231.5387 239.9183  339.7270   100
       wmapply(x = x, width = width, FUN = prod, na.rm = T)  98.42682 105.2641 117.4923 183.4472  245.4577   100
    forloopply(x = x, width = width, FUN = prod, na.rm = T) 533.95641 602.0652 646.7420 672.7483  922.3317   100

# 2b. length(x) = 10000, window = 50-200
x <- runif(10000,0.5,1.5)
width <- rep(seq(from = 50, to = 200, by = 50), length(x)/4)
microbenchmark(
  rollapplyr(data = x, width = width, FUN = prod, fill = NA),
  base_mapply(x = x, width = width, FUN = prod, na.rm=T),
  wmapply(x = x, width = width, FUN = prod, na.rm=T),
  forloopply(x = x, width = width, FUN = prod, na.rm=T),
  times=100L
)
Unit: milliseconds
                                                       expr      min       lq    median        uq       max neval
 rollapplyr(data = x, width = width, FUN = prod, fill = NA) 912.5829 946.2971 1024.7245 1071.5599 1431.5289   100
   base_mapply(x = x, width = width, FUN = prod, na.rm = T) 171.3189 180.6014  260.8817  269.5672  344.4500   100
       wmapply(x = x, width = width, FUN = prod, na.rm = T) 123.1964 131.1663  204.6064  221.1004  484.3636   100
    forloopply(x = x, width = width, FUN = prod, na.rm = T) 561.2993 696.5583  800.9197  959.6298 1273.5350   100

Answer 2

作为参考，如果您只有一个窗口长度可以“翻转”，您绝对应该查看RcppRoll：

library(RcppRoll) ## install.packages("RcppRoll")
library(microbenchmark)
x <- runif(1E5)
all.equal( rollapplyr(x, 10, FUN=prod), roll_prod(x, 10) )
microbenchmark( times=5,
  rollapplyr(x, 10, FUN=prod),
  roll_prod(x, 10)
)

给我

> library(RcppRoll)
> library(microbenchmark)
> x <- runif(1E5)
> all.equal( rollapplyr(x, 10, FUN=prod), roll_prod(x, 10) )
[1] TRUE
> microbenchmark( times=5,
+   zoo=rollapplyr(x, 10, FUN=prod),
+   RcppRoll=roll_prod(x, 10)
+ )
Unit: milliseconds
     expr        min         lq     median         uq         max neval
      zoo 924.894069 968.467299 997.134932 1029.10883 1079.613569     5
 RcppRoll   1.509155   1.553062   1.760739    1.90061    1.944999     5

它有点快 ;) 并且包足够灵活，用户可以定义和使用他们自己的滚动函数（使用 C++）。将来我可能会扩展该包以允许多个窗口宽度，但我确信要正确处理会很棘手。

如果您想自己定义prod，您可以这样做——RcppRoll 允许您定义自己的 C++ 函数来传递并根据需要生成“滚动”函数。 rollit 提供了一个更好的界面，而 rollit_raw 只是让您自己编写一个 C++ 函数，有点像您对 Rcpp::cppFunction 所做的那样。哲学是，您应该只需要表达您希望在特定窗口上执行的计算，RcppRoll 可以负责迭代某些大小的窗口。

library(RcppRoll)
library(microbenchmark)
x <- runif(1E5)
my_rolling_prod <- rollit(combine="*")
my_rolling_prod2 <- rollit_raw("
double output = 1;
for (int i=0; i < n; ++i) {
  output *= X(i);
}
return output;
")
all.equal( roll_prod(x, 10), my_rolling_prod(x, 10) )
all.equal( roll_prod(x, 10), my_rolling_prod2(x, 10) )
microbenchmark( times=5,
  rollapplyr(x, 10, FUN=prod),
  roll_prod(x, 10),
  my_rolling_prod(x, 10),
  my_rolling_prod2(x, 10)
)

给我

> library(RcppRoll)
> library(microbenchmark)
> # 1a. length(x) = 1000, window = 5-20
> x <- runif(1E5)
> my_rolling_prod <- rollit(combine="*")
C++ source file written to /var/folders/m7/_xnnz_b53kjgggkb1drc1f8c0000gn/T//RtmpcFMJEV/file80263aa7cca2.cpp .
Compiling...
Done!
> my_rolling_prod2 <- rollit_raw("
+ double output = 1;
+ for (int i=0; i < n; ++i) {
+   output *= X(i);
+ }
+ return output;
+ ")
C++ source file written to /var/folders/m7/_xnnz_b53kjgggkb1drc1f8c0000gn/T//RtmpcFMJEV/file802673777da2.cpp .
Compiling...
Done!
> all.equal( roll_prod(x, 10), my_rolling_prod(x, 10) )
[1] TRUE
> all.equal( roll_prod(x, 10), my_rolling_prod2(x, 10) )
[1] TRUE
> microbenchmark(
+   rollapplyr(x, 10, FUN=prod),
+   roll_prod(x, 10),
+   my_rolling_prod(x, 10),
+   my_rolling_prod2(x, 10)
+ )

> microbenchmark( times=5,
+   rollapplyr(x, 10, FUN=prod),
+   roll_prod(x, 10),
+   my_rolling_prod(x, 10),
+   my_rolling_prod2(x, 10)
+ )
Unit: microseconds
                          expr        min          lq      median          uq         max neval
 rollapplyr(x, 10, FUN = prod) 979710.368 1115931.323 1117375.922 1120085.250 1149117.854     5
              roll_prod(x, 10)   1504.377    1635.749    1638.943    1815.344    2053.997     5
        my_rolling_prod(x, 10)   1507.687    1572.046    1648.031    2103.355    7192.493     5
       my_rolling_prod2(x, 10)    774.381     786.750     884.951    1052.508    1434.660     5

真的，只要您能够通过rollit 接口或通过rollit_raw 传递的C++ 函数（其接口有点死板，但仍然有效），你的状态很好。

Answer 3

不知何故，人们错过了 base R（统计包）中的超快runmed()。据我了解原始问题，它不是自适应的，但是对于滚动中位数，它很快！在这里与来自 RcppRoll 的 roll_median() 进行比较。

> microbenchmark(
+   runmed(x = x, k = 3),
+   roll_median(x, 3),
+   times=1000L
+ )
Unit: microseconds
                 expr     min      lq      mean  median      uq     max neval
 runmed(x = x, k = 3)  41.053  44.854  47.60973  46.755  49.795 117.838  1000
    roll_median(x, 3) 101.872 105.293 108.72840 107.574 111.375 178.657  1000