Ziggurat algorithm

Last updated

The ziggurat algorithm is an algorithm for pseudo-random number sampling. Belonging to the class of rejection sampling algorithms, it relies on an underlying source of uniformly-distributed random numbers, typically from a pseudo-random number generator, as well as precomputed tables. The algorithm is used to generate values from a monotonically decreasing probability distribution. It can also be applied to symmetric unimodal distributions, such as the normal distribution, by choosing a value from one half of the distribution and then randomly choosing which half the value is considered to have been drawn from. It was developed by George Marsaglia and others in the 1960s.

Contents

A typical value produced by the algorithm only requires the generation of one random floating-point value and one random table index, followed by one table lookup, one multiply operation and one comparison. Sometimes (2.5% of the time, in the case of a normal or exponential distribution when using typical table sizes)[ citation needed ] more computations are required. Nevertheless, the algorithm is computationally much faster[ citation needed ] than the two most commonly used methods of generating normally distributed random numbers, the Marsaglia polar method and the Box–Muller transform, which require at least one logarithm and one square root calculation for each pair of generated values. However, since the ziggurat algorithm is more complex to implement it is best used when large quantities of random numbers are required.

The term ziggurat algorithm dates from Marsaglia's paper with Wai Wan Tsang in 2000; it is so named because it is conceptually based on covering the probability distribution with rectangular segments stacked in decreasing order of size, resulting in a figure that resembles a ziggurat.

The Ziggurat algorithm used to generate sample values with a normal distribution. (Only positive values are shown for simplicity.) The pink dots are initially uniform-distributed random numbers. The desired distribution function is first segmented into equal areas "A". One layer i is selected at random by the uniform source at the left. Then a random value from the top source is multiplied by the width of the chosen layer, and the result is x tested to see which region of the layer it falls into with 3 possible outcomes: 1) (left, solid black region) the sample is clearly under the curve and may be output immediately, 2) (right, vertically striped region) the sample value may lie under the curve, and must be tested further. In that case, a random y value within the chosen layer is generated and compared to f(x). If less, the point is under the curve and the value x is output. If not, (the third case), the chosen point x is rejected and the algorithm is restarted from the beginning. Ziggurat method.gif
The Ziggurat algorithm used to generate sample values with a normal distribution. (Only positive values are shown for simplicity.) The pink dots are initially uniform-distributed random numbers. The desired distribution function is first segmented into equal areas "A". One layer i is selected at random by the uniform source at the left. Then a random value from the top source is multiplied by the width of the chosen layer, and the result is x tested to see which region of the layer it falls into with 3 possible outcomes: 1) (left, solid black region) the sample is clearly under the curve and may be output immediately, 2) (right, vertically striped region) the sample value may lie under the curve, and must be tested further. In that case, a random y value within the chosen layer is generated and compared to f(x). If less, the point is under the curve and the value x is output. If not, (the third case), the chosen point x is rejected and the algorithm is restarted from the beginning.

Theory of operation

The ziggurat algorithm is a rejection sampling algorithm; it randomly generates a point in a distribution slightly larger than the desired distribution, then tests whether the generated point is inside the desired distribution. If not, it tries again. Given a random point underneath a probability density curve, its x coordinate is a random number with the desired distribution.

The distribution the ziggurat algorithm chooses from is made up of n equal-area regions; n  1 rectangles that cover the bulk of the desired distribution, on top of a non-rectangular base that includes the tail of the distribution.

Given a monotone decreasing probability density function f(x), defined for all x ≥ 0, the base of the ziggurat is defined as all points inside the distribution and below y1 = f(x1). This consists of a rectangular region from (0, 0) to (x1, y1), and the (typically infinite) tail of the distribution, where x > x1 (and y < y1).

This layer (call it layer 0) has area A. On top of this, add a rectangular layer of width x1 and height A/x1, so it also has area A. The top of this layer is at height y2 = y1 + A/x1, and intersects the density function at a point (x2, y2), where y2 = f(x2). This layer includes every point in the density function between y1 and y2, but (unlike the base layer) also includes points such as (x1, y2) which are not in the desired distribution.

Further layers are then stacked on top. To use a precomputed table of size n (n = 256 is typical), one chooses x1 such that xn = 0, meaning that the top box, layer n − 1, reaches the distribution's peak at (0, f(0)) exactly.

Layer i extends vertically from yi to yi+1, and can be divided into two regions horizontally: the (generally larger) portion from 0 to xi+1 which is entirely contained within the desired distribution, and the (small) portion from xi+1 to xi, which is only partially contained.

