加快R中的循环操作

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

有什么办法可以加快这次行动吗?

106667 次浏览
小开

这可以通过使用索引或嵌套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."
小开

在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]

向量化算术需要更多的时间和思考问题，但有时可以节省几个数量级的执行时间。

小开
最大的问题和无效的根源是索引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中，我不需要打乱名称。 那么它有多好呢?< / p >

我用nrow从1000到10,000 × 1000运行data.frame的每个函数，并用system.time测量时间

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

结果是

performance

你可以看到你的版本指数依赖于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。

performance

调谐调谐

另一个调整是在循环中将temp[i,9]的索引改为res[i](这在第i次循环迭代中完全相同)。 这也是索引向量和索引data.frame的区别。
第二件事:当你查看循环时，你可以看到没有必要遍历所有i，而只遍历符合条件的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行小于1秒的计算时间

performance

如果你想更进一步，我认为至少有两件事是可以做到的:

  • 编写C代码来执行条件cumsum

  • 如果你知道你的数据Max序列不是很大，那么你可以将循环改为向量化，而

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

Code used for simulations and figures is available on GitHub.

小开

如果你正在使用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
小开

加速R代码的一般策略

首先，找出在哪里慢的部分。没有必要优化运行速度不慢的代码。对于少量的代码，简单地思考一下就可以了。如果失败了，RProf和类似的分析工具可能会有所帮助。

一旦你找到了瓶颈，考虑更高效的算法来做你想做的事情。如果可能，计算应该只运行一次，因此:

使用更多的高效的功能可以产生中等或较大的速度增益。例如，paste0产生了一个小的效率增益，而.colSums()和它的亲戚产生了一些更明显的增益。mean是特别慢。

然后你可以避免一些特别的常见的问题:

尝试更好的向量化，这通常可以但不总是有帮助。在这方面，像ifelse、diff之类的固有向量化命令将比apply系列命令提供更多的改进(在编写良好的循环中，apply几乎没有提供速度提升)。

你也可以尝试为R函数提供更多的信息。例如，使用__ABC0而不是sapply，并指定当读入基于文本的数据时colClasses。速度的增加取决于你消除了多少猜测。

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

接下来，尝试通过更有效地调用R来提高速度:

  • 编译你的R脚本。或者同时使用Ra和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)。

其他指引的连结

小开

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

#includeing也有类似的过程，你只需要传递一个参数

inc <- '#include <header.h>

到cxxfunction，如include=inc。最酷的是它为你做了所有的链接和编译，所以原型制作非常快。

免责声明:我不完全确定tmp的类应该是数字而不是数字矩阵或其他东西。但我基本确定。

编辑:如果在此之后仍然需要更快的速度，OpenMP是一个适合C++的并行化工具。我还没有尝试从inline使用它，但它应该工作。这个想法是，在n内核的情况下，让循环迭代k由k % n执行。一个合适的介绍可以在Matloff的R编程的艺术中找到，可用在这里，在第16章，求助于C中。

小开

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

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

小开

我不喜欢重写代码……当然，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)
小开

这里的答案很好。一个没有涉及到的小方面是，问题声明“我的电脑还在工作(大约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)
}
小开

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

提交答案