Параллельные выборки из случайного нормального распределения - не быстрее?

Question

Я использую R для создания симуляции, которая берет выборки из случайного нормального распределения, и неудивительно, что это довольно медленно.Итак, я искал несколько способов ускорить его с помощью Rcpp и наткнулся на пакет RcppZiggurat для более быстрых случайных нормальных выборок и пакет RcppParallel для многопоточных вычислений, и подумалпочему бы не использовать более быстрый алгоритм и рисовать сэмплы параллельно?

Итак, я приступил к созданию прототипов и в итоге получил три метода сравнения:

1. Сэмплы с использованием RcppParallel и RcppZigguratвместе
2. Сэмплы с использованием только RcppZiggurat
3. Сэмплы с использованием старого доброго rnorm

Ниже приведены мои реализации с использованием RcppParallel + RcppZiggurat (функция parallelDraws), ипросто RcppZiggurat (функция serialDraws):

#include <Rcpp.h>
// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
// [[Rcpp::depends(RcppZiggurat)]]
#include <Ziggurat.h>

static Ziggurat::Ziggurat::Ziggurat zigg;

using namespace RcppParallel;

struct Norm : public Worker
{   
  int input;

  // saved draws
  RVector<double> draws;

  // constructors
  Norm(const int input, Rcpp::NumericVector draws)
    : input(input), draws(draws) {}

  void operator()(std::size_t begin, std::size_t end) {
    for (std::size_t i = begin; i < end; i++) {
      draws[i] = zigg.norm();
    }
  }
};

// [[Rcpp::export]]
Rcpp::NumericVector parallelDraws(int x) {

  // allocate the output vector
  Rcpp::NumericVector draws(x);

  // declare the Norm instance 
  Norm norm(x, draws);

  // call parallelFor to start the work
  parallelFor(0, x, norm);

  // return the draws
  return draws;
};

// [[Rcpp::export]]
Rcpp::NumericVector serialDraws(int x) {

  // allocate the output vector
  Rcpp::NumericVector draws(x);

  for (int i = 0; i < x; i++) {
    draws[i] = zigg.norm();
  }

  // return the draws
  return draws;
};

Когда я сравнил их, я обнаружил некоторые удивительные результаты:

library(microbenchmark)
microbenchmark(parallelDraws(1e5), serialDraws(1e5), rnorm(1e5))

Unit: microseconds
                 expr      min       lq     mean    median       uq        max neval
 parallelDraws(1e+05) 3113.752 3539.686 3687.794 3599.1540 3943.282   5058.376   100
   serialDraws(1e+05)  695.501  734.593 2536.940  757.2325  806.135 175712.496   100
         rnorm(1e+05) 6072.043 6264.030 6655.835 6424.0195 6661.739  18578.669   100

Использование только RcppZiggurat было примерно в 8 раз быстрее, чемrnorm, но использование RcppParallel и RcppZiggurat вместе было медленнее, чем использование только RcppZiggurat!Я попытался поиграться с размером зерна для функции RcppParallel ParallelFor, но это не привело к какому-либо заметному улучшению.

Мой вопрос: что может быть причиной, по которой добавлениепараллелизм на самом деле здесь хуже?Я знаю, что «накладные расходы» в параллельных вычислениях могут перевесить преимущества в зависимости от различных факторов.Это то, что здесь происходит?Или я совершенно не понимаю, как эффективно использовать пакет RcppParallel?

Answer 1

Как упомянуто в комментариях, издержки могут быть проблематичными, особенно когда общее время выполнения короткое, лучше не инициализировать выходные векторы с нулем и использовать потоки локальных RNG. Пример реализации:

#include <Rcpp.h>
// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
// [[Rcpp::depends(RcppZiggurat)]]
#include <Ziggurat.h>


using namespace RcppParallel;

struct Norm : public Worker
{   
  // saved draws
  RVector<double> draws;

  // constructors
  Norm(Rcpp::NumericVector draws)
    : draws(draws) {}

  void operator()(std::size_t begin, std::size_t end) {
    Ziggurat::Ziggurat::Ziggurat zigg(end);
    for (std::size_t i = begin; i < end; i++) {
      draws[i] = zigg.norm();
    }
  }
};

// [[Rcpp::export]]
Rcpp::NumericVector parallelDraws(int x) {
  // allocate the output vector
  Rcpp::NumericVector draws(Rcpp::no_init(x));
  Norm norm(draws);
  parallelFor(0, x, norm);
  return draws;
}

// [[Rcpp::export]]
Rcpp::NumericVector serialDraws(int x) {
  // allocate the output vector
  Rcpp::NumericVector draws(Rcpp::no_init(x));
  Ziggurat::Ziggurat::Ziggurat zigg(42);
  for (int i = 0; i < x; i++) {
    draws[i] = zigg.norm();
  }
  return draws;
}

Обратите внимание, что я использую "параллельный RNG бедняка", то есть различные семена для разных потоков, и надеюсь на лучшее. Я использую end в качестве начального числа, поскольку begin может быть нулем, и я не уверен, что RNG в RcppZiggurat это нравится. Поскольку для создания Ziggurat объекта требуется некоторое время (и память), я также использую локальный объект для последовательных вычислений, чтобы быть справедливым.

Для 10 ^ 5 случайных дро все еще нет выгоды от использования параллельных вычислений:

> bench::mark(parallelDraws(1e5), serialDraws(1e5), check = FALSE, min_iterations = 10)[,1:5]
# A tibble: 2 x 5
  expression                min   median `itr/sec` mem_alloc
  <bch:expr>           <bch:tm> <bch:tm>     <dbl> <bch:byt>
1 parallelDraws(1e+05)   1.08ms   1.78ms      558.     784KB
2 serialDraws(1e+05)   624.16µs  758.6µs     1315.     784KB

Но за 10 ^ 8 ничьих я получаю хорошее ускорение на моем двухъядерном ноутбуке:

> bench::mark(parallelDraws(1e8), serialDraws(1e8), check = FALSE, min_iterations = 10)[,1:5]
# A tibble: 2 x 5
  expression                min   median `itr/sec` mem_alloc
  <bch:expr>           <bch:tm> <bch:tm>     <dbl> <bch:byt>
1 parallelDraws(1e+08)    326ms    343ms      2.91     763MB
2 serialDraws(1e+08)      757ms    770ms      1.30     763MB

Так что имеет ли смысл использовать параллельные вычисления, во многом зависит от того, сколько случайных дро вам нужно.

Кстати, в комментариях упоминается мой пакет dqrng . В этом пакете также используется метод Ziggurat для обычных (и экспоненциальных) отрисовок в сочетании с очень быстрыми 64-битными ГСЧ, что придает ему сопоставимую последовательную скорость с RcppZiggurat для обычных отрисовок. Кроме того, используемые RNG * предназначены для параллельных вычислений , т. Е. Нет необходимости надеяться получить непересекающиеся случайные потоки с использованием различных начальных чисел.