How do I deal with Random Number generation when multithreading

parallel
#1

If I am using Threads.@threads to run my Monte Carlo simulation in paralell, how do I set up Mersenne Twister RNGs for each thread so that I am both able to set a seed, in order to be able to reproduce the results, but also set a random seed when not wanting reproducability?

Now I am using a “home cooked” skip-ahead solution (which I am not even sure is correct/at least not robust) .

  function set_rngs(num_draws_per_path,num_paths)
      seed_1 = 1234;
      seed_2 = seed1+num_draws_per_path*num_paths/4;
      seed_3 = seed2+num_draws_per_path*num_paths/4;
      seed_4 = seed3+num_draws_per_path*num_paths/4;
      return [MersenneTwister(seed_1),MersenneTwister(seed_2),MersenneTwister(seed_3),MersenneTwister(seed_4)];
  end

then setting rngs = set_rngs(s,num_paths) and

Threads.@threads for p=1:num_paths
      rand_nums = randn(rngs[Threads.threadid()],num_draws_per_path);

      do stuff.....
end
``
#2

You can use randjump to get independent streams (by default randjump jumps ahead 10^20 steps)

EDITED

num_paths = 10
num_draws_per_path = 1000

# set random seed and output with results for later reproducibility
rng_seed = Base.Random.make_seed()

# to reproduce results, set explicit seed from previous run
# rng_seed = UInt32[0xbd8e7dc4,0x0d00f383,0xdccee4c9,0x553afb99]

rng = MersenneTwister(rng_seed)

# RNG per thread
rngs = randjump(rng, Threads.nthreads());
rand_nums = Vector{Vector{Float64}}(num_paths);

Threads.@threads for p = 1:num_paths
    rand_nums[p] = randn(rngs[Threads.threadid()], num_draws_per_path);
end

Reproducibility using a separate RNG per thread relies on loop iterations being allocated to threads deterministically. I think this is true with current implementation. However changing the number of threads available will break reproducibility.

In that case, a safer option might be to use a separate RNG per path:

# RNG per path
rngs = randjump(rng, num_paths);
rand_nums = Vector{Vector{Float64}}(num_paths);

Threads.@threads for p = 1:num_paths
    rand_nums[p] = randn(rngs[p], num_draws_per_path);
end
#3

Thanks! Here’s my set_rngs function. Will only run with 4 threads, so do not need to worry about reproducibility with varying number of threads.

  function set_rngs(s,num_paths,random_seed)

    if random_seed == true
      rng_seed = Base.Random.make_seed()
      rng = MersenneTwister(rng_seed)
      return randjump(rng, Threads.nthreads());
    else
      rng_seed =UInt32[0xbd8e7dc4,0x0d00f383,0xdccee4c9,0x553afb99];
      rng = MersenneTwister(rng_seed)
      return randjump(rng, Threads.nthreads());
    end
  end
#4

@vgdev I presume that you’re aware that, while this approach will give separated sequences of random numbers there is no guarantee of statistical independence. RND generators produce a cyclical sequence, the seed just starts you at a different point in the sequence, and it’s not possible to decide how far apart two seeds place you in that sequence.

If you want truly independent streams, you need to either, (i) generate all random deviates using a single generator instance (single seed), and assign the individual deviates to separate tasks, or (ii) use a different generator, one suited to parallel generation, (take a look at Random123).

I hope this helps.

#5

Create one RNG object per thread?

#6

If you do (i) as you suggest above, don’t you also just get two points on the same (cyclical) sequence?

AFAIK randjump just codes best practices, cf this implementation and the reference there.

#7

Hi @Tamas_Papp,

Sort of. If you use a single random generator, you get a sequence of random deviates which (depending on the quality of the generator) are statistically independent from one another. My suggestion is, if you know in advance how many deviates (total) you will need, you simply generate the entire sequence then slice it into sections of length N/p where p is the number of parallel threads, which consume the deviates, and N the total number of deviates.

randjump takes advantage of the fact that we typically don’t need full independence in our random deviates. It is often sufficient if each of the parallel threads generates its own sequence which is independent with respect to itself. If there are spatial correlations (across threads) or temporo-spatial correlates (between different thread instances at different points in time) this only leads to problems in certain applications. randjump is better than hardcoding the jump since we do not know the effect of different seed values on these problems.

Mersenne Twister is beautiful because it requires so few lines of code. It is also incredibly well tested as a generator of random sequences. But this makes it very tempting to pop it inside of a parallelised instance of code and just use each instance with a different seed. I just wanted to point out that this can lead to issues in certain applications. (And they may be difficult to identify - The best debug would be to do a test with an external single instance of Mersenne Twister, much slower, but good for testing.)