Функция адаптивного поворотного окна - верхняя производительность в R

Question

Я ищу некоторые улучшения производительности с точки зрения функций качения/скольжения в R. Это довольно распространенная задача, которая может быть использована в любом упорядоченном наборе данных наблюдений. Я хотел бы поделиться некоторыми своими выводами, возможно, кто-то сможет предоставить обратную связь, чтобы сделать это еще быстрее.
Важным примечанием является то, что я фокусируюсь на корпусе align="right" и прокатного окна width в виде вектора (такой же длины, что и наш вектор наблюдения). В случае, если у нас есть width в качестве скаляра, есть уже очень хорошо разработанные функции в пакетах zoo и TTR, которые очень сложно обыграть, поскольку некоторые из них даже используют Fortran (но все же пользовательские FUN могут быть быстрее, используя приведенные ниже wapply) ,
RcppRoll Пакет стоит упомянуть из-за его высокой производительности, но до сих пор нет функции, которая отвечает на этот вопрос. Было бы здорово, если бы кто-то мог продлить его, чтобы ответить на вопрос.

Рассмотрим мы имеем следующие данные:

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

И мы хотим применить функцию прокатки над x вектора с окном переменной качению width.

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

Я выбрал 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 или с помощью функции C++ прошел через rollit_raw (интерфейс которого немного жесткий, но функциональный), вы в хорошей форме.

Answer 3

Как-то люди пропустили сверхбыстрые runmed() в базе R (пакет статистики). Это не адаптивный, насколько я понимаю исходный вопрос, но для скользящей медианы это БЫСТРО! Сравнивая здесь с roll_median() от RcppRoll.

> 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