«Эти сэмплеры нельзя использовать в распараллеленном коде»

Question

Я читал виньетка для пакета rgen , который предоставляет заголовки для выборки из некоторых распространенных дистрибутивов. В первом абзаце говорится, что:

Обратите внимание, что эти сэмплеры, так же как и в armadillo, не могут использоваться в распараллеленном коде, так как подпрограммы генерации основаны на однопоточных вызовах R.

Это было для меня новостью, и я уже довольно давно использую RcppArmadillo. Мне было интересно, если кто-то может уточнить этот вопрос (или предоставить ссылки, где я могу прочитать о проблеме). Мне особенно интересно узнать, что здесь означает «нельзя использовать»; результаты будут неправильными или просто не будут распараллеливаться?

Answer 1

В этих функциях используется генератор случайных чисел R, который нельзя использовать в распараллеленном коде, поскольку это приводит к неопределенному поведению. Неопределенное поведение может привести практически к чему угодно. С моей точки зрения, вам повезло в случае сбоя программы, поскольку это ясно говорит о том, что что-то идет не так.

В представлении задачи HPC перечислены некоторые ГСЧ, которые подходят для параллельных вычислений. Но вы не можете легко использовать их с дистрибутивами, предоставленными rgen или RcppDist . Вместо этого можно сделать следующее:

• Функция копирования для многомерного нормального распределения из rgen настраивает его сигнатуру так, чтобы он принимал std::function<double()> в качестве источника для N(0, 1) распределенных случайных чисел.
• Используйте быстрый RNG вместо R's 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 

Здесь реальное улучшение производительности достигается за счет использования более быстрого ГСЧ, что понятно, поскольку здесь основное внимание уделяется множеству случайных чисел, а не матричным операциям. Если я перейду больше к матричным операциям, используя 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). Таким образом, есть немало факторов, которые необходимо учитывать, прежде чем выбирать метод.