加速R循环

Question

使用apply系列中的函数可以轻松地在R中加速循环.如何在下面的代码中使用apply函数来加快速度？注意,在循环内,在每次迭代时,置换一列并将函数应用于新数据帧(即,具有一列置换的初始数据帧).我似乎无法申请工作,因为必须在循环内构建新的数据框.

#x <- data.frame(a=1:10,b=11:20,c=21:30) #small example
x <- data.frame(matrix(runif(50*100),nrow=50,ncol=100)) #larger example
y <- rowMeans(x)

start <- Sys.time()  

totaldiff <- numeric()

for (i in 1:ncol(x)){
    x.after <- x

    x.after[,i] <- sample(x[,i])

    diff <- abs(y-rowMeans(x.after))

    totaldiff[i] <- sum(diff)

}

colnames(x)[which.max(totaldiff)]

Sys.time() - start

Answer 1

在完成此回复和其他回复之后,这里的优化策略(以及近似加速)似乎就是这样

• (30x)选择适当的数据表示 - 矩阵,而不是data.frame
• (1.5x)减少不必要的数据副本 - 列的差异,而不是rowMeans
• 循环结构作为*apply函数(强调代码结构,简化内存管理,并提供类型一致性)
• (2x)循环外的提升矢量运算 - 列上的abs和sum变为矩阵上的abs和colSums

总体速度提升约100倍.对于此代码的大小和复杂性,编译器或并行包的使用将无效.

我把你的代码放到一个函数中

f0 <- function(x) {
    y <- rowMeans(x)
    totaldiff <- numeric()
    for (i in 1:ncol(x)){
        x.after <- x
        x.after[,i] <- sample(x[,i])
        diff <- abs(y-rowMeans(x.after))
        totaldiff[i] <- sum(diff)
    }
    which.max(totaldiff)
}

我们在这里

x <- data.frame(matrix(runif(50*100),nrow=50,ncol=100)) #larger example
set.seed(123)
system.time(res0 <- f0(x))
##   user  system elapsed 
##  1.065   0.000   1.066 

您的数据可以表示为矩阵,R矩阵上的操作比data.frames上的操作更快.

m <- matrix(runif(50*100),nrow=50,ncol=100)
set.seed(123)
system.time(res0.m <- f0(m))
##   user  system elapsed 
##  0.036   0.000   0.037 
identical(res0, res0.m)
##[1] TRUE

这可能是最大的加速.但是对于这里的具体操作,我们不需要计算更新矩阵的行平均值,只需要改变平均值来改变一列

f1 <- function(x) {
     y <- rowMeans(x)
    totaldiff <- numeric()
    for (i in 1:ncol(x)){
        diff <- abs(sample(x[,i]) - x[,i]) / ncol(x)
        totaldiff[i] <- sum(diff)
    }
    which.max(totaldiff)
}

该for循环不遵循正确的模式填补了结果向量totaldiff(要"预分配,并填写",这样totaldiff <- numeric(ncol(x))),但我们可以使用一个sapply,让有关的R忧(此内存管理的优势之一使用apply系列函数)

f2 <- function(x) {
    totaldiff <- sapply(seq_len(ncol(x)), function(i, x) {
        sum(abs(sample(x[,i]) - x[,i]) / ncol(x))
    }, x)
    which.max(totaldiff)
}
set.seed(123); identical(res0, f1(m))
set.seed(123); identical(res0, f2(m))

时间是

> library(microbenchmark)
> microbenchmark(f0(m), f1(m), f2(m))
Unit: milliseconds
  expr      min       lq   median       uq      max neval
 f0(m) 32.45073 33.07804 33.16851 33.26364 33.81924   100
 f1(m) 22.20913 23.87784 23.96915 24.06216 24.66042   100
 f2(m) 21.02474 22.60745 22.70042 22.80080 23.19030   100

@flodel指出vapply可以更快(并提供类型安全)

f3 <- function(x) {
    totaldiff <- vapply(seq_len(ncol(x)), function(i, x) {
        sum(abs(sample(x[,i]) - x[,i]) / ncol(x))
    }, numeric(1), x)
    which.max(totaldiff)
}

然后

f4 <- function(x)
    which.max(colSums(abs((apply(x, 2, sample) - x))))

仍然更快(ncol(x)是一个恒定的因素,所以删除) - abs并且sum被提升在外面sapply,可能以额外的内存使用为代价.评论中对编译函数的建议总的来说是好的; 这里有一些进一步的时间安排

>     microbenchmark(f0(m), f1(m), f1.c(m), f2(m), f2.c(m), f3(m), f4(m))
Unit: milliseconds
    expr      min       lq   median       uq       max neval
   f0(m) 32.35600 32.88326 33.12274 33.25946  34.49003   100
   f1(m) 22.21964 23.41500 23.96087 24.06587  24.49663   100
 f1.c(m) 20.69856 21.20862 22.20771 22.32653 213.26667   100
   f2(m) 20.76128 21.52786 22.66352 22.79101  69.49891   100
 f2.c(m) 21.16423 21.57205 22.94157 23.06497  23.35764   100
   f3(m) 20.17755 21.41369 21.99292 22.10814  22.36987   100
   f4(m) 10.10816 10.47535 10.56790 10.61938  10.83338   100

其中".c"是编译版本和

编译在用for循环编写的代码中特别有用,但对矢量化代码没有太大作用; 这里显示的是编译f1 for循环的一个小但一致的改进,但不是f2的sapply.