在迭代算法中使用 Rcpp 加速替换列表和向量的元素是不是合法？

Posted 2023-04-18

【中文标题】在迭代算法中使用 Rcpp 加速替换列表和向量的元素是不是合法？

【英文标题】：Is it legitimate to use Rcpp to speed up replacing elements of lists and vectors in iterative algorithm?在迭代算法中使用 Rcpp 加速替换列表和向量的元素是否合法？

【发布时间】：2018-09-04 12:36:41

【问题描述】：

上下文

我最近一直在研究迭代算法，其中每次迭代 n 取决于迭代 n-1。在每次迭代期间，大部分计算时间是通过子设置和/或替换向量、列表或 data.tables (N > 10^6) 的元素来花费的。

我最近遇到了 Rcpp 并用它玩了一下，我发现替换向量或列表的元素 k 可以加快两个或三个数量级（下面很少有基准测试）。

但是，当在 for 和 while 循环中使用 Rcpp 子集代码时，R 似乎变得不稳定，会话在随机点中止，没有任何问题的提示。

问题

我的问题： 这种使用 Rcpp 是否合法，或者它会导致我不知道的问题？

示例

下面是我正在使用的 Rcpp 代码和一些基准测试。总体而言，该算法应该调用替换函数约 55 亿次和子集函数约 500 亿次。

请注意，使用 Rcpp 替换列表和双向量的元素会更快，而对于整数向量，基本 R 解决方案是首选（基准 1）；数据表是替换元素的好选择，但如果您必须重复子集化以访问其元素，则向量方法要快得多（benchmark 2）。

功能：

#include <Rcpp.h>
using namespace Rcpp;


// [[Rcpp::export]]

void assign_list(List x, int k, NumericVector new_element)
  x[k-1] = new_element;


// [[Rcpp::export]]
void assign_dbl(NumericVector x, int k, double y)
  x[k-1] = y;


// [[Rcpp::export]]
void assign_int(IntegerVector x, int k, int y)
  x[k-1] = y;

基准测试：

输入

set.seed(747474)

int <- 1:10^7
dou <- rnorm(10^7, 1000, 300)
l   <- lapply(sample(5:20, 10^7,  replace = T), rnorm, mean = 1000, sd = 300)
dt  <- data.table(int = int, dou = dou, l = l)

i <- 999999
z <- 2222
k <- 30000
s <- 552877

1)

Unit: nanoseconds
                                     expr      min       lq        mean      median        uq       max neval
                             int[i] <- -1L     488     2439  36938108.9      4146.0  15651119 799720107    30
                             dou[i] <- -1      732     3170  19101960.4   6609193.5  16187500 212369197    30
                             l[i]   <- -1      489     3902 159442538.1 186035314.5 227131872 618326686    30
                               assign_int 19853910 22028692  81055363.5  24665494.0  39352345 872241539    30
                               assign_dbl     1220     5852     48023.2      8534.5     13167   1158828    30
                              assign_list     1464     6828     52866.9     10850.5     13411   1243430    30
 dt[k, ':=' (int = -1, dou = -1, l = -1)]   206020   340116    481850.0    425326.5    529312   1414341    30

2)

microbenchmark(times = 30L,

               "subset vector + list" = int[s]; dou[s]; l[s],
               "subset datatable"     = dt[s, int]; dt[s, dou]; dt[s, l])

Unit: nanoseconds
                 expr    min     lq       mean   median     uq     max neval
 subset vector + list    488    488   1715.533   1585.5   2926    4389    30
     subset datatable 563688 574417 719304.467 600138.0 875765 1308040    30

【问题讨论】：

从[] 切换到() 以启用边界检查。任何“不稳定”都可能是触发未定义行为 (UB) 导致崩溃的“越界”错误。确保您考虑到 C++ 索引从 0 开始而不是 1 的事实。

@coatless，谢谢。一旦我更改为x(k-1) = new_element，如果出现“越界”错误会发生什么？

代码将出错，您将收到一条漂亮的新 Rcpp 错误消息，说明已超出范围。

【参考方案1】：

这是非常危险的，因为这里显示了副作用，即使您只将 x 传递给 Rcpp 函数，x 和 y 也会发生变化

> x <- y <- 1:10
> assign_int(x, 1, 2)
> y
 [1]  2  2  3  4  5  6  7  8  9 10

似乎没有更快；对于这些功能

f0 <- function(x) 
    for (i in seq_along(x))
        x[i] = -i


f1 <- function(x) 
    for (i in seq_along(x))
        assign_int(x, i, -i)

我有

> int <- 1:10^5
> microbenchmark(f0(int), f1(int), times=5)
Unit: milliseconds
    expr       min        lq      mean    median        uq       max neval
 f0(int)  14.78777  14.80264  15.05683  15.03138  15.17678  15.48556     5
 f1(int) 659.67346 669.00095 672.93343 670.48917 676.16930 689.33429     5

