循环和引导脚本运行时间过长

Question

我有下面的R 脚本，它需要超过 24 小时才能运行，但最终在 10-gigabyte ram 和 core M7 的 Windows 10 上运行。该脚本执行以下操作：

这就是我想用R做的事情

• A.我已经生成了 50 个时间序列数据集。

• B.我将相同的时间序列数据集分割成以下大小的块：2,3,...,48,49 使我从上面的步骤 1 中形成了 48 个不同的时间序列。

• C.我将每个 48 时间序列数据集划分为 train 和 test 集，因此我可以使用 Metrics 包中的 rmse 函数来获取步骤 2 中形成的 48 个子序列的均方根误差 (RMSE)。

• D.然后根据其块大小将每个系列的 RMSE 制成表格

• E.我为每个 48 个不同的时间序列数据集获得了最好的 ARIMA 模型。

我的 R 脚本

# simulate arima(1,0,0)
library(forecast)
library(Metrics)

n=50
phi <- 0.5
set.seed(1)

wn <- rnorm(n, mean=0, sd=1)
ar1 <- sqrt((wn[1])^2/(1-phi^2))

for(i in 2:n){
  ar1[i] <- ar1[i - 1] * phi + wn[i]
}
ts <- ar1

t <- length(ts)    # the length of the time series
li <- seq(n-2)+1   # vector of block sizes to be 1 < l < n (i.e to be between 1 and n exclusively)

# vector to store block means
RMSEblk <- matrix(nrow = 1, ncol = length(li))
colnames(RMSEblk) <-li

for (b in 1:length(li)){
    l <- li[b]# block size
    m <- ceiling(t / l)                                 # number of blocks
    blk <- split(ts, rep(1:m, each=l, length.out = t))  # divides the series into blocks

    # initialize vector to receive result from for loop
    singleblock <- vector()                     
    for(i in 1:1000){
        res<-sample(blk, replace=T, 10000)        # resamples the blocks
        res.unlist<-unlist(res, use.names = F)    # unlist the bootstrap series
        # Split the series into train and test set
        train <- head(res.unlist, round(length(res.unlist) * 0.6))
        h <- length(res.unlist) - length(train)
        test <- tail(res.unlist, h)

        # Forecast for train set
        model <- auto.arima(train)
        future <- forecast(test, model=model,h=h)
        nfuture <- as.numeric(future$mean)        # makes the `future` object a vector            
        RMSE <- rmse(test, nfuture)               # use the `rmse` function from `Metrics` package

        singleblock[i] <- RMSE # Assign RMSE value to final result vector element i
    }

    RMSEblk[b] <- mean(singleblock) # store into matrix
}

RMSEblk

R 脚本实际运行，但需要超过 24 小时才能完成。 loops 中的运行次数（10000 和 1000）是使任务完美所需的最小值。

请问我该怎么做才能在更短的时间内完成脚本？

Answer 1

tl;dr你可能不得不以某种方式并行化它。

---

一个问题是您正在增长一个对象；也就是说，首先分配一个长度为零的向量 (singleblock <- vector())，然后一次将其递增一个元素 (singleblock[i] <- RMSE)。正如R Inferno 第 2 章中所讨论的，这是非常低效的。对于这个示例，它慢了 5 倍。

f1 <- function(x) { p <- numeric(0); for (i in 1:1000) p[i] <- 0 }
f2 <- function(x) { p <- numeric(1000); for (i in 1:1000) p[i] <- 0 }
microbenchmark(f1(),f2())
## Unit: microseconds
##  expr     min       lq      mean  median      uq     max neval cld
##  f1() 202.519 207.2105 249.84095 210.574 221.340 3504.95   100   b
##  f2()  40.274  40.6710  69.83741  40.9615  42.8275 2811.779   100  a 

但是：这并不重要。低效版本（增长向量）需要 210 微秒的中位时间。

microbenchmark(auto.arima(train),times=20L)
## Unit: milliseconds
##               expr      min       lq     mean   median       uq      max neval
##  auto.arima(train) 630.7335 648.3471 679.2703 657.6697 668.0563 829.1648    20

您的auto.arima() 调用大约需要 660 毫秒 - 大约长 3000 倍。对预测步骤使用类似的microbenchmark 调用会产生大约 20 毫秒的中值时间。

你可以做更正式的profiling，或者像这里显示的那样继续点点滴滴，但我没有在你的代码中看到任何看起来需要很长时间的东西（我可能会检查sample()下一个，但我怀疑它可以与auto.arima()相媲美。）

除非您能找到更快的 auto.arima() 版本（我对此表示怀疑），或者精简一些内容（例如限制搜索空间），否则您唯一剩下的选择就是并行化。您可以使用许多不同的工具在许多不同的级别上执行此操作，但首先要查看的是parallel option to auto.arima。您可能会选择并行化循环（在“R 中的并行计算”上进行网络搜索会提供大量资源）；请注意，尝试在多个级别上进行并行化可能会给您带来麻烦。

PS 粗略计算（48000 * 660 毫秒）大约需要 9 小时 - 只占大约 1/3 的时间（我原本预计它会达到 80% 左右）；也许你的处理器比我的慢？

Answer 2

为了演示，为避免循环中的对象增长，请考虑apply family 解决方案，例如vapply。请注意 RMSEblk 和 singleblock 现在是如何直接分配 vapply 的结果，而不需要按索引分配元素的簿记。

...

# DEFINED METHOD
proc_bootstrap <- function(b) {
    l <- li[b]                                          # block size
    m <- ceiling(t / l)                                 # number of blocks
    blk <- split(ts, rep(1:m, each=l, length.out = t))  # divides the series into blocks

    # initialize vector to receive result from for loop
    singleblock <- vapply(1:1000, function(i) {
      res <- sample(blk, replace=TRUE, 10000)        # resamples the blocks
      res.unlist <- unlist(res, use.names = FALSE)   # unlist the bootstrap series

      # Split the series into train and test set
      train <- head(res.unlist, round(length(res.unlist) * 0.6))
      h <- length(res.unlist) - length(train)
      test <- tail(res.unlist, h)

      # Forecast for train set
      model <- auto.arima(train)
      future <- forecast(test, model=model,h=h)
      nfuture <- as.numeric(future$mean)        # makes the `future` object a vector

      RMSE <- Metrics::rmse(test, nfuture)      # RETURN RMSE
    }, numeric(1))

    mean(singleblock)                           # RETURN MEAN
  }

# VAPPLY CALL
RMSEblk <- vapply(1:length(li), proc_bootstrap, numeric(1))

或者，填充您最初定义的单行矩阵（也许作为命名向量更好？）：

# MATRIX to store block means
RMSEblk <- matrix(nrow = 1, ncol = length(li))
colnames(RMSEblk) <-li

RMSEblk[] <- vapply(1:length(li), proc_bootstrap, numeric(1))

注意：上面的时间可能与嵌套的 for 循环没有本质区别，因为您仍会迭代 48,000 个模型调用。不过，这个解决方案可能在更大的迭代中可以更好地扩展。但正如所讨论的，请查看并行处理（参见parallel、doParallel、foreach 包），它可以从for 或apply 解决方案翻译。

---

请务必同时 profile 显示（在建模调用之外）unlist、head、tail 有时间问题：

utils::Rprof(tmp <- tempfile(), memory.profiling = TRUE)
RMSEblk <- vapply(1:length(li), proc_bootstrap, numeric(1))
utils::Rprof(NULL)
summaryRprof(tmp, memory="both")
unlink(tmp)