How to use different RNGs in different parts of a simulation avoiding correlation?-CodePudding

This question arises from my attempt to mix two different RNGs. I'd like to mix them choosing the best of the two according to the operations I need to carry out to achieve better performance. More concretely, the two RNGs are:

The Mersenne Twister (MT19937) coming from the random module of Python;
Another RNG that can be any of the Numpy RNGs, listed here.

I'm restricted to use the Mersenne Twister of the random module since in some situations it is better than any of the others in Numpy. Also, a requirement of my project is that I need to have reproducible results so unpredictable entropy should be avoided to set the seeds, but I can set one of the seed with a pseudorandom number generated with one of the RNG if needed.

So far I've only been able to implement the "safer" solution which is to use a MT19937 also in Numpy so that each time I need to use it, I pass the state from the Mersenne Twister of the random module (here for example you can obtain it with random.getstate()), do some operations with it and then pass the state back to the other. The problem with this solution is that passing the state creates a relevant overhead.

I'm unsure if other more performant solutions, like initializing two different RNGs at the start of the simulation and using them at will, can be problematic in terms of the quality/correlation of the sequence generated since I read here that using two differently seeded Mersenne Twister is not very good because the two sequences of pseudo-random numbers can be more correlated than one generated from a single one. However, in my situation I can use any of the Numpy RNG (a PCG-64 generator for example) in combination with the Python Mersenne Twister from the random module, so this is what I'd like to ask: is initializing a different RNG (with a different seed if useful) good enough? And also, what would be the best choice in Numpy to mix with a MT19937? Thank you in advance.

EDIT:

These are some timings to show that some operations can be made more performant with one RNG instead of the other:

import timeit

python_shuffle = timeit.timeit("random.shuffle(a)", "import random; a = list(range(500))", number=10000)
numpy_shuffle = timeit.timeit("rng.shuffle(a)", "from numpy.random import default_rng; a = list(range(500)); rng = default_rng()", number=10000)

python_gen1rand = timeit.timeit("random.random()", "import random", number=10000)
numpy_gen1rand = timeit.timeit("rng.random()", "from numpy.random import default_rng; rng = default_rng()", number=10000)

print(python_shuffle/numpy_shuffle)
print(python_gen1rand/numpy_gen1rand)

which prints something like:

14.992528343855014
0.15127849559797613

so the default Numpy PCG64 RNG, which should be one of the fastest, is 15 times better than the other one in shuffling but 7 times worse than the one in the random module to generate a single random number.

CodePudding user response：

I probably found an answer myself.

The way to do this is to use the same seeded MT19937 generator in Numpy (actually I need to pass the initial state of the MT19937 of the random module to the one in numpy since seeds are computed differently) and then jump ahead the state! This is possible with the jumped function which accept a jumps parameter in the Numpy library where each jump makes the sequence shift as if 2^128 random numbers had been generated. Doing this, it is possible to use the two RNGs at will because it is like we are creating random numbers from a single RNG. The method used to jump ahead is described in this paper, as a support for what I claim it is said

Random number generators (RNGs) with multiple disjoint streams and substreams are an important component of any good general-purpose simulation or statistical software. they are very handy, for example, to obtain parallel RNGs and to support the implementation of variance reduction techniques (Hellekalek 1998, Kelton 2006, Law and Kelton 2000, L’Ecuyer et al. 2002). The most convenient way of getting these streams and substreams is to start with a backbone RNG having a huge period and partition its output sequence into long disjoint subse quences and subsubsequences whose starting points are at equidistant lags (Law and Kelton 2000, L’Ecuyer 1990, L’Ecuyer and Côté 1991). When a new stream is needed, we find its starting point by jumping ahead from the starting point of the current subsequence to the starting point of the next one.

As a verification of the cost of the jumping, this shows that jumping ahead by 1 is cheap

numpy_jump = timeit.timeit("rng.bit_generator.jumped(1)", "from numpy.random import default_rng, MT19937; rng = default_rng(MT19937())", number=1)
print(numpy_jump) # 0.0042395820000820095 on my computer