1

有没有人对如何提高以下代码示例的效率有什么想法或建议？

首先，我通过一个简单的示例来定义问题，在这篇文章的底部有一个完整的 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]

接下来，我要减去test_1从train_1每个可能的连续位置train_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分别为run_sum120.04 秒和 28.69 秒。

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

4

2 回答 2

3

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

  • 少做；
  • 使用矢量化。

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

问题中的关键功能单元可以表述如下：给定train_start（ 的子集的起始索引train。我们将为此子集使用单词“train”），test_start（的起始索引test）和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）。我们现在应该减去所有过度使用的元素，即那些我们不应该在总和中使用的元素，因为其他集合中的相应元素是NAor NaN。在计算这些信息时，我们还将计算分母。请参阅 中的“有关额外元素的信息” run_mean_diff。

这段代码的美妙之处在于train_start,test_start并且test_length现在可以是向量：i每个向量的第一个元素被视为我们任务的单个元素。这就是矢量化。我们现在的工作是构建适合我们任务的这些向量。见功能generate_run_data。

呈现的代码使用更少的 RAM，不需要额外的zoo依赖，并且在小型train_length和test_length. 在 big *_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
)
于 2017-06-11T16:35:50.807 回答
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 方法的时间安排：

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
于 2017-06-12T13:25:00.257 回答