r - 如何增强 for 循环

Question

我在 R 中遇到 for 循环执行缓慢的问题。这里我提供了一部分产生延迟的代码。

## subsitutes for original data
DC <- matrix(rnorm(10), ncol=101, nrow=6400)
C <- matrix(rnorm(20), ncol=101, nrow=6400)


N <- 80
Vcut <- ncol(DC) 
V <- seq(-2.9,2.5,length=Vcut)
fNC <- matrix(NA, nrow=(N*N), ncol=Vcut)
fNDC <- matrix(NA, nrow=(N*N), ncol=Vcut)


Arbfunc <- function(dV){

b <- matrix(NA, nrow=1, ncol=Vcut)

  for(i in 1:(N*N)) {
    for (n in 1:Vcut) {
      for (k in 1:Vcut) {
        b[k] = (V[2]-V[1])*(exp((-1)*abs(V[k])))*exp(abs(V[n]-V[k])/dV)*(C[i,k]/V[k])
      }
      fNC[i,n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum(b[]))
      fNDC[i,n] = DC[i,n]/fNC[i,n]
    }
  }   
}

Arbfunc(0.5)

因为我需要比较 dV's 的不同值之间的结果，这段代码至少应该在几秒钟内运行。但是结果是

user   system  elapsed
40.15   0.03   40.24

这对于足够的比较来说太慢了。我尝试了几种并行化方法，但结果并不令人满意(40 -> 25 秒，尽管我在我的电脑中使用了 11 个线程)。

因此，我的猜测是瓶颈是这个 for 循环本身，而不是非并行代码。你能给我一些建议来改进这个 for 循环或并行化提示吗？只需简短的评论，我们将不胜感激。

Answer 1

非常感谢@Mikko Marttila 更正功能 3 和 4 并提供功能 5 的想法。

R 最好使用矢量化选项而不是显式循环。例如，带有 k 的内部循环:

for (k in 1:Vcut) {
  b[k] = (V[2]-V[1])*(exp((-1)*abs(V[k])))*exp(abs(V[n]-V[k])/dV)*(C[i,k]/V[k])
}

这和说是一样的