Ignoring for a moment the problem of layer 0, and given uniform random variables U0 and U1 ∈ [0,1), the ziggurat algorithm can be described as:

  1. Choose a random layer 0 ≤ i < n.
  2. Let x = U0xi.
  3. If x<xi+1, return x.
  4. Let y = yi + U1(yi+1yi).
  5. Compute f(x). If y<f(x), return x.
  6. Otherwise, choose new random numbers and go back to step 1.

Step 1 amounts to choosing a low-resolution y coordinate. Step 3 tests if the x coordinate is clearly within the desired density function without knowing more about the y coordinate. If it is not, step 4 chooses a high-resolution y coordinate, and step 5 does the rejection test.

With closely spaced layers, the algorithm terminates at step 3 a very large fraction of the time. For the top layer n − 1, however, this test always fails, because xn = 0.

Layer 0 can also be divided into a central region and an edge, but the edge is an infinite tail. To use the same algorithm to check if the point is in the central region, generate a fictitious x0 = A/y1. This will generate points with x < x1 with the correct frequency, and in the rare case that layer 0 is selected and xx1, use a special fallback algorithm to select a point at random from the tail. Because the fallback algorithm is used less than one time in a thousand, speed is not essential.

Thus, the full ziggurat algorithm for one-sided distributions is:

  1. Choose a random layer 0 ≤ i < n.
  2. Let x = U0xi
  3. If x<xi+1, return x.
  4. If i = 0, generate a point from the tail using the fallback algorithm.
  5. Let y = yi + U1(yi+1yi).
  6. Compute f(x). If y<f(x), return x.
  7. Otherwise, choose new random numbers and go back to step 1.

For a two-sided distribution, the result must be negated 50% of the time. This can often be done conveniently by choosing U0 ∈ (−1,1) and, in step 3, testing if |x| <xi+1.

Fallback algorithms for the tail

Because the ziggurat algorithm only generates most outputs very rapidly, and requires a fallback algorithm whenever x > x1, it is always more complex than a more direct implementation. The specific fallback algorithm depends on the distribution.

For an exponential distribution, the tail looks just like the body of the distribution. One way is to fall back to the most elementary algorithm E = −ln(U1) and let x = x1  ln(U1). Another is to call the ziggurat algorithm recursively and add x1 to the result.

For a normal distribution, Marsaglia suggests a compact algorithm:

  1. Let x = −ln(U1)/x1.
  2. Let y = −ln(U2).
  3. If 2y > x2, return x + x1.
  4. Otherwise, go back to step 1.

Since x1 ≈ 3.5 for typical table sizes, the test in step 3 is almost always successful. Since −ln(U1) is an exponentially distributed variate, an implementation of the exponential distribution may be used.

Optimizations

The algorithm can be performed efficiently with precomputed tables of xi and yi = f(xi), but there are some modifications to make it even faster:

Generating the tables

It is possible to store the entire table precomputed, or just include the values n, y1, A, and an implementation of f−1(y) in the source code, and compute the remaining values when initializing the random number generator.

As previously described, you can find xi = f−1(yi) and yi+1 = yi + A/xi. Repeat n  1 times for the layers of the ziggurat. At the end, you should have yn = f(0). There will be some round-off error, but it is a useful sanity test to see that it is acceptably small.

When actually filling in the table values, just assume that xn = 0 and yn = f(0), and accept the slight difference in layer n  1's area as rounding error.

Finding x1 and A

Given an initial (guess at) x1, you need a way to compute the area t of the tail for which x > x1. For the exponential distribution, this is just ex1, while for the normal distribution, assuming you are using the unnormalized f(x) = ex2/2, this is π/2 erfc(x/2). For more awkward distributions, numerical integration may be required.

With this in hand, from x1, you can find y1 = f(x1), the area t in the tail, and the area of the base layer A = x1y1 + t.

Then compute the series yi and xi as above. If yi > f(0) for any i<n, then the initial estimate x1 was too low, leading to too large an area A. If yn<f(0), then the initial estimate x1 was too high.

Given this, use a root-finding algorithm (such as the bisection method) to find the value x1 which produces yn−1 as close to f(0) as possible. Alternatively, look for the value which makes the area of the topmost layer, xn−1(f(0)  yn−1), as close to the desired value A as possible. This saves one evaluation of f−1(x) and is actually the condition of greatest interest.

McFarland's variation

Christopher D. McFarland has proposed a further-optimized version. [1] This applies three algorithmic changes, at the expense of slightly larger tables.

First, the common case considers only the rectangular portions, from (0, yi−1) to (xi, yi) The odd-shaped regions to the right of these (mostly almost triangular, plus the tail) are handled separately. This simplifes and speeds up the algorithm's fast path.

Second, the exact area of the odd-shaped regions is used; they are not rounded up to include the entire rectangle to (xi−1, yi). This increases the probability that the fast path will be used.

