numpy performance and random numbers

116 views
Skip to first unread message

Carl Johan Rehn

unread,
Dec 19, 2009, 5:05:17 AM12/19/09
to
Dear friends,

I plan to port a Monte Carlo engine from Matlab to Python. However,
when I timed randn(N1, N2) in Python and compared it with Matlab's
randn, Matlab came out as a clear winner with a speedup of 3-4 times.
This was truly disappointing. I ran tthis test on a Win32 machine and
without the Atlas library.

Why is there such a large difference in speed and how can I improve
the speed of randn in Python! Any help with this matter is truly
appreciated since I really would like to get away from Matlab and move
over to Python instead.

Yours Carl

Steven D'Aprano

unread,
Dec 19, 2009, 6:29:18 AM12/19/09
to
On Sat, 19 Dec 2009 02:05:17 -0800, Carl Johan Rehn wrote:

> Dear friends,
>
> I plan to port a Monte Carlo engine from Matlab to Python. However, when
> I timed randn(N1, N2) in Python and compared it with Matlab's randn,

What's randn? I don't know that function. I know the randint, random, and
randrange functions, but not randn. Where does it come from?


> Matlab came out as a clear winner with a speedup of 3-4 times. This was
> truly disappointing. I ran tthis test on a Win32 machine and without the
> Atlas library.
>
> Why is there such a large difference in speed and how can I improve the
> speed of randn in Python! Any help with this matter is truly appreciated
> since I really would like to get away from Matlab and move over to
> Python instead.

Could be many reasons. Python could be generally slower than Matlab. Your
timing code might have been faulty and you weren't comparing equal
amounts of work (if your Python code was doing four times as much work as
the Matlab code, then naturally it will be four times slower). Perhaps
the Matlab random number generator is a low-quality generator which is
fast but not very random. Python uses a very high quality RNG which is
not cheap.

But does it really matter if the RNG is slower? Your Monte Carlo engine
is a lot more than just a RNG. What matters is whether the engine as a
whole is faster or slower, not whether one small component is slower.


--
Steven

Carl Johan Rehn

unread,
Dec 19, 2009, 8:06:53 AM12/19/09
to
On Dec 19, 12:29 pm, Steven D'Aprano <st...@REMOVE-THIS-

randn is given by

>> import numpy
>>> numpy.random.randn(2,3)
array([[-2.66435181, -0.32486419, 0.12742156],
[-0.2387061 , -0.55894044, 1.20750493]])

Generally, at least in my MC application, I need a large number of
random numbers. Usually I execute, for example, r = randn(100, 10000)
sequentially a relatively large number of times until sufficient
accuracy has been reached. Thus, randn is in my case a mission
critical component for obtaining an acceptable overall run time.
Matlab and numpy have (by chance?) the exact names for the same
functionality, so I was very happy with numpy's implementation until I
timed it. So the basioc question is, how can I speed up random number
generation?

Carl

sturlamolden

unread,
Dec 19, 2009, 8:49:17 AM12/19/09
to
On 19 Des, 11:05, Carl Johan Rehn <car...@gmail.com> wrote:

> I plan to port a Monte Carlo engine from Matlab to Python. However,
> when I timed randn(N1, N2) in Python and compared it with Matlab's
> randn, Matlab came out as a clear winner with a speedup of 3-4 times.
> This was truly disappointing. I ran tthis test on a Win32 machine and
> without the Atlas library.

This is due to the algorithm. Matlab is using Marsaglia's ziggurat
method. Is is the fastest there is for normal and gamma random
variates. NumPy uses the Mersenne Twister to produce uniform random
deviates, and then applies trancendental functions to transform to the
normal distribution. Marsaglia's C code for ziggurat is freely
available, so you can compile it yourself and call from ctypes, Cython
or f2py.

The PRNG does not use BLAS/ATLAS.

sturlamolden

unread,
Dec 19, 2009, 8:53:52 AM12/19/09
to
On 19 Des, 12:29, Steven D'Aprano <st...@REMOVE-THIS-
cybersource.com.au> wrote:

> Perhaps
> the Matlab random number generator is a low-quality generator which is
> fast but not very random. Python uses a very high quality RNG which is
> not cheap.

Marsaglia and Matlab's implementation of ziggurat uses a slightly
lower quality RNG for uniform deviates than NumPy's Mersenne Twister.
But the real speed advantage comes from avoiding trancendental
functions. I have for some time thought of contributing a ziggurat
generator to NumPy, while retaining the Mersenne Twister internally,
but I have not got around to do it.


sturlamolden

unread,
Dec 19, 2009, 9:16:16 AM12/19/09
to
On 19 Des, 14:06, Carl Johan Rehn <car...@gmail.com> wrote:

> Matlab and numpy have (by chance?) the exact names for the same
> functionality,

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


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


> So the basioc question is, how can I speed up random number
> generation?

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.


Carl Johan Rehn

unread,
Dec 19, 2009, 10:14:16 AM12/19/09
to

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.

Well, I guess I knew that random numbers were not handled by BLAS/
ATLAS, but wasn't 100% sure.

Carl

Carl Johan Rehn

unread,
Dec 19, 2009, 10:20:10 AM12/19/09
to
On Dec 19, 3:16 pm, sturlamolden <sturlamol...@yahoo.no> wrote:
> On 19 Des, 14:06, Carl Johan Rehn <car...@gmail.com> wrote:
>
> > Matlab and numpy have (by chance?) the exact names for the same
> > functionality,
>
> Common ancenstry, NumPy and Matlab borrowed the name from IDL.
>
> LabView, Octave and SciLab uses the name randn as well.
>
> > So the basioc question is, how can I speed up random number
> > generation?
>
> 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.

How about mulit-core or (perhaps more exciting) GPU and CUDA? I must
admit that I am extremely interested in trying the CUDA-alternative.

Obviously, cuBLAS is not an option here, so what is the safest route
for a novice parallel-programmer?

Carl


Carl

Steven D'Aprano

unread,
Dec 19, 2009, 10:22:45 AM12/19/09
to
On Sat, 19 Dec 2009 05:06:53 -0800, Carl Johan Rehn wrote:

> so I was very happy with numpy's implementation until I timed it.

How did you time it?


--
Steven

sturlamolden

unread,
Dec 19, 2009, 10:47:48 AM12/19/09
to
On 19 Des, 16:20, Carl Johan Rehn <car...@gmail.com> wrote:

> How about mulit-core or (perhaps more exciting) GPU and CUDA? I must
> admit that I am extremely interested in trying the CUDA-alternative.
>
> Obviously, cuBLAS is not an option here, so what is the safest route
> for a novice parallel-programmer?

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.

Carl Johan Rehn

unread,
Dec 19, 2009, 12:02:38 PM12/19/09
to

>>>How did you time it?

