Randomness in NEST Simulations

osx-arm64: missing random number generators

Due to a cross-compiling issue in the conda NEST package, some random number generators are not available if you are using macOS arm64 architecture. The available generators are the Mersenne Twister generators mt19937 and mt19937_64.

Random numbers in network simulations

Random numbers are used for a variety of purposes in neuronal network simulations, e.g.,

  1. during node creation to randomize node parameters;

  2. during connection creation

    1. to select connections to be created,

    2. to parameterize connections;

  3. during simulation

    1. to generate stochastic input (randomized spike trains, noisy currents),

    2. to support stochastic spike generation.

This document discusses how NEST provides random numbers for these purposes, how you can choose which random number generator (RNG) to use, and how to set the seed of RNGs in NEST. We use the term “random number” here for ease of writing, even though we are always talking about pseudorandom numbers generated by some algorithm.

Compiler dependency

NEST uses random generators provided by the C++ Standard Library to obtain random deviates, i.e., random numbers following different distributions such as normal, binomial, or Poissonian. Different versions of the C++ Standard Library use different algorithms to generate these deviates. Therefore, simulation results can differ in detail depending on the compiler and the C++ Standard Library implementation used. Note in particular that GCC and Clang come with different library implementations.

How to work with randomness in NEST

NEST makes randomness easy: Setting a single seed gives you full control over randomness throughout a simulation, and you can choose seed values at liberty. NEST’s powerful mechanism to parameterize nodes and synapses and NEST’s probabilistic connection rules allow you to construct your network with parameters with a wide range of random properties, from uniformly distributed numbers drawn from a fixed interval to mathematical expressions combining spatial information with probability distributions. Neuron, device, and synapse models apply randomness during simulation as described in the model documentation.

This section provides a practical introduction into managing random number generators and seeds in NEST, followed by an overview of the nest.random module and examples of using randomness. A separate section provides suggestions for additional randomization from the Python level.

Programs that are entirely serial typically use a single source of random numbers. In parallel simulations, such a single source of random numbers would form an unacceptable bottleneck. To avoid this, NEST uses parallel streams of random numbers. See section Random number internals if you are interested in the technical details.

Select the type of random number generator

NEST provides different random number generators, presently taken from the C++ Standard Library and the Random123 library. To see all available generators, run

nest.rng_types

This currently gives the following generators:

'Philox_32', 'Philox_64',
'Threefry_32', 'Threefry_64',
'mt19937', 'mt19937_64'

The Philox and Threefry generators are cryptographic generators from Random123, while mt19937 is the Mersenne Twister generator. Each generator comes in a 32- and a 64-bit implementation. The number of bits has no consequence for the quality of randomness, but depending on computer and compiler, one or the other may show better performance. More generators may be added in the future. Currently, only a small number of generators is available because we require generators with a sufficiently long period as discussed in Random number internals.

The default random number generator set in NEST is mt19937_64. To choose a different generator, set it in the following way:

nest.rng_type = 'Philox_32'

It is a good idea to cross-check your simulation results using a different random number generator type. Even though generators and our understanding of them has become much better in recent years, there always remains a risk of RNG artifacts affecting simulations.

Seed the random number generator

NEST uses a built-in default seed if you do not specify one. This means that unless you explicitly set the seed, all simulations will be run with the same sequence of random numbers. So do not forget to seed your simulations!

You can use any number \(s\) with \(1\leq s \leq 2^{31}-1\) as seed:

nest.rng_seed = 12345

As long as you use different seed values, NEST will ensure that all random number streams in a simulation are seeded properly; see Random number internals for details.

You can inspect the RNG type and seed value used with

print(nest.rng_type, nest.rng_seed)

Any simulation run with the same seed shall return identical results (provided the same compiler/C++ Standard Library was used).

The NEST random module

The nest.random module provides a range of random distributions that can be used to specify parameters for neurons, synapses, and connection rules. See below for examples on how to use them in practice and some details on randomizing delays.

nest.random.hl_api_random.exponential(beta=1.0)

Draws samples from an exponential distribution.

Parameters:

beta (float, optional) – Scale parameter the distribution. Default value is 1.0.

Returns:

Object yielding values drawn from the distribution.

Return type:

Parameter

nest.random.hl_api_random.lognormal(mean=0.0, std=1.0)

Draws samples from a log-normal distribution.