One major consequence of this is that the number of layers is slightly less than n. Even though the area of the odd-shaped portions is taken exactly, the total adds up to more than one layer's worth. The area per layer is adjusted so that the number of rectangular layers is an integer. If the initial 0  i < n exceeds the number of rectangular layers, phase 2 proceeds.

If the value sought lies in any of the odd-shaped regions, the alias method is used to choose one, based on its true area. This is a small amount of additional work, and requires additional alias tables, but chooses one of the layers' right-hand sides.

The chosen odd-shaped region is rejection sampled, but if a sample is rejected, the algorithm does not return to the beginning. The true area of each odd-shaped region was used to choose a layer, so the rejection-sampling loop stays in that layer until a point is chosen.

Third, the almost-triangular shape of most odd-shaped portions is exploited, although this must be divided into three cases depending on the second derivative of the probability distribution function in the selected layer.

If the function is convex (as the exponential distribution is everywhere, and the normal distribution is for |x| > 1), then the function is strictly contained within the lower triangle. Two unit uniform deviates U1 and U2 are chosen, and before they are scaled to the rectangle enclosing the odd-shaped region, their sum is tested. If U1 + U2 > 1, the point is in the upper triangle and can be reflected to (1−U1, 1−U2). Then, if U1 + U2 < 1−ε, for some suitable tolerance ε, the point is definitely below the curve and can immediately be accepted. Only for points very close to the diagonal is it necessary to compute the distribution function f(x) to perform an exact rejection test. (The tolerance ε should in theory depend on the layer, but a single maximum value can be used on all layers with little loss.)

If the function is concave (as the normal distribution is for |x| < 1), it includes a small portion of the upper triangle so reflection is impossible, but points whose normalized coordinates satisfy U1 + U2  1 can be immediately accepted, and points for which U1 + U2 > 1+ε can be immediately rejected.

In the one layer which straddles |x| = 1, the normal distribution has an inflection point, and the exact rejection test must be applied if 1−ε <U1 + U2 < 1+ε.

The tail is handled as in the original Ziggurat algorithm, and can be thought of as a fourth case for the shape of the odd-shaped region to the right.

Related Research Articles

A pseudorandom number generator (PRNG), also known as a deterministic random bit generator (DRBG), is an algorithm for generating a sequence of numbers whose properties approximate the properties of sequences of random numbers. The PRNG-generated sequence is not truly random, because it is completely determined by an initial value, called the PRNG's seed. Although sequences that are closer to truly random can be generated using hardware random number generators, pseudorandom number generators are important in practice for their speed in number generation and their reproducibility.

<span class="mw-page-title-main">Box–Muller transform</span> Statistical transform

The Box–Muller transform, by George Edward Pelham Box and Mervin Edgar Muller, is a random number sampling method for generating pairs of independent, standard, normally distributed random numbers, given a source of uniformly distributed random numbers. The method was first mentioned explicitly by Raymond E. A. C. Paley and Norbert Wiener in their 1934 treatise on Fourier transforms in the complex domain. Given the status of these latter authors and the widespread availability and use of their treatise, it is almost certain that Box and Muller were well aware of its contents.

Bresenham's line algorithm is a line drawing algorithm that determines the points of an n-dimensional raster that should be selected in order to form a close approximation to a straight line between two points. It is commonly used to draw line primitives in a bitmap image, as it uses only integer addition, subtraction, and bit shifting, all of which are very cheap operations in historically common computer architectures. It is an incremental error algorithm, and one of the earliest algorithms developed in the field of computer graphics. An extension to the original algorithm called the midpoint circle algorithm may be used for drawing circles.

In mathematics, a Diophantine equation is an equation of the form P(x1, ..., xj, y1, ..., yk) = 0 (usually abbreviated P(x, y) = 0) where P(x, y) is a polynomial with integer coefficients, where x1, ..., xj indicate parameters and y1, ..., yk indicate unknowns.

<span class="mw-page-title-main">Gamma distribution</span> Probability distribution

In probability theory and statistics, the gamma distribution is a versatile two-parameter family of continuous probability distributions. The exponential distribution, Erlang distribution, and chi-squared distribution are special cases of the gamma distribution. There are two equivalent parameterizations in common use:

  1. With a shape parameter k and a scale parameter θ
  2. With a shape parameter and an inverse scale parameter , called a rate parameter.
<span class="mw-page-title-main">Perlin noise</span> Type of gradient noise in computer graphics

Perlin noise is a type of gradient noise developed by Ken Perlin in 1983. It has many uses, including but not limited to: procedurally generating terrain, applying pseudo-random changes to a variable, and assisting in the creation of image textures. It is most commonly implemented in two, three, or four dimensions, but can be defined for any number of dimensions.

In mathematics, a low-discrepancy sequence is a sequence with the property that for all values of N, its subsequence x1, ..., xN has a low discrepancy.

