subject:"Numpy Performance"


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


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.
--
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers


sturlamolden wrote:

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.


My parallel version of ziggurat works fine, with Fortran/OpenMP.
--
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers


sturlamolden wrote:

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.



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.

--
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

2009-12-21 Thread Robert Kern


On 2009-12-21 16:57 PM, r0g wrote:

sturlamolden wrote:

On 19 Des, 16:20, Carl Johan Rehncar...@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.


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

--
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers


David Cournapeau wrote:

On Sun, Dec 20, 2009 at 6:47 PM, Lie Ryan lie.1...@gmail.com wrote:

On 12/20/2009 2:53 PM, sturlamolden wrote:

On 20 Des, 01:46, Lie Ryanlie.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


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.

--
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

2009-12-20 Thread Lie Ryan


On 12/20/2009 2:53 PM, sturlamolden wrote:

On 20 Des, 01:46, Lie Ryanlie.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.


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.

--
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

2009-12-20 Thread David Cournapeau

On Sun, Dec 20, 2009 at 6:47 PM, Lie Ryan lie.1...@gmail.com wrote:
 On 12/20/2009 2:53 PM, sturlamolden wrote:

 On 20 Des, 01:46, Lie Ryanlie.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
-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

2009-12-20 Thread Lie Ryan


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.

--
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

2009-12-20 Thread Peter Pearson

On Sun, 20 Dec 2009 07:26:03 +1100, Lie Ryan lie.1...@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?  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.
-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

2009-12-20 Thread Rami Chowdhury


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

 On Sun, 20 Dec 2009 07:26:03 +1100, Lie Ryan lie.1...@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)

-- 
http://mail.python.org/mailman/listinfo/python-list

numpy performance and random numbers

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
-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

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
-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

On Dec 19, 12:29 pm, Steven D'Aprano st...@remove-this-
cybersource.com.au wrote:
 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

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



-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

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.
-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

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.


-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

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.












-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

On Dec 19, 2:49 pm, sturlamolden sturlamol...@yahoo.no wrote:
 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.

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
-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

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
-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

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
-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

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.





-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

On Dec 19, 4:47 pm, sturlamolden sturlamol...@yahoo.no wrote:
 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.


How did you time it?

Well, in Matlab I used tic; for i = 1:1000, randn(100, 1), 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
-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

2009-12-19 Thread Lie Ryan


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

--
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

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.

--
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

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.

-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

On 19 Dec, 23:09, sturlamolden sturlamol...@yahoo.no wrote:
 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.

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
-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

2009-12-19 Thread Gregory Ewing


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
--
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

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, 1), 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, 1)', '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
-- 
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers

2009-12-19 Thread Lie Ryan


On 12/20/2009 8:58 AM, sturlamolden wrote:

On 19 Des, 21:26, Lie Ryanlie.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.

--
http://mail.python.org/mailman/listinfo/python-list

Re: numpy performance and random numbers