加快R中的循环操作速度

Question

我在R中遇到了一个严重的性能问题。我编写了一个函数，该函数对data.frame对象进行迭代操作。它只是向data.frame添加一个新列并累加一些东西（简单操作）。该data.frame大约有85万行。我的电脑仍在运行（现在已经约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)
}

有什么想法可以加速这个操作吗？

- Kay

---

在测试函数时，考虑添加类似于if(i%%1000) {print(i)}的内容，以获取大致的运行时间。 - David

Answer 1

最大的问题和无效性的根源是数据框的索引，我的意思是你使用 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 的函数运行测试，并使用 system.time 测量从 1,000 到 10,000 的 nrow，间隔为 1,000。

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

结果是

您可以看到，您的版本与nrow(X)成指数关系。修改后的版本具有线性关系，并且简单的lm模型预测，对于850,000行计算需要6分10秒。

向量化的优势

正如Shane和Calimo在他们的答案中所述，向量化是提高性能的关键。 从您的代码中，您可以移出循环:

• 条件语句
• 结果的初始化（即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个循环迭代中完全相同）中进行更改。 这再次是向量索引和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 代码来执行条件累加。

• 如果您知道您的数据中最大序列不是很大，那么您可以将循环改为向量化 while 循环，类似这样：

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

---

用于模拟和图表的代码可以在 GitHub 上获取。

- Marek

---

由于我找不到私下询问Marek的方法，那些图表是如何生成的？ - carbontwelve

@carbontwelve 你是在询问数据还是图表？这些图表是使用lattice包制作的。如果我有时间，我会将代码放在网上并通知您。 - Marek

@carbontwelve 哎呀，我错了 :) 这些是标准图（来自基本的R）。 - Marek

@Gregor 很遗憾，不能进行向量化。简单的例子：res = c(1,2,3,4)，cond 全部为 TRUE，那么最终结果应该是：1，3（因为 1+2），6（因为第二个现在是 3，第三个也是 3），10（6+4）。如果只做简单的求和，你会得到 1，3，5，7。 - Marek

哇..这个答案可能是该网站上最有帮助的答案之一。将我的运行时间从几分钟缩短到了几秒钟。 - neuron

Answer 2

加速 R 代码的一般策略

首先，找出真正缓慢的部分在哪里。没有必要优化运行速度不慢的代码。对于少量的代码，简单思考即可解决。如果失败，可以使用 RProf 等性能分析工具。

一旦您找出瓶颈，考虑使用更高效的算法来实现您想要的功能。如果可能，应该只运行计算一次：

• 存储结果并访问它们而不是反复重新计算
• 将非循环相关计算从循环中提取出来
• 避免不必要的计算（例如 不要使用带有固定搜索的正则表达式）

使用更高效的函数可以产生适度或大幅的速度增益。例如，paste0 可以产生小的效率收益，但 .colSums() 及其相关函数可以产生更为明显的增益。 mean 特别慢。

接下来，您可以避免一些特别常见的问题:

• cbind 会很快减慢你的速度。
• 初始化数据结构，然后填充它们，而不是每次扩展它们。
• 即使进行了预分配，您也可以切换到传递引用而不是传递值的方法，但这可能不值得麻烦。
• 查看R Inferno以避免更多陷阱。

尝试更好的向量化，这通常可以但并非总是有帮助。在这方面，本质上向量化的命令（如ifelse、diff等）将提供比apply系列命令更大的改进（后者在编写良好的循环时几乎没有速度提升）。

您还可以尝试为R函数提供更多信息。例如，使用vapply而不是sapply, 并在读取基于文本的数据时指定colClasses。根据您消除了多少猜测，速度增益将有所不同。

接下来，考虑使用优化的软件包：data.table软件包可以在数据操作和读取大量数据(fread)时产生巨大的速度提升。

接下来，尝试通过更有效的调用R的方法来获得速度提升：

• 编译您的R脚本。或者使用Ra和jit软件包进行即时编译（Dirk在this presentation中有一个示例）。
• 确保您正在使用优化的BLAS。这些提供全面的速度提升。老实说，R没有在安装时自动使用最有效的库，这真是一件遗憾的事情。希望Revolution R将他们在这里所做的工作贡献给整个社区。
• Radford Neal进行了大量的优化，其中一些被采纳到了R Core中，许多其他优化则被分叉到了pqR中。

最后，如果以上方法仍无法满足您的速度需求，您可能需要转向一种更快的语言来处理缓慢的代码段。在此处，Rcpp和inline的组合使得用C++代码替换算法中仅有的最慢部分变得特别容易。例如，这里是我第一次尝试这样做，它甚至能够击败高度优化的R解决方案。

如果所有这些都无法解决您的问题，您只需要更多的计算能力。研究一下并行化（http://cran.r-project.org/web/views/HighPerformanceComputing.html）甚至是基于GPU的解决方案（gpu-tools）。 其他指南链接

• http://www.noamross.net/blog/2013/4/25/faster-talk.html

- Ari B. Friedman

Answer 3

