Random matrix generation

[option72]

I would like to generate a random matrix of uniform distributed numbers in a "sequential" way.
The parallel version, Ralf wrote, works very well because it is fast and it has good statistical properties.
In particular the cor() function gives almost negligible values even for small matrices (eg 10 x 10).
I cannot say the same thing for the "sequential" version, poor cor() results. I think I should set
the RNG (eg auto rng = dqrng::generatordqrng::xoshiro256plusplus();) but I am unable.
The following code, I modified, doesn't give the expected result: too high correlation. Thanks for you help!

#include <Rcpp.h>
// [[Rcpp::depends(dqrng, BH)]]
#include <dqrng_distribution.h>
using namespace Rcpp;
// [[Rcpp::export]]
Rcpp::NumericMatrix random_mtx(int n, int m) {
  dqrng::uniform_distribution dist(0.0, 1.0);               // Uniform distribution [0,1)
  auto rng = dqrng::generator<dqrng::xoshiro256plusplus>();
   // seeded from R's RNG
  Rcpp::NumericMatrix res(n,m);
  for (int i = 0; i < m; ++i) {
    
    for (int j = 0; j < n; ++j) {
        
        res(j,i) = dist(*rng);
      
    }
  }
  return res;
}

/*** R

 vrng <- random_mtx(100, 10)
print(symnum(x = cor(vrng), cutpoints = c(0.001, 0.003, 0.999),
             symbols = c(" ", "?", "!", "1"),
             abbr.colnames = FALSE, corr = TRUE))

*/

and and even the following has the same problem:

#include <Rcpp.h>
// [[Rcpp::depends(dqrng, BH)]]
#include <dqrng_distribution.h>
#include <xoshiro.h>

template<typename RNG>
Rcpp::NumericVector runif_rng(int n, int seme, int procx) {
  auto rng = dqrng::generator<RNG>(seme, procx);
  dqrng::uniform_distribution dist;
  Rcpp::NumericVector result(Rcpp::no_init(n));
  std::generate(result.begin(), result.end(), [rng, dist] () {return dist(*rng);});
  
  return result;
}

// [[Rcpp::export]]
Rcpp::NumericMatrix runif_256plusplus(int n, int m, int seme) {
  
  Rcpp::NumericMatrix v0(n,m);
  for (int i = 0; i < m; ++i) {
  Rcpp::NumericVector v1 = runif_rng<dqrng::xoshiro256plusplus>(n, seme * i +1,i);
  for (int j = 0; j < n; ++j){
    
     v0(j,i)=v1(j);
    
  }
  } 
  
  return v0; 
 
  
}



/*** R
library(bench)

N <- 100
bm <- bench::mark(
  runif_256plusplus(N, 10, 2025),
  check = FALSE,
  min_time = 1
)
print(bm[, 1:6])

vec1 <- runif_256plusplus(N, 10, 2025)
print(symnum(x = cor(vec1), cutpoints = c(0.001, 0.003, 0.999),
             symbols = c(" ", "?", "!", "1"),
             abbr.colnames = FALSE, corr = TRUE))

*/

[rstub]

The cut points c(0.001, 0.003, 0.999) for c(" ", "?", "!", "1") are very aggressive, i.e. anything above 0.003 is marked as !. I am not sure what the distribution function of correlation coefficients is, but it seems for vectors of 100 random numbers, it is unlikely to go below 0.003. Things look much better if you use 1e6 instead of 100, like it is done in the parallel vignette:

> vrng <- random_mtx(1e6, 10)

