[Rcpp-devel] Parallel random numbers using Rcpp and OpenMP

Matteo Fasiolo matteo.fasiolo at gmail.com
Wed Apr 30 11:15:58 CEST 2014


Hi Mark,

many thanks for all the info, I will certainly have a detailed look at what
you are doing.

I think it would be nice to have a package that uses C++ level parallelism
(OpenMP or not)
to speed up random number generation at R level. For example:

> library("microbenchmark")
> library("mvtnorm")
> library("mvnfast")
> library("MASS")

> # Generating 1e4 50-dimensional multivariate normal r.v.
> N <- 10000
> d <- 50
> mu <- 1:d
> tmp <- matrix(rnorm(d^2), d, d)
> mcov <- tcrossprod(tmp, tmp)
> myChol <- chol(mcov)
>
> microbenchmark(rmvn(N, mu, mcov, ncores = 4),
+                rmvn(N, mu, mcov),
+                rmvnorm(N, mu, mcov, method = "chol"),
+                mvrnorm(N, mu, mcov))

Unit: milliseconds
                                  expr       min       lq   median
  uq      max neval
         rmvn(N, mu, mcov, ncores = 4)  8.339534 12.85869 13.23915
13.82640 38.04688   100
                     rmvn(N, mu, mcov) 31.645092 32.32220 33.70345
34.03561 57.78102   100
 rmvnorm(N, mu, mcov, method = "chol") 57.956183 58.99730 59.84668
82.11610 90.71133   100
                  mvrnorm(N, mu, mcov) 60.621391 62.06400 82.86829
85.42350 90.50662   100


 Here I followed Matt's suggestion and used C++11 mt19937_64 rng, with
seeds generated by R::runif.

Obviously what you have done by using RngStreams with Boost's
distributions is much more rigorous, so probably

that the approach rights to provide faster version of rnorm(), rpois() etc.


Thanks,


Matteo



On Wed, Apr 30, 2014 at 9:10 AM, Mark Clements <mark.clements at ki.se> wrote:

> My two cents worth:
>
> For the microsimulation package, we needed uniform random number streams
> and sub-streams at the C++ level, while supporting R's non-uniform
> random number distributions[*]. For this, we used the C++ RngStreams
> library and provided "double *user_unif_rand ()" for user-defined RNGs.
> Essentially, this reproduces the "L'Ecuyer-CMRG" RNG provided by R, with
> seed manipulation at the C++ level (cf. using .Random.seed and R's C
> machinery). We also allow for the C++ seeds to be set to and from R.
>
> See:
>
> https://github.com/mclements/microsimulation/blob/master/src/microsimulation.cc
>
> https://github.com/mclements/microsimulation/blob/master/src/microsimulation.h
> (circa lines 290-355)
>
> https://github.com/mclements/microsimulation/blob/master/R/rcpp_hello_world.R
> (circa lines 6-36)
>
> This is not implemented using OpenMP. For our purpose, where we are
> doing 10^7 small simulations, we chunk the simulations at the R level
> and use the parallel package to call the C++ code on each chunk. This
> approaches scales well with more processors.
>
> We have also looked at using the C++ RngStream library with
> Boost.Random's and C++11's non-uniform random number distributions.
> Again, this has not been implemented using OpenMP (todo?). For a simple
> wrapper for Boost, see:
>
> https://github.com/mclements/microsimulation/blob/master/src/rngstream-boost.hpp
>
> with an example:
>
> https://github.com/mclements/microsimulation/blob/master/src/rngstream-example.cpp
>
> Two brief notes: first, we have put the RngStream library in the ssim
> namespace for use with the microsimulation package. Other uses of the
> Boost wrapper would probably want to change the namespace.
>
> Second, the implementation for C++11 random number generators and
> distributions required a small change to the RngStream library -
> specifically, to output the random uniform as an unsigned long.
> Interestingly, the resulting C++11 random numbers drop every second
> value compared with R and Boost.Random.
>
> For the C+++11 wrapper, see the code RngStream* and rngstream* in:
> https://github.com/mclements/microsimulation/tree/master/test
>
> Sincerely, Mark.
>
> [*] The BH package was not available, C++11 compiler requirements were
> not accepted on CRAN,  reimplementing non-uniform random number
> distributions using UNURAN would have taken too long, and the R
> non-uniform random number distributions are extensive and well tested.
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20140430/0b994c7c/attachment.html>


More information about the Rcpp-devel mailing list