How can I seed a parallel linear congruential pseudo-random number generator for maximal period?

Question

Normally when I seed a sequential random number generator in C, I use the call

srand(time(NULL))

then use

rand() mod N

to obtain a random integer between 0 and N-1. However, when I do this in parallel, the calls to time(NULL) are so close to each other that they end up being exactly the same number.

I have tried using a linear congruential random number generator:

$x_{n+1} := a x_n + c\;\;\; (\bmod m) $ for some integers $a,c,$ and $m$.

I know that choosing $m=2^k$ for some large integer $k$ produces fast results because the modulus operator can be calculated by truncating digits. However, I find it difficult to establish seeds that produce random sequences with a large period in parallel. I know that a period length is maximal if

1. c and m are primes with respect to eachother
2. a-1 is divisible by all prime factors of m
3. if m is a multiple of 4, a-1 must also be a multiple of 4.

(source: wikipedia)

But how do I ensure that all random number streams have this maximal property? In terms of MPI, how do I incorporate the rank and size to produce maximal periods using the linear congruential method? Would it be easier to use a Lagged Fibonacci or Mersenne Twister to produce longer parallel random streams?

Answer 1

The trick is to interleave each process's LCG stream: for $p$ processes, we modify the LCG

$$ x_{n+1} := a x_n + c\;\;\; (\bmod m),$$

to be

$$ x_{n+p} := A_p x_n + C_p\;\;\; (\bmod m),$$

where $A_p$ and $C_p$ effectively step forward $p$ steps. We can quickly derive them by expanding the original LCG step:

$$ x_{n+2} = a (a x_n + c) + c\;\;\; (\bmod m) = a^2 x_n + (a+1)c\;\;\; (\bmod m), $$

and the pattern is that $A_p \equiv a^p \mod m,$ and $C_p$ is the result of, starting with the number $0$, successively multiplying by $a$ and adding $1$ $p$ times, then multiplying by $c$, all $\bmod m$.

The last step is to ensure that each process's $p$-stride LCG does not overlap: simply initialize the process with rank $r$ with $x_r$ and a parallel LCG with individual periods of $N/p$ is ready, where $N$ is the original period and $p$ is the number of processes. If each process's modified LCG is used equally, the full period is recovered in parallel.

I implemented this about six months ago (perhaps naively), and you can find the code here.

Answer 2

There's a very nice tutorial overview paper by Katzgrabber, Random Numbers In Scientific Computing: An Introduction, which I point people to who want to be a user of PRNGs for scientific computation. Linear congruential generators are fast, but that's about all they have going for them; they have short periods, and they can very easily go wrong; perfectly reasonable looking combinations of a, c, and m can end up with horrifically correlated outputs, even if you satisfy the usual requirements between a, c, and m.

Worse, in one common case where m is a power of two (so the mod operation is fast), the lower order bits have a much shorter period than the sequence as a whole, so if you're doing rand() % N, you have an even shorter period than you would expect.

As a general rule of thumb, lagged-fibbonacci, MT, and WELL generators have much better properties and they're still pretty fast.

In terms of seeding in parallel, Jack Poulson's method is nice because it gives a well-defined sequence of numbers divided evenly between the processors. If that doesn't matter, you can do anything reasonable to seed the different PRNGs; the same paper suggests something lots of people have come up with independantly, hashing the PID or MPI task number with the time. The particular formula suggested there is

long seedgen(void)  {
    long s, seed, pid;

    pid = getpid();
    s = time ( &seconds );

    seed = abs(((s*181)*((pid-83)*359))%104729); 
    return seed;
}

I have no particular opinions about that specific implementation, but the general approach is certainly reasonable.

Answer 3

A simple idea for spreading a typical sequential RNG over a decent number of threads is to have a single thread advance the seed as fast as possible and send only every thousandth or so seed out to memory. Then have each of your other threads pick up one of these spaced out reference seeds and process the 1000 values in that block, i.e. regenerate again the 1000 seeds in the block, generate their psuedorandom draws, and then do whatever other processing your task involves.

This works because for RNGs that don't compute all that much (LCG is certainly one, but many others should be in this category) the real bottleneck is sending the seeds out to memory (and also probably subsequent processing). If you run a LCG without sending anything out to memory, the whole thing ought to stay in CPU registers and be extremely fast. Even for a more complicated RNG, you should stay in L1 cache and be very fast.

I've used this very simple approach with an LCG that for legacy reasons we have to keep. We basically get linear speedup up to the 4-8 threads in a typical multicore workstation. But now I will try the method from Jack Poulson's answer and hopefully get even faster :).

OTOH, I believe this simple trick should work for other inherently sequential RNGs.

Answer 4

This answer is late in coming, but you should have a look at SPRNG. It's specifically designed for scalability in parallel and supports a handful of types of PRNGs.