> print(symnum(x = cor(vrng), cutpoints = c(0.001, 0.003, 0.999),
+              symbols = c(" ", "?", "!", "1"),
+              abbr.colnames = FALSE .... [TRUNCATED] 
                         
 [1,] 1                  
 [2,]   1                
 [3,] ?   1              
 [4,] ? ?   1            
 [5,] ?       1          
 [6,]   ?     ? 1        
 [7,]       ? ?   1      
 [8,] ?           ? 1    
 [9,]     ?   ?   ?   1  
[10,]   ? ? ?   ? ?     1
attr(,"legend")
[1] 0 ‘ ’ 0.001 ‘?’ 0.003 ‘!’ 0.999 ‘1’ 1

BTW, since it is possible to iterate over a Rcpp::NumericMatrix, you could also simplify your code by using the generate function supported by the RNG:

#include <Rcpp.h>
// [[Rcpp::depends(dqrng, BH)]]
#include <dqrng_distribution.h>

// [[Rcpp::export]]
Rcpp::NumericMatrix random_mtx(int n, int m) {
  auto rng = dqrng::generator<dqrng::xoshiro256plusplus>();
  // seeded from R's RNG
  Rcpp::NumericMatrix res(n,m);
  rng->generate<dqrng::uniform_distribution>(res, 0.0, 1.0);
  return res;
}

[option72]

Thank you very much for your precious help!

Forgive me for some further questions. I would like to build a stochastic model for simulating the evolution of a population of insured people. The simulation process implies a huge number of random drawings at different steps (eg a time orizon of 10 years and 100 replications need 1000 random drawings), so I am wondering whether there are Rccp instructions to set the seed of the random number generator at each iteration.

Morover I hope it is possible to get control over rng and to manage all the process by Rcpp, speed is crucial.

[rstub]

You are welcome!

While it is possible to set the seed at each iteration, it is not considered best practice. Where would that seed come from? How would you make sure that there is no overlap between the different series of random numbers you are getting?

For serial work, it is best to set the RNG's state once in the beginning of your code and then forget about that. This is also how RNGs are being tested. In the case of dqrng, you can use dqset.seed() from the R side to set the state of the global RNG. Within C++, you can then either use the functions from the compiled library or access the global RNG directly. The latter is only necessary if you need distribution functions beyond the currently supported uniform, normal and exponential. In both cases you are using the global RNG and don't have to think about its state any further.

But maybe I am misunderstanding your use-case. In that case it would be good to explain why you think you need to set the seed regularly.

[rstub]

I forgot one other use-case for accessing the global RNG: Many calls to obtain a small number of random numbers within the same function. Using the functions from the compiled library is possible but inefficient in that case, c.f. #40.

[option72]

Sorry for taking some days before answering. I think you solved my problem, I had no dubt about!
I needed a way to build random matrices from a uniform distribution after having set a seed. The code below (unfortunately not well written) works as desidered; nRip is the number of iterations of the process, where random_mtx is called 5 times, and each times 5 different random matrices are returned.

#include <Rcpp.h>
#include <random>
//#include <algorithm>
// [[Rcpp::depends(dqrng, BH)]]
#include <dqrng_distribution.h>

using namespace Rcpp;

inline Rcpp::NumericMatrix random_mtx(int n, int m) {
  auto rng = dqrng::generator<dqrng::xoshiro256plusplus>();
  // seeded from R's RNG
  Rcpp::NumericMatrix res(n,m);
  rng->generate<dqrng::uniform_distribution>(res, 0.0, 1.0);
  return res;
}

inline Rcpp::IntegerMatrix permutazioni(int n,Rcpp::IntegerVector x)  {
  Rcpp::IntegerMatrix perm(n,x.size());
  for(int i = 0; i < n; ++i) {
    std::random_device rd;
    std::mt19937 g(rd());
    std::shuffle(x.begin(), x.end(), g);
    for(int j = 0; j <x.size();j++){
      perm(i,j)=x[j];
    }
  }
  return perm;
}

// [[Rcpp::export]]
Rcpp::NumericMatrix scenari_seq(int nR, int nC, IntegerVector mEventiAtt, IntegerVector mEventiCes, int nAnni, int nRip) {
  
  Rcpp::NumericMatrix mSim(nR, nC + 3 * nAnni);
  double mpcM=0.0;
  double mpcU =0.0;
  int nEA = mEventiAtt.size();
  int nEC = mEventiCes.size();
  
  for(int r = 1; r <= nRip; r++) {
  
  Rcpp::NumericMatrix mpcMorte(nR, nAnni);
  Rcpp::NumericMatrix mpcPensione(nR, nAnni);
  Rcpp::NumericMatrix mpcRiscatto(nR, nAnni);
  Rcpp::NumericMatrix mpcUscita(nR, nAnni);
  
  Rcpp::IntegerMatrix mEventiAtt1(nR, nEA);
  Rcpp::IntegerMatrix mEventiCes1(nR, nEC);
  
  mpcMorte =    random_mtx(nR, nAnni);
  mpcPensione = random_mtx(nR, nAnni);
  mpcRiscatto = random_mtx(nR, nAnni);
  mpcUscita =   random_mtx(nR, nAnni);
  
  mEventiAtt1 = permutazioni(nR, mEventiAtt);
  mEventiCes1 = permutazioni(nR, mEventiCes);
  
 
   for(int t = 1; t <= nAnni; t++) {
     
       for(int i = 0; i < nR; i++) {
       
       } //next i
       mpcM = mpcMorte(0,t-1);
       mpcU = mpcUscita(0,t-1);
       Rcout << "iterazione: "<<  r << ", epoca: "<< t <<": "  << mpcU << "; " << mpcM << std::endl;     
   }  //next t
   
  }  //next r
   
  
  return mSim;
}


/*** R
mEventiAtt <- as.matrix(c(1,2,3,4))
mEventiCes <- as.matrix(c(1,2,3,4))
set.seed(20250420)
scenari_seq(100, 8, mEventiAtt, mEventiCes, 10,100)

*/

[rstub]

I see two problems in your approach:

1. Due to using std::random_device, the results are no reproducible. It would be cleaner to use dqrng here as well.
2. I don't like the pattern of using dqrng::generator<dqrng::xoshiro256plusplus>() multiple times in a loop. This way one is seeding one RNG with the consecutive output from another RNG. It is not likely that you will get any overlaps, but still it would be cleaner to seed the RNG from dqrng only once.

The second issue can be addressed using the global RNG that dqrng provides, e.g. by accessing it via dqrng::random_64bit_accessor:

#include <Rcpp.h>
// [[Rcpp::depends(dqrng, BH)]]
#include <dqrng_distribution.h>

inline Rcpp::NumericMatrix random_mtx(int n, int m) {
  auto rng = dqrng::random_64bit_accessor{};
  Rcpp::NumericMatrix res(n,m);
  rng.generate<dqrng::uniform_distribution>(res, 0.0, 1.0);
  return res;
}

Alternatively, one could also use the functions from the compiled library, which automatically make use if the global RNG:

#include <Rcpp.h>
// [[Rcpp::depends(dqrng)]]
#include <dqrng.h>

inline Rcpp::NumericMatrix random_mtx(int n, int m) {
  Rcpp::NumericVector res = dqrng::dqrunif(n * m, 0.0, 1.0);
  res.attr("dim") = Rcpp::Dimension(n, m);
  return Rcpp::as<Rcpp::NumericMatrix>(res);
}

Using the dqrng.h header also gives you access to the sampling functions, which can be used to implement permutazioni:

#include <Rcpp.h>
// [[Rcpp::depends(dqrng)]]
#include <dqrng.h>

inline Rcpp::IntegerMatrix permutazioni(int n,Rcpp::IntegerVector x)  {
  Rcpp::IntegerMatrix perm(n,x.size());
  for(int i = 0; i < n; ++i) {
    Rcpp::IntegerVector indices = dqrng::dqsample_int(x.size(), x.size());
    for(int j = 0; j <x.size();j++){
      perm(i,j)=x[indices[j]];
    }
  }
  return perm;
}

(This assumes that x.size() < 2^31.)

Now that we know how to use the global RNG instead of a local instance: How do we control the RNG kind? How to set the seed properly?

There are two ways to control the type:

1. Use dqrng::dqRNGkind("Xoshiro256++") early in your controlling R code.
2. Use dqrng::dqRNGkind("Xoshiro256++") early in your R-facing C++ code.

There are two ways to set the seed based on R's RNG:

1. Use dqrng::dqset.seed(NULL) early in your controlling R code.
2. Use dqrng::dqset_seed(R_NilValue) early in your R-facing C++ code.

Here the question is to what extend you want to use code from dqrng within your R code, or if the C++ functions should be self-contained. With self contained C++ code, i.e. using the second options above, you will still be seeding the RNG every time you are calling the function from R. This is still not ideal, but much better than doing so inside a loop in C++.

Controlling the RNG from R via dqset.seed solves that properly, but has the disadvantage that a later set.seed has no effect on dqrng's RNG. Two possibilities:

1. Every set.seed is followed by dqset.seed(NULL).
2. Use dqrng::register_methods() to use dqrng's RNG as replacement for R's build in RNG. This way, set.seed affects the global RNG and it is enough to that once.

dqrng::dqRNGkind("Xoshiro256++")
dqrng::register_methods()
set.seed(20250428)

[arbitrary R code]

[option72]

Great, it is even faster! I implemented the new version of permutazioni and for random_mtx I chose to access the global RNG from the compiled library. Finally I added the last three lines of code to my R program.
I would like to thank you again, without your precious advice and patience I would never have been able to write an efficient and correct code.