1Parallel Random Number Generation 2================================= 3 4There are three strategies implemented that can be used to produce 5repeatable pseudo-random numbers across multiple processes (local 6or distributed). 7 8.. currentmodule:: numpy.random 9 10.. _seedsequence-spawn: 11 12`~SeedSequence` spawning 13------------------------ 14 15`~SeedSequence` `implements an algorithm`_ to process a user-provided seed, 16typically as an integer of some size, and to convert it into an initial state for 17a `~BitGenerator`. It uses hashing techniques to ensure that low-quality seeds 18are turned into high quality initial states (at least, with very high 19probability). 20 21For example, `MT19937` has a state consisting of 624 22`uint32` integers. A naive way to take a 32-bit integer seed would be to just set 23the last element of the state to the 32-bit seed and leave the rest 0s. This is 24a valid state for `MT19937`, but not a good one. The Mersenne Twister 25algorithm `suffers if there are too many 0s`_. Similarly, two adjacent 32-bit 26integer seeds (i.e. ``12345`` and ``12346``) would produce very similar 27streams. 28 29`~SeedSequence` avoids these problems by using successions of integer hashes 30with good `avalanche properties`_ to ensure that flipping any bit in the input 31input has about a 50% chance of flipping any bit in the output. Two input seeds 32that are very close to each other will produce initial states that are very far 33from each other (with very high probability). It is also constructed in such 34a way that you can provide arbitrary-sized integers or lists of integers. 35`~SeedSequence` will take all of the bits that you provide and mix them 36together to produce however many bits the consuming `~BitGenerator` needs to 37initialize itself. 38 39These properties together mean that we can safely mix together the usual 40user-provided seed with simple incrementing counters to get `~BitGenerator` 41states that are (to very high probability) independent of each other. We can 42wrap this together into an API that is easy to use and difficult to misuse. 43 44.. code-block:: python 45 46 from numpy.random import SeedSequence, default_rng 47 48 ss = SeedSequence(12345) 49 50 # Spawn off 10 child SeedSequences to pass to child processes. 51 child_seeds = ss.spawn(10) 52 streams = [default_rng(s) for s in child_seeds] 53 54.. end_block 55 56Child `~SeedSequence` objects can also spawn to make grandchildren, and so on. 57Each `~SeedSequence` has its position in the tree of spawned `~SeedSequence` 58objects mixed in with the user-provided seed to generate independent (with very 59high probability) streams. 60 61.. code-block:: python 62 63 grandchildren = child_seeds[0].spawn(4) 64 grand_streams = [default_rng(s) for s in grandchildren] 65 66.. end_block 67 68This feature lets you make local decisions about when and how to split up 69streams without coordination between processes. You do not have to preallocate 70space to avoid overlapping or request streams from a common global service. This 71general "tree-hashing" scheme is `not unique to numpy`_ but not yet widespread. 72Python has increasingly-flexible mechanisms for parallelization available, and 73this scheme fits in very well with that kind of use. 74 75Using this scheme, an upper bound on the probability of a collision can be 76estimated if one knows the number of streams that you derive. `~SeedSequence` 77hashes its inputs, both the seed and the spawn-tree-path, down to a 128-bit 78pool by default. The probability that there is a collision in 79that pool, pessimistically-estimated ([1]_), will be about :math:`n^2*2^{-128}` where 80`n` is the number of streams spawned. If a program uses an aggressive million 81streams, about :math:`2^{20}`, then the probability that at least one pair of 82them are identical is about :math:`2^{-88}`, which is in solidly-ignorable 83territory ([2]_). 84 85.. [1] The algorithm is carefully designed to eliminate a number of possible 86 ways to collide. For example, if one only does one level of spawning, it 87 is guaranteed that all states will be unique. But it's easier to 88 estimate the naive upper bound on a napkin and take comfort knowing 89 that the probability is actually lower. 90 91.. [2] In this calculation, we can ignore the amount of numbers drawn from each 92 stream. Each of the PRNGs we provide has some extra protection built in 93 that avoids overlaps if the `~SeedSequence` pools differ in the 94 slightest bit. `PCG64` has :math:`2^{127}` separate cycles 95 determined by the seed in addition to the position in the 96 :math:`2^{128}` long period for each cycle, so one has to both get on or 97 near the same cycle *and* seed a nearby position in the cycle. 98 `Philox` has completely independent cycles determined by the seed. 99 `SFC64` incorporates a 64-bit counter so every unique seed is at 100 least :math:`2^{64}` iterations away from any other seed. And 101 finally, `MT19937` has just an unimaginably huge period. Getting 102 a collision internal to `SeedSequence` is the way a failure would be 103 observed. 104 105.. _`implements an algorithm`: http://www.pcg-random.org/posts/developing-a-seed_seq-alternative.html 106.. _`suffers if there are too many 0s`: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html 107.. _`avalanche properties`: https://en.wikipedia.org/wiki/Avalanche_effect 108.. _`not unique to numpy`: https://www.iro.umontreal.ca/~lecuyer/myftp/papers/parallel-rng-imacs.pdf 109 110 111.. _independent-streams: 112 113Independent Streams 114------------------- 115 116`Philox` is a counter-based RNG based which generates values by 117encrypting an incrementing counter using weak cryptographic primitives. The 118seed determines the key that is used for the encryption. Unique keys create 119unique, independent streams. `Philox` lets you bypass the 120seeding algorithm to directly set the 128-bit key. Similar, but different, keys 121will still create independent streams. 122 123.. code-block:: python 124 125 import secrets 126 from numpy.random import Philox 127 128 # 128-bit number as a seed 129 root_seed = secrets.getrandbits(128) 130 streams = [Philox(key=root_seed + stream_id) for stream_id in range(10)] 131 132.. end_block 133 134This scheme does require that you avoid reusing stream IDs. This may require 135coordination between the parallel processes. 136 137 138.. _parallel-jumped: 139 140Jumping the BitGenerator state 141------------------------------ 142 143``jumped`` advances the state of the BitGenerator *as-if* a large number of 144random numbers have been drawn, and returns a new instance with this state. 145The specific number of draws varies by BitGenerator, and ranges from 146:math:`2^{64}` to :math:`2^{128}`. Additionally, the *as-if* draws also depend 147on the size of the default random number produced by the specific BitGenerator. 148The BitGenerators that support ``jumped``, along with the period of the 149BitGenerator, the size of the jump and the bits in the default unsigned random 150are listed below. 151 152+-----------------+-------------------------+-------------------------+-------------------------+ 153| BitGenerator | Period | Jump Size | Bits | 154+=================+=========================+=========================+=========================+ 155| MT19937 | :math:`2^{19937}` | :math:`2^{128}` | 32 | 156+-----------------+-------------------------+-------------------------+-------------------------+ 157| PCG64 | :math:`2^{128}` | :math:`~2^{127}` ([3]_) | 64 | 158+-----------------+-------------------------+-------------------------+-------------------------+ 159| Philox | :math:`2^{256}` | :math:`2^{128}` | 64 | 160+-----------------+-------------------------+-------------------------+-------------------------+ 161 162.. [3] The jump size is :math:`(\phi-1)*2^{128}` where :math:`\phi` is the 163 golden ratio. As the jumps wrap around the period, the actual distances 164 between neighboring streams will slowly grow smaller than the jump size, 165 but using the golden ratio this way is a classic method of constructing 166 a low-discrepancy sequence that spreads out the states around the period 167 optimally. You will not be able to jump enough to make those distances 168 small enough to overlap in your lifetime. 169 170``jumped`` can be used to produce long blocks which should be long enough to not 171overlap. 172 173.. code-block:: python 174 175 import secrets 176 from numpy.random import PCG64 177 178 seed = secrets.getrandbits(128) 179 blocked_rng = [] 180 rng = PCG64(seed) 181 for i in range(10): 182 blocked_rng.append(rng.jumped(i)) 183 184.. end_block 185 186When using ``jumped``, one does have to take care not to jump to a stream that 187was already used. In the above example, one could not later use 188``blocked_rng[0].jumped()`` as it would overlap with ``blocked_rng[1]``. Like 189with the independent streams, if the main process here wants to split off 10 190more streams by jumping, then it needs to start with ``range(10, 20)``, 191otherwise it would recreate the same streams. On the other hand, if you 192carefully construct the streams, then you are guaranteed to have streams that 193do not overlap. 194