Random seed is replication across child processes

[mosco]

When spawning child processes using the multiprocessing module, it appears that all child processes share the parent's random seed.

This creates a subtle and difficult to detect bug in the common use case of embarrassingly parallel multicore monte-carlo simulations. In such simulations the different processes will generate the same data, thus the final averaged result will be much less accurate (i.e. have higher variance) than expected.

Many people seem to have been bitten by this issue (and surely, many more are unaware of a silent bug in their simulation). e.g.
https://stackoverflow.com/questions/12915177/same-output-in-different-workers-in-multiprocessing
http://forum.cogsci.nl/index.php?p=/discussion/1441/solved-numpy-random-state-seems-to-repeat-across-multiple-os-runs

Note that Python's random module does not suffer from this problem.

[bashtage]

If you copy a RandomState you get that RandomState. That means the state -- not the seed -- is the same. The state is what matters for determining the sequence of random numbers.

Any correct method requires you to initialize a RandomState within your child processes. Better is to use the improved RandomState here which explicitly supports generating 1000s or guaranteed distinct streams using either a stream-based generator (PCG) or jump-ahead to advance to state as-if many draws have been made (2**128).

[bashtage]

Worth noting that the random module also doesn't automatically reseed.

import random
from joblib import Parallel, delayed
def gs(v):
    return random.getstate()

res = Parallel()(delayed(gs)(i) for i in range(10))

for r in res:
    print(r==res[0])
    
True
True
True
True
True
True
True
True
True
True

[njsmith]

Huh, I thought we had an issue open for this already, but I can't find it.

Related: #9248

This is very difficult to fix in Python versions before 3.7. In 3.7 it's not too bad using os.register_at_fork, but no-one has done the work yet.

[mosco]

For some reason, when using multiprocessing.Pool, one gets the desired behavior using Python's random but not using NumPy's!

This behaves differently than the joblib example of @bashtage. Here is my code:

import multiprocessing
import random
import numpy.random

def python_random(whatever):
    return random.random()

def numpy_random(whatever):
    return numpy.random.random()

pool = multiprocessing.Pool(10)
print('Python random:', pool.map(python_random, [0]*10))
print('Numpy random:', pool.map(numpy_random, [0]*10))

Output on Python 3.6.2, Numpy 1.13.1:

Python random: [0.3232979890381721, 0.3440537969129853, 0.8325653603083111, 0.6378306617495579, 0.5425156816989956, 0.2300538134522857, 0.27904636944778316, 0.45440811632035105, 0.6899393498838224, 0.9521238318896579]

Numpy random: [0.7768863343872144, 0.7768863343872144, 0.7768863343872144, 0.7768863343872144, 0.7768863343872144, 0.7768863343872144, 0.7768863343872144, 0.7768863343872144, 0.7768863343872144, 0.7768863343872144]

[bashtage]

The difference is that NumPy is using n copies of a RandomState while random is using a single instance. This means that when you are drawing randoms from random you will be effectively in single threaded mode since they are all using the same state. This will limit multiprocessing if random number generation is an important part of your time. NumPy should, OTOH, have no restrictions across workers and so you will get the maximum benefit.

It is very difficult to enforce the behavior you are seeing from random using NumPy AFAICT since the sharing a single complex object in multiprocessing is challenging (see remote manager in the multiprocessing docs)

[bashtage]

See this answer for details and the complexity it introduces.

[bashtage]

Here is a minimal example using a shared RandomState

from joblib import Parallel, delayed
import numpy as np
from multiprocessing.managers import BaseManager


class SharedRandomState(object):
    random_state = np.random.RandomState(1234)
    def random_sample(self):
        return self.random_state.random_sample()

class MyManager(BaseManager):
    pass

MyManager.register('SharedRandomState', SharedRandomState)

def worker(rs):
    return rs.random_sample()

