r语言 - 为什么循环代码的速度与 apply() 一样快?

Question

>已编辑以提供可重复的结果

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

#

我最初使用 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 (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
}

应用(( 版本：

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 函数;使用泛函的优点是简单明了。下面是一个将 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.

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

也就是说，对于您的特定问题，很可能还有其他优化。您可能希望提供一个可重现的示例，并在代码审查中提出有关性能改进的新问题。

^{创建于 2019-09-02 由 reprex 软件包 (v0.3.0.9000(}

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
)