r - "These samplers cannot be used in parallelized code"

Question

我正在阅读 vignette为 rgen包提供了从一些常见分布中采样的头文件。在第一段中，它说:

Please note, these samplers, just like the ones in armadillo cannot be used in parallelized code as the underlying generation routines rely upon R calls that are single-threaded.

这对我来说是个新闻，而且我已经使用 RcppArmadillo 有一段时间了。我想知道是否有人可以详细说明这一点(或提供有关我可以在哪里阅读该问题的引用)。我特别有兴趣了解这里的“不能使用”是什么意思；结果会是错误的，还是不会并行化？

Answer 1

这些函数使用 R 的随机数生成器，不能在并行代码中使用，因为这会导致未定义的行为。未定义的行为几乎可以导致任何事情。从我的角度来看，如果程序崩溃，你很幸运，因为这清楚地告诉你出了问题。

HPC task view列出了一些适合并行计算的 RNG。但是您不能通过 rgen 或 RcppDist 提供的发行版轻松使用它们。 .相反，可以执行以下操作:

来自 rgen 的多元正态分布的复制函数调整它的签名，使其需要 std::function<double()>作为 N(0, 1) 的来源分布式随机数。

使用快速 RNG 而不是 R 的 RNG。

在并行模式下使用相同的快速 RNG。

在代码中作为一个快速黑客:

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(dqrng)]]
#include <xoshiro.h>
#include <dqrng_distribution.h>
// [[Rcpp::plugins(openmp)]]
#include <omp.h>

inline arma::mat rmvnorm(unsigned int n, const arma::vec& mu, const arma::mat& S,
                         std::function<double()> rnorm = norm_rand){
  unsigned int ncols = S.n_cols;
  arma::mat Y(n, ncols);
  Y.imbue( rnorm ) ;
  return arma::repmat(mu, 1, n).t() + Y * arma::chol(S);
}

// [[Rcpp::export]]
arma::mat defaultRNG(unsigned int n, const arma::vec& mu, const arma::mat& S) {
  return rmvnorm(n, mu, S);
}

// [[Rcpp::export]]
arma::mat serial(unsigned int n, const arma::vec& mu, const arma::mat& S) {
  dqrng::normal_distribution dist(0.0, 1.0);
  dqrng::xoshiro256plus rng(42);
  return rmvnorm(n, mu, S, [&](){return dist(rng);});
}

// [[Rcpp::export]]
std::vector<arma::mat> parallel(unsigned int n, const arma::vec& mu, const arma::mat& S, unsigned int ncores = 1) {
  dqrng::normal_distribution dist(0.0, 1.0);
  dqrng::xoshiro256plus rng(42);
  std::vector<arma::mat> res(ncores);

  #pragma omp parallel num_threads(ncores)
  {
    dqrng::xoshiro256plus lrng(rng);      // make thread local copy of rng 
    lrng.jump(omp_get_thread_num() + 1);  // advance rng by 1 ... ncores jumps 
    res[omp_get_thread_num()] = rmvnorm(n, mu, S, [&](){return dist(lrng);});
  }
  return res;
}


/*** R
set.seed(42)
N <- 1000000
M <- 100
mu <- rnorm(M)
S <- matrix(rnorm(M*M), M, M)
S <- S %*% t(S)
system.time(defaultRNG(N, mu, S))
system.time(serial(N, mu, S))
system.time(parallel(N/2, mu, S, 2))
*/

结果:

> system.time(defaultRNG(N, mu, S))
   user  system elapsed 
  6.984   1.380   6.881 

> system.time(serial(N, mu, S))
   user  system elapsed 
  4.008   1.448   3.971 

> system.time(parallel(N/2, mu, S, 2))
   user  system elapsed 
  4.824   2.096   3.080 

这里真正的性能改进来自使用更快的 RNG，这是可以理解的，因为这里的重点在于许多随机数而不是矩阵运算。如果我通过使用 N <- 100000 更多地转向矩阵运算和 M <- 1000我得到:

> system.time(defaultRNG(N, mu, S))
   user  system elapsed 
 16.740   1.768   9.725 

> system.time(serial(N, mu, S))
   user  system elapsed 
 13.792   1.864   6.792 

> system.time(parallel(N/2, mu, S, 2))
   user  system elapsed 
 14.112   3.900   5.859 

在这里我们清楚地看到，在 全部 如果用户时间大于耗时。这样做的原因是我使用的并行 BLAS 实现 (OpenBLAS)。因此，在决定一种方法之前，需要考虑很多因素。