Random() Adventures

Better to remain silent and be thought a fool than to speak and to remove all doubt.

Often, but probably incorrectly, attributed to Mark Twain

The idea of writing about random number generation (I mean, pseudo-random number generation, but that term involves way too much typing) makes me extremely nervous.  When it comes to that topic, I think of two groups of people:

1. The High Priesthood. The individuals (one or two dozen, in my imagination) who develop, test and understand these generators
2. The thousands who make fools of themselves writing about random numbers on the Internet.

I am definitely not in the first category; I’ll keep my fingers crossed that this post doesn’t put me solidly into the second.

Back in the day, when dinosaurs roamed the earth and computers were Big Machines in Large Rooms with Raised Floors, there was a popular expression in the industry – No one ever got fired for buying IBM. The Mersenne Twister is the IBM of pseudo-random number generators.  It’s not the fastest or most efficient generator.  It may or may not deliver the “best” results (a debate I’ll leave to the aforementioned priesthood). But it is fast enough, good enough and best-known of the modern vintage of random number generators – so it has become a de facto standard.  If you were to recommend a different generator that I haven’t heard of, I’d probably ask how it compares to the Mersenne Twister.  Python’s random module uses MT19937, the most commonly used (and heavily tested) version of Mersenne Twister.  The standard library designers definitely could have chosen worse.  But is it good enough for our discrete event simulator?

Law (Ch. 7) lists five properties of a “good” random number generator for simulation purposes:

1. The output should be uniformly distributed on the interval [0, 1] and uncorrelated.
2. It should be efficient “enough”, in terms of both speed and memory requirements.
3. The output should be reproducible – i.e., given an initial state (such as seed), it should always generate the same output.
4. To facilitate simulation output analysis, it should be able to generate multiple independent (uncorrelated) streams of sufficient length.
5. It should be portable, ensuring that given an initial state, the output is the same on every target computing platform.

The generator provided by the random module fulfills properties 1, 2, 3 and 5.  The fourth is, however, a bit more problematic.

The canonical way to generate multiple independent streams is to divide a single stream into multiple non-overlapping substreams, each of which should meet the same statistical randomness tests as the main stream and be independent of each other.  Many random number generators include functions that can be used to initialize substreams from a single stream, either through a jump-ahead or other means.  The Python 3.x random module does not.

First, a bit of history is in order…

• The random module initially (or at least as of Python 2.0) included a jumpahead() method, which was supported by what was then the default random number generator – Wichmann-Hill. The Mersenne Twister replaced Wichmann-Hill in version 2.3, released in 2003.
• The original Mersenne Twister publications (1998) did not include a true jump ahead algorithm. As a result, the behavior of Python jumpahead() was changed in version 2.3 as well. Quoting the documentation: Instead of jumping to a specific state, n steps ahead, jumpahead(n) jumps to another state likely to be separated by many steps.
• The jumpahead() function was eventually removed from Python as an API defect, as the implementation no longer really matched the function’s implied semantics.
• In 2013, Mersenne Twister jumpahead algorithms were published (code here). Unfortunately, this functionality has not yet been incorporated into the Python random module.

So where does that leave us? Before I dive in, a couple of questions to consider:

• How many streams do I need?
• How long is a “sufficiently long” stream?

How many streams do I need? In a perfect world, the answer is obvious – as many as I want! Failing that, we need to think about the two ways that multiple random number streams are used:

1. Dedicating streams to particular components of the simulation model, which is usually done to facilitate variance reduction during output analysis.
2. Running replications of the simulation in order to (for example) generate confidence intervals on the output measures. This requires an independent stream (or set of streams) for each replication; the more replications, the lower the variance (and the tighter the confidence intervals).

The numbers required is, of course, model and context-dependent.  In my prototype simulator, I’ve assumed up to 100 streams per model, and up to 100 replications – which implies 10,000 independent random number streams.  That’s probably the right order of magnitude for the models I’m thinking of.

How long is “sufficiently long”?  For a given model, this depends on the number of random variates used (model size, essentially) and the desired length of the simulation run.  What might this be for some mythical “maximal” model?  It’s worth a quick back-of-the-envelope estimation…

Let’s consider the M/M/1 queuing model described here.  It takes roughly two seconds to execute a run of one million simulated time units, with an interarrival mean of ten time units.  That translates into roughly 100,000 processes serviced, requiring 200,000 random values (an interarrival time and a service time for each process serviced).  So roughly speaking, this model requires 100,000 random numbers per second, 360 million per hour, or roughly 8.6 billion per day (all on my laptop, of course).  I’m going to go out on a limb and say that if a simulation takes days to execute, than my simulator is probably not the right tool for the job.  So I’ll guess that streams of 100 billion are almost certainly more than sufficient for virtually all cases – though being conservative on these sorts of things, I’d really like to add some zeros to that number.

Say we’d like 10,000 independent random number streams, each with a period of at least 100 billion (say roughly 2³⁷). Given the current limitations of Python’s random module, how should we proceed?

