加速R中的循环操作

Question

我在R中遇到了很大的性能问题.我编写了一个迭代data.frame对象的函数.它只是添加一个新列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)
}

有什么想法如何加快这个操作？

Answer 1

最大的问题和无效的根源是索引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秒.

矢量化的力量

正如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)
}

比较此功能的结果,此时间为nrow10,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上获得.

Answer 2

加速R代码的一般策略

首先,弄清楚其中的慢部分确实是.没有必要优化运行缓慢的代码.对于少量代码,只需通过思考即可.如果失败,RProf和类似的分析工具可能会有所帮助.

一旦找出瓶颈,就要考虑更有效的算法来做你想做的事情.如果可能,计算应仅运行一次,因此:

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

使用更高效的功能可以产生中等或大的速度增益.例如,paste0产生小的效率增益,但.colSums()其亲属产生一些更明显的收益. mean是特别慢.

然后你可以避免一些特别常见的麻烦:

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

尝试更好的矢量化,这通常但不总是有帮助.在这一点上,本质上向量化的命令一样ifelse,diff等会提供比更改善apply家庭命令(其中没有速度提升提供有点过了精心编写的循环)的.

您还可以尝试向R函数提供更多信息.例如,使用vapply而不是sapply,并colClasses在读取基于文本的数据时指定.速度增益将根据您消除的猜测量而变化.

接下来,考虑优化的软件包:data.table软件包可以在可能的情况下,在数据操作和读取大量数据时产生大量的速度增益(fread).

接下来,通过更有效的方式调用R来尝试提高速度:

• 编译你的R脚本.或者使用Raand和jit音乐包进行即时编译(Dirk在本演示文稿中有一个例子).
• 确保您使用的是优化的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

Answer 3

如果您正在使用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 

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."

Answer 5

我不喜欢重写代码......当然,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 

使用矢量列表比使用data.frame快17倍.

关于为什么内部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)

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,as include=inc.真正酷的是它为你完成了所有的链接和编译,因此原型设计非常快.

免责声明:我不完全确定tmp类应该是数字而不是数字矩阵或其他东西.但我很肯定.

编辑:如果此后仍需要更高的速度,OpenMP是一个适合的并行化工具C++.我没有尝试使用它inline,但它应该工作.在n核心的情况下,想法是k通过执行循环迭代k % n.合适的引入在Matloff的发现[R编程的艺术,可在这里,在第16章,诉诸到C.

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)
}