Ускорить цикл работы в R

193

У меня большая проблема с производительностью в R. Я написал функцию, которая перебирает data.frameобъект. Он просто добавляет новый столбец в a data.frameи что-то накапливает. (простая операция). data.frameИмеет примерно 850K строк. Мой компьютер все еще работает (около 10 часов), и я понятия не имею о времени выполнения.

dayloop2 <- function(temp){
    for (i in 1:nrow(temp)){    
        temp[i,10] <- i
        if (i > 1) {             
            if ((temp[i,6] == temp[i-1,6]) & (temp[i,3] == temp[i-1,3])) { 
                temp[i,10] <- temp[i,9] + temp[i-1,10]                    
            } else {
                temp[i,10] <- temp[i,9]                                    
            }
        } else {
            temp[i,10] <- temp[i,9]
        }
    }
    names(temp)[names(temp) == "V10"] <- "Kumm."
    return(temp)
}

Есть идеи как ускорить эту операцию?

Кей
источник

Ответы:

435

Самая большая проблема и корень неэффективности заключается в индексации data.frame, я имею в виду все эти строки, где вы используете temp[,].
Старайтесь избегать этого как можно больше. Я взял твою функцию, поменяй индексацию и вот version_A

dayloop2_A <- function(temp){
    res <- numeric(nrow(temp))
    for (i in 1:nrow(temp)){    
        res[i] <- i
        if (i > 1) {             
            if ((temp[i,6] == temp[i-1,6]) & (temp[i,3] == temp[i-1,3])) { 
                res[i] <- temp[i,9] + res[i-1]                   
            } else {
                res[i] <- temp[i,9]                                    
            }
        } else {
            res[i] <- temp[i,9]
        }
    }
    temp$`Kumm.` <- res
    return(temp)
}

Как видите, я создаю вектор, resкоторый собирает результаты. В конце я добавляю это вdata.frame и мне не нужно связываться с именами. Так как же лучше?

Я запускаю каждую функцию data.frameс nrowот 1000 до 10000 на 1000 и измеряю время сsystem.time

X <- as.data.frame(matrix(sample(1:10, n*9, TRUE), n, 9))
system.time(dayloop2(X))

Результат

производительность

Вы можете видеть, что ваша версия экспоненциально зависит от nrow(X). Модифицированная версия имеет линейную зависимость и простуюlm модель предсказывает, что для 850 000 строк вычисление занимает 6 минут и 10 секунд.

Сила векторизации

Как утверждает Шейн и Калимо в своих ответах, векторизация является ключом к повышению производительности. Из вашего кода вы можете выйти за пределы цикла:

  • кондиционирование
  • инициализация результатов (которые есть temp[i,9])

Это приводит к этому коду

dayloop2_B <- function(temp){
    cond <- c(FALSE, (temp[-nrow(temp),6] == temp[-1,6]) & (temp[-nrow(temp),3] == temp[-1,3]))
    res <- temp[,9]
    for (i in 1:nrow(temp)) {
        if (cond[i]) res[i] <- temp[i,9] + res[i-1]
    }
    temp$`Kumm.` <- res
    return(temp)
}

Сравните результат для этих функций, на этот раз nrowс 10 000 до 100 000 на 10 000.

производительность

Тюнинг настроенного

Другой твик заключается в изменении индексации цикла temp[i,9]на res[i](что в точности повторяется в итерации i-го цикла). Это опять разница между индексированием вектора и индексированием a data.frame.
Второе: когда вы смотрите на цикл, вы видите, что нет необходимости циклически повторять все i, а только те, которые соответствуют условию.
Итак, поехали

dayloop2_D <- function(temp){
    cond <- c(FALSE, (temp[-nrow(temp),6] == temp[-1,6]) & (temp[-nrow(temp),3] == temp[-1,3]))
    res <- temp[,9]
    for (i in (1:nrow(temp))[cond]) {
        res[i] <- res[i] + res[i-1]
    }
    temp$`Kumm.` <- res
    return(temp)
}

Производительность, которую вы получаете, сильно зависит от структуры данных. Точно - на процент TRUEзначений в состоянии. Для моих смоделированных данных требуется время вычисления на 850 000 строк ниже одной секунды.

производительность