如果你正在使用for循环，那么你很可能像编写C、Java或其他语言一样编写R。适当向量化的 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)}) 有类似的结果。 - Henry

@Henry，你说得很对，但正如你所指出的，结果是相同的。我已经修改了示例，将a初始化为rep(1, 1e5) - 时间相同。 - Andrie

矢量化是尽可能采用的方法，但有些循环无法以这种方式重新排列，这是真实的。 - David

Answer 4

通过使用索引或嵌套的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."

- Shane

---

谢谢您的回答。我会努力理解您的陈述。第4行：“temp[idx1,10] <- temp[which(idx1)-1,10] + temp[idx1,9]”引起了错误，因为较长对象的长度不是较短对象长度的倍数。“temp[idx1,9] = num [1:11496]”和“temp[which(idx1)-1,10] = int [1:11494]”，所以缺少了2行。 - Kay

如果您提供一个数据样本（使用dput()并包含几行），那么我会为您修复它。由于which()-1位，索引是不相等的。但是从这里开始，您应该看到它是如何工作的：没有任何循环或应用的必要；只需使用向量化函数即可。 - Shane

哇！我刚刚将一个嵌套的if..else函数块和mapply改为了一个嵌套的ifelse函数，速度提升了200倍！ - James

你的一般建议是正确的，但在代码中你忽略了一个事实，即第i个值取决于第i-1个值，因此它们不能以你所做的方式设置（使用which()-1）。 - Marek

Answer 5

我不喜欢重写代码……当然ifelse和lapply是更好的选择，但有时很难让它们适应。

经常我使用数据框作为像列表一样的工具来使用，例如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
})

数据框版本：

   user  system elapsed 
   0.53    0.00    0.53

版本列表：

   user  system elapsed 
   0.04    0.00    0.03 

使用向量列表比使用数据框快17倍。

有没有关于为什么内部数据框在这方面如此缓慢的评论？人们会认为它们的操作方式类似于列表...

要编写更快的代码，请使用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)

- Chris

---

可能是由于[<-.data.frame的开销，当您执行d$foo[i] = mark时，它会被某种方式调用，并且可能会在每个<-修改上制作一个新的向量副本或整个数据框。这将成为SO上一个有趣的问题。 - Frank

@Frank，它必须确保修改后的对象仍然是有效的数据框，并且据我所知，它至少会进行一次复制，可能会进行多次复制。众所周知，数据框子赋值操作很慢，如果你看一下长长的源代码，这并不奇怪。 - Roland

@Frank，@Roland：df$var[i]这种表示法是否也会经过[<-.data.frame函数？我发现它确实很长。如果不是，它使用的是哪个函数？ - Chris

@Chris 我相信 d$foo[i]=mark 大致上可以翻译成 d <- `$<-`(d, 'foo', `[<-`(d$foo, i, mark))，但需要使用一些临时变量。 - Tim Goodman

Answer 6

正如Ari在他的回答中提到的那样，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>

要转换为cxxfunction，需要使用include=inc。这真的很酷，因为它可以为您完成所有的链接和编译，所以原型制作非常快速。

免责声明：我不确定tmp的类是否应该是numeric而不是numeric matrix或其他什么。但我大部分是确定的。

编辑：如果您在此之后仍需要更快的速度，则OpenMP是一个适用于C++的并行化工具。我没有尝试从inline中使用它，但它应该可以工作。其思想是，在n个核心的情况下，使循环迭代k由k％n执行。Matloff的《R编程艺术》一书中的第16章“求助于C”中提供了一个适当的介绍，可在此处找到。

- jclancy

Answer 7

这里的答案很棒。但是有一个小细节没有被提及，就是问题陈述中提到了“我的电脑仍在运行（大约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)
}

- rookie

---

类似的选项，打印分数 i/n。我通常会使用 cat(sprintf("\nNow running... %40s, %s/%s \n", nm[i], i, n)) 这样的语句，因为我通常会循环遍历具有名称（名称在 nm 中）的事物。 - Frank

Answer 8

在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

---

你说向量函数比循环或apply（）更快是正确的，但不是所有情况下apply（）都比循环快。在许多情况下，apply（）只是将循环从用户中抽象出来，但仍然进行循环。请参阅此前的问题：https://dev59.com/73E95IYBdhLWcg3wheNR - JD Long

Answer 9

看一下来自{purrr}的accumulate()函数：

dayloop_accumulate <- function(temp) {
  temp %>%
    as_tibble() %>%
     mutate(cond = c(FALSE, (V6 == lag(V6) & V3 == lag(V3))[-1])) %>%
    mutate(V10 = V9 %>% 
             purrr::accumulate2(.y = cond[-1], .f = function(.i_1, .i, .y) {
               if(.y) {
                 .i_1 + .i
               } else {
                 .i
               }
             }) %>% unlist()) %>%
    select(-cond)
}

- Elena Salette

Answer 10

使用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

如果忽略条件过滤可能带来的收益，这非常快。显然，如果你可以在数据子集上进行计算，它会有所帮助。

- Bulat

---

你为什么要重复建议使用data.table呢？早在之前的回答中已经提到过多次了。 - IRTFM