For starters, there’s always the brute force method; call the built-in random() function one quadrillion (10¹⁵) times while saving the state (obtained via Random.getstate())  every 100 billion calls.  When we’re done, we’ll have the states for 10,000 substreams of length 100 billion.  That should only take about four years of runtime on my laptop 🙂 .

A few other options worth considering, if we don’t have the patience to wait that long:

1. Create a separate generator for each required stream using different seeds. This is certainly the easiest solution; unfortunately, there are no guarantees that these streams are independent. While there are differences of opinion, this is not an approach I would use for anything other than demoware (until or unless someone publishes a table of seeds that have been shown to have good results in this regard)
2. Use a different random number generator. We could use a Python implementation (or wrapper for) a generator that does support multiple independent streams.  There is a package wrapping DCMT; unfortunately, the most recent version is pretty old (2011), and there appears to be no support for Python 3.x.  There is another package implementing the Well1024a generator, but it does not include jump-ahead functionality.  In other words, if I want a Python package that provides multiple independent streams by replacing or augmenting the existing random module, it looks like I’ll have to build at least some of it myself.
3. Generate (and save) the MT19937 states for 10,000 streams offline using the jump-ahead algorithm. If I don’t want to go to the trouble of creating, maintaining and delivering a Python wrapper, I could write, cobble together or (better yet) find a non-Python MT19937 implementation that incorporates the jump-ahead algorithm noted above.  We could then run that code offline to generate the state vectors for all 10,000 substreams, and save them to a file of some sort.  At simulation run time, we should then be able to simply create generators from the standard random module and initialize their states using values from this file.  Using jump-ahead allows us to make the substream length considerably longer, perhaps 2⁵⁰.  (The MT19937 generator has a period of 2¹⁹⁹³⁷ – 1.)

Being the intellectually lazy person that I am, my current vote is for option (3). The C++ 2011 standard library includes a <random> module that includes an MT19937 generator.  The API includes a discard() function.  While the default implementation is typically dumb/brute-force (discard(n) invokes the generator n times), the Boost version of the library implements discard() via the aforementioned jump-ahead algorithm (for large values of n).  So it’s time to dust off my rusty C++ skills:

#include <iostream>
#include <sstream>
#include <fstream>
#include <boost/random/mersenne_twister.hpp>
#include <math.h>


using namespace std;

int main(int argc, char** argv)
{
    unsigned long long substream_sz = pow(2,50);
    int nstreams = 10000;
    const char* outfilename = “mt19937_states.dat”;

    boost::random::mt19937 rng;
    ofstream fout(outfilename, ios::out|ios::binary);

    for (int i = 0; i < nstreams; ++i) {
        rng.discard(substream_sz);
        std::stringstream input;
        input << rng;
 
        std::vector<unsigned> state;
        unsigned p;
        while (input >> p) {
            state.push_back(p);
        }

        fout.write(reinterpret_cast<const char*>(&state[0]), 
                   state.size() * sizeof(unsigned));
    }
}

This small program generates the 10,000 substream state vectors and writes them to a binary file.  It runs in about 18 minutes on my laptop.

At this point, I’d like to run that binary file through a short Python script in order to turn it into a portable (machine-independent) NPY format:

state_array = numpy.fromfile(“mt19937_states.dat”)
state_array.reshape(10000, 624)
state_array.save(“mt19937_states.npy”)

And finally, just to close the loop, a snippet of Python code demonstrating how a simulator might initialize random number generator instances for each of these streams:

# Load the .npy file, which contains an nstreams x 624 array 
# of unsigned ints (type numpy.uint32)
nparray = numpy.load(“mt19937_states.npy”)
nstreams = nparray.shape[0]

# Create and initialize a generator for each defined stream
rng = []
for i in range(nstreams):
    # Convert the array of numpy.uint32s to a list of Python ints
    state_list = [int(x) for x in nparray[i]]
    
    # Append the required value of 624 to the list
    state_list.append(624)
    
    # Define the required internalstate tuple and use it to 
    # set the state for a new RNG instance.
    rng_state = (3, tuple(state_list), None)
    rng.append(random.Random())
    rng[-1].setstate(rng_state)

A few notes on the code above:

• The argument to Random.setstate() is the tuple (version, internalstate, gauss_next). version is 3 for the Python 3.x distributions that I use, gauss_next is always None when initializing a Random instance, as near as I can tell.
• internalstate is a tuple containing 635 Python ints – the 624 elements of the MT19937 state vector plus a value of 624 at the end. Random is very particular about types here; internalstate must be a tuple (not a list!), and the elements must be ints (not NumPy types such as uint32). Hence the extra data massaging in the code above.

All in all, I’m pretty comfortable with this approach, at least until the Python standard library (or, more likely, NumPy) incorporates jump-ahead functionality and/or an additional, pluggable random number generator such as DCMT.  And hopefully I haven’t made an utter fool of myself writing about it!