# 8. Pseudo-Random Numbers and Monte Carlo Sampling Methods¶

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.

## 8.1. Generating Pseudo-Random Numbers¶

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.

### 8.1.1. Parallel random numbers¶

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_rngs 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_rngs`

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.

## 8.2. Low-Discrepancy Sequences¶

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:

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 : i32
val norm : f64
val independent : i32 -> [D]u32
val recurrent : i32 -> [D]u32 -> [D]u32
val sobol : (n: i32) -> [n][D]f64
}
module Sobol : (DM : sobol_dir) -> (X : { val D : i32 }) -> 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:

```
import "lib/github.com/diku-dk/sobol/sobol-dir-50"
import "lib/github.com/diku-dk/sobol/sobol"
module sobol = Sobol sobol_dir { let D = 2 }
```

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:

```
let sqr (x:f64) = x * x
let in_circle (p:[2]f64) : bool =
sqr p[0] + sqr p[1] < 1.0f64
let pi_arr [n] (arr: [n][2]f64) : f64 =
let bs = map (i32.bool <-< in_circle) arr
let sum = reduce (+) 0 bs
in 4f64 * r64 sum / f64.i32 n
let main (n:i32) : f64 =
sobol.sobol n |> pi_arr
```

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 [HEO18]. In general, there are other good reasons to avoid uniform sampling in relation to Monte-Carlo methods.