numpy performance and random numbers

[Lie Ryan]

Except there is no way to find two very distant states and prove they
are distant enough.

Except only theoretical scientist feel the need to prove it and perhaps
perhaps for cryptographic-level security. Random number for games,
random number for tmp files, and 99.99% random number users doesn't
really need such proves.

And don't forget the effect of the very long period of PRNG like
Mersenne Twister (2**19937 − 1, according to Wikipedia) makes it very
unlikely to choose two different seeds and ended up in nearby entry
point. Let's just assume we're using a 2**64-bit integer as the seeds
and let's assume the entry point defined by these seeds are uniformly
distributed you would to generate (on average) 2.34E+5982 numbers before
you clashes with the nearest entry point. Such amount of numbers would
require decades to generate. Your seed generator guard would only need
to prevent seeding parallel generators with the same seed (which is
disastrous), and that's it, the big period covers everything else.

 

[Lie Ryan]

But the OP case mostly like falls in your estimated 0.01% case. PRNG
quality is essential for reliable Monte Carlo procedures. I don't
think long period is enough to guarantee those good properties for //
random generators - at least it is not obvious to me.

Now it's not, long periods are not indicator of quality. I was
responding to the chance of unexpected repetition of sequence because of
collision of entry points. Long periods is an indicator that the chance
of entry point collision should be low enough. Long periods (alone)
doesn't mean anything to the quality of the randomness itself.

 

[Peter Pearson]

If you don't care about "repeatability" (which is already extremely
difficult in parallel processing even without random number generators),
you can just start two PRNG at two distinct states (and probably from
two different algorithms) and they will each spews out two independent
streams of random numbers. What was "unsolved" was the "pseudo-" part of
the random number generation, which guarantee perfect replayability in
all conditions.

Why not use a good cipher, such as AES, to generate a pseudorandom
bit stream by encrypting successive integers? If you use a different
key for each instance, you'll have exactly the independence you want.
And if you can detect any statistical anomaly in the output, you
automatically have the most exciting paper to be published in the next
issue of the Journal of Cryptology.

Minor caveats:

1. This might be slower than another approach, but maybe not:
AES is much faster than the ciphers of the old days.

2. Since AES(key,i) != AES(key,j) if i != j, there will be a
dearth of duplicates, which will become statistically
detectable around the time i has been incremented
2**64 times. There are many reasons why this might not
bother you (one: you don't plan to use so many values;
two: you might use the 128-bit AES output in pieces, rather
than as a single value, in which case duplicates will
appear among the pieces at the right rate), but if it *does*
bother you, it can be fixed by using i^AES(key,i) instead
of just AES(key,i).

 

[Rami Chowdhury]

Why not use a good cipher, such as AES, to generate a pseudorandom
bit stream by encrypting successive integers?

Isn't the Fortuna PRNG based around that approximate concept?

 

[Robert Kern]

Thank you, this was very informative. I know about the Mersenne
Twister, but had no idea about Marsaglia's ziggurat method. I will
definitely try f2py or cython.

It's worth noting that the ziggurat method can be implemented incorrectly, and
requires some testing before I will accept such in numpy. That's why we don't
use the ziggurat method currently. C.f.

http://www.doornik.com/research/ziggurat.pdf

Requests for the ziggurat method come up occasionally on the numpy list, but so
far no one has shown code or test results. Charles Harris and Bruce Carneal seem
to have gotten closest. You can search numpy-discussion's archive for their
email addresses and previous threads. Additionally, you should ask future numpy
questions on the numpy-discussion mailing list.

http://www.scipy.org/Mailing_Lists

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco

 

[r0g]

The problem with PRNG is that they are iterative in nature, and
maintain global states. They are therefore very hard to vectorize. A
GPU will not help. The GPU has hundreds of computational cores that
can run kernels, but you only get to utilize one.

Parallel PRNGs are an unsolved problem in computer science.

