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

У меня есть функция rollapply, которая делает что-то очень простое, однако более миллиона точек данных эта простая функция работает довольно медленно. Я хотел бы знать, можно ли предоставить информацию для 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. Каков следующий логический шаг для дальнейшей оптимизации?


r optimization rollapply
person channon    schedule 07.02.2019    source источник
comment
Если посты, на которые вы ссылаетесь (предлагающие циклы lapply), написаны несколько лет назад, то они могут говорить об этом, потому что циклы for раньше работали медленнее. С тех пор это было исправлено (R-3.0?), поэтому предположение о том, что циклы for всегда медленнее, возможно, опрометчиво. Бывают случаи, когда имеет смысл не использовать цикл for (когда вы не хотите побочных эффектов и хотите зафиксировать результат каждой итерации в цикле), и в этом случае один из функции *apply могут иметь больше смысла.   -  person r2evans    schedule 07.02.2019
comment
При использовании параллельных функций, таких как mclapply, меня больше всего беспокоит объем данных, передаваемых между ядрами. Например, если вы переходите от mc.cores=8 к (скажем) =4, удваивается ли время выполнения? А как насчет =12 (если у вас есть больше)? Если вы не видите явного повышения производительности при увеличении количества ядер, возможно, вы испытываете отставание из-за копирования данных (в зависимости от вашей параллельной настройки). Так что, возможно, ускорение связано не только с вашим UDF (который настолько эффективен, насколько вы можете получить в необработанном виде).   -  person r2evans    schedule 07.02.2019
comment
Поскольку я новичок в r, я не знал об изменениях производительности в циклах for. Я посмотрю, как окно, отсортированное по циклу for, работает с w.r.t. роллприменить. 8 ядер был выбран осознанно, т.к. 8 роллов (и 8 ядер). Я собираюсь перейти на новую машину с 36 ядрами и, возможно, тогда я смогу исследовать эффекты закона Амдалса. После того, как я разберусь с удаленным r в ess. Спасибо за совет   -  person channon    schedule 07.02.2019
comment
Remote R в ESS: зависит от того, что вам нужно. У меня на удаленной системе установлен emacs/ess, поэтому ssh и tmux (screen) достаточно; когда мне нужно увидеть графики, я использую rmote. Я читал о ESS-Remote, но никогда не пробовал.   -  person r2evans    schedule 07.02.2019

Ответы (2)


arrow_upward
3
arrow_downward

Не уверен, что это будет применяться в среде mclapply, но вы можете немного ускориться, используя оптимизированную функцию rollmax zoo. Поскольку у них нет дополнения 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 каждый раз.)

person r2evans    schedule 07.02.2019
comment
Я не знал rollmax, я не могу себе представить, что было бы слишком много работы, чтобы реализовать оптимизированный rollmin. - person channon; 07.02.2019
comment
Это просто: rollmin <- function(x, ...) -zoo::rollmax(-x, ...). - person r2evans; 07.02.2019
comment
время работы снижено с 124 до 65 с... Браво! - person channon; 07.02.2019

arrow_upward
3
arrow_downward

Если вы действительно хотите максимально снизить производительность, используйте 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
)
person thc    schedule 07.02.2019
comment
Чтобы действительно снизить производительность настолько, насколько это возможно, поможет ли удаление привязок std:: и выполнение грубой силы на низком уровне? (Это может показаться смешным... ваше ускорение здесь впечатляет.) - person r2evans; 08.02.2019
comment
Я не гуру С++, так что это было не столько предложение, сколько еще один вопрос, @thc - person r2evans; 08.02.2019
comment
Нет, std означает стандартную библиотеку С++. Предполагается, что он очень хорошо оптимизирован, поэтому используйте его, когда сможете;) - person thc; 08.02.2019
comment
Я был достаточно осведомлен о том, что это было, но не знал, насколько это хорошо в общем смысле, так что спасибо за это утешение. Мой опыт работы с RcppSugar показывает, что это хорошо и делает код более удобным для сопровождения и чтения, но имеет небольшое снижение производительности, поэтому я спроецировал это ухудшение и на std::. Спасибо! - person r2evans; 08.02.2019