# Multithreading and random number generators

**URL:** <https://discourse.julialang.org/t/multithreading-and-random-number-generators/49777>\
**Category:** General Usage\
**Tags:** multithreading\
**Created:** [November 8, 2020, 8:47am UTC](https://discourse.julialang.org/t/multithreading-and-random-number-generators/49777 "2020-11-08T08:47:55Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Skoffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skoffer/32/378_2.png) [@Skoffer](https://discourse.julialang.org/u/Skoffer)\
**Post date:** [November 8, 2020, 8:47am UTC](https://discourse.julialang.org/t/multithreading-and-random-number-generators/49777/1 "2020-11-08T08:47:55Z")

</div>

How one can approach multithreading simulations, which uses random numbers generators? Especially if one wants to use non Base generators, like `StableRNGs.jl`, `RandomNumbers.jl` and others.

Here is a simple example of a random walk, which I want to make parallel

```julia
using Random
using StableRNGs

function randomwalk!(tbl, i, rng)
    len = size(tbl, 1)
    tbl[1, i] = rand(rng)
    @inbounds for j in 2:len
        tbl[j, i] = (tbl[j - 1, i]) + rand(rng)
    end
end

function randomwalk(nparticles, t, rng)
    tbl = Array{Float64, 2}(undef, t, nparticles)
    map(i -> randomwalk!(tbl, i, rng), 1:nparticles)

    return tbl
end

rng = StableRNG(2020)
randomwalk(10, 1000, rng)

```

Here is what I want to achieve

1. I want to be able to use different number generators (MersenneTwister, StableRNG, etc)
2. I want calculations to be run on multiple threads
3. Calculations should be random in a sense, that there shouldn’t be a correlation between random numbers on different threads. Well, as much as we can talk of non-correlation of pseudo-random numbers anyway.
4. It would be good for this algorithm to be not too complicated, so it can be easily used in other similar calculations.
5. As a bonus it would be good to have these calculations reproducible. But this is not a very strict requirement, i.e. if one thread is running ahead of another and grab RNG first, it’s ok, as long as the random sequence is not correlated with other threads.

I guess I can do some bookkeeping in `randomwalk` function and generate some amount of RNGs, but which way is more appropriate;

1. Generate new RNG on-the-fly for each new index `i`, where RNG seed should be generate from the main generator.
2. Pregenerate RNGs for all elements, and then feed them to `randomwalk!` function.

```julia
rngs = StableRNG.(rand(rng, 1:typemax(Int), length(nparticles))
map(i -> randomwalk!(tbl, i, rngs[i]), 1:nparticles)

```

by the way, is it correct to use `1:typemax(Int)` or range should be made even larger?  
3. Generate only `Threads.nthreads()` number of new RNGs and determine which one should be used inside `randomwalk!` function?

In this discussion [Reproducible multithreaded monte-carlo / Task local random](https://discourse.julialang.org/t/reproducible-multithreaded-monte-carlo-task-local-random/35269) if I understand it correctly, it is claimed, that `MersenneTwister` is thread-safe now. Is it possible to somehow “upgrade” other RNG to make them thread-safe too? And on the other hand, wouldn’t thread safety of RNG slow down multithreading calculations?

---

<div class="post-metadata">

**Author:** ![rfourquet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rfourquet/32/3610_2.png) [@rfourquet](https://discourse.julialang.org/u/rfourquet)\
**Post date:** [November 8, 2020, 9:18am UTC](https://discourse.julialang.org/t/multithreading-and-random-number-generators/49777/2 "2020-11-08T09:18:48Z")

</div>

I don’t really have an answer, but just few comments.

> [@Skoffer](#):
>
> Generate new RNG on-the-fly for each new index `i` , where RNG seed should be generate from the main generator.

That would really depend on the RNG type. Generating a new `MersenneTwister` on each iteration might be very expensive, whereas `StableRNG` is very cheap to create.

> [@Skoffer](#):
>
> 1. Generate only `Threads.nthreads()` number of new RNGs and determine which one should be used inside `randomwalk!` function?

So that might be the most general solution, as it works also for expensive-to-create RNGs.

> [@Skoffer](#):
>
> it is claimed, that `MersenneTwister` is thread-safe now.

What is thread-safe is to use the “global RNG” concurrently, as in fact there a multiple “global” RNGs, one for each thread. But it’s unsafe to use one given RNG concurrently, as is the case for most RNGs.

> [@Skoffer](#):
>
> Is it possible to somehow “upgrade” other RNG to make them thread-safe too?

I think your solution 3. above is the one generally shared on Julia forums for that.

> [@Skoffer](#):
>
> is it correct to use `1:typemax(Int)` or range should be made even larger?

I think that really depends on the RNG type, i.e. how it does initialization. Note also that `StableRNG` was not really intended for this use case; it currently accepts a seed in `0 <= seed <= typemax(UInt)` (but it could be made to accept a seed almost equal to `typemax(UInt128)`), but no facility is provided to produce two `StableRNGs` which are as much not correlated as possible. But maybe random seeds would be good enough.

---

<div class="post-metadata">

**Author:** ![Skoffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skoffer/32/378_2.png) [@Skoffer](https://discourse.julialang.org/u/Skoffer)\
**Post date:** [November 8, 2020, 9:31am UTC](https://discourse.julialang.org/t/multithreading-and-random-number-generators/49777/3 "2020-11-08T09:31:09Z")

</div>

So, what is happening if I am using RNG naively in multithreading context?

Something like

```julia
using Random
using StableRNGs

function f(rng)
  res = Vector{Float64}(undef, 100_000)
  Threads.@threads for i in 1:100_000
     res[i] = rand(rng)
  end
  return res
end

f(Random.GLOBAL_RNG)
f(StableRNG(2020))

```

They are definitely use all cores in both cases and produce some results, but what is going on with RNG at this moment? Are values duplicated somehow or `rng` is locked between threads or it is something else entirely?

---

<div class="post-metadata">

**Author:** ![rfourquet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rfourquet/32/3610_2.png) [@rfourquet](https://discourse.julialang.org/u/rfourquet)\
**Post date:** [November 8, 2020, 9:45am UTC](https://discourse.julialang.org/t/multithreading-and-random-number-generators/49777/4 "2020-11-08T09:45:22Z")

</div>

> [@Skoffer](#):
>
> but what is going on with RNG at this moment? Are values duplicated somehow or `rng` is locked between threads or it is something else entirely?

For `f(Random.GLOBAL_RNG)`, something special happens: `GLOBAL_RNG` is _not_ a `MersenneTwister`, since Julia 1.3 it’s a special type which was created to not break compatibility (nice trick from Jameson), so that you can still call `rand(GLOBAL_RNG)` while being thread-safe; but it’s equivalent to calling `rand(Random.default_rng())`, where `default_rng()` returns the “thread-local global RNG”. So `rand(GLOBAL_RNG)` is thread-safe. But `rand(::MersenneTwister)` is not, and your `f(StableRNG(2020))` is not, and will result in corrupted state.

---

<div class="post-metadata">

**Author:** ![Skoffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skoffer/32/378_2.png) [@Skoffer](https://discourse.julialang.org/u/Skoffer)\
**Post date:** [November 8, 2020, 9:53am UTC](https://discourse.julialang.org/t/multithreading-and-random-number-generators/49777/5 "2020-11-08T09:53:50Z")

</div>

So, I shouldn’t use them naively, because state is corrupted and results are unreliable. Good to know, thank you.

One last question, is it a viable approach to extend `rand` function to accept vector of `RNG` in this manner

```julia
rand(rngs::AbstractVector) = rand(rngs[Threads.threadid()])

# Use it like this.
# single thread
rng = StableRNG(2020)
randomwalk(10, 1000, rng)

# multi-thread (of course `@Threads.threads` or something similar should be used inside `randomwalk` to actually use multithreading
rngs = StableRNG.(2020 .+ collect(1:Threads.nthreads()))
randomwalk(10, 1000, rngs)

```

Of course, other versions of `rand` should be extended too if needed, may be something like

```julia
# to avoid type piracy
struct VectorRNG{T}
  rngs::Vector{T}
end

rand(rngs::VectorRNG, args...) = rand(rngs.rngs[Threads.threadid()], args...)

```

---

<div class="post-metadata">

**Author:** ![rfourquet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rfourquet/32/3610_2.png) [@rfourquet](https://discourse.julialang.org/u/rfourquet)\
**Post date:** [November 8, 2020, 10:42am UTC](https://discourse.julialang.org/t/multithreading-and-random-number-generators/49777/6 "2020-11-08T10:42:25Z")

</div>

> [@Skoffer](#):
>
> One last question, is it a viable approach to extend `rand` function to accept vector of `RNG` in this manner

I would use your second solution, building on it here is how I would go about it:

```julia

struct VectorRNG{T} <: AbstractRNG
    rngs::Vector{T}
end

struct VectorRNGSampler{T,X} <: Random.Sampler{X}
    sp::T

    VectorRNGSampler(sp::T) where T = new{T,Random.gentype(sp)}(sp)
end

for N = (Val{1}, Val{Inf})
    for S = (Any, AbstractFloat)
        @eval Random.Sampler(::Type{VectorRNG{T}}, ::Type{X}, n::$N) where {T,X<:$S} =
            VectorRNGSampler(Random.Sampler(T, X, n))
    end
    for S = (Any, AbstractUnitRange{<:Base.BitInteger64})
        @eval Random.Sampler(::Type{VectorRNG{T}}, x::X, n::$N) where {T,X<:$S} =
            VectorRNGSampler(Random.Sampler(T, x, n))
    end
end

Random.rand(rngs::VectorRNG, t::VectorRNGSampler) = rand(rngs.rngs[Threads.threadid()], t.sp)

```

There has to be loops to avoid ambiguities, maybe it’s possible to improve `Random` to get less of these.

[DISCLAIMER: I almost didn’t use multi-threading in Julia!]

---

<div class="post-metadata">

**Author:** ![palday](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palday/32/12640_2.png) [@palday](https://discourse.julialang.org/u/palday)\
**Post date:** [December 4, 2020, 5:00pm UTC](https://discourse.julialang.org/t/multithreading-and-random-number-generators/49777/7 "2020-12-04T17:00:26Z")

</div>

I unfortunately missed this discussion, but my [JuliaCon talk](https://youtu.be/qhfQKIDR7Ms) is very relevant for this and discusses how we solved this problem in [MixedModels.jl](https://github.com/JuliaStats/MixedModels.jl/). There are various things we could do to make it more performant, but we have a strong reproducibility guarantee that is independent of the number of threads and doesn’t require the big allocation I think is inherent in your (2).

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [March 15, 2021, 2:41pm UTC](https://discourse.julialang.org/t/multithreading-and-random-number-generators/49777/8 "2021-03-15T14:41:01Z")

</div>

I found an other solution that doesn’t use locks, and I can have, depending of the needs, (a) different results, (b) different results but the same “sequence” whenever I reset the random number generator or (c) the same result at each call. However it works only with `MersenneTwister`. I can’t get it working with `StableRNG`, that is indeed what I use in the rest of my library (above all, to guarantee the same random number independently than the specific Julia version).

```julia
using Statistics, Random,StableRNGs

x = rand(100)

function generateParallelRngs(rng::Union{Random._GLOBAL_RNG,MersenneTwister,StableRNGs.LehmerRNG}, n::Integer)
    step = rand(rng,big(10)^20:big(10)^40) # making the step random too !
    rngs = Vector{Union{MersenneTwister,Random._GLOBAL_RNG,StableRNGs.LehmerRNG}}(undef, n)
    rngs[1] = copy(rng)
    for i = 2:n
        rngs[i] = Future.randjump(rngs[i-1], step)
    end
    return rngs
end

function innerFunction(bootstrappedx; rng=MersenneTwister())
     sum(bootstrappedx .* rand(rng) ./ 0.5)
end

function outerFunction(x;rng = MersenneTwister())
    rngs = generateParallelRngs(rng,Threads.nthreads())
    results = Array{Float64,1}(undef,30)
    Threads.@threads for i in 1:30
        tsrng = rngs[Threads.threadid()] # Thread safe random number generator
        toSample = rand(tsrng, 1:100,100)
        bootstrappedx = x[toSample]
        innerResult = innerFunction(bootstrappedx, rng=tsrng)
        results[i] = innerResult
    end
    overallResult = mean(results)
    return overallResult
end

# Different sequences..
outerFunction(x)
outerFunction(x)
outerFunction(x)

# Different values, but same sequence
mainRng = MersenneTwister(123)
outerFunction(x, rng=mainRng)
outerFunction(x, rng=mainRng)
outerFunction(x, rng=mainRng)

# Same value at each call
outerFunction(x,rng=MersenneTwister(123))
outerFunction(x,rng=MersenneTwister(123))
outerFunction(x,rng=MersenneTwister(123))

```

If I try `outerFunction(x,rng=StableRNG(123))` I have however that `randJump` is not implemented for `StableRNGs`. Any alternative ?