Если вы хотите, вы можете пойти дальше, я вижу, по крайней мере, две вещи, которые можно сделать:

  • Напиши C код, чтобы сделать условное cumsum
  • если вы знаете, что в ваших данных максимальная последовательность не велика, вы можете изменить цикл на векторизованное время, что-то вроде

    while (any(cond)) {
        indx <- c(FALSE, cond[-1] & !cond[-n])
        res[indx] <- res[indx] + res[which(indx)-1]
        cond[indx] <- FALSE
    }

Код, используемый для моделирования и рисунков, доступен на GitHub .

Marek
источник
2
Поскольку я не могу найти способ частно спросить Марека, как генерировались эти графики?
carbontwelve
@carbontwelve Вы спрашиваете о данных или графиках? Участки были сделаны с решеткой пакета. Если у меня есть время, я размещаю код где-то в сети и сообщаю вам.
Марек
@carbontwelve Упс, я был неправ :) Это стандартные графики (из базы R).
Марек
@ Грегор К сожалению нет. Это накопительное, поэтому вы не можете векторизовать его. Простой пример: res = c(1,2,3,4)и condэто все TRUE, то конечный результат должен быть: 1, 3(причина 1+2), 6(причина второй теперь 3, и третий 3и), 10( 6+4). Делая простое суммирование вы получили 1, 3, 5, 7.
Марек
Ах, я должен был продумать это более тщательно. Спасибо, что показали мне ошибку.
Грегор Томас
132

Общие стратегии для ускорения кода R

Во-первых, выяснить, где медленная часть на самом деле. Нет необходимости оптимизировать код, который не работает медленно. Для небольшого количества кода, просто продумывая его, может работать. Если это не помогло, RProf и подобные инструменты профилирования могут быть полезны.

Как только вы обнаружите узкое место, подумайте о более эффективных алгоритмах для выполнения того, что вы хотите. Расчеты должны выполняться только один раз, если это возможно, поэтому:

Использование более эффективных функций может привести к умеренному или большому увеличению скорости. Например, paste0дает небольшой прирост эффективности, но .colSums()и его родственники дают несколько более выраженный прирост. meanявляется особенно медленным .

Тогда вы можете избежать некоторых наиболее распространенных проблем :

  • cbind замедлит вас очень быстро.
  • Инициализируйте ваши структуры данных, затем заполните их, вместо того, чтобы каждый раз расширять их .
  • Даже с предварительным распределением вы можете переключиться на подход с передачей по ссылке, а не с передачей по значению, но это может не стоить хлопот.
  • Взгляните на R Inferno, чтобы избежать других подводных камней.

Попробуйте улучшить векторизацию , которая часто, но не всегда, помогает. В связи с этим, по своей сути векторизованных команды типа ifelse, diffи тому подобное обеспечит больше , чем улучшениеapply семейство команд (которые обеспечивают практически полное повышение скорости по сравнению с хорошо написанным циклом).

Вы также можете попытаться предоставить дополнительную информацию для функций R . Например, используйте vapplyвместоsapply , и укажите colClassesпри чтении в текстовых данных . Прирост скорости будет изменяться в зависимости от того, сколько угадывания вы устраняете.

Далее рассмотрим оптимизированные пакеты . data.tableПакет может значительно увеличить скорость, где его использование возможно, при манипулировании данными и при чтении больших объемов данных ( fread).

Затем попробуйте увеличить скорость за счет более эффективных способов вызова R :

  • Скомпилируйте ваш R-скрипт. Или использовать Raи jitпакеты концерта для точно в момент компиляции (Dirk есть пример в данной презентации ).
  • Убедитесь, что вы используете оптимизированный BLAS. Они обеспечивают повсеместный прирост скорости. Честно говоря, это позор, что R не использует автоматически самую эффективную библиотеку при установке. Надеюсь, что Revolution R внесет свой вклад в общее сообщество.
  • Рэдфорд Нил провел ряд оптимизаций, некоторые из которых были перенесены в R Core, а многие другие были разветвлены в pqR .

И, наконец, если все вышеперечисленное все еще не дает вам такой скорости, как вам нужно, вам может потребоваться перейти на более быстрый язык для медленного фрагмента кода . Комбинация Rcppи inlineздесь делает замену только самой медленной части алгоритма кодом C ++ особенно легкой. Вот, например, моя первая попытка сделать это , и она поражает даже высоко оптимизированные R-решения.