Parameters:
  • mean (float, optional) – Mean value of the underlying normal distribution. Default value is 0.

  • std (float, optional) – Standard deviation of the underlying normal distribution. Default value is 1.0.

Returns:

Object yielding values drawn from the distribution.

Return type:

Parameter

nest.random.hl_api_random.normal(mean=0.0, std=1.0)

Draws samples from a normal distribution.

Parameters:
  • mean (float, optional) – Mean of the distribution. Default value is 0.

  • std (float, optional) – Standard deviation of the distribution. Default value is 1.0.

Returns:

Object yielding values drawn from the distribution.

Return type:

Parameter

nest.random.hl_api_random.uniform(min=0.0, max=1.0)

Draws samples from a uniform distribution.

Samples are distributed uniformly in [min, max) (includes min, but excludes max).

Note

See this documentation for details on the effect of time discretization on delays drawn from a uniform distribution.

Parameters:
  • min (float, optional) – Lower boundary of the sample interval. Default value is 0.

  • max (float, optional) – Upper boundary of the sample interval. Default value is 1.0.

Returns:

Object yielding values drawn from the distribution.

Return type:

Parameter

nest.random.hl_api_random.uniform_int(max)

Draws integer samples from a uniform distribution.

Samples are distributed uniformly in [0, max) (includes 0, but excludes max).

Parameters:

max (integer) – Upper boundary of the sample interval.

Returns:

Object yielding values drawn from the distribution.

Return type:

Parameter

Examples of using randomness

Randomize the membrane potential

To set the membrane potential at creation, just pass a random distribution as the V_m value.

nest.Create('iaf_psc_alpha', 10000, {'V_m': nest.random.normal(mean=-60.0, std=10.0)})

You may also set the membrane potential after creation.

n = nest.Create('iaf_psc_alpha', 10000)
n.V_m = nest.random.normal(mean=-60.0, std=10.0)

Randomize other node parameters

Other parameters may be randomized in the same way as the membrane potential.

nest.Create('iaf_psc_alpha', 10000,
            {'C_m': nest.random.uniform(min=240.0, max=260.0),
             'I_e': nest.random.uniform(min=0.0, max=5.0)})

Randomize connection parameters

Likewise, synapse parameters can be specified using the random distributions.

nest.Connect(n, n, syn_spec={'weight': nest.random.normal(mean=0., std=1.),
                             'delay': nest.random.uniform(min=0.5, max=1.5)})

Rounding effects when randomizing delays

Connection delays in NEST are rounded to the nearest multiple of the simulation resolution, even if delays are drawn from a continuous distribution. This will work as expected if the distribution has infinite support, for example, the normal distribution. For the uniform distribution, though, this rounding will usually lead to lower probabilities for the delays at the edges of the distribution. This then also applies to a truncated normal distribution. Consider the following case:

nest.resolution = 0.1
n = nest.Create('iaf_psc_alpha', 100)
nest.Connect(n, n, syn_spec={'delay': nest.random.uniform(min=1, max=2)})

This will create 10000 connections in total with delay values 1.0, 1.1, …, 2.0. But while the interior delay values 1.1, …, 1.9 will occur approximately 1000 times each, the first and last cases, 1.0 and 2.0, will occur only approximately 500 times each. This happens because NEST first draws the delay uniformly from \([1, 2)\) and then rounds to a fixed delay. Thus, any number from \([1.05, 1.15)\) will be rounded to 1.1, but only numbers in \([1.0, 1.05)\) will be rounded to 1.0.

To achieve equal probabilities of the first and last values, you can either extend the interval of the uniform distribution by half the resolution:

nest.Connect(n, n, syn_spec={'delay': nest.random.uniform(min=1 - 0.5 * nest.resolution,
                                                          max=2 + 0.5 * nest.resolution)})

or use the uniform integer distribution. Since it always draws numbers beginning with zero, this is slightly more cumbersome:

nest.Connect(n, n, syn_spec={'delay': 1 + 0.1 * nest.random.uniform_int(11)})

An in-depth analysis of delay rounding is available in a master thesis.

Randomize spatial positions

Spatial positions of neurons may be randomized using the same random distributions. Note that we have to specify the number of dimensions we want when specifying a single distribution.

nest.Create('iaf_psc_alpha', 1000,
            positions=nest.spatial.free(nest.random.uniform(min=0.0, max=10.0),
                                        num_dimensions=2))

