使用嵌套 for 循环改进 R 的 运行 时间

Improving run time for R with nested for loops

我的可重现 R 示例:

f = runif(1500,10,50)
p = matrix(0, nrow=1250, ncol=250)
count = rep(0, 1250)
for(i in 1:1250) {
    ref=f[i]
    for(j in 1:250) {
        p[i,j] = f[i + j - 1] / ref-1
        if(p[i,j] == "NaN") {
           count[i] = count[i]
           }
        else if(p[i,j] > (0.026)) {
                count[i] = (count[i] + 1) 
                ref = f[i + j - 1] 
                } 
        } 
    }

更准确地说,我有一组 600 个 f 系列,每个 f 系列的代码 运行s 200 次。目前我正在循环中进行迭代,并且大多数操作都是按元素进行的。我的随机变量是 f、条件 if(p[i,j] > (0.026)) 和数字 0.026 本身。

可以通过矢量化我的代码和使用函数(特别是 apply 系列)来显着减少 运行 时间,但我对 apply 感到生疏,正在寻找一些建议以朝着正确的方向前进。

在Rcpp中放入for loop。我只是将您的代码复制粘贴到 Rcpp 并且没有检查有效性。如有差异,请告诉我。 fCpp returns pc 的列表。

cppFunction('List fCpp(NumericVector f) {
    const int n=1250; 
            const int k=250;
            NumericMatrix p(n, k);
            NumericVector c(n);

            for(int i = 0; i < n; i++) {
            double ref=f[i];
            for(int j = 0; j < k; j++) {
            p(i,j) = f[i+j+1]/ref-1;
            if(p(i,j) == NAN){
            c[i]=c[i];
            }
            else if(p(i,j) > 0.026){
            c[i] = c[i]+1; 
            ref = f[i+j+1]; 
            } 
            }
            }
            return List::create(p, c);
            }')

基准

set.seed(1)
f = runif(1500,10,50)

f1 <- function(f){
    p = matrix(0, nrow=1250, ncol=250)
    count = rep(0, 1250)
    for(i in 1:1250) {
        ref=f[i]
        for(j in 1:250) {
            p[i,j] = f[i + j - 1] / ref-1
            if(p[i,j] == "NaN") {
                count[i] = count[i]
            }
            else if(p[i,j] > (0.026)) {
                count[i] = (count[i] + 1) 
                ref = f[i + j - 1] 
            } 
        } 
    }
    list(p, count)
}


microbenchmark::microbenchmark(fCpp(f), f1(f), times=10L, unit="relative")
Unit: relative
    expr      min       lq     mean   median       uq      max neval
 fCpp(f)   1.0000   1.0000   1.0000   1.0000   1.0000   1.0000    10
   f1(f) 785.8484 753.7044 734.4243 764.5883 718.0868 644.9022    10

fCpp(f)f1(f) 返回的值基本相同,除了 f1 返回的 p 矩阵的第 1 列填充了 0。

system.time(a <- f1(f))[3]
#elapsed 
#    2.8 
system.time(a1 <- fCpp(f))[3]
#elapsed 
#      0 
all.equal( a[[1]], a1[[1]])
#[1] "Mean relative difference: 0.7019406"
all.equal( a[[2]], a1[[2]])
#[1] TRUE

这是一个使用 while 的实现,虽然它比嵌套 for 循环花费的时间长得多,这有点违反直觉。

f1 <- function() {
    n <- 1500
    d <- 250
    f = runif(n,1,5)
    f = embed(f, d)
    f = f[-(n-d+1),]
    count = rep(0, n-d)
    for(i in 1:(n-d)) {
        tem <- f[i,]/f[i,1] - 1
        ti <- which(t[-d] > 0.026)[1]
        while(ti < d & !is.na(ti)) {
            ti.plus = ti+1
            tem[ti.plus:d] = f[i, ti.plus:d] / tem[ti]
            count[i] = count[i] + 1
            ti <- ti + which(tem[ti.plus:d-1] > 0.026)[1]
        }
        f[i] = tem
    }
    list(f, count)
}

system.time(f1())

#elapsed 
#6.365

@ajmartin,你的逻辑更好,减少了我尝试的迭代次数。这是 R 中代码的改进版本:

f1 <- function() {
  n <- 1500
  d <- 250
  f = runif(n,1,5)
  count = rep(0, n-d)
  for(i in 1:(n-d)) {
    tem <- f[i:(i+d-1)] / f[i] - 1
    ind = which(tem>0.026)[1]
    while(length(which(tem>0.026))){
      count[i] = count[i] + 1
      tem[ind:d] = f[ind:d] / tem[ind] - 1
      ind = ind - 1 + (which(tem[ind:d] > 0.026)[1])
    }
  }
  list(f, count)
}

system.time(f1())[3]
# elapsed 
#    0.09 

在 Rcpp 中实现此功能将进一步减少系统时间,但我无法安装 Rtools,因为我当前的计算机没有管理员权限。同时这有帮助。