Parallel pseudorandom number generators

Question

This question is primarily related to a practical software-engineering problem, but I would be curious to hear if theoreticians could provide more insight in it.

---

Put simply, I have a Monte Carlo simulation that uses a pseudorandom number generator, and I would like to parallelise it so that there are 1000 computers running the same simulation in parallel. Therefore I need 1000 independent streams of pseudorandom numbers.

Can we have 1000 parallel streams with the following properties? Here $X$ should be a very well-known and widely-studied PRNG with all kinds of nice theoretical and empirical properties.

1. The streams are provably as good as what I would get if I simply used $X$ and split the stream generated by $X$ into 1000 streams.

2. Generating the next number in any stream is (almost) as fast as generating the next number with $X$.

Put otherwise: can we get multiple independent streams "for free"?

Of course if we simply used $X$, always discarding 999 numbers and picking 1, then we certainly would have property 1, but we would lose in the running time by factor 1000.

A simple idea would be to use 1000 copies of $X$, with seeds 1, 2, ..., 1000. This certainly would be fast, but it is not obvious if the streams have good statistical properties.

---

After some Googling, I have found, for example, the following:

• The SPRNG library seems to be designed for exactly this purpose, and it supports multiple PRNGs.

• Mersenne twister seems to be a popular PRNG nowadays, and I found some references to a variant that is able to produce multiple streams in parallel.

But all this is so far from my own research areas, that I couldn't figure out what is really the state-of-the-art, and which constructions work well not only in theory but also in practice.

---

Some clarifications: I do not need any kind of cryptographic properties; this is for scientific computation. I will need billions of random numbers, so we can forget any generator with a period of $< 2^{32}$.

Edit: I cannot use a true RNG; I need a deterministic PRNG. Firstly, it helps a lot with debugging and makes everything repeatable. Secondly, it allows me to do, e.g., median-finding very efficiently by exploiting the fact that I can use the multi-pass model (see this question).

Edit 2: There is a closely related question @ StackOverflow: Pseudo-random number generator for cluster environment.

Answer 1

You can use an evolution of the Mersenne Twister algorithm developed by Saito and Matsumoto:

SIMD-oriented Fast Mersenne Twister (SFMT)

SFMT is a Linear Feedbacked Shift Register (LFSR) generator that generates a 128-bit pseudorandom integer at one step. SFMT is designed with recent parallelism of modern CPUs, such as multi-stage pipelining and SIMD (e.g. 128-bit integer) instructions. It supports 32-bit and 64-bit integers, as well as double precision floating point as output. SFMT is much faster than MT, in most platforms. Not only the speed, but also the dimensions of equidistributions at v-bit precision are improved. In addition, recovery from 0-excess initial state is much faster. See Master's Thesis of Mutsuo Saito for detail.

The period varies from $2^{607}-1$ to $2^{216091}-1$.

Using one same pesudorandom number generator for generating multiple independent streams by changing the initial values may cause a problem (with negligibly small probability). To avoid the problem, using different parameters for each generation is preferred. This technique is called dynamic creation of the MT parameters.

In the SFMT source code you can find some examples of parameter sets (of variable periods) and an awk script to convert a CSV file to a compilable parameter set. There is also a tool called "Dynamic Creation of Mersenne Twister generators".

The authors recently developed another modified version of the Mersenne Twister - Mersenne Twister for Graphic Processors - designed to run in GPUs and take advantage of their native parallel execution threads. The key feature is speed: $5 \times 10^7$ random integers every 4.6ms on a GeForce GTX 260.

The periods of generated sequence are $2^{11213}-1$ , $2^{23209}-1$ and $2^{44497}-1$ for 32-bit version, and $2^{23209}-1$, $2^{44497}-1$, $2^{110503}-1$ for 64-bit version. It It support 128 parameter sets for each period, in other words, it can generate 128 independent pseudorandom number sequences for each period. We have developed Dynamic Creator for MTGP, which generates more parameter sets

