Patch with Ziggurat method for Normal distribution (Trac #1450)

[thouis]

Original ticket http://projects.scipy.org/numpy/ticket/1450 on 2010-04-10 by trac user ilan, assigned to @charris.

I've written a patch which replaces the Box-Muller transform
in numpy/random/mtrand/randomkit.c with the faster Ziggurat method.
The patch also includes updates to doc-strings which include the
relevant references:

Doornik, J.A. (2005), "An Improved Ziggurat Method to Generate Normal
Random Samples", mimeo, Nuffield College, University of Oxford, and
www.doornik.com/research.

Marsaglia, G. and Tsang, W. W. (2000), "The Ziggurat Method for
Generating Random Variables", Journal of Statistcal Software 5,
http://www.jstatsoft.org/v05/i08/

Marsaglia, G. (1964) "Generating a variable from the tail of the
normal distribution", Technometrics 6, 101-102.

My patch is based on the support code which J. Doornik has made available
on his web page (first reference). I've contacted him, and he has given
permission to include his code in numpy on the condition that his paper
is referenced (which the patch does).

The patch was written against http://svn.scipy.org/svn/numpy/tags/1.4.0,
and I've just checked to make sure it can still be applied to the
trunk (revision 8324). I've tested the patch on Windows, MacOSX, Linux
and Solaris.

The speedup is a bit over 2 times, when creating a large number of Standard
Normal distributed random numbers. To verify that the method works and
produces good random numbers, I have:

• created a large array, using randn()
• transformed (using the cumulative distribution function) it back to
  uniform distributed random numbers
• written the data to a large file
• performed the diehard tests (I used the dieharder program) on the data

Feel free to hammer the patch. I hope this patch will make its way
into numpy soon.

[thouis]

Attachment added by trac user ilan on 2010-04-10: ziggurat.patch

[thouis]

@charris wrote on 2010-04-10

We've had better ziggurat methods than those in Doornik for some years now, but there has been little push in getting them into numpy.

Hi there.

I've completed timing and robustness tests that should
allow apples to apples comparisons of the mtrand normal distribution
possibilities.  The short story goes like this:

 0) One of the 'Harris' ziggurats, zig7, passes all TestU01 crush tests
   when running atop each of five PRNGS.

 1) Depending on how fast the underlying PRNG is, zig7
   outperforms Box_Muller by ~1.3X to ~8.5X.

 2) using gcc, mwc8222+zig7 outperforms mt19937+Box_Muller by 7X.

See you both at SciPy.

Bruce


And here are the details:

In what follows the uniform variate generators are:
 lcg64
 mwc8222
 mt19937
 mt19937_64
 yarn5

And the normal distribution codes are:
 trng  - default normal distribution code in TRNG
 boxm  - Box-Muller, mtrand lookalike, remembers/uses 2nd value
 zig7  - a 'Harris' ziggurat indexed by 7 bits
 zig8  - a 'Harris' ziggurat indexed by 8 bits
 zig9  - a 'Harris' ziggurat indexed by 9 bits

Here are the numbers in more detail:

# Timings from icc -O2 running on 2.4GhZ Core-2

lcg64 trng: 6.52459e+06 ops per second
lcg64 boxm: 2.18453e+07 ops per second
lcg64 zig7: 1.80616e+08 ops per second
lcg64 zig8: 2.01865e+08 ops per second
lcg64 zig9: 2.05156e+08 ops per second

mwc8222 trng: 6.52459e+06 ops per second
mwc8222 boxm: 2.08787e+07 ops per second
mwc8222 zig7: 9.44663e+07 ops per second
mwc8222 zig8: 1.05326e+08 ops per second
mwc8222 zig9: 1.03478e+08 ops per second

mt19937 trng: 6.41112e+06 ops per second
mt19937 boxm: 1.64986e+07 ops per second
mt19937 zig7: 4.23762e+07 ops per second
mt19937 zig8: 4.52623e+07 ops per second
mt19937 zig9: 4.52623e+07 ops per second