In probability theory, a Chernoff bound is an exponentially decreasing upper bound on the tail of a random variable based on its moment generating function. The minimum of all such exponential bounds forms the Chernoff or Chernoff-Cramér bound, which may decay faster than exponential. It is especially useful for sums of independent random variables, such as sums of Bernoulli random variables.

A random permutation is a random ordering of a set of objects, that is, a permutation-valued random variable. The use of random permutations is often fundamental to fields that use randomized algorithms such as coding theory, cryptography, and simulation. A good example of a random permutation is the shuffling of a deck of cards: this is ideally a random permutation of the 52 cards.

The diehard tests are a battery of statistical tests for measuring the quality of a random number generator. They were developed by George Marsaglia over several years and first published in 1995 on a CD-ROM of random numbers. In 2006, the original diehard tests were extended into the dieharder tests.

George Marsaglia was an American mathematician and computer scientist. He is best known for creating the diehard tests, a suite of software for measuring statistical randomness.

<span class="mw-page-title-main">Random number generation</span> Producing a sequence that cannot be predicted better than by random chance

Random number generation is a process by which, often by means of a random number generator (RNG), a sequence of numbers or symbols that cannot be reasonably predicted better than by random chance is generated. This means that the particular outcome sequence will contain some patterns detectable in hindsight but impossible to foresee. True random number generators can be hardware random-number generators (HRNGs), wherein each generation is a function of the current value of a physical environment's attribute that is constantly changing in a manner that is practically impossible to model. This would be in contrast to so-called "random number generations" done by pseudorandom number generators (PRNGs), which generate numbers that only look random but are in fact predetermined—these generations can be reproduced simply by knowing the state of the PRNG.

Random self-reducibility (RSR) is the rule that a good algorithm for the average case implies a good algorithm for the worst case. RSR is the ability to solve all instances of a problem by solving a large fraction of the instances.

In computer science, multiply-with-carry (MWC) is a method invented by George Marsaglia for generating sequences of random integers based on an initial set from two to many thousands of randomly chosen seed values. The main advantages of the MWC method are that it invokes simple computer integer arithmetic and leads to very fast generation of sequences of random numbers with immense periods, ranging from around to .

<span class="mw-page-title-main">Truncated normal distribution</span> Type of probability distribution

In probability and statistics, the truncated normal distribution is the probability distribution derived from that of a normally distributed random variable by bounding the random variable from either below or above. The truncated normal distribution has wide applications in statistics and econometrics.

In computer graphics, a digital differential analyzer (DDA) is hardware or software used for interpolation of variables over an interval between start and end point. DDAs are used for rasterization of lines, triangles and polygons. They can be extended to non linear functions, such as perspective correct texture mapping, quadratic curves, and traversing voxels.

<span class="mw-page-title-main">Relationships among probability distributions</span> Topic in probability theory and statistics

In probability theory and statistics, there are several relationships among probability distributions. These relations can be categorized in the following groups:

In sequence design, a Feedback with Carry Shift Register is the arithmetic or with carry analog of a linear-feedback shift register (LFSR). If is an integer, then an N-ary FCSR of length is a finite state device with a state consisting of a vector of elements in and an integer . The state change operation is determined by a set of coefficients and is defined as follows: compute . Express s as with in . Then the new state is . By iterating the state change an FCSR generates an infinite, eventually periodic sequence of numbers in .

Non-uniform random variate generation or pseudo-random number sampling is the numerical practice of generating pseudo-random numbers (PRN) that follow a given probability distribution. Methods are typically based on the availability of a uniformly distributed PRN generator. Computational algorithms are then used to manipulate a single random variate, X, or often several such variates, into a new random variate Y such that these values have the required distribution. The first methods were developed for Monte-Carlo simulations in the Manhattan project, published by John von Neumann in the early 1950s.

<span class="mw-page-title-main">Alias method</span> Family of algorithms for sampling from discrete probability distributions

In computing, the alias method is a family of efficient algorithms for sampling from a discrete probability distribution, published in 1974 by Alastair J. Walker. That is, it returns integer values 1 ≤ in according to some arbitrary discrete probability distribution pi. The algorithms typically use O(n log n) or O(n) preprocessing time, after which random values can be drawn from the distribution in O(1) time.

References

  1. McFarland, Christopher D. (24 June 2015). "A modified ziggurat algorithm for generating exponentially and normally distributed pseudorandom numbers". Journal of Statistical Computation and Simulation. 86 (7): 1281–1294. arXiv: 1403.6870 . doi:10.1080/00949655.2015.1060234. Note that the Bitbucket repository mentioned in the paper is no longer available and the code is now at https://github.com/cd-mcfarland/fast_prng