Parallel Random Number Generation¶
There are five strategies implemented that can be used to produce repeatable pseudo-random numbers across multiple processes (local or distributed). The final method is suitable for parallel applications where reproducibility is not a concern.
Using a SeedSequence¶
Using a SeedSequence is a universal way to initialize many streams.
When using with a bit generator with a moderately large state (say 256-bits), the probability
of a collision is negligible. This said, some generators can have difficult to detect
correlation when producing many streams in corner cases.
import numpy as np
from randomgen import SFC64
NUM_STREAMS = 8196
seed_seq = SeedSequence(603484028490308141)
children = seed_seq.spawn(NUM_STREAMS)
streams = [SFC64(child) for child in children]
Distinct Keys in Cryptographic Generators¶
The cryptographic pseudo-random number generators (PRNGs) support using distinct
keys to produce distinct sequence.
Generators that support this form of parallelization include AESCounter,
ChaCha, HC128,
ThreeFry, Philox, and
SPECK128.
This example shows how many streams can be created by passing in different index values in the second input while using the same seed in the first.
import numpy as np
from randomgen import AESCounter, random_entropy
NUM_STREAMS = 8196
keys = random_entropy(4 * NUM_STREAMS).view(np.uint64)
keys = set([tuple(key.tolist()) for key in keys.reshape((-1, 2))])
# Essentially 0 probability this is needed
while len(keys) < NUM_STREAMS:
new_keys = random_entropy(4 * (NUM_STREAMS - len(keys))).view(np.uint8)
new_keys = set([tuple(key.tolist()) for key in new_keys.reshape((-1, 2))])
keys.update(new_keys)
# Distinct 128-bit numbers as a key encoded as 2 uint64 values
streams = [AESCounter(key=np.array(key,dtype=np.uint64)) for key in keys]
Advancing the PRNG’s state¶
Most of the cryptographic PRNGs are counter-based, and so support advancing which increments the counter. Advancing a PRNG updates the underlying PRNG state as if a given number of calls to the underlying PRNG have been made. In general there is not a one-to-one relationship between the number output random values from a particular distribution and the number of draws from the core PRNG. This occurs for two reasons:
The random values are simulated using a rejection-based method and so, on average, more than one value from the underlying PRNG is required to generate a single draw.
The number of bits required to generate a simulated value differs from the number of bits generated by the underlying PRNG. For example, two 16-bit integer values can be simulated from a single draw of a 32-bit PRNG.
Advancing the PRNG state resets any pre-computed random numbers. This is required to ensure exact reproducibility.
import numpy as np
from randomgen import SPECK128, SeedSequence
PHI = (np.sqrt(5) - 1) / 2
STEP = int(PHI * 2**96)
NUM_STREAMS = 8196
seed_seq = SeedSequence(603484028490308141)
base = SPECK128(seed_seq)
streams = [base]
for i in range(1, NUM_STREAMS):
next_gen = SPECK128(seed_seq)
streams.append(next_gen.advance(i * STEP))
In addition to the cryptographic PRNGs,
the PCG-based generators also support advance: PCG64,
PCG64DXSM, LCG128Mix, and
PCG32.
Note that HC128 is based on a stream cipher and so
does not support advancing a counter.
Warning
Do not split a single stream of a PCG-based generator by advancing it by a
large power of two, such as \(2^{127}\) (\(2^{63}\) for
PCG32). These generators are built on an LCG
modulo a power of two, whose lowest \(e\) bits repeat every
\(2^e\) steps, so the streams share the lowest bits of their states and
their outputs are strongly related. For example, after advancing
PCG64 by \(2^{127}\), each output is the
other with its 32-bit halves swapped and one bit flipped. The advance
method of each generator describes the effect on its output functions.
Instead, advance stream i by i * step with an odd step. With
any odd step the states of two streams differ as if at random, apart from
their lowest few bits. The size of the step matters for overlap: it must be
at least the number of draws that each stream makes, and n streams are
only guaranteed not to overlap when n * step is no larger than the
period. A convenient step is the odd integer closest to
\(2^{XX}/\phi\), where \(\phi\) is the golden ratio, for example
0x9e3779b97f4a7c15 (\(XX = 64\), about \(1.1 \times 10^{19}\)
draws per stream) or, for PCG32,
0x9e3779b9 (\(XX = 32\), about \(2.7 \times 10^{9}\) draws per
stream). Each generator’s advance_step method, such as
PCG64.advance_step, takes the
largest number of draws that any stream makes and returns this step for the
smallest \(XX\) that works. jumped uses a golden-ratio step and is
also a good choice.
The period of PCG32 is only \(2^{64}\), which
limits it to, for example, about \(1.4 \times 10^{7}\) streams of
\(10^{12}\) draws. Use PCG64, or the related
PCG64DXSM or LCG128Mix,
whose period is \(2^{128}\), when many streams or many draws per stream
are needed.
from randomgen import PCG64
# Each stream makes at most 10**12 draws
STEP = PCG64.advance_step(10**12)
streams = [PCG64(1234).advance(i * STEP) for i in range(16)]
Jumping the PRNG state¶
jumped advances the state of the PRNG as if a large number of random
numbers have been drawn, and returns a new instance with this state. Jumping
is more universal than advance since the multiplier needed to jump some PRNGs
is expensive to compute. However, this multiplier can be pre-computed for a fixed
step size when the PRNG does not support an arbitrary advance which enables the jump.
The specific number of as if draws varies by PRNG, and ranges from around \(2^{64}\) to
\(2^{512}\). Additionally, the as if draws also depend on the size of
the default random number produced by the specific PRNG. The PRNGs that
support jumped, along with the period of the PRNG, the size of the jump
and the bits in the default unsigned random are listed below.
PRNG |
Period |
Jump Size |
Bits |
|---|---|---|---|
AESCounter |
\(2^{128}\) |
\(2^{64}\) |
64 |
ChaCha |
\(2^{128}\) |
\(2^{64}\) |
64 |
LCG128Mix |
\(2^{128}\) |
\(\phi\) |
64 |
DSFMT |
\(2^{19937}\) |
\(2^{128}\) |
53 |
LXM |
\(2^{256}\) |
\(2^{128}\) |
64 |
MT19937 |
\(2^{19937}\) |
\(2^{128}\) |
32 |
PCG32 |
\(2^{64}\) |
\(\phi\) |
32 |
PCG64 |
\(2^{128}\) |
\(\phi\) |
64 |
Philox |
\(2^{256}\) |
\(2^{128}\) |
64 |
SFMT |
\(2^{19937}\) |
\(2^{128}\) |
64 |
SPECK128 |
\(2^{128}\) |
\(2^{64}\) |
64 |
ThreeFry |
\(2^{256}\) |
\(2^{128}\) |
64 |
Xoroshiro128 |
\(2^{128}\) |
\(2^{64}\) |
64 |
Xorshift1024 |
\(2^{1024}\) |
\(2^{512}\) |
64 |
Xoshiro256 |
\(2^{256}\) |
\(2^{128}\) |
64 |
Xoshiro512 |
\(2^{512}\) |
\(2^{256}\) |
64 |
jumped can be used to produce long blocks that are long enough to not
overlap. A jump size of \(\phi\) is the integer value of \(\sqrt{5}/2 - 1\), the
golden ratio, times the full period.
from randomgen.entropy import random_entropy
from randomgen import Xoshiro512, SeedSequence
NUM_STREAMS = 8196
seed = SeedSequence()
blocked_rng = []
last_rng = rng = Xoshiro512(seed)
for i in range(NUM_STREAMS):
blocked_rng.append(last_rng)
last_rng = last_rng.jumped()
Weyl Sequences¶
SFC64 uniquely supports generating streams using distinct
Weyl sequences. Distinct sequences are produced by setting unique values in k`
when initializing SFC64. These increments must be odd.
There is some indication that sparse bit patterns in the Weyl sequence with 50%
or fewer non-zero bits perform best.
import numpy as np
from randomgen.entropy import random_entropy
from randomgen import SFC64, SeedSequence
NUM_STREAMS = 8196
seed_seq = SeedSequence()
# Helper function to generate values with 32 non-zero bits
# The number of non-zero bits is configurable
weyl_inc = SFC64(seed_seq).weyl_increments(1000)
streams = [SFC64(seed_seq, k=k) for k in weyl_inc]
Non-reproducible Sequences¶
randomgen.rdrand.RDRAND uses a hardware based random number generator to
produce non-reproducible sequences on random numbers. These can be used on any
system that supports the RDRAND instruction, subject to the key limitation that
results cannot be reproduced without storing the values generated.
from randomgen import RDRAND
NUM_STREAMS = 8196
streams = [RDRAND() for _ in range(NUM_STREAMS)]