使用 openMP 和 Rcpp 设置 RNG 状态

Question

我有一个澄清问题。据我了解，sourceCpp 会自动传递 RNG 状态，因此 set.seed(123( 在调用 Rcpp 代码时为我提供了可重现的随机数。编译包时，我必须添加一个集合 RNG 语句。现在，这一切如何在源代码Cpp或软件包中与openMP一起工作？

考虑以下 Rcpp 代码

#include <Rcpp.h>
#include <omp.h>
// [[Rcpp::depends("RcppArmadillo")]]
// [[Rcpp::export]]
Rcpp::NumericVector rnormrcpp1(int n, double mu, double sigma  ){
  Rcpp::NumericVector out(n);
        for (int i=0; i < n; i++) {
          out(i) =R::rnorm(mu,sigma);
        }
  return(out);
}

// [[Rcpp::export]]
Rcpp::NumericVector rnormrcpp2(int n, double mu, double sigma, int cores=1  ){
  omp_set_num_threads(cores);
  Rcpp::NumericVector out(n);
    #pragma omp parallel for schedule(dynamic)   
        for (int i=0; i < n; i++) {
          out(i) =R::rnorm(mu,sigma);
        }
  return(out);
}

然后运行

   set.seed(123)
    a1=rnormrcpp1(100,2,3,2)
    set.seed(123)
    a2=rnormrcpp1(100,2,3,2)
    set.seed(123)
    a3=rnormrcpp2(100,2,3,2)
    set.seed(123)
    a4=rnormrcpp2(100,2,3,2)
    all.equal(a1,a2)
    all.equal(a3,a4)

虽然 a1 和 a2 相同，但 a3 和 a4 不是。如何使用 openMP 循环调整 RNG 状态？我可以吗？

Answer 1

为了扩展Dirk Eddelbuettel已经说过的内容，几乎不可能并行生成相同的PRN序列并具有所需的加速。其根源是PRN序列的生成本质上是一个顺序过程，其中每个状态都依赖于前一个状态，这创建了一个向后依赖链，该链可以追溯到初始种子状态。

此问题有两种基本解决方案。其中一个需要大量内存，另一个需要大量 CPU 时间，两者实际上更像是解决方法，而不是真正的解决方案：

预生成的 PRN 序列：一个线程按顺序生成一个巨大的 PRN 数组，然后所有线程都以与顺序情况一致的方式访问此数组。此方法需要大量内存才能存储序列。另一种选择是将序列存储到稍后进行内存映射的磁盘文件中。后一种方法的优点是可以节省一些计算时间，但通常 I/O 操作很慢，因此它只在处理能力有限或 RAM 量较小的计算机上有意义。

预缠 PRNG：当工作在线程之间静态分布时，例如使用 schedule(static) .每个线程都有自己的 PRNG，并且所有 PRNG 都使用相同的初始种子进行种子设定。然后，每个线程绘制与其起始迭代一样多的虚拟 PRN，基本上将其 PRNG 预卷到正确的位置。例如：

• 线程 0：绘制 0 个虚拟 PRN，然后绘制 100 个 PRN 并填充out(0:99)
线程 1：绘制 100 个
• 虚拟 PRN，然后绘制 100 个 PRN 并填充out(100:199)
• 线程 2：绘制 200 个虚拟 PRN，然后绘制 100 个 PRN 并填充out(200:299)

等等。当每个线程除了绘制 PRN 之外还要进行大量计算时，此方法效果很好，因为在某些情况下(例如，多次迭代(，预缠 PRNG 的时间可能很长。

对于除了绘制 PRN 之外还有大量数据处理的情况，存在第三个选项。这个使用 OpenMP 有序循环(请注意，迭代块大小设置为 1(：

#pragma omp parallel for ordered schedule(static,1)
for (int i=0; i < n; i++) {
  #pragma omp ordered
  {
     rnum = R::rnorm(mu,sigma);
  }
  out(i) = lots of processing on rnum
}

尽管循环排序本质上序列化了计算，但它仍然允许并行执行lots of processing on rnum，因此可以观察到并行加速。请参阅此答案以更好地解释为什么会这样。

Answer 2

是的，sourceCpp() 等，并且是RNGScope的实例化，因此 RNG 保持正确的状态。

是的，人们可以做OpenMP。 但是在 OpenMP 段中，您无法控制线程的执行顺序 - 因此您延长了相同的序列。我对正在开发的软件包有同样的问题，我希望在使用 OpenMP 的同时进行可重现的绘图。 但似乎你不能。