Если после всего этого у вас все еще есть проблемы, вам просто нужно больше вычислительной мощности. Посмотрите на распараллеливание ( http://cran.r-project.org/web/views/HighPerformanceComputing.html ) или даже решения на основе GPU ( gpu-tools).

Ссылки на другие руководства

Ari B. Friedman
источник
36

Если вы используете forциклы, вы, скорее всего, кодируете R, как если бы это был C, Java или что-то еще. Правильно векторизованный R-код очень быстр.

Возьмем, к примеру, эти два простых фрагмента кода для генерации списка из 10 000 целых чисел в последовательности:

Первый пример кода - это то, как можно было бы закодировать цикл, используя традиционную парадигму кодирования. Это займет 28 секунд, чтобы завершить

system.time({
    a <- NULL
    for(i in 1:1e5)a[i] <- i
})
   user  system elapsed 
  28.36    0.07   28.61 

Вы можете получить улучшение почти в 100 раз простым действием предварительного выделения памяти:

system.time({
    a <- rep(1, 1e5)
    for(i in 1:1e5)a[i] <- i
})

   user  system elapsed 
   0.30    0.00    0.29 

Но используя базовую векторную операцию R с использованием оператора двоеточия, :эта операция практически мгновенная:

system.time(a <- 1:1e5)

   user  system elapsed 
      0       0       0 
Andrie
источник
+1 хотя я бы расценил ваш второй пример как неубедительный, так как a[i]не меняется. Но system.time({a <- NULL; for(i in 1:1e5){a[i] <- 2*i} }); system.time({a <- 1:1e5; for(i in 1:1e5){a[i] <- 2*i} }); system.time({a <- NULL; a <- 2*(1:1e5)})имеет аналогичный результат.
Генри
@ Генри, честный комментарий, но, как вы заметили, результаты такие же. Я изменил пример для инициализации a rep(1, 1e5)- время идентично.
Андри
17

Это можно сделать намного быстрее, пропустив циклы с помощью индексов или вложенных ifelse()операторов.

idx <- 1:nrow(temp)
temp[,10] <- idx
idx1 <- c(FALSE, (temp[-nrow(temp),6] == temp[-1,6]) & (temp[-nrow(temp),3] == temp[-1,3]))
temp[idx1,10] <- temp[idx1,9] + temp[which(idx1)-1,10] 
temp[!idx1,10] <- temp[!idx1,9]    
temp[1,10] <- temp[1,9]
names(temp)[names(temp) == "V10"] <- "Kumm."
Шейн
источник
Спасибо за ответ. Я пытаюсь понять ваши заявления. Строка 4: «temp [idx1,10] <- temp [idx1,9] + temp [which (idx1) -1,10]» вызвала ошибку, поскольку длина более длинного объекта не кратна длине короче объект. «temp [idx1,9] = num [1: 11496]» и «temp [which (idx1) -1,10] = int [1: 11494]», поэтому отсутствуют 2 строки.
Кей
Если вы предоставите образец данных (используйте dput () с несколькими строками), я исправлю это для вас. Из-за which () - 1 бит, индексы неравны. Но вы должны увидеть, как это работает отсюда: нет необходимости в зацикливании или применении; просто используйте векторизованные функции.
Шейн
1
Вот Это Да! Я только что изменил вложенный функциональный блок if..else и mapply на вложенную функцию ifelse и получил ускорение в 200 раз!
Джеймс
Ваш общий совет верен, но в коде вы пропустили тот факт, что i-th значение зависит от i-1-th, поэтому они не могут быть установлены так, как вы это делаете (используя which()-1).
Марек
8

Мне не нравится переписывать код ... Также, конечно, ifelse и lapply - лучшие варианты, но иногда это трудно сделать.

Часто я использую data.frames, как можно использовать списки, такие как df$var[i]

Вот вымышленный пример:

nrow=function(x){ ##required as I use nrow at times.
  if(class(x)=='list') {
    length(x[[names(x)[1]]])
  }else{
    base::nrow(x)
  }
}

system.time({
  d=data.frame(seq=1:10000,r=rnorm(10000))
  d$foo=d$r
  d$seq=1:5
  mark=NA
  for(i in 1:nrow(d)){
    if(d$seq[i]==1) mark=d$r[i]
    d$foo[i]=mark
  }
})

system.time({
  d=data.frame(seq=1:10000,r=rnorm(10000))
  d$foo=d$r
  d$seq=1:5
  d=as.list(d) #become a list
  mark=NA
  for(i in 1:nrow(d)){
    if(d$seq[i]==1) mark=d$r[i]
    d$foo[i]=mark
  }
  d=as.data.frame(d) #revert back to data.frame
})

версия data.frame:

   user  system elapsed 
   0.53    0.00    0.53

версия списка:

   user  system elapsed 
   0.04    0.00    0.03 

В 17 раз быстрее использовать список векторов, чем в data.frame.

Любые комментарии о том, почему внутренне data.frames так медленно в этом отношении? Казалось бы, они работают как списки ...

Для еще более быстрого кода сделайте это class(d)='list'вместо d=as.list(d)иclass(d)='data.frame'

system.time({
  d=data.frame(seq=1:10000,r=rnorm(10000))
  d$foo=d$r
  d$seq=1:5
  class(d)='list'
  mark=NA
  for(i in 1:nrow(d)){
    if(d$seq[i]==1) mark=d$r[i]
    d$foo[i]=mark
  }
  class(d)='data.frame'
})
head(d)
Крис
источник
1
Вероятно, это связано с накладными расходами [<-.data.frame, которые как-то вызываются, когда вы это делаете, d$foo[i] = markи могут в конечном итоге создать новую копию вектора, возможно, всего data.frame в каждой <-модификации. Это сделало бы интересный вопрос о SO.
Фрэнк
2
@Frank Он (i) должен убедиться, что измененный объект по-прежнему является действительным data.frame, и (ii) afaik делает как минимум одну копию, возможно, более одной. Известно, что подчинение подкадров данных происходит медленно, и если вы посмотрите на длинный исходный код, это не удивительно.
Роланд
@Frank, @Roland: df$var[i]Нотация проходит через ту же [<-.data.frameфункцию? Я заметил, что это действительно довольно долго. Если нет, то какую функцию он использует?
Крис
@ Крис, я думаю, d$foo[i]=markпереводится примерно d <- `$<-`(d, 'foo', `[<-`(d$foo, i, mark)), но с некоторым использованием временных переменных.
Тим Гудман
7

Как Ари уже упоминался в конце своего ответа, Rcppи inlineпакеты делают его невероятно легко сделать вещи быстро. Как пример, попробуйте этот inlineкод (предупреждение: не проверено):

body <- 'Rcpp::NumericMatrix nm(temp);
         int nrtemp = Rccp::as<int>(nrt);
         for (int i = 0; i < nrtemp; ++i) {
             temp(i, 9) = i
             if (i > 1) {
                 if ((temp(i, 5) == temp(i - 1, 5) && temp(i, 2) == temp(i - 1, 2) {
                     temp(i, 9) = temp(i, 8) + temp(i - 1, 9)
                 } else {
                     temp(i, 9) = temp(i, 8)
                 }
             } else {
                 temp(i, 9) = temp(i, 8)
             }
         return Rcpp::wrap(nm);
        '

settings <- getPlugin("Rcpp")
# settings$env$PKG_CXXFLAGS <- paste("-I", getwd(), sep="") if you want to inc files in wd
dayloop <- cxxfunction(signature(nrt="numeric", temp="numeric"), body-body,
    plugin="Rcpp", settings=settings, cppargs="-I/usr/include")

dayloop2 <- function(temp) {
    # extract a numeric matrix from temp, put it in tmp
    nc <- ncol(temp)
    nm <- dayloop(nc, temp)
    names(temp)[names(temp) == "V10"] <- "Kumm."
    return(temp)
}

Есть аналогичная процедура для #includeвещей, где вы просто передаете параметр

inc <- '#include <header.h>

к функции, как include=inc . Что действительно круто в этом, так это то, что он выполняет всю компоновку и компиляцию за вас, поэтому прототипирование действительно быстрое.

Отказ от ответственности: я не совсем уверен, что класс tmp должен быть числовым, а не числовой матрицей или чем-то еще. Но я в основном уверен.

Изменить: если вам все еще нужно больше скорости после этого, OpenMP является средством параллелизации хорошо для C++. Я не пробовал использовать его inline, но он должен работать. Идея была бы, в случае nядер, есть цикл итерация kбыть осуществлена путем k % n. Подходящее введение можно найти в Matloff's Art of R Programming , доступной здесь , в главе 16, « Использование C» .

jclancy
источник
3

Ответы здесь великолепны. Один незначительный аспект, который не был затронут, заключается в том, что вопрос гласит: « Мой ПК все еще работает (около 10 часов), и я понятия не имею о времени выполнения ». Я всегда вставляю следующий код в циклы при разработке, чтобы понять, как изменения влияют на скорость, а также для отслеживания того, сколько времени потребуется для завершения.

dayloop2 <- function(temp){
  for (i in 1:nrow(temp)){
    cat(round(i/nrow(temp)*100,2),"%    \r") # prints the percentage complete in realtime.
    # do stuff
  }
  return(blah)
}

Работает с lapply также.

dayloop2 <- function(temp){
  temp <- lapply(1:nrow(temp), function(i) {
    cat(round(i/nrow(temp)*100,2),"%    \r")
    #do stuff
  })
  return(temp)
}

Если функция в цикле довольно быстрая, но число циклов велико, рассмотрите возможность периодической печати, так как печать на самой консоли имеет накладные расходы. например

dayloop2 <- function(temp){
  for (i in 1:nrow(temp)){
    if(i %% 100 == 0) cat(round(i/nrow(temp)*100,2),"%    \r") # prints every 100 times through the loop
    # do stuff
  }
  return(temp)
}
новобранец
источник
Похожая опция, выведите дробь i / n. У меня всегда есть что-то вроде, cat(sprintf("\nNow running... %40s, %s/%s \n", nm[i], i, n))так как я обычно перебираю именованные вещи (с именами в nm).
Фрэнк
2

В R вы часто можете ускорить циклическую обработку, используя applyсемейные функции (в вашем случае, вероятно, так и будет replicate). Посмотрите наplyr пакет, который предоставляет индикаторы выполнения.

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

temp[1:nrow(temp), 10] <- temp[1:nrow(temp), 9] + temp[0:(nrow(temp)-1), 10]

Это будет намного быстрее, и тогда вы сможете отфильтровать строки с вашим условием:

cond.i <- (temp[i, 6] == temp[i-1, 6]) & (temp[i, 3] == temp[i-1, 3])
temp[cond.i, 10] <- temp[cond.i, 9]

Векторизованная арифметика требует больше времени и размышлений о проблеме, но иногда вы можете сэкономить несколько порядков времени выполнения.

Calimo
источник
14
вы заметили, что векторные функции будут быстрее, чем циклы или apply (), но это неправда, что apply () работает быстрее, чем циклы. Во многих случаях apply () просто абстрагирует цикл от пользователя, но все еще делает цикл. Смотрите предыдущий вопрос: stackoverflow.com/questions/2275896/…
JD Long
0

Обработка с помощью data.tableявляется жизнеспособным вариантом:

n <- 1000000
df <- as.data.frame(matrix(sample(1:10, n*9, TRUE), n, 9))
colnames(df) <- paste("col", 1:9, sep = "")

library(data.table)

dayloop2.dt <- function(df) {
  dt <- data.table(df)
  dt[, Kumm. := {
    res <- .I;
    ifelse (res > 1,             
      ifelse ((col6 == shift(col6, fill = 0)) & (col3 == shift(col3, fill = 0)) , 
        res <- col9 + shift(res)                   
      , # else
        res <- col9                                 
      )
     , # else
      res <- col9
    )
  }
  ,]
  res <- data.frame(dt)
  return (res)
}

res <- dayloop2.dt(df)

m <- microbenchmark(dayloop2.dt(df), times = 100)
#Unit: milliseconds
#       expr      min        lq     mean   median       uq      max neval
#dayloop2.dt(df) 436.4467 441.02076 578.7126 503.9874 575.9534 966.1042    10

Если вы игнорируете возможные выгоды от фильтрации условий, это очень быстро. Очевидно, что если вы можете сделать расчет на подмножестве данных, это помогает.

Булат
источник
2
Почему вы повторяете предложение использовать data.table? Это уже было сделано несколько раз в предыдущих ответах.
IRTFM