嵌套for循环的向量化

Question

我有以下函数：

g.Bn = function(n) {
  Bn = 0
  for(k in 0:n) {
    res.of.loop = 0
    for(j in 0:k) {
      res.of.loop = res.of.loop + (-1)^j * (j + 1)^n * choose(k, j)
    }
    Bn = res.of.loop * 1/(k+1) + Bn
  }
  return(Bn)
}

这里有一种向量化的方法可以代替使用for循环吗？

- JWallmark

---

你可以尝试用j<-0:k; res.of.loop<- sum((-1)^j * (j + 1)^n * choose(k, j))替换内部循环。 - Dave2e

修改后的 g.Bn 应该是一个向量化函数吗？意思是当提供输入向量时返回一组答案的函数？ - slava-kohut

@slava-kohut 或许我的问题表述不够清晰。我想使用向量操作，但返回值仍应为实数。 - JWallmark

哦等等，你可以用C语言实现。这可能会大大提高速度。 - waferthin

Answer 1

你可以将内部循环向量化（按照@DaveT所述），并使用sapply：

g.Bn2 = function(n) {
  sum(sapply(0:n, function(k) {
    sum((-1)^(0:k) * (0:k + 1)^n * choose(k, 0:k)) * 1/(k+1)
  }))
}

或者另一种将外层循环向量化的可能性：

g.Bn3 = function(n) {
  f <- function(k, n) sum((-1)^(0:k) * (0:k + 1)^n * choose(k, 0:k)) * 1/(k+1)
  sum(Vectorize(f, vectorize.args = "k")(0:n, n))
}

> microbenchmark(g.Bn(100), g.Bn2(100), g.Bn3(100))
       expr      min        lq      mean   median        uq      max neval
  g.Bn(100) 1493.086 1533.9280 1841.3455 1585.354 1675.3575 9023.316   100
 g.Bn2(100)  617.063  650.7850  905.6899  738.230  788.7305 9224.460   100
 g.Bn3(100)  685.094  772.3785 1015.9182  816.945  860.1775 8213.777   100

- waferthin

---

尝试了你的代码，发现在 n 超过某个阈值后，g.Bn2 (g.Bn3) 会产生不同的结果。请尝试使用 g.Bn(15) 和 g.Bn2(15) 进行测试。这是为什么呢？ - slava-kohut

老实说，我不确定，但我认为精度在某种程度上很重要。在这方面，我无法说哪种方法更好。 - waferthin

我猜这可能与浮点算术有关...尝试在对数刻度上进行操作？ - user20650

@waferthin 感谢回复。我发现只使用 DaveT 的方法来处理内部循环，保持外部循环不变，可以给我比你的任何解决方案都更好的微基准测试结果。你认为这是为什么呢？sapply 是否只是以类似于循环的方式迭代？ g.Bn = function(n) { Bn = 0 for(k in 0:n) { j = c(0:k) res.of.loop = sum((-1)^j * (j + 1)^n * choose(k, j)) Bn = Bn + res.of.loop * 1/(k+1) } return(Bn) } - JWallmark

是的，真正的速度提升可以通过向量化内部循环来实现。所有其他更改实际上只是不同形式的循环（可以说是具有更好的语法，但没有真正的速度增益）。我认为进一步加速的唯一方法是考虑方程是否可以以任何方式简化。 - waferthin

Answer 2

你可以将 for-loop 转换为 map 和 reduce。

在下面的示例中，purrr::map 迭代所有数据，而 sum 将数值向量缩减为一个标量（长度为一的数值向量）。

g.Bn = function(n) {
    sum(
        purrr::map_dbl(0:n,function(k){
            sum(
                purrr::map_dbl(0:k,function(j){
                    (-1)^j * (j + 1)^n * choose(k, j)
                })
            ) * 1/(k+1)

        })
    )
}

看起来内部循环中的所有j都可以替换为0:k

g.Bn = function(n) {
    sum(
        purrr::map_dbl(0:n,function(k){
            sum((-1)^(0:k) * (0:k + 1)^n * choose(k, 0:k)) * 1/(k+1)
        })
    )
}

- yusuzech

---

两个 for 被两个 purrr:: 替换了。并且两者都使用迭代。这使它更加强大的原因是什么？我不是在批评，我是认真询问。 - maydin

使用 purrr:: 可以避免至少分配初始化变量的麻烦。对于我来说，使用嵌套的匿名函数可以使每个层级的逻辑非常清晰。 - yusuzech