mt19937_64 trng: 6.42509e+06 ops per second
mt19937_64 boxm: 1.93226e+07 ops per second
mt19937_64 zig7: 5.8762e+07 ops per second
mt19937_64 zig8: 6.17213e+07 ops per second
mt19937_64 zig9: 6.29146e+07 ops per second

yarn5 trng: 5.95781e+06 ops per second
yarn5 boxm: 1.19156e+07 ops per second
yarn5 zig7: 1.48945e+07 ops per second
yarn5 zig8: 1.54809e+07 ops per second
yarn5 zig9: 1.53201e+07 ops per second


# Timings from g++ -O2 running on a 2.4GhZ Core-2

lcg64 trng: 6.72163e+06 ops per second
lcg64 boxm: 1.50465e+07 ops per second
lcg64 zig7: 1.31072e+08 ops per second
lcg64 zig8: 1.48383e+08 ops per second
lcg64 zig9: 1.6036e+08 ops per second

mwc8222 trng: 6.64215e+06 ops per second
mwc8222 boxm: 1.44299e+07 ops per second
mwc8222 zig7: 8.903e+07 ops per second
mwc8222 zig8: 1.00825e+08 ops per second
mwc8222 zig9: 1.03478e+08 ops per second

mt19937 trng: 6.52459e+06 ops per second
mt19937 boxm: 1.28223e+07 ops per second
mt19937 zig7: 5.00116e+07 ops per second
mt19937 zig8: 5.41123e+07 ops per second
mt19937 zig9: 5.47083e+07 ops per second

mt19937_64 trng: 6.58285e+06 ops per second
mt19937_64 boxm: 1.42988e+07 ops per second
mt19937_64 zig7: 6.72164e+07 ops per second
mt19937_64 zig8: 7.39591e+07 ops per second
mt19937_64 zig9: 7.46022e+07 ops per second

yarn5 trng: 6.25144e+06 ops per second
yarn5 boxm: 8.93672e+06 ops per second
yarn5 zig7: 1.50465e+07 ops per second
yarn5 zig8: 1.57496e+07 ops per second
yarn5 zig9: 1.56038e+07 ops per second

You should contact Bruce Carneal bcarneal@gmail.com and see if you can get his code. It is in c++ but can probably be translated to c without too much trouble. If we are going the ziggurat direction, it would be nice to implement the exponential distribution also. As you can see, the overall speedup also depends on the underlying rng, so if you can figure out a way of making choices available there that will also help. Of course, for small sets of values call overhead will dominate everything.

Thanks for looking into this problem, all it needs is for someone to take it in hand and push it through.

[thouis]

trac user ilan wrote on 2010-04-10

Thanks for your quick response. Yes, the speedup depends on the underlying RNG, and the patch I propose is just using what is already there. The method is quite simple and already a more than 2 times speedup. My felling is that by using faster underlying RNGs, a lot more could be done, but this is really a separate issue.

Writing the Ziggurat method more general, such that exponential (and possibly other) distributions can be done with the same code, might be possible but I don't think its worth the effort, as the code is very small, and also the tail distribution is different in other cases.

[thouis]

@charris wrote on 2010-04-10

The 'Harris' ziggurat is mine and the code is actually simpler than the Marsaglia version (Doornik). I also did the exponential (as did Marsaglia), which is even easier. Notice that the Marsaglia version actually uses an exponential to bound the tail. So do I, but I make it part of the ziggurat, which simplifies matters. It takes maybe ten lines of python to compute the needed constants for the various ziggurats. Note that in general you can get about 4x speedup vs Box-Muller for the MT rng.

[thouis]

@charris wrote on 2010-04-10

A few additional comments.

1. The random variable

u = 2.0 * rk_double(state) - 1.0