Well, in Matlab I used "tic; for i = 1:1000, randn(100, 10000), end;
toc" and in IPython i used a similar construct but with "time" instead
of tic/(toc.


>>> Parallel PRNGs are an unsolved problem in computer science.

Thanks again for sharing your knowledge. I had no idea. This means
that if I want to speed up my application I have to go for the fastest
random generator and focus on other parts of my code that can be
vectorized.

Carl

Lie Ryan

unread,
Dec 19, 2009, 3:26:03 PM12/19/09
to
On 12/20/2009 4:02 AM, Carl Johan Rehn wrote:
>>>> How did you time it?
>
> Well, in Matlab I used "tic; for i = 1:1000, randn(100, 10000), end;
> toc" and in IPython i used a similar construct but with "time" instead
> of tic/(toc.

Code?

>>>> Parallel PRNGs are an unsolved problem in computer science.
>
> Thanks again for sharing your knowledge. I had no idea. This means
> that if I want to speed up my application I have to go for the fastest
> random generator and focus on other parts of my code that can be
> vectorized.

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.

sturlamolden

unread,
Dec 19, 2009, 4:58:37 PM12/19/09
to
On 19 Des, 21:26, Lie Ryan <lie.1...@gmail.com> wrote:

> you can just start two PRNG at two distinct states

No. You have to know for certain that the outputs don't overlap.

If you pick two random states (using any PRNG), you need error-
checking that states are always unique, i.e. that each PRNG never
reaches the starting state of the other(s). If I had to do this, I
would e.g. hash the state to a 32 bit integer. You can deal with
errors by jumping to a different state or simply aborting the
simulation. The problem is that the book-keeping will be costly
compared to the work involved in generating a single pseudorandom
deviate. So while we can construct a parallel PRNG this way, it will
likely run slower than one unique PRNG.

It has been suggested to use short-period MTs with different
characteristic polynomials. This is also available for the nVidia GPU
using CUDA. This is probably the generator I would pick if I wanted to
produce random numbers in parallel. However, here we should be aware
that the math has not been thoroughly proven.

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/DC/dgene.pdf
http://developer.download.nvidia.com/compute/cuda/2_2/sdk/website/projects/MersenneTwister/doc/MersenneTwister.pdf

There is also a version of MT that use SIMD instructions internally,
which gives a factor-of-two speedup.

I would never use different algorithms. It will introduce a systematic
difference in the simulation.

sturlamolden

unread,
Dec 19, 2009, 5:09:54 PM12/19/09
to
On 19 Des, 22:58, sturlamolden <sturlamol...@yahoo.no> wrote:

> If you pick two random states (using any PRNG), you need error-
> checking that states are always unique, i.e. that each PRNG never
> reaches the starting state of the other(s).

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.

Carl Johan Rehn

unread,
Dec 19, 2009, 5:57:15 PM12/19/09
to

Thank you for pointing me to the short-period MT reference and
especially the reference on the CUDA-version of parallel MT (even
though I would have wished the author had included a benchmark
comparison in the report). This is a very interesting topic. I agree
that it may work to start PRNGs at distinct and different states, but
that bookkeeping may slow down the algorithm so that it is not worth
the effort. However, the CUDA-version sounds interesting and should be
easy enough to use in a practical application.

Carl

Gregory Ewing

unread,
Dec 19, 2009, 7:04:10 PM12/19/09
to
Lie Ryan wrote:

> 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)

There's no need for different algorithms, and no problem
with repeatability. There exist algorithms with a very
long period (on the order of 2*100 or more) for which it's
easy to calculate seeds for a set of guaranteed non-
overlapping subsequences.

http://or.journal.informs.org/cgi/content/abstract/44/5/816
http://or.journal.informs.org/cgi/content/abstract/47/1/159

In these kinds of generator, moving from one state to
the next involves multiplying a state vector by a matrix
using modulo arithmetic. So you can jump ahead N steps
by raising the matrix to the power of N.

--
Greg

Steven D'Aprano

unread,
Dec 19, 2009, 7:23:46 PM12/19/09
to
On Sat, 19 Dec 2009 09:02:38 -0800, Carl Johan Rehn wrote:

> Well, in Matlab I used "tic; for i = 1:1000, randn(100, 10000), end;
> toc" and in IPython i used a similar construct but with "time" instead
> of tic/(toc.


I don't know if this will make any significant difference, but for the
record that is not the optimal way to time a function in Python due to
the overhead of creating the loop sequence.

In Python, the typical for-loop construct is something like this:

for i in range(1, 1000)

but range isn't a funny way of writing "loop over these values", it
actually creates a list of integers, which has a measurable cost. In
Python 2.x you can reduce that cost by using xrange instead of range, but
it's even better to pre-compute the list outside of the timing code. Even
better still is to use a loop sequence like [None]*1000 (precomputed
outside of the timing code naturally!) to minimize the cost of memory
accesses.

The best way to perform timings of small code snippets in Python is using
the timeit module, which does all these things for you. Something like
this:

from timeit import Timer
t = Timer('randn(100, 10000)', 'from numpy import randn')
print min(t.repeat())

This will return the time taken by one million calls to randn. Because of
the nature of modern operating systems, any one individual timing can be
seriously impacted by other processes, so it does three independent
timings and returns the lowest. And it will automatically pick the best
timer to use according to your operating system (either clock, which is
best on Windows, or time, which is better on Posix systems).

--
Steven

Lie Ryan

unread,
Dec 19, 2009, 7:46:37 PM12/19/09
to
On 12/20/2009 8:58 AM, sturlamolden wrote:
> On 19 Des, 21:26, Lie Ryan<lie.1...@gmail.com> wrote:
>
>> you can just start two PRNG at two distinct states
>
> No. You have to know for certain that the outputs don't overlap.

Not necessarily, you only need to be certain that the two streams don't
overlap in any reasonable amount of time. For that purpose, you can use
a PRNG that have extremely high period like Mersenne Twister and puts
the generators to very distant states.

> If you pick two random states (using any PRNG), you need error-
> checking that states are always unique, i.e. that each PRNG never
> reaches the starting state of the other(s). If I had to do this, I
> would e.g. hash the state to a 32 bit integer. You can deal with
> errors by jumping to a different state or simply aborting the
> simulation. The problem is that the book-keeping will be costly
> compared to the work involved in generating a single pseudorandom
> deviate. So while we can construct a parallel PRNG this way, it will
> likely run slower than one unique PRNG.

You will need some logic to ensure that your generator's states are
distant enough, but that won't incur any generation overhead. Of course
this relies on the very large period of the PRNG algorithm.

Steven D'Aprano

unread,
Dec 19, 2009, 7:55:48 PM12/19/09
to
On Sat, 19 Dec 2009 13:58:37 -0800, sturlamolden wrote:

> On 19 Des, 21:26, Lie Ryan <lie.1...@gmail.com> wrote:
>
>> you can just start two PRNG at two distinct states
>
> No. You have to know for certain that the outputs don't overlap.

"For certain"? Why?

Presumably you never do a Monte Carlo simulation once, you always repeat
the simulation. Does it matter if (say) one trial in fifty is lower-
quality than the rest because the outputs happened to overlap? What if
it's one trial in fifty thousand? So long as the probability of overlap
is sufficiently small, you might not even care about doing repeated
simulations.


> If you pick two random states (using any PRNG), you need error- checking
> that states are always unique, i.e. that each PRNG never reaches the
> starting state of the other(s). If I had to do this, I would e.g. hash
> the state to a 32 bit integer. You can deal with errors by jumping to a
> different state or simply aborting the simulation. The problem is that
> the book-keeping will be costly compared to the work involved in
> generating a single pseudorandom deviate. So while we can construct a
> parallel PRNG this way, it will likely run slower than one unique PRNG.

Since truly random sequences sometimes produce repeating sequences, you
should be able to tolerate short periods of overlap. I'd consider the
following strategy: for each of your parallel PRNGs other than the first,
periodically jump to another state unconditionally. The periods should be
relatively prime to each other, e.g. the second PRNG jumps-ahead after
(say) 51 calls, the third after 37 calls, etc. (the exact periods
shouldn't be important). Use a second, cheap, PRNG to specify how far to
jump ahead. The overhead should be quite small: a simple counter and test
for each PRNG, plus a small number of calls to a cheap PRNG, and
statistically you should expect no significant overlap.

Is this workable?

--
Steven

sturlamolden

unread,
Dec 19, 2009, 10:53:28 PM12/19/09
to
On 20 Des, 01:46, Lie Ryan <lie.1...@gmail.com> wrote:

> Not necessarily, you only need to be certain that the two streams don't
> overlap in any reasonable amount of time. For that purpose, you can use
> a PRNG that have extremely high period like Mersenne Twister and puts
> the generators to very distant states.

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

Lie Ryan

unread,
Dec 20, 2009, 4:47:39 AM12/20/09
to
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.

David Cournapeau

unread,
Dec 20, 2009, 9:13:57 AM12/20/09
to Lie Ryan, pytho...@python.org
On Sun, Dec 20, 2009 at 6:47 PM, Lie Ryan <lie....@gmail.com> wrote:
> On 12/20/2009 2:53 PM, sturlamolden wrote:
>>
>> On 20 Des, 01:46, Lie Ryan<lie.1...@gmail.com>  wrote:
>>
>>> Not necessarily, you only need to be certain that the two streams don't
>>> overlap in any reasonable amount of time. For that purpose, you can use
>>> a PRNG that have extremely high period like Mersenne Twister and puts
>>> the generators to very distant states.
>>
>> 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.

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

Lie Ryan

unread,
Dec 20, 2009, 7:18:38 PM12/20/09
to
On 12/21/2009 1:13 AM, David Cournapeau wrote:
> 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

unread,
Dec 20, 2009, 8:41:17 PM12/20/09
to

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).

--
To email me, substitute nowhere->spamcop, invalid->net.

Rami Chowdhury

unread,
Dec 21, 2009, 12:02:02 AM12/21/09
to Peter Pearson, pytho...@python.org

On Dec 20, 2009, at 17:41 , Peter Pearson wrote:

> On Sun, 20 Dec 2009 07:26:03 +1100, Lie Ryan <lie....@gmail.com> wrote:
>> On 12/20/2009 4:02 AM, Carl Johan Rehn wrote:
>>
>>>>>> Parallel PRNGs are an unsolved problem in computer science.
>>>
>>> Thanks again for sharing your knowledge. I had no idea. This means
>>> that if I want to speed up my application I have to go for the fastest
>>> random generator and focus on other parts of my code that can be
>>> vectorized.
>>
>> 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?

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

-------------
Rami Chowdhury
"Never assume malice when stupidity will suffice." -- Hanlon's Razor
408-597-7068 (US) / 07875-841-046 (UK) / 0189-245544 (BD)

Robert Kern

unread,
Dec 21, 2009, 11:40:08 AM12/21/09
to pytho...@python.org

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

unread,
Dec 21, 2009, 5:57:32 PM12/21/09
to

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

unread,
Dec 21, 2009, 6:22:52 PM12/21/09
to

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

unread,
Dec 21, 2009, 6:23:56 PM12/21/09
to

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

Gib Bogle

unread,
Dec 21, 2009, 6:26:55 PM12/21/09
to

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

unread,
Dec 21, 2009, 6:40:24 PM12/21/09
to

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

unread,
Dec 21, 2009, 6:46:14 PM12/21/09
to pytho...@python.org

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.

Peter Pearson

unread,
Dec 22, 2009, 9:15:06 PM12/22/09
to
On Sun, 20 Dec 2009 21:02:02 -0800, Rami Chowdhury wrote:
>
> On Dec 20, 2009, at 17:41 , Peter Pearson wrote:
>
>> 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?

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.

Reply all
Reply to author
Forward
0 new messages