Pseudo-random number generation and Monte Carlo sampling are concepts
that apply to a large number of application areas. In a data-parallel
setting, these concepts require special treatment beyond the usual
sequential methods. In this chapter, we first present a Futhark
package, called cpprandom for generating pseudo-random numbers in
parallel. We then present a Futhark package, called sobol, for
generating Sobol sequences, which are examples of so-called
low-discrepancy sequences, sequences that make numerical
multi-dimensional integration converge faster than if pseudo-random
numbers were used.
The cpprandom package is inspired by the C++ library <random>,
which is very elaborate, but also very flexible. Due to Futhark's
purity, it is up to the programmer to explicitly manage the state of
the pseudo-random number engine (the RNG state). In particular, it is
the programmer's responsibility to ensure that the same state is not
used more than once (unless that is what is desired).
The following program constructs a uniform distribution of single precision
floats using minstd_rand as the underlying RNG engine.
module dist = uniform_real_distribution f32 minstd_rand let rng = minstd_rand.rng_from_seed [123] let (rng, x) = dist.rand (1,6) rng
The dist module is constructed at the program top level, while we
use it at the expression level. We use the minstd_rand module for
initialising the random number state using a seed, and then we pass
that state to the rand function in the generated distribution
module, along with a description of the distribution we desire. We
get back not just the random number, but also the new state of the
engine.
The dist.rand function, coming from uniform_real_distribution,
simply takes a pair of numbers describing the range. Consider instead the following code:
module norm_dist = normal_distribution f32 minstd_rand
let (rng, y) = norm_dist.rand {mean=50, stddev=25} rng
In contrast to dist.rand, the norm_dist.rand function, coming
from normal_distribution takes a record specifying the mean and
the standard deviation. Since both dist and norm_dist have
been initialised with the same underlying rng_engine, we can reuse
the same RNG state. Such reuse is often convenient when a program
needs to generate random numbers from several different distributions,
as we still only have to manage a single RNG state.
Random number generation is inherently sequential. The rand
functions take an RNG state as input and produce a new RNG state.
This dependence creates challenges when we wish to map a function
f across some array xs, and each application of the function
must produce some random numbers. We generally don't want to pass the
exact same state to every application, as that means each element will
see the exact same stream of random numbers. The common procedure is
to use split_rng, which creates any number of RNG states from one,
and then pass one to each application of f:
let rngs = minstd_rand.split_rng n rng let (rngs, ys) = unzip (map2 f rngs xs) let rng = minstd_rand.join_rng rngs
We assume here that the function f returns not just the result,
but also the new RNG state. Generally, all functions that accept
random number states should behave like this. We subsequently use
join_rng to combine all resulting states back into a single state.
Thus, parallel programming with random numbers involves frequently
splitting and rejoining RNG states. For most RNG engines, these
operations are generally very cheap.
The Futhark package sobol is a package for generating Sobol
sequences, which are examples of so-called low-discrepancy
sequences, sequences that, when combined with Monte-Carlo methods,
make numeric integration converge faster than if ordinary
pseudo-random numbers are used and are more flexible than if uniform
sampling techniques are used. Sobol sequences may be multi-dimensional
and a key property of using Sobol sequences is that we can freely
choose the number of points that should span the multi-dimensional
space. In contrast, if we set out to use a simpler uniform sampling
technique for spanning two dimensions, we can only span the space
properly if we choose the number of points to be on the form
x^2, for some natural number x. This spanning problem
becomes worse for higher dimensions.
As an example, we shall see how we can use Sobol sequences together with Monte-Carlo simulation to compute the value of \pi. We shall also see that doing so will result in faster convergence towards the true value of \pi compared to if pseudo-random numbers are used.
To calculate an approximation to the value of \pi, we will use a simple dart-throwing approach. We will throw darts at a 2 by 2 square, centered around the origin, and then establish the ratio between the number of darts hitting within the unit circle with the number of darts hitting the square. This ratio multiplied with 4 will be our approximation of \pi. The more darts we throw, the better our approximation, assuming that the darts we throw hit the board somewhat evenly. To calculate whether a particular dart, thrown at the point (x,y), is within the unit circle, we can apply the standard Pythagoras formula:
\pi ~~\approx~~ \frac{4}{N} \sum_{i=1}^N \left \{ \begin{array}{ll} 1 & \mbox{if} ~ x_i^2 + y_i^2 < 1 \\ 0 & \mbox{otherwise} \end{array} \right .
For the actual throwing of darts, we need to establish N pairs of numbers, each in the interval [-1;1]. Now, it turns out that it matters significantly how we choose to throw the darts. Some obvious choices would be to throw the darts in a regular grid (i.e., uniform sampling), or to choose points using a pseudo-random number generator.
The Futhark package makes essential use of an independent formula for calculating, independently, the n'th Sobol number. However, even though such a formula is essential for achieving parallelism, it performs poorly compared to the more efficient recurrent formula, which makes it possible to calculate the n'th Sobol number if we know the previous Sobol number. The Futhark package makes essential use of both formulas. The calculation of a sequence of Sobol numbers depends on a set of direction vectors, which are also provided by the package.
The key functionality of the package comes in the form of a higher-order module Sobol, which takes as arguments a direction vector module and a module specifying the dimensionality of the generated Sobol numbers:
module type sobol_dir = { ... }
module sobol_dir : sobol_dir -- file sobol-dir-50, e.g.
module type sobol = {
val D : i64
val norm : f64
val independent : i32 -> [D]u32
val recurrent : i32 -> [D]u32 -> [D]u32
val sobol : (n: i64) -> [n][D]f64
}
module Sobol : (DM : sobol_dir) -> (X : { val D : i64 }) -> sobol
For estimating the value of \pi, we will need a two-dimensional Sobol sequence, thus we apply the Sobol higher-order module to the direction vector module that works for up-to 50 dimensions and a module specifying a dimensionality of two:
.. literalinclude:: src/pi.fut :lines: 1-4
We can now complete the program by writing a main function that computes an array of Sobol numbers of a size given by the parameter given to main and feed this array into a function that will compute the estimation of \pi using the function shown above:
.. literalinclude:: src/pi.fut :lines: 6-17
The use of Sobol numbers for estimating \pi turns out to be about three times slower than using a uniform grid on a standard GPU. However, it converges towards \pi equally well (with increasing N) and is superior for larger dimensions :cite:`futhark:fhpc18`. In general, there are other good reasons to avoid uniform sampling in relation to Monte-Carlo methods.