Pass random state to setRandom with RcppEigen

Question

Is there a way to pass the random state to Eigen's setRandom with RcppEigen or do I need to use runif?

Here is an example:

// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
using namespace Rcpp;
using  Eigen::MatrixXd;
using  Eigen::VectorXd;

// [[Rcpp::export]]
NumericVector fx() {
  RNGScope   scope;
  MatrixXd   x(3,2);
  x=x.setRandom();
  x.col(1)=as<VectorXd>(runif(3,0,1));

return wrap(x);
}

Testing it:

set.seed(42); fx()
#           [,1]      [,2]
#[1,] -0.8105760 0.9148060
#[2,]  0.6498853 0.9370754
#[3,]  0.6221027 0.2861395

set.seed(42); fx()
#           [,1]      [,2]
#[1,] -0.9449154 0.9148060
#[2,]  0.8063267 0.9370754
#[3,] -0.0673205 0.2861395

Note how column 2 (i.e., runif) is reproducible, but column 1 (i.e., setRandom) is not.

Answer 1

The RNG in Eigen is orthogonal to R's RNGs. We have RNGScope to deal with R and allow you to set.seed() as you know; what you do with Eigen's RNG is separate.

In particular, you would have to add the glue to register Eigen's RNG as a user-provided RNG for R. By default, R is unaware of Eigen, which makes sense as R users expect R's RNGs.

Your issue here seems to be a vanilla 'programming with Eigen' issue where you fail to initialize the Eigen RNG. It has nothing to do with Rcpp.

Edit: Here is your example, corrected. Turns out Eigen uses the system RNG, so you need srand() to seed it.

R> sourceCpp("/tmp/roland.cpp")
R> set.seed(42); fx(42)
          [,1]     [,2]
[1,] -0.933060 0.914806
[2,] -0.340072 0.937075
[3,]  0.381271 0.286140
R> set.seed(42); fx(42)
          [,1]     [,2]
[1,] -0.933060 0.914806
[2,] -0.340072 0.937075
[3,]  0.381271 0.286140
R>

It uses this modified version of your code:

// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
using namespace Rcpp;
using  Eigen::MatrixXd;
using  Eigen::VectorXd;

// [[Rcpp::export]]

  RNGScope   scope;
  MatrixXd   x(3,2);
  srand(seed);
  x=x.setRandom();
  x.col(1)=as<VectorXd>(runif(3,0,1));

  return wrap(x);
}

Edit 2 in response to OP's comments:

I would not expect much of a speed difference between Rcpp::runif() and the srand() calls made by Eigen, so you are still stuck with the fact that srand() has been having issues, and is possibly behaving differently between systems.

Quick demo script:

#include <RcppEigen.h> 

// [[Rcpp::depends(RcppEigen)]]

// [[Rcpp::export]]
Rcpp::NumericVector v1(int n) {
  return Rcpp::runif(n);
}

// [[Rcpp::export]]
Rcpp::NumericVector v2(int n) {
  Eigen::VectorXd   x(n);
  x = x.setRandom();
  return Rcpp::wrap(x);
}

/*** R
library(rbenchmark)
N <- 1e7
benchmark(v1(N), v2(N))
*/

which produces

R> sourceCpp("/tmp/roland.cpp")

R> library(rbenchmark)

R> N <- 1e7

R> benchmark(v1(N), v2(N))
   test replications elapsed relative user.self sys.self user.child sys.child
1 v1(N)          100  12.633    1.000    11.356    1.261          0         0
2 v2(N)          100  17.222    1.363    13.981    3.198          0         0
R> 

and note that RcppEigen is slower here even in a simpler setup of just creating vectors. But we are talking microseconds here and that is probably not something I'd worry about either way in a real application which will most likely have other bottlenecks.