isn't evenly distributed between positive and negative. I don't think that makes any difference in practice given the other roundoff errors, but is worth noting because the steps in the ziggurat are symmetric about zero. Accepting the roundoff errors, it can be fixed by shifting and scaling the ziggurat layer to [0, 1-eps], using the positive double, and shifting and scaling the return value. This also avoids the call to fabs and saves one floating multiplication.

       /* first try the rectangular boxes */
       if (fabs(u) < s_adZigR[i])
           return u * s_adZigX[i];

uses a double in the comparison. One of the virtues of the original ziggurat was using integers here, which is most likely faster even on modern hardware.

3. It would be nice to move the initialization to load time so that there is no need to check it for every call. That isn't likely to make much difference given the overhead of the call to the rng but it looks a bit out of place.

4. The function zig_NormalTail can be removed if an improved ziggurat is used.

5. You can avoid the occasional extra call to rk_random if you handle the double conversion yourself. There are enough bits left over from the double to deal with the ziggurat layers. Having the pre-conversion integer form will also allow you to use integer comparisons in 2).

Chuck

[charris]

IIRC, this was in the Enthought edition for a while, but messed up results since variates were not produced in the same order as the current algorithm. Adding something like this either needs a new name, or a method argument for the normal generators.

[bashtage]

@thouis @charris

Is this patch still active? NumPy's random number generation performance is somewhat behind other choices (Julia or MATLAB, around a factor of 4), and this seems like a simple method to improve the performance for Gaussians which are both directly interesting and important inputs for random numbers from other distributions.

[charris]

No, but we could take it up again. The current generator needs to be kept for backward compatibility, so a new one would need a new name.

[argriffing]

It seems overzealous to require backward compatibility for the specific numbers sampled by the random number generators. Everything else equal of course it's better, but if sampling can be sped up 2x or 4x then there is a tradeoff.

[rgommers]

@argriffing it's not unimportant for reproducibility of unit tests using random number generators (same reason that a seed needs to be used in such tests).

[njsmith]

If anyone really wants the speedup, it's pretty easy to satisfy both
impulses. Save the old generator as a subclass, and write a function that
takes a numpy version number and fetches a RandomState instance with the
corresponding implementations. Then people who need to reproduce old
results have an easy way to do it, and people who don't care can get the
fastest version.
On 19 Jul 2014 15:29, "Ralf Gommers" notifications@github.com wrote:

@argriffing https://github.com/argriffing it's not unimportant for
reproducibility of unit tests using random number generators (same reason
that a seed needs to be used in such tests).

—
Reply to this email directly or view it on GitHub
#2047 (comment).

[charris]

We do require that the same sequences be generated by the same functions, it is even tested. Reproducibility requires that. But there is no reason not to add new functions, say zignormal. Bruce disappeared after all that work he did to test and verify the algorithms. I still have my code sitting around somewhere, although it is pretty ancient now. IIRC, there was a small bug that Bruce found, overflow or some such, which I think I fixed in my code. Hmm, I've still got all the emails from the summer of 2008. Looks like Bruce also had an SIMD version which sped things up even more.

PS: As before, the scaling, relative to a single serial thread, is 16X
(4 wide SIMD X 4 cores).

So if we really wanted fast, that would be one way to go. Multiple cores may make preserving the sequence difficult. Super high performance can probably be pushed off to a specialized package.

[seberg]

I also don't like changing random state for some speed up (no matter how awesome). Now personally, I doubt that will happen, but I used to write scripts with a fixed seed for reproducibility. Some of that reproducibility is gone anyway because even tiny float rounding changes will break it, but still. I don't mind having a method keyword. I also don't mind adding logic such as "If the seed was fixed, warn the user of a future change and use the old version, otherwise just use the new version", then keep the old version around and after a while make the faster one default. And maybe right away or a little later remove the warning.

[charris]

Bruce is on Linkedin. I never got any of my code back, so it may be worth contacting him first if we decide to pursue this.