如何增强循环

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的各个值之间比较结果，因此该代码至少应在几秒钟内运行。但是结果是

user   system  elapsed
40.15   0.03   40.24

这对于足够的比较来说太慢了。我尝试了几种并行化方法，但结果并不令人满意（40-> 25秒，尽管我在PC中使用了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的消除）。优化的解决方案比原始解决方案快2600倍。

> 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 ms的命中，但是如果此强制执行处于循环中，则它实际上可能加起来。

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