如何提高双滚动window操作的效率？
How to increase efficiency of double rolling window operation?
有没有人对如何提高以下使用 "kind-of" 双滚动 window 占用我所有 ram 的代码示例的效率有想法或建议？
首先，我通过一个简单的例子来定义问题，在这个 post 的底部有一个完整的 MWE（实现）。
首先，考虑以下"random"测试向量（通常长度>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 被分成 "kind-of" 训练集和测试集，两者都滚动 windows。在此 MWE 中，考虑了长度为 4 的训练集开始和 2 的测试集长度（通常长度 >200）。所以最初，以下值是训练和测试集的一部分：
train_1 <- A[1:4]
test_1 <- A[5:6]
接下来，我想在 train_1 的每个可能的连续位置（因此第一次滚动 window）从 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,]))
})
在下一步中，"kind-of" 训练集和测试集通过将它们的顺序从 A 增加一个来更新（因此第二次滚动 window）：
train_2 <- A[2:5]
test_2 <- A[6:7]
如前所述，在 train_2 中的每个可能位置减去 test_2，然后计算 run_2_sub 和 run_2_sum。这个过程一直持续到测试集代表 A 的最后两个值，最后我（在这个 MWE 中）以 6 run_sum 矩阵结束。但是，我的实现速度很慢，我想知道是否有人可以帮助我提高它的效率？
这是我的实现：
# 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。
关于如何提高和改进速度和代码有什么建议吗？
通常R中代码优化的前两步是：
- 少做点；
- 使用矢量化。
我们将完成这两个步骤。让我们同意将 x 作为输入向量（在您的示例中为 A）。
您的问题中的关键功能单元可以表述如下：给定 train_start（[=15= 的子集的起始索引]。我们将为该子集使用词 '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 中的 'Preparatory computations'。
res 现在包含 x_mod 所需的差异总和（它是 x 的副本，但用 0 而不是 NAs 和 NaN秒）。我们现在应该减去所有过度使用的元素，即那些我们不应该在求和中使用的元素，因为其他集合中的相应元素是 NA 或 NaN。在计算此信息时，我们还将计算分母。请参阅 run_mean_diff 中的 'Info about extra elements'。
这段代码的美妙之处在于 train_start、test_start 和 test_length 现在可以是向量：每个向量的第 i 个元素被视为单个元素我们的任务。这就是矢量化。我们现在的工作是构建适合我们任务的这些向量。参见函数 generate_run_data.
提供的代码使用更少的 RAM，不需要额外的 zoo 依赖，并且在小型 train_length 和 test_length 上比原始代码快得多。在大 *_length 上也更快，但不是很多。
接下来的步骤之一可能是使用 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
)
您的代码使用如此多 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
下面是 的时间：
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
有没有人对如何提高以下使用 "kind-of" 双滚动 window 占用我所有 ram 的代码示例的效率有想法或建议？
首先，我通过一个简单的例子来定义问题，在这个 post 的底部有一个完整的 MWE（实现）。
首先，考虑以下"random"测试向量（通常长度>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 被分成 "kind-of" 训练集和测试集，两者都滚动 windows。在此 MWE 中，考虑了长度为 4 的训练集开始和 2 的测试集长度（通常长度 >200）。所以最初，以下值是训练和测试集的一部分：
train_1 <- A[1:4]
test_1 <- A[5:6]
接下来，我想在 train_1 的每个可能的连续位置（因此第一次滚动 window）从 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,]))
})
在下一步中，"kind-of" 训练集和测试集通过将它们的顺序从 A 增加一个来更新（因此第二次滚动 window）：
train_2 <- A[2:5]
test_2 <- A[6:7]
如前所述，在 train_2 中的每个可能位置减去 test_2，然后计算 run_2_sub 和 run_2_sum。这个过程一直持续到测试集代表 A 的最后两个值，最后我（在这个 MWE 中）以 6 run_sum 矩阵结束。但是，我的实现速度很慢，我想知道是否有人可以帮助我提高它的效率？
这是我的实现：
# 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。
关于如何提高和改进速度和代码有什么建议吗？
通常R中代码优化的前两步是：
- 少做点；
- 使用矢量化。
我们将完成这两个步骤。让我们同意将 x 作为输入向量（在您的示例中为 A）。
您的问题中的关键功能单元可以表述如下：给定 train_start（[=15= 的子集的起始索引]。我们将为该子集使用词 '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 中的 'Preparatory computations'。
res 现在包含 x_mod 所需的差异总和（它是 x 的副本，但用 0 而不是 NAs 和 NaN秒）。我们现在应该减去所有过度使用的元素，即那些我们不应该在求和中使用的元素，因为其他集合中的相应元素是 NA 或 NaN。在计算此信息时，我们还将计算分母。请参阅 run_mean_diff 中的 'Info about extra elements'。
这段代码的美妙之处在于 train_start、test_start 和 test_length 现在可以是向量：每个向量的第 i 个元素被视为单个元素我们的任务。这就是矢量化。我们现在的工作是构建适合我们任务的这些向量。参见函数 generate_run_data.
提供的代码使用更少的 RAM，不需要额外的 zoo 依赖，并且在小型 train_length 和 test_length 上比原始代码快得多。在大 *_length 上也更快，但不是很多。
接下来的步骤之一可能是使用 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
)
您的代码使用如此多 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
下面是
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