如何提高双滚动窗口操作的效率？

Question

对于如何提高以下使用“某种”双滚动窗口耗尽我所有内存的代码示例的效率，有人有什么想法或建议吗？

首先，我通过一个简单的示例来定义问题，在这篇文章的底部有一个完整的 MWE（实现）。

---

首先，考虑以下“随机”测试向量（通常长度 >25000）：

A <- c(1.23,5.44,6.3,8.45,NaN,3.663,2.63,1.32,6.623,234.6,252.36)

A 被分割成“种类”的训练和测试集，都带有滚动窗口。在这个 MWE 中，考虑了长度为 4 的训练集起点和长度为 2 的测试集长度（通常长度 >200）。所以最初，以下值是训练集和测试集的一部分：

train_1 <- A[1:4]
test_1 <- A[5:6]

接下来，我想在train_1（因此是第一个滚动窗口）的每个可能连续位置从train_1 中减去test_1，生成run_1_sub 矩阵。

run_1_sub <- matrix(NaN,3,2)
run_1_sub[1,] <- train_1[1:2] - test_1
run_1_sub[2,] <- train_1[2:3] - test_1
run_1_sub[3,] <- train_1[3:4] - test_1

之后，我想在run_1_sub 的每一行上找到每一行的总和除以每行中不是NaN 的条目数。

run_1_sum <-
    sapply(1:3, function(x) {
       sum(run_1_sub[x,], na.rm = T) / sum(!is.na(run_1_sub[x,]))
})

在下一步中，“种类”的训练和测试集通过将它们从 A 的顺序增加一来更新（因此是第二个滚动窗口）：

train_2 <- A[2:5] 
test_2 <- A[6:7]  

如前所述，test_2 在train_2 中的每个可能位置减去run_2_sub 和run_2_sum。这个过程一直持续到测试集代表 A 的最后两个值，最后我以 6 个run_sum 矩阵结束（在这个 MWE 中）。但是，我的实现非常缓慢，我想知道是否有人可以帮助我提高它的效率？

---

这是我的实现：

# Initialization
library(zoo) 
#rm(list = ls())
A <- c(1.23, 5.44, 6.3, 8.45, NaN, 3.663, 2.63, 1.32, 6.623, 234.6, 252.36) # test vector
train.length <- 4
test.length <- 2
run.length <- length(A) - train.length - test.length + 1
# Form test sets
test.sets <- sapply(1:run.length, function(x) {
A[(train.length + x):(train.length + test.length + x - 1)]
})
# Generate run_sub_matrices
run_matrix <- lapply(1:run.length, function(x) {
    rollapply(A[x:(train.length + x - 1)], width = test.length, by = 1,
        function(y) {
            y - test.sets[, x]
            })
})
# Genereate run_sum_matrices
run_sum <- sapply(1:length(run_matrix), function(x) {
rowSums(run_matrix[[x]], na.rm = T) / apply(run_matrix[[x]], 1,  function(y) {
sum(!is.na(y))})
})

当然，以下初始化设置会显着减慢 run_sum 和 run_sub 的生成速度：

A <- runif(25000)*400
train.length <- 400
test.length <- 200

这里，生成run_sub的时间分别为120.04s和run_sum 28.69s。

关于如何提高和改进速度和代码的任何建议？

Answer 1

通常R中代码优化的前两步是：

• 少做事；
• 使用矢量化。

我们将完成这两个步骤。让我们同意将x 记为输入向量（在您的示例中为A）。

您的问题中的关键功能单元可以表述如下：给定train_start（train 子集的起始索引。我们将为此子集使用单词“train”），test_start（@ 的起始索引987654326@) 和test_length（test 的长度）计算：

train_inds <- train_start + 0:(test_length-1)
test_inds <- test_start + 0:(test_length-1)
run_diff <- x[train_inds] - x[test_inds]
sum(run_diff, na.rm = TRUE) / sum(!is.na(run_diff))

这个单元被多次调用，总和和!is.na 的计算也是如此。我们将少做：我们预先计算累积和并使用这些数据，而不是用它们的总和计算多次差。请参阅run_mean_diff 中的“准备计算”。

res 现在包含所需的x_mod 差异总和（这是x 的副本，但用0 代替NAs 和NaNs）。我们现在应该减去所有过度使用的元素，即那些我们不应该在总和中使用的元素，因为其他集合中的相应元素是NA 或NaN。在计算这些信息时，我们还将计算分母。请参阅run_mean_diff 中的“有关额外元素的信息”。

这段代码的美妙之处在于train_start、test_start 和test_length 现在可以成为向量：每个向量的ith 元素被视为我们任务的单个元素。这就是矢量化。我们现在的工作是构建适合我们任务的这些向量。见函数generate_run_data。

呈现的代码使用更少的 RAM，不需要额外的 zoo 依赖，并且在小型 train_length 和 test_length 上的原始代码要快得多。在大 *_lengths 上也更快，但不是很多。

接下来的步骤之一可能是使用 Rcpp 编写此代码。

代码：

