为什么这个 for 循环代码和 apply() 一样快？

Question

经过编辑以提供可重现的结果

最初并不清楚，但我需要将结果考虑到原始数据中的 NA (df)

#

我最初使用 for 循环编写代码以使概念证明有效，但现在我需要加快速度。它使用 for 循环在大约 2-3 分钟内运行，但是当我使用 apply() 重新编写它时，它并没有更快。我认为 apply() 应该是矢量化解决方案，因此速度更快，但也许我的整个前提不正确？ （我对 R 并不陌生，但计算速度对我来说通常不是问题。）

我正在处理 1000 多个案例和约 100 个变量，需要对数据执行 5000 多个模拟（打开和关闭不同的条件）。

起始定义和示例数据：

cases = 1000
variables = 100
simulations = 5000

df <- as.data.frame(array(rnorm(cases * variables, 0, 5), dim=c(cases, variables)))

montecarlo <- matrix(rbinom(simulations * variables, 1, 20/variables), simulations, variables)
montecarlo[montecarlo==0] <- NA

calc <- array(,dim=c(cases, variables, simulations))
interim <- array(,dim=c(cases, variables, simulations))
results <- array(,dim=c(variables, simulations))

for (j in 1:simulations) {
  calc[,,j] <- exp(t(t(df) * as.numeric(montecarlo[j,])))
}

For循环版本：

for (j in 1:simulations) {
  interim[,,j] <- t(apply(calc[,,j], 1, function(x) x/sum(x, na.rm = TRUE))) # re-share
  results[,j] <- apply(interim[,,j], 2, sum) # aggregates results
}

Apply() 版本：

interim <- apply(calc, c(1,3), function(x) x/sum(x, na.rm = TRUE)) # re-share
results <- as.data.frame(t(apply(interim, c(1,3), sum))) # aggregates results

我愿意接受任何加快速度的建议和/或 apply() 版本没有更快的原因。谢谢！

Answer 1

一般来说：for 循环本身并不慢。如果您正在预分配输出（即不增长向量，导致多个副本），它们在速度上几乎可以与 *apply() 函数相媲美。迭代的开销来自于从 C 代码中重复调用 R 函数；而使用泛函的好处就是clarity。这是一个包含在函数中的 for 循环的示例：

foo <- function(x, f, ...) {
  out <- vector("list", length(x))
  for (i in seq_along(x)) {
    out[[i]] <- f(x[[i]], ...)
  }
  out
}

x <- replicate(10000, rnorm(30), simplify = FALSE)
bench::mark(foo(x, mean), lapply(x, mean))
#> # A tibble: 2 x 6
#>   expression           min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr>      <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 foo(x, mean)      36.9ms   38.4ms      26.1   157.3KB     52.3
#> 2 lapply(x, mean)   42.6ms   44.9ms      22.3    78.2KB    100.

在这些情况下提高速度的方法是将所有计算转移到编译代码中。

也就是说，对于您的特定问题，很可能还有其他优化。您可能想提供一个可重现的示例，并在 Code Review 上提出一个关于性能改进的新问题。

^{由reprex package (v0.3.0.9000) 于 2019 年 9 月 2 日创建}

Answer 2

正如@Mikko Marttila 所指出的，apply() 系列不保证更快的代码。使用 sweep 和 aperm()，对于 1,000 x 100 x 70 数组（即只有 700 万个元素），下面的代码大约快 3 倍。

results4 <- colSums(sweep(calc, c(1,3), colSums(aperm(calc, c(2,1,3)), na.rm = T), FUN = '/'), na.rm = T)

或者，性能稍差，但与您原来的性能更相似：

interim3 <- sweep(calc, c(1,3), apply(calc, 3, rowSums, na.rm = T), FUN = '/')
results3 <- apply(interim3, c(2,3), sum, na.rm = T)

性能：

Unit: milliseconds
          expr      min       lq     mean   median       uq      max neval
      for_loop 510.9131 514.9030 537.0344 518.2491 524.5709 705.4087    10
      apply_OP 446.0352 458.4940 491.6710 500.1995 523.1843 533.9654    10
 sweep_rowSums 225.5855 233.2632 252.6149 240.7245 284.1517 292.3476    10
   sweep_aperm 136.2519 140.8912 163.7498 154.6984 191.5337 217.8015    10

数据

cases = 1000
variables = 100
simulations = 70

set.seed(123)
calc <- array(sample(cases *variables * simulations),dim=c(cases, variables, simulations))

interim <- array(,dim=c(cases, variables, simulations))
results <- array(,dim=c(variables, simulations))

# Original Loop
for (j in seq_len(simulations)) {
  interim[,,j] <- t(apply(calc[,,j], 1, function(x) x/sum(x, na.rm = TRUE))) # re-share
  results[,j] <- apply(interim[,,j], 2, sum) # aggregates results
}

# original apply

interim2 <- apply(calc, c(1,3), function(x) x/sum(x, na.rm = TRUE)) # re-share
results2 <- apply(interim2, c(1,3), sum) # aggregates results

# using sweep

interim3 <- sweep(calc, c(1,3), apply(calc, 3, rowSums, na.rm = T), FUN = '/')
results3 <- apply(interim3, c(2,3), sum, na.rm = T)

#using sweep and aperm
# interim4 <- sweep(calc, c(1,3), colSums(aperm(calc, c(2,1,3)), na.rm = T), FUN = '/')
results4 <- colSums(sweep(calc, c(1,3), colSums(aperm(calc, c(2,1,3)), na.rm = T), FUN = '/'), na.rm = T)

all.equal(results4, results3, results2, results)



library(microbenchmark)

microbenchmark(
  for_loop = {
    for (j in seq_len(simulations)) {
      interim[,,j] <- t(apply(calc[,,j], 1, function(x) x/sum(x, na.rm = TRUE))) # re-share
      results[,j] <- apply(interim[,,j], 2, sum) # aggregates results
    }
  }
  ,
  apply_OP = {
    interim2 <- apply(calc, c(1,3), function(x) x/sum(x, na.rm = TRUE)) # re-share
    results2 <- apply(interim2, c(1,3), sum) # aggregates results
  }
  ,
  sweep_rowSums = {
    interim3 <- sweep(calc, c(1,3), apply(calc, 3, rowSums, na.rm = T), FUN = '/')
    results3 <- apply(interim3, c(2,3), sum, na.rm = T)
  }
  ,
  sweep_aperm = {
    results4 <- colSums(sweep(calc, c(1,3), colSums(aperm(calc, c(2,1,3)), na.rm = T), FUN = '/'), na.rm = T)
  }
  , times = 10
)