1

我有一个澄清问题。据我了解，sourceCpp 会自动传递 RNG 状态，因此 set.seed(123) 在调用 Rcpp 代码时会为我提供可重现的随机数。编译包时，我必须添加一个 set RNG 语句。现在，这一切如何在 sourceCpp 或包中与 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 状态？我可以吗？

4

2 回答 2

2

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

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

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

prewound 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并行执行，因此可以观察到并行加速。请参阅此答案以更好地解释为什么会这样。

于 2013-09-24T09:08:20.203 回答
1

是的，sourceCpp()等等，以及RNGScopeRNG 的一个实例，它们都处于适当的状态。

是的，可以做 OpenMP。但是在 OpenMP 段内部，您无法控制线程执行的顺序——因此您需要更长的相同顺序。我对正在开发的包有同样的问题，我希望有可重现的绘图但使用 OpenMP。但似乎你不能。

于 2013-09-23T20:17:46.560 回答