"parallel" random number generators

[Bob Carpenter]

Suppose we're in the MCMC case and want to generate
four parallel Markov chains.

The theoretical concern is that if you just take four
"random" seeds, the sequence of numbers they generate
might overlap. If you're running a whole lot of parallel
MCMC chains and the period of the random number generator (the
number of random numbers before it repeats) is small, then
there's a non-negligible chance of overlap.

The "parallel" random number generators aren't parallel
in the usual sense of thread-safe, etc. Instead, what
they do (or at least L'Ecuyer's linear congruential generator
does) is let you efficiently skip ahead any number of
elements in the random-number sequence. This allows you to
get "independent" samples for multiple Markov chains. With
four chains and a total chain length of M, we can use:

chain 1: seed N
chain 2: seed N, skip ahead M/4 items
chain 3: seed N, skip ahead 2M/4 items
chain 4: seed N, skip ahead 3M/4 items

Then, each chain draws from a different part of the underlying
random number sequence starting from seed N, and are
effectively "independent".

I don't think that's much of a concern for us.
And the solutions to overcome it are relatively slow
(around 8% of the speed of faster generators). So, what
I'm proposing is:

PROPOSED DESIGN
---------------
Mersenne Twister RNG (Boost mt19937)

Specify for each chain the base seed and chain ID
and use seed base seed + ID for the chain.

The alternative is to go with:

FANCY PANTS RNG DESIGN
----------------------
L'Ecuyer RNG (Boost ecuyer1988)

Specify for each chain the base seed, the chain ID and
the total number of chains. Asisgn seeds to chains
as above.

- Bob

[Andrew Gelman]

I don't understand why these methods aren't all equally fast. Aren't they all equivalent to just picking 4 seeds? I mean, I know there's some subtlety here I'm not guessing, but I don't have a sense of what it is that I'm not knowing here.
A

[Bob Carpenter]

There are a variety of configurable random number
generators with different properties. Sort of like
MCMC algorithms.

The two we're looking at are:

Mersenne Twister:
reasonably fast samples
very long cycle before repeats (2**19000)
slow skip-ahead

L'Ecuyer:
slower samples (8% speed of Mersenne)
shorter cycle before repeats (2**61)
fast skip-ahead

To absolutely guarantee non-overlap of
the sequences, you need something like L'Ecuyer.

Matt, Daniel, and I discussed it today and thought
it wasn't necessary. But then Ben and Malecki have
been urging us to use a parallel version all along,
so we thought we might still be missing something,
hence the e-mail.

It'll have the biggest impact on performance in
simple examples like RATS where the random number
generation is a bigger proportion of the compute
time.

I could also make it configurable, but that makes
the software more complicated to write, document,
and maintain.

- Bob

[Ben Goodrich]

On Feb 3, 10:59 pm, Bob Carpenter <c...@alias-i.com> wrote:
> Matt, Daniel, and I discussed it today and thought
> it wasn't necessary.  But then Ben and Malecki have
> been urging us to use a parallel version all along,
> so we thought we might still be missing something,
> hence the e-mail.

> >> PROPOSED DESIGN
> >> ---------------
> >> Mersenne Twister RNG (Boost mt19937)
>
> >> Specify for each chain the base seed and chain ID
> >> and use seed base seed + ID for the chain.

My impression is that this is insufficient, although this is not my
area of expertise and I couldn't quickly demonstrate why.

I do know that the people who invented Mersenne Twister have proposed
a different way of doing it in parallel

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/DC/dc.html

that is similar in that the ID is a parameter but more complicated
than just ++ing the base seed.

Ben

[Matt Hoffman]

This looks interesting, might be a good compromise if it's not too huge of a pain to implement.

[Ben Goodrich]

I think it is the case that you can "ask" their DC library for the
right seeds to use given your job requirements and then just pass the
result to boost's PRNGs.