run_mean_diff <- function(x, train_start, test_start, test_length) {
  # Preparatory computations
  x_isna <- is.na(x)
  x_mod <- ifelse(x_isna, 0, x)
  x_cumsum <- c(0, cumsum(x_mod))

  res <- x_cumsum[train_start + test_length] - x_cumsum[train_start] -
    (x_cumsum[test_start + test_length] - x_cumsum[test_start])

  # Info about extra elements
  extra <- mapply(
    function(cur_train_start, cur_test_start, cur_test_length) {
      train_inds <- cur_train_start + 0:(cur_test_length-1)
      test_inds <- cur_test_start + 0:(cur_test_length-1)

      train_isna <- x_isna[train_inds]
      test_isna <- x_isna[test_inds]

      c(
        # Correction for extra elements
        sum(x_mod[train_inds][test_isna]) -
              sum(x_mod[test_inds][train_isna]),
        # Number of extra elements
        sum(train_isna | test_isna)
      )
    },
    train_start, test_start, test_length, SIMPLIFY = TRUE
  )

  (res - extra[1, ]) / (test_length - extra[2, ])
}

generate_run_data <- function(n, train_length, test_length) {
  run_length <- n - train_length - test_length + 1
  num_per_run <- train_length - test_length + 1

  train_start <- rep(1:num_per_run, run_length) +
    rep(0:(run_length - 1), each = num_per_run)
  test_start <- rep((train_length + 1):(n - test_length + 1),
                    each = num_per_run)

  data.frame(train_start = train_start,
             test_start = test_start,
             test_length = rep(test_length, length(train_start)))
}

A <- c(1.23, 5.44, 6.3, 8.45, NaN, 3.663,
       2.63, 1.32, 6.623, 234.6, 252.36)
train_length <- 4
test_length <- 2
run_data <- generate_run_data(length(A), train_length, test_length)

run_sum_new <- matrix(
  run_mean_diff(A, run_data$train_start, run_data$test_start,
                run_data$test_length),
  nrow = train_length - test_length + 1
)

Answer 2

您的代码使用这么多 RAM 的原因是您保留了很多中间对象，主要是 run_matrix 中的所有元素。通过Rprof 进行的分析表明，大部分时间都花在了rollapply。

避免所有中间对象的最简单和最简单的方法是使用 for 循环。它还使代码清晰。然后你只需要用更快的方式替换对rollapply 的调用。

您要应用于每个滚动子集的函数很简单：减去测试集。您可以使用stats::embed 函数创建滞后矩阵，然后利用 R 的回收规则从每一列中减去测试向量。我创建的函数是：

calc_run_sum <- function(A, train_length, test_length) {
  run_length <- length(A) - train_length - test_length + 1L
  window_size <- train_length - test_length + 1L

  # Essentially what embed() does, but with column order reversed
  # (part of my adaptation of echasnovski's correction)
  train_lags <- 1L:test_length +
                rep.int(1L:window_size, rep.int(test_length, window_size)) - 1L
  dims <- c(test_length, window_size)  # lag matrix dims are always the same

  # pre-allocate result matrix
  run_sum <- matrix(NA, window_size, run_length)

  # loop over each run length
  for (i in seq_len(run_length)) {
    # test set indices and vector
    test_beg <- (train_length + i)
    test_end <- (train_length + test_length + i - 1)

    # echasnovski's correction
    #test_set <- rep(test_set, each = train_length - test_length + 1)
    #lag_matrix <- embed(A[i:(test_beg - 1)], test_length)
    #run_sum[,i] <- rowMeans(lag_matrix - test_set, na.rm = TRUE)

    # My adaptation of echasnovski's correction
    # (requires train_lags object created outside the loop)
    test_set <- A[test_beg:test_end]
    train_set <- A[i:(test_beg - 1L)]
    lag_matrix <- train_set[train_lags]
    dim(lag_matrix) <- dims
    run_sum[,i] <- colMeans(lag_matrix - test_set, na.rm = TRUE)
  }
  run_sum
}

现在，进行一些基准测试。我使用了以下输入数据：

library(zoo) 
set.seed(21)
A <- runif(10000)*200
train.length <- 200
test.length <- 100

以下是您最初方法的时间安排：

system.time({
  run.length <- length(A) - train.length - test.length + 1
  # Form test sets
  test.sets <- sapply(1:run.length, function(x) {
    A[(train.length + x):(train.length + test.length + x - 1)]
  })
  # Generate run_sub_matrices
  run_matrix <- lapply(1:run.length, function(x) {
    rm <- rollapply(A[x:(train.length + x - 1)], width = test.length, by = 1,
                    FUN = function(y) { y - test.sets[, x] })
  })
  # Genereate run_sum_matrices
  run_sum <- sapply(run_matrix, function(x) {
    rowSums(x, na.rm = T) / apply(x, 1,  function(y) {
  sum(!is.na(y))})
  })
})
#    user  system elapsed 
#  19.868   0.104  19.974 

这是echasnovski's approach的时间安排：

system.time({
  run_data <- generate_run_data(length(A), train.length, test.length)

  run_sum_new <- matrix(
    run_mean_diff(A, run_data$train_start, run_data$test_start,
                  run_data$test_length),
    nrow = train.length - test.length + 1
  )
})
#    user  system elapsed 
#  10.552   0.048  10.602 

以及我的方法的时间安排：

system.time(run_sum_jmu <- calc_run_sum(A, train.length, test.length))
#    user  system elapsed 
#   1.544   0.000   1.548 

所有 3 种方法的输出都是相同的。

identical(run_sum, run_sum_new)
# [1] TRUE
identical(run_sum, run_sum_jmu)
# [1] TRUE