Методы ускорения роллаплинга с пользовательской функцией (r)

Question

У меня есть функция rollapply, которая делает что-то очень простое, однако более миллиона точек данных эта простая функция довольно медленная.Я хотел бы знать, возможно ли предоставить информацию для выполнения, как сделать следующий переход, а не определять саму функцию.

Конкретно, я выполняю скользящее окно для базового статистического обнаружения аномалий.

Функция применения кувырка:

minmax <- function(x) { max(x) - min(x) }

, вызываемая:

mclapply(data[,eval(vars),with=F], 
         function(x) rollapply(x,width=winSize,FUN=minmax,fill=NA),
         mc.cores=8)

Где data - таблица данных из 8 столбцов, а winsize - 300

Этот вызов занимает около 2 минут на 8 ядрах.Это одно из главных узких мест в вычислительной технике.Однако я могу представить, что мы можем сохранить их отсортированными (по значению и индексу), а затем делать сравнения Olog (n) при каждом скольжении.

Однако я часто вижу сообщения, предлагающие отойти от циклов for и использовать lapply.Каков следующий логический шаг для дальнейшей оптимизации?

Answer 1

Если вы действительно хотите снизить производительность, используйте Rcpp.Пользовательские циклы - отличный пример использования C ++, особенно когда ваша функция довольно проста.

Первые результаты, затем код:

microbenchmark::microbenchmark(
  minmax = zoo::rollapply(aa, width=100, FUN=minmax, fill=NA),
  dblmax = zoo::rollmax(aa, k=100, fill=NA) + zoo::rollmax(-aa, k=100, fill=NA),
  cminmax = crollapply(aa, width=width), times = 10
)
    Unit: milliseconds
    expr       min         lq       mean    median         uq        max neval cld
  minmax 154.04630 162.728871 188.198416 173.13427 200.928005 298.568673    10   c
  dblmax  37.38127  38.541603  44.818505  41.42796  50.001888  61.024250    10  b 
 cminmax   2.31766   2.363676   2.406835   2.39237   2.438109   2.512162    10 a  

C ++ / Rcpp код:

#include <Rcpp.h>
#include <algorithm>
using namespace Rcpp;

// [[Rcpp::export]]
std::vector<double> crollapply(std::vector<double> aa, int width) {
  if(width > aa.size()) throw exception("width too large :(");
  int start_offset = (width-1) / 2;
  int back_offset = width / 2;
  std::vector<double> results(aa.size());
  int i=0;
  for(; i < start_offset; i++) {
    results[i] = NA_REAL;
  }
  for(; i < results.size() - back_offset; i++) {
    double min = *std::min_element(&aa[i - start_offset], &aa[i + back_offset + 1]);
    double max = *std::max_element(&aa[i - start_offset], &aa[i + back_offset + 1]);
    results[i] = max - min;
  }
  for(; i < results.size(); i++) {
    results[i] = NA_REAL;
  }
  return results;
}

R код:

library(dplyr)
library(zoo)
library(microbenchmark)
library(Rcpp)

sourceCpp("~/Desktop/temp.cpp")

minmax <- function(x) max(x) - min(x)
aa <- runif(1e4)
width <- 100
x1 <- zoo::rollapply(aa, width=width, FUN=minmax, fill=NA)
x3 <- crollapply(aa, width=width)
identical(x1,x3)

width <- 101
x1 <- zoo::rollapply(aa, width=width, FUN=minmax, fill=NA)
x3 <- crollapply(aa, width=width)
identical(x1,x3)

microbenchmark::microbenchmark(
  minmax = zoo::rollapply(aa, width=100, FUN=minmax, fill=NA),
  dblmax = zoo::rollmax(aa, k=100, fill=NA) + zoo::rollmax(-aa, k=100, fill=NA),
  cminmax = crollapply(aa, width=width), times = 10
)

Answer 2

Не уверен, что / как это будет применяться в среде mclapply, но вы можете немного ускориться, используя оптимизированную функцию zoo rollmax.Поскольку они не имеют дополняющего rollmin, вам нужно будет адаптироваться.

minmax <- function(x) max(x) - min(x)
aa <- runif(1e4)
identical(
  zoo::rollapply(aa, width=100, FUN=minmax, fill=NA),
  zoo::rollmax(aa, k=100, fill=NA) + zoo::rollmax(-aa, k=100, fill=NA)
)
# [1] TRUE

microbenchmark::microbenchmark(
  minmax = zoo::rollapply(aa, width=100, FUN=minmax, fill=NA),
  dblmax = zoo::rollmax(aa, k=100, fill=NA) + zoo::rollmax(-aa, k=100, fill=NA)
)
# Unit: milliseconds
#    expr     min      lq     mean   median      uq      max neval
#  minmax 70.7426 76.0469 84.81481 77.99565 81.8047 148.8431   100
#  dblmax 15.6755 17.4501 19.09820 17.93665 18.8650  52.4849   100

(Улучшение будет зависеть от размера окна, поэтому ваши результаты могут отличаться, но я думаю, что использование оптимизированногофункция zoo::rollmax почти всегда превосходит вызов UDF каждый раз.)