116 views

Skip to first unread message

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

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

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

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.

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:

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.

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.

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

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

>

> 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

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

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.

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

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.

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

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.

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.

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

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

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

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.

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

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

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.

Dec 20, 2009, 4:47:39 AM12/20/09

to

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.

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.

> 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

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.

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

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.

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)

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

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.

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.

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

to

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

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.

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.

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.

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?

>

> 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

Search

Clear search

Close search

Google apps

Main menu