(V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(C[i,]/V)

这个小改动让我们这部分函数的性能提升了 500 倍:

Unit: microseconds
         expr     min      lq      mean  median      uq     max neval
       k_loop 13186.7 13603.2 14605.471 13832.9 14517.8 41935.1   100
 k_vectorized    16.4    17.6    25.559    28.8    32.0    52.7   100

现在，如果我们查看带有 i 的外层循环，我们会发现实际上没有必要按每一行进行循环。我们可以为 sum(b[k]) 语句创建一个矩阵，将此转换为:

(V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(C[i,]/V)

进入这个:

(V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(t(C)/V)

这刚刚为我们节省了 N*N*k 循环。在您的情况下，这是 646,400 个循环。

总而言之，我们有:

Arbfunc3 <- function(dV){
    for (n in 1:Vcut) {
      sum_b = colSums((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(t(C)/V))
      fNC[, n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum_b)
      fNDC[, n] = DC[,n]/fNC[,n]
    }
}

对于此替代方案，我的微基准测试平均时间为 750 毫秒。

为了进一步提高性能，我们需要处理 V[n] - V。值得庆幸的是，R 有一个函数 - outer(V, V, '-')，这将生成一个包含我们需要的所有组合的矩阵。

Arbfunc4 <- function(dV) {
  sum_b = apply((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(outer(V, V, '-')) / dV) / V, 2, function(x) colSums(x * t(C)))

  fNC = exp(1*abs(V))*(1/(2*dV))*t(sum_b)
  fNDC= DC/t(fNC)
  fNDC
}

感谢@Mikko Marttila 提出的摆脱点积应用的建议。

Arbfunc5 <- function(dV) {
  a = (V[2] - V[1]) * exp(-abs(V)) * t(C) / V
  b = exp(abs(outer(V, V, "-")) / dV) %*% a

  fNC = exp(1*abs(V))*(1/(2*dV))*(b)
  fNDC= DC/t(fNC)
  fNDC
}

这是每个解决方案的 system.time(Arbfunc2 是 k_loop 的消除)。优化后的解决方案比原始解决方案快 2,600 倍。

> system.time(Arbfunc(0.5))
   user  system elapsed 
  78.03    0.39   79.72 
> system.time(Arbfunc2(0.5))
   user  system elapsed 
  10.41    0.03   10.46 
> system.time(Arbfunc3(0.5))
   user  system elapsed 
   0.69    0.13    0.81 
> system.time(Arbfunc4(0.5))
   user  system elapsed 
   0.43    0.05    0.47 
> system.time(Arbfunc5(0.5))
   user  system elapsed 
   0.03    0.00    0.03 

最终编辑:这是我在重新启动 R 并清空我的环境后运行的完整代码。没有错误:

## subsitutes for original data
DC <- matrix(rnorm(10), ncol=101, nrow=6400)
C <- matrix(rnorm(20), ncol=101, nrow=6400)

N <- 80
Vcut <- ncol(DC) 
V <- seq(-2.9,2.5,length=Vcut)

# Unneeded for Arbfunc4 adn Arbfunc5
# Corrected from NA to NA_real_ to prevent coercion from logical to numeric
# h/t to @HenrikB
fNC <- matrix(NA_real_, nrow=(N*N), ncol=Vcut)
fNDC <- matrix(NA_real_, nrow=(N*N), ncol=Vcut)

Arbfunc <- function(dV){
  b <- matrix(NA, nrow=1, ncol=Vcut)

  for(i in 1:(N*N)) {
    for (n in 1:Vcut) {
      for (k in 1:Vcut) {
        b[k] = (V[2]-V[1])*(exp((-1)*abs(V[k])))*exp(abs(V[n]-V[k])/dV)*(C[i,k]/V[k])
      }
      fNC[i,n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum(b[]))
      fNDC[i,n] = DC[i,n]/fNC[i,n]
    }
  }
  fNDC
}

Arbfunc2 <- function(dV){
  b <- matrix(NA, nrow=1, ncol=Vcut)

  for(i in 1:(N*N)) {
    for (n in 1:Vcut) {
      sum_b = sum((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(C[i,]/V))
      fNC[i,n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum_b)
      fNDC[i,n] = DC[i,n]/fNC[i,n]
    }
  }
  fNDC
}

Arbfunc3 <- function(dV){
  for (n in 1:Vcut) {
    sum_b = colSums((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(t(C)/V))
    fNC[, n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum_b)
    fNDC[, n] = DC[,n]/fNC[,n]
  }
  fNDC
}

Arbfunc4 <- function(dV) {
  sum_b = apply((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(outer(V, V, '-')) / dV) / V, 2, function(x) colSums(x * t(C)))

  fNC = exp(1*abs(V))*(1/(2*dV))*t(sum_b)
  DC/t(fNC)
}

Arbfunc5 <- function(dV) {
#h/t to Mikko Marttila for dot product
  a = (V[2] - V[1]) * exp(-abs(V)) * t(C) / V
  b = exp(abs(outer(V, V, "-")) / dV) %*% a

  fNC = exp(1*abs(V))*(1/(2*dV))*(b)
  DC/t(fNC)
}

#system.time(res <- Arbfunc(0.5))
system.time(res2 <- Arbfunc2(0.5))
system.time(res3 <- Arbfunc3(0.5))
system.time(res4 <- Arbfunc4(0.5))
system.time(res5 <- Arbfunc5(0.5))

all.equal(res2,res3,res4,res5)

正如@HenrikB 提到的，fNC 和fNDC 初始化为逻辑矩阵。这意味着我们在将它们强制转换为 real 矩阵时会受到性能影响。不正确的做法是该数据集的一次性命中 1 毫秒，但如果这种强制转换处于循环中，它可能真的会累加起来。

mat_NA_real_ <- function() {
  mat = matrix(NA_real_, nrow = 6400, ncol = 101)
  mat[1,1] = 1
}

mat_NA <- function() {
  mat = matrix(NA, nrow = 6400, ncol = 101)
  mat[1,1] = 1
}
microbenchmark(mat_NA_real_(), mat_NA())

Unit: microseconds
           expr    min      lq     mean  median     uq     max neval
 mat_NA_real_()  979.5  992.25 1490.081  998.65 1021.1  7612.5   100
       mat_NA() 1865.8 1883.30 3793.119 1911.30 5335.4 53635.2   100