Indeed they provide a MTGPDC tool to create up to $2^{32}$ parameter sets (i.e. independent streams).

The algorithm passes the main randomness tests like Diehard and NIST.

A preliminary paper is also availbale on arXiv: A Variant of Mersenne Twister Suitable for Graphic Processors

Answer 2

There seem to be many ways to tackle this problem, but one simple way would be to use the Blum Blum Shub PRNG. This PRNG is defined by the recurrance relation $x_{i+1} = x_i^2 \mbox{ mod }N$, where $N$ is a semiprime. To get a random bit out of this you can simply take the bit parity of $x_i$. What's nice about this is that since $x_{i+k} = x_i^{2^k}\mbox{ mod }N = x_i^{2^k \mbox{ mod } \lambda(N)}\mbox{mod }N$ you can directly calculate any step in time constant in $k$ (i.e. $O(\log(N)^3)$ or faster depending on which multiplication algorithm you use for the modular exponential). Thus is you have $M$ machines, then for the machine indexed by $y$ you can use the generator $x_{i+1,y} = x_i^{2^M \mbox{mod }\lambda(N)}\mbox{ mod }N$, where $x_{0,y} = x_0^{2^y \mbox{ mod }\lambda(N)}\mbox{ mod }N$, where $x_0$ is your seed. Conveniently this generates exactly the same stream of numbers as if you used a single stream and distributed it's output to each of the machines in turn.

This isn't the fastest of PRNGs, though, so it will only be useful if the overhead of whatever you are doing in the simulation is significantly more than the cost of the PRNG. However it is worth pointing out that it will be much faster for certain combinations of $M$ and $N$ than others, particularly if the binary representation of $2^M \mbox{ mod }\lambda(N)$ contains few 1s or is small.

Answer 3

How about a preprocessing phase? Given a random seed $s$ (of size $n$), run $X$ to get a pseudorandom stream of size $1000n$. Denote this stream by $s_1,s_2,\ldots,s_{1000}$, where for $1\le i \le 1000$, $s_i$ is a contiguous portion of the stream of size $n$.

This preprocessing phase can be done with a very low overhead, given the fact that $X$ is an efficient PRNG (today, we have very fast PRNG's).

Now, give $s_i$ as the seed to the $i$th machine, which uses $X$ to generate its own pseudorandom stream.

Given the nice properties of $X$, unless $s$ is known, for any $1 \le i < j \le 1000$, the seeds $s_i$ and $s_j$ are computationally independent. Moreover, you only have to generate and save one small seed (i.e. $s$); therefore, this approach does not need a great deal of true randomness or storage.

Answer 4

You could use a pseudorandom function $f$ such as AES or ChaCha with a single random key, encrypting a counter. Assign each of the $M = 1000$ parallel processes a unique starting value in $\{ 0, 1, \ldots, M - 1 \}$, and then compute the $j$th random block of bits for process $i$ as $f(i + jM)$, i.e. increment the counter in each process by $M$ for every subsequent block of random bits.

This will give you a cryptographic RNG on every process, but it does not necessarily come with a performance cost. AES is fast if you have hardware that supports it, and ChaCha is fast regardless. Of course, you'll want to measure this in your specific setting to be sure.

Both desired properties 1 and 2 are directly satisfied by this. It's moreover convenient that the behavior of the entire system of parallel tasks is controlled by a single "seed" (the key for $f$).

Answer 5

There is now a jump function for SFMT (a fast Mersenne Twister implementation).

This allows me to initialise 1000 MTs so that there is no cycle overlap. And SFMT should be faster than MTGP. Almost perfect for my purposes.

Answer 6

You can just use 1000 instances of the Mersenne Twister initialized with different seeds.

You can sample the seeds from another Mersenne Twister, or, to be surer of their independence, from the OS cryptographic pseudorandom number generator (/dev/urandom in Linux).

The Mersenne Twister always operates on the same cyclic sequence, the seed controls where you start generating it. With indepenently sampled seeds, each generator will start at different, typically very far points, with a very small probability of intersection.