if __name__ == '__main__':
    manager = MyManager()
    manager.start()
    shared = manager.SharedRandomState()
    res = Parallel(n_jobs=4)(delayed(worker)(shared) for i in range(10))
    print('From multiple workers')
    print(np.array(res))
    
    # Verify
    random_state = np.random.RandomState(1234)
    print('From a single instance')
    print(random_state.random_sample(size=10))

From multiple workers
[ 0.19151945  0.62210877  0.43772774  0.78535858  0.77997581  0.27259261
  0.27646426  0.80187218  0.95813935  0.87593263]
From a single instance
[ 0.19151945  0.62210877  0.43772774  0.78535858  0.77997581  0.27259261
  0.27646426  0.80187218  0.95813935  0.87593263]

[mosco]

Upon further examination, Python 3.6's default behavior of reseeding in the child processes was hard-coded in the _launch method of multiprocessing.popen_fork.Popen and is specific to the built-in random module. Starting from Python 3.7 this was changed to use the new os.register_at_fork mechanism.

So it may be difficult to fix on older versions, but surely we can reseed on Python 3.7 using the same mechanism as @njsmith suggested. I'm willing to do the work if one of the core developers is willing to accept the pull request.

[bashtage]

The hard question for numpy.random is whether this is a behavior break since the stream of generated randoms will differ.

Reseeding in child processes is not a good way to generate multiple streams of random numbers since it is not easy to verify the statistical properties of the generators. It would be intellectually better to use a generator that supports multiple streams -- this would probably require a new version of RandomState due to the compat guarantee.

[bashtage]

In 3.5 random is not reseeding but is using a shared state across workers. This can be seen below since the only thing that changes across draws is the position within the 624 element state vector that underlies the MT.

import random
from joblib import Parallel, delayed

def f(v):
    random.random()
    return random.getstate()
    
    res = Parallel(n_jobs=4)(delayed(f)(i) for i in range(10))
    
for i in range(len(res)-1):
    print('624 element state identical? {0}'.format(res[i][1][:-1] == res[i+1][1][:-1]))
    print('Position identical? {0}'.format(res[i][1][-1] == res[i+1][1][-1]))
    print('Positions: {0}, {1}'.format(res[i][1][-1],res[i+1][1][-1]))

624 element state identical? True
Position identical? False
Positions: 44, 46
624 element state identical? True
Position identical? False
Positions: 46, 48
624 element state identical? True
Position identical? False
Positions: 48, 50
624 element state identical? True
Position identical? False
Positions: 50, 52
624 element state identical? True
Position identical? False
Positions: 52, 54
624 element state identical? True
Position identical? False
Positions: 54, 56
624 element state identical? True
Position identical? False
Positions: 56, 58
624 element state identical? True
Position identical? False
Positions: 58, 60
624 element state identical? True
Position identical? False
Positions: 60, 62

[bashtage]

@mosco You are right -- my examples were too trivial. joblib is finishing 10 before spawning a second process and so they were all in the same process. Running more shows clearly many different states after seed has been called.

[njsmith]

The hard question for numpy.random is whether this is a behavior break since the stream of generated randoms will differ.

We should re-seed each RandomState object that had seed=None, and only those.

# no initial seed, re-seeded on fork
r1 = np.random.RandomState()
# explicit seed given, NOT re-seeded on fork
r2 = np.random.RandomState(0)
# seed was given, but then randomized; re-seeded on fork
r3 = np.random.RandomState(0); r3.seed()

[bashtage]

That is probably reasonable since when you use a default seed you don't really know what you get anyway.

I took a look at random and it ignores provided seed values when using in multiprocessing so that values produced in multiple processes after random.seed(123) is called in the main process are not reproducible.

[pv]

If I write ``` r = np.random.RandomState() z = r.random() pid = os.fork() ``` then why should I expect the internal state of r changes in the child process? It also sounds a bit dangerous that it would get reseeded based on the *time*.

[bashtage]

It is normally reseeded with /dev/urandom or the Windows equiv. In theory if one tried to spawn a huge number it might be possible to exhaust the pool, in which case there is a has of time + pid which should be safe (but not necessarily very random, even after hashing).