如何矢量化三重嵌套循环？

Question

我已经搜索过类似的问题，但我对我应该做什么有一个模糊的想法：将所有内容矢量化或使用apply()family。但是我是R编程的初学者，以上两种方法都很混乱。

这是我的源代码：

x<-rlnorm(100,0,1.6)
j=0
k=0
i=0
h=0
lambda<-rep(0,200)
sum1<-rep(0,200)
constjk=0
wj=0
wk=0
for (h in 1:200)
{
   lambda[h]=2+h/12.5
   N=ceiling(lambda[h]*max(x))
   for (j in 0:N)
   {
      wj=(sum(x<=(j+1)/lambda[h])-sum(x<=j/lambda[h]))/100
      for (k in 0:N)
      {
         constjk=dbinom(k, j + k, 0.5)
         wk=(sum(x<=(k+1)/lambda[h])-sum(x<=k/lambda[h]))/100
         sum1[h]=sum1[h]+(lambda[h]/2)*constjk*wk*wj
      }
   }
}

让我解释一下。我想收集 200 个 sum1 值（这是第一个循环），对于每个 sum1 值，它是 (lambda[h]/2)*constjk*wk*wj 的总和，因此是其他两个循环的总和。最繁琐的是N随h变化，所以我不知道如何向量化j-loop和k-loop。但是我当然可以用lambda<-seq() 和N<-ceiling() 对h-loop 进行矢量化，这是我能做的最好的。有没有办法进一步简化代码？

Answer 1

您的代码可以通过 3 个嵌套的 sapply 调用完美验证。对于未经训练的人来说可能有点难以阅读，但它的本质是，我们不是一次将一个值添加到sum1[h]，而是一次性计算出最内层循环产生的所有项并将它们相加。

虽然这个矢量化解决方案比您的三倍 for 循环更快，但改进并不显着。如果您打算多次使用它，我建议您在 C 或 Fortran 中实现它（使用常规的 for 循环），这会提高速度很多。请注意，尽管它具有很高的时间复杂度，并且会随着 lambda 值的增加而严重扩展，最终达到一个点，即无论实现如何都无法在合理的时间内进行计算。

lambda <- 2 + 1:200/12.5
sum1 <- sapply(lambda, function(l){
    N <- ceiling(l*max(x))
    sum(sapply(0:N, function(j){
        wj <- (sum(x <= (j+1)/l) - sum(x <= j/l))/100
        sum(sapply(0:N, function(k){
            constjk <- dbinom(k, j + k, 0.5)
            wk <- (sum(x <= (k+1)/l) - sum(x <= k/l))/100
            l/2*constjk*wk*wj
        }))
    }))
})

顺便说一句，您不需要预定义变量，例如 h、j、k、wj 和 wk。特别是因为不是在矢量化时，因为在馈送到sapply 的函数中对它们的赋值将创建具有相同名称的重叠局部变量（即忽略您预先定义的那些）。

Answer 2

让我们将您的模拟封装在一个函数中并对其计时：

sim1 <- function(num=20){
  set.seed(42)
  x<-rlnorm(100,0,1.6)
  j=0
  k=0
  i=0
  h=0
  lambda<-rep(0,num)
  sum1<-rep(0,num)
  constjk=0
  wj=0
  wk=0

  for (h in 1:num)
  {
    lambda[h]=2+h/12.5
    N=ceiling(lambda[h]*max(x))
    for (j in 0:N)
    {
      wj=(sum(x<=(j+1)/lambda[h])-sum(x<=j/lambda[h]))/100
      for (k in 0:N)
      {
        set.seed(42)
        constjk=dbinom(k, j + k, 0.5)
        wk=(sum(x<=(k+1)/lambda[h])-sum(x<=k/lambda[h]))/100
        sum1[h]=sum1[h]+(lambda[h]/2)*constjk*wk*wj
      }
    }
  }

  sum1
}

system.time(res1 <- sim1())
#   user  system elapsed 
#    5.4     0.0     5.4

现在让我们加快速度：

sim2 <- function(num=20){
  set.seed(42) #to make it reproducible
  x <- rlnorm(100,0,1.6)

  h <- 1:num
  sum1 <- numeric(num)
  lambda <- 2+1:num/12.5
  N <- ceiling(lambda*max(x))

  #functions for wj and wk
  wjfun <- function(x,j,lambda,h){
    (sum(x<=(j+1)/lambda[h])-sum(x<=j/lambda[h]))/100
  }
  wkfun <- function(x,k,lambda,h){
    (sum(x<=(k+1)/lambda[h])-sum(x<=k/lambda[h]))/100
  }

  #function to calculate values of sum1
  fun1 <- function(N,h,x,lambda) {
    sum1 <- 0
    set.seed(42) #to make it reproducible
    #calculate constants using outer
    const <- outer(0:N[h],0:N[h],FUN=function(j,k) dbinom(k, j + k, 0.5))
    wk <- numeric(N[h]+1)
    #loop only once to calculate wk
    for (k in 0:N[h]){
      wk[k+1] <- (sum(x<=(k+1)/lambda[h])-sum(x<=k/lambda[h]))/100 
    }

    for (j in 0:N[h])
    {
      wj <- (sum(x<=(j+1)/lambda[h])-sum(x<=j/lambda[h]))/100
      for (k in 0:N[h])
      {
        sum1 <- sum1+(lambda[h]/2)*const[j+1,k+1]*wk[k+1]*wj
      }
    }
    sum1
  }

  for (h in 1:num)
  {
    sum1[h] <- fun1(N,h,x,lambda)
  }  
  sum1
}

system.time(res2 <- sim2())
#user  system elapsed 
#1.25    0.00    1.25 

all.equal(res1,res2)
#[1] TRUE

@Backlin 代码的时间安排（有 20 次交互）以进行比较：

   user  system elapsed 
   3.30    0.00    3.29 

如果这仍然太慢并且您不能或不想使用另一种语言，那么并行化也是可能的。据我所知，外循环是令人尴尬的平行。有一些很好的和简单的并行化包。