We can also specify a separate distribution for each dimension.

nest.Create('iaf_psc_alpha', 1000,
            positions=nest.spatial.free([nest.random.uniform(min=0.0, max=10.0),
                                         nest.random.normal(mean=5.0, std=10.0)])

Combine random distributions

Random distributions are NEST parameter objects which support arithmetic with constants.

vm = -50.0 + nest.random.exponential(beta=2.0)
nest.Create('iaf_psc_alpha', 100, {'V_m': vm})

They also support arithmetic with other NEST parameters.

n = nest.Create('iaf_psc_alpha', 1000,
                positions=nest.spatial.free(nest.random.uniform(min=0.0, max=10.0),
                                            num_dimensions=2))

nest.Connect(n, n,
             conn_spec={'rule': 'fixed_indegree', 'indegree': 1},
             syn_spec={'delay': nest.spatial.distance + nest.random.normal(mean=1.0, std=0.5)})

Randomization from the Python level

In rare occasions, the capabilities provided by NEST for parameterizing your network will not cover your needs and require additional randomization at the Python level. This requires some additional consideration when performing simulations in parallel, especially when using MPI parallelization. You may want to consult the section on Random number internals before continuing.

A key principle of parallel simulation in NEST is that a simulation performed with a fixed number of virtual processes \(N_{\text{vp}} = M \times T\) shall produce identical results independent of the number of MPI processes \(M\) and threads \(T\) that go into each virtual process. To observe this principle when also randomizing from the Python level, it is essential to create one Python random number generator per virtual process and use the random number generator for the virtual process to which a node belongs (for synapses: the VP of the target node).

We consider first an example setting a random membrane potential. We use the modern random package introduced with NumPy 1.17.

Warning

Don’t do this in Python!

The randomization example below is shown only to demonstrate how you in principle could use random values generated at the Python level in a NEST simulation in a way consistent with NEST’s approach to random numbers in parallel simulations. You should only need to use this approach under very particular circumstances (e.g., if the distribution you want to use is not available in nest.random).

The proper way to randomize the membrane potential in NEST with uniform random numbers as in the example below is

nest.Create('iaf_psc_alpha', 12, params={'V_m': nest.random.uniform()})
import numpy as np

n_vp = 4
py_seed = 987654

nest.total_num_virtual_procs = n_vp
nrns = nest.Create('iaf_psc_alpha', 12)

rngs = {vp: np.random.default_rng(py_seed + vp) for vp in nest.GetLocalVPs()}

for n in nest.GetLocalNodeCollection(nrns):
    n.set({'V_m': rngs[n.get('vp')].uniform()})

After creating the neurons, we create n_vp random number generators. We then use nest.GetLocalNodeCollection() to obtain all those neurons that belong to the local MPI rank, since each rank can only access local nodes. We then use n.get('vp') to obtain the virtual processes responsible for the node and use it to pick out the correct RNG, from which we draw the random membrane potential.

To parameterize a connection, we need to use the random generator for the virtual process of the target neuron. Continuing from the example above, we can randomize the connection weight as follows:

nest.Connect(nrns, nrns, 'one_to_one')

for c in nest.GetConnections():
    c.set({'weight': rngs[c.get('vp')].uniform()})

Random number internals

Global and local generators

As described above, random numbers are consumed in NEST

  1. during node creation to randomize node parameters;

  2. during connection creation

    1. to select connections to be created,

    2. to parameterize connections;

  3. during simulation

    1. to generate stochastic input (randomized spike trains, noisy currents),

    2. to support stochastic spike generation.

In all cases except certain sub-cases of 2.a, randomization is tied to a node, which in turn is assigned to a specific virtual process (VP), either because the node is concerned (cases 1, 3) or because the node is the target of a connection (case 2). NEST guarantees that nodes assigned to different VPs can be updated independently of each other, and the same is true for connections if created or parameterized from a target-node perspective. Therefore, for all these purposes, one random number stream per VP is required and sufficient to ensure that different VPs can operate independently.

The only exception are certain sub-cases of 2.a and 3 which require identical random number streams either per MPI process or per virtual process. This pertains in particular to fixed_outdegree and fixed_total_number connection rules and the mip_generator.

NEST therefore provides three kinds of random number streams

  1. a rank-synchronized (global) random number stream to be used only by the main thread of each MPI rank; all MPI ranks must draw equally many numbers from this stream and discard those they do not need.

  2. a VP-synchronized (global) random number stream to be used in parallel by all VPs.

  3. one asynchronous random number stream per virtual process.

This results in a total of \(N_{\text{vp}}+2\) random number streams. To avoid unnecessary complications in the code using random numbers, serial simulations also use all three kinds. The generators for all streams are of the same type. If the RNG type is changed, the change applies to all generators.

NEST regularly checks during a simulation that the rank- and VP-synchronized random number streams indeed stay in sync.

Quantity of random numbers required

Considering the future with exascale computers, we need to prepare for \(N_{\text{vp}}=\mathcal{O}(10^7)\) random streams. We can estimate the number of random numbers required per stream as follows:

  • Based on current practice, we assume \(N = 1000\) neurons per stream; this results in a total of \(\mathcal{O}(10^{11})\) neurons, on the scale of the human brain and thus a sensible upper limit.

  • Each neuron receives excitatory and inhibitory Poisson background input, i.e., two Poisson distributed random numbers per time step per neuron, each of which typically requires multiple uniform random numbers; we assume for simplicity two uniform numbers per Poisson number, i.e., four plain random numbers per time step per neuron.

  • A time step is 0.1ms and a typical simulation duration 100s simulated time resulting in \(\mathcal{O}(10^6)\) time steps.

  • In total, a simulation will consume \(\sim 1000 \times 4 \times 10^6 = 4 \times 10^{10}\) random numbers per stream (consumption during network construction will usually be negligible).

  • To allow for a margin of error, we assume that \(\mathcal{O}(10^{11})\) random numbers are consumed per stream during a single simulation.

Collision risk of parallel random number streams

L’Ecuyer et al (2017) provide a recent account of random number generation in highly parallel settings. The main approach to providing independent random number streams in parallel simulations is to use a high-quality random number generator with long period and to seed it with a different seed for each stream. They then argue that if we have \(s\) streams of length \(l\) and a generator with period \(r\), the probability that two streams seeded with different seeds will overlap is given by

\[p_{\text{overlap}} \sim s^2\frac{l}{r}\; .\]

We have

\[s = N_{\text{vp}} = 10^7 \sim 2^{23}\quad\text{and}\quad l = 10^{11} \sim 2^{37}\;,\]

so assuming an RNG period of \(r = 2^{128}\) we obtain

\[p_{\text{overlap}} \sim \left(2^{23}\right)^2 \times 2^{37} / 2^{128} = 2^{-45} \sim 3 \times 10^{-14}\]

while for a period of \(r = 2^{256}\) we obtain

\[p_{\text{overlap}} \sim 2^{-173} \sim 10^{-52}\; .\]

The probability of stream collisions is thus negligibly small, provided we use random number generators with a period of at least \(2^{128}\).

L’Ecuyer (2012, Ch 3.6) points out that certain random number generator classes will show too much regularity in their output if more than \(r^{1/3}\) numbers are used (this relation is based on exercises in Knuth, TAOCP vol 2, see L’Ecuyer and Simard (2001)). While it is not certain that that analysis also applies to (short) subsequences of the period of an RNG, one might still re-consider the analysis above, by using \(l^3 \sim 2^{111}\) instead of \(l\) when computing \(p_{\text{overlap}}\). We then obtain

Period \(r\)

\(p_{\text{overlap}}\)

\(2^{128} \sim 10^{38}\)

\(>> 1\)

\(2^{256} \sim 10^{77}\)

\(2^{-99} \sim 10^{-30}\)

\(2^{512} \sim 10^{154}\)

\(2^{-355} \sim 10^{-107}\)

Thus, for generators with periods of at least \(2^{256}\), we have negligible collision probability also under this tightened criterium.

Among the generators offered by the C++ Standard Library, only the Mersenne Twister generator has a sufficiently long period, so only this generator is included. The generators from the Random123 library are cryptographic generators and thus no period.

Seed parallel random number streams

In NEST, the user provides a single seed \(1 \leq S \leq 2^{31}-1\) to initialize all \(N_{\text{vp}}+2\) random number generators in a parallel simulation. NEST uses a seed sequence generator to derive seeds for the individual RNGs from the user-provided seed \(S\).

Since std::seed_seq from the C++ Standard Library has weaknesses, we use Melissa O’Neill’s improved seed sequence generator.

Further background

For more background on the NEST random number generation architecture, see also NEST Github issue #1440.