Surely you could have as many totally independent cores as you like
independently spitting out random bits whenever they feel like it to
make an equally random bitstream? would have thought the only issue
would be ensuring high quality bitstream was used to seed each thread no?

Roger.

 

[Gib Bogle]

Common ancenstry, NumPy and Matlab borrowed the name from IDL.


LabView, Octave and SciLab uses the name randn as well.



The obvious thing would be to compile ziggurat yourself, and turn on
optimization flags for your hardware.
http://www.jstatsoft.org/v05/i08/


P.S. Be careful if you consider using more than one processor.
Multithreading is a very difficult issue with PRNGs, becuase it is
difficult to guarrantee they are truely independent. But you can use a
producer-consumer pattern, though: one thread constantly producing
random numbers (writing into a buffer or pipe) and another thread(s)
consuming them.

In case you're interested, I've made a multi-thread version of ziggurat
(actually in Fortran for use with OpenMP). It is a simple enough mod, there is
an additional argument (thread number) in each call to the ziggurat functions,
and ziggurat maintains separate variables for each thread (not just the seed).
There was one non-trivial issue. To avoid cache collisions, the seed values (an
array corresponding to jsr in the original code) need to be spaced sufficiently
far apart. Without this measure the performance was disappointingly slow.

I agree with your recommendation - ziggurat is really fast.

 

[Gib Bogle]

The problem with PRNG is that they are iterative in nature, and
maintain global states. They are therefore very hard to vectorize. A
GPU will not help. The GPU has hundreds of computational cores that
can run kernels, but you only get to utilize one.

Parallel PRNGs are an unsolved problem in computer science.

My parallel version of ziggurat works fine, with Fortran/OpenMP.

 

[Gib Bogle]

Another note on this:

Ideally, we would e.g. know how to find (analytically) MT states that
are very far apart. But to my knowledge no such equation has been
derived.

But often in Monte Carlo simulations, the PRNG is not the dominant
computational bottleneck. So we can simply start N PRNGs from N
consequtive states, and for each PRNG only use every N-th pseudorandom
deviate.

OK, now I understand. In my application I don't care about a thread's PRNG
reaching the start of another threads. I do understand that this is not the
case in all applications.

 

[Gib Bogle]

But the OP case mostly like falls in your estimated 0.01% case. PRNG
quality is essential for reliable Monte Carlo procedures. I don't
think long period is enough to guarantee those good properties for //
random generators - at least it is not obvious to me.

David

My simulation program is Monte Carlo, but the complexity and variety of all the
interactions ensure that PRNG sequence overlap will have no discernible effect.

 

[Robert Kern]

Surely you could have as many totally independent cores as you like
independently spitting out random bits whenever they feel like it to
make an equally random bitstream? would have thought the only issue
would be ensuring high quality bitstream was used to seed each thread no?

No. For most quality PRNGs, all seeds are equal (with a few minor exceptions.
E.g. for the Mersenne Twister, a state vector of all zeros will yield only
zeros, but any nontrivial state vector puts the PRNG into the same orbit, just
at different places). There is no notion of a "high quality seed". The problem
is that during the run, the separate PRNGs may overlap, which will reduce the
independence of the samples.

That said, the enormous length of the Mersenne Twister's period helps a lot. You
can use an ungodly number of streams and run length without having a physically
realizable chance of overlap. The chance of having a bug in your simulation code
is overwhelmingly greater.

There are also algorithms that can initialize a given number of PRNGs with
different parameters (*not* seeds) that are guaranteed not to overlap. No one
has implemented this for numpy, yet.

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco

 

[Peter Pearson]

Isn't the Fortuna PRNG based around that approximate concept?

Yes, although the interesting part of Fortuna is its strategy for
entropy accumulation, which is (I think) not of interest to the
Original Poster. In cryptology, one finds many places where a
statistically significant departure from random behavior would
be a publishable discovery.