在您的基准测试中，int[i] <- 1, '1' 是一个数字（双精度）值，因此您将强制转换为双精度向量（在赋值后查看 class(int)），从而触发完整副本。使用int[i] <- 1L 强制右侧为整数。

> int0 <- int1 <- 1:10^7
> microbenchmark(int0[1] <- 1, int1[1] <- 1L)
Unit: microseconds
          expr   min    lq      mean median     uq       max neval
  int0[1] <- 1 1.047 1.102 1770.9911  1.143 1.2650 176960.52   100
 int1[1] <- 1L 1.105 1.176  339.4264  1.213 1.2655  33815.97   100
> class(int0)
[1] "numeric"
> class(int1)
[1] "integer"

在基础 R 中仅更新单个元素作为基准是昂贵的，因为它会在每次分配时触发向量的副本；在f0() 中，副本只发生一次。在第一次迭代中，R 制作了一个副本，因为它知道整数值的向量至少被两个符号（函数 int 的参数和函数 x 中使用的符号）引用，所以它制作了一个内存的副本并将其分配给函数内部的x。它这样做是为了避免在您的 Rcpp 代码中看到的副作用（即，避免修改 int）。下一次循环时，R 会识别出只有一个符号引用了该向量，因此在不复制的情况下进行替换。

注意

> int <- 1:10^5
> f1(int)
> head(int)
[1] -1 -2 -3 -4 -5 -6

说明 Rcpp 代码的副作用可能会产生意想不到的结果。

此外，有几种方法可以减慢迭代循环，例如，

f2 <- function(x) 
    for (i in seq_along(x)) 
        x[i] = -i
        y <- x
    


f3 <- function(x) 
    result <- integer()
    for (i in seq_along(x))
        result <- c(result, -i)

与

> int <- 1:10^3
> microbenchmark(f0(int), f2(int), f3(int), times = 1)
Unit: microseconds
    expr      min       lq     mean   median       uq      max neval
 f0(int)  150.507  150.507  150.507  150.507  150.507  150.507     1
 f2(int)  667.201  667.201  667.201  667.201  667.201  667.201     1
 f3(int) 4379.005 4379.005 4379.005 4379.005 4379.005 4379.005     1

f2() 导致 R 每次通过循环复制x（以避免修改y 的副作用）。 f3() 在连续迭代中复制长度为 0, 1, 2, 3, ... n - 1（其中n = length(x)）的向量，导致复制n * (n - 1) / 2 元素，并且算法缩放为x 的长度。

一般原则也适用于其他类型，包括带有

的列表

f0l <- function(x) 
    for (i in seq_along(x))
        x[[i]] <- i
    x


f1l <- function(x) 
    for (i in seq_along(x))
        assign_list(x, i, i)

导致

> set.seed(123)
> l0   <- lapply(sample(5:20, 10^6,  replace = T), rnorm, mean = 1000, sd = 300)
> set.seed(123)
> l1   <- lapply(sample(5:20, 10^6,  replace = T), rnorm, mean = 1000, sd = 300)
> microbenchmark(f0l(l0), f1l(l1), times=1)
Unit: milliseconds
    expr       min        lq      mean    median        uq       max neval
 f0l(l0)  239.9865  239.9865  239.9865  239.9865  239.9865  239.9865     1
 f1l(l1) 6767.9172 6767.9172 6767.9172 6767.9172 6767.9172 6767.9172     1

【讨论】：

亲爱的马丁非常感谢您的回复。事实上，assign_int 并不比 base R 快（事实上我并没有在我的实际脚本中使用它——我发布它是为了完整性）。我确实在基准测试中使用了 -1L，尽管标签具有误导性（我现在已经更正了）。 您所说的副作用正是我所寻找的，而且，特别是在列表情况下以及双向量情况下，使代码更快，因为据我了解没有复制。重点是我必须在每次迭代时访问和更新单个元素。

我希望我们意见一致； f0() 每次迭代更新单个元素，比 Rcpp 解决方案更快。

确实，assign_int 并不比 base R 快（事实上，我并没有在我的实际脚本中使用它——为了完整起见，我发布了它）。对于列表和双向量，Rcpp 确实加快了我的代码速度，并且由于我的意图是修改 dou 和 l，因此您描述的副作用实际上是这些函数的唯一意图。

再一次，在 R 不会进行过多复制的情况下编写代码（f0()，现在也是 f0l()）也很简单，例如，请参阅我回复末尾的添加。此外，副作用也不是一件好事。您可能认为您已经了解它何时重要，但您几乎肯定会在您的代码中引入难以理解的行为。

澄清一下，我提到的“副作用”是在我回复的第一行中对y 的意外修改，而不是您对x 的预期修改。