# Random numbers and threads

**URL:** <https://discourse.julialang.org/t/random-numbers-and-threads/77364>\
**Category:** General Usage\
**Tags:** question, multithreading, random\
**Created:** [March 3, 2022, 3:31pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364 "2022-03-03T15:31:24Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![jw3126](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jw3126/32/3086_2.png) [@jw3126](https://discourse.julialang.org/u/jw3126)\
**Post date:** [March 3, 2022, 3:31pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/1 "2022-03-03T15:31:24Z")

</div>

How to use random numbers in conjunction with threads? In the past I used the following pattern:

```julia
using Random
function simulation(rngs)
    Threads.@threads for i in 1:1000000
        rng = rngs[Threads.threadid()]
        x = rand(rng)
        # do more stuff with rng
    end
end

baserng = MersenneTwister(0)
rngs = [Random.jump(baserng, i*10^15) for i in 1:Threads.nthreads()]
simulation(rngs)

```

However it is my understanding, that this is unsafe. As I understand julia is free start executing a task on one thread and continue execution on another thread. For instance julia could start two tasks on thread1, pause and finish them on thread2 and thread3. In this case thread2 and thread3 would both access `rngs[1]` and cause havoc.

- Is my understanding correct?
- As of julia 1.7 can the described corruption actually occur? Or may it only happen with future julia versions?
- What would be a better pattern here?

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [March 3, 2022, 3:48pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/2 "2022-03-03T15:48:01Z")

</div>

I seem to remember that

```julia
help?> rand()
  rand([rng=GLOBAL_RNG], [S], [dims...])

```

has a generator argument. I believe you should be on the safe side using thread local generator objects then?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 3, 2022, 4:01pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/3 "2022-03-03T16:01:41Z")

</div>

You can just use the default RNG and it is thread safe. In Julia 1.7, it is not only thread-safe but is reproducible independent of the thread scheduling: [Julia 1.7 Highlights](https://julialang.org/blog/2021/11/julia-1.7-highlights/#new_rng_reproducible_rng_in_tasks)

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [March 3, 2022, 4:04pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/4 "2022-03-03T16:04:47Z")

</div>

But what about the overhead associated to the concurrent access to the global random state (if that matters in the specific application). For example:

```julia
julia> import Random
       function simulation_default(nthreads)
           s = zeros(nthreads)
           Threads.@threads for it in 1:nthreads
               for j in it:nthreads:10^6 # easy splitter
                   s[it] += rand()
               end
           end
           sum(s)
       end
simulation_default (generic function with 1 method)

julia> import Random
       function simulation_perthread(nthreads)
           s = zeros(nthreads)
           Threads.@threads for it in 1:nthreads
               rng = Random.MersenneTwister()
               for j in it:nthreads:10^6 # easy splitter
                   s[it] += rand(rng)
               end
           end
           sum(s)
       end
simulation_perthread (generic function with 1 method)

julia> @btime simulation_default(4);
  3.945 ms (22 allocations: 2.06 KiB)

julia> @btime simulation_perthread(4);
  2.405 ms (66 allocations: 80.25 KiB)

```

(ps: note that in this particular example running without threading at all is faster than any of the alternatives)

---

<div class="post-metadata">

**Author:** ![jw3126](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jw3126/32/3086_2.png) [@jw3126](https://discourse.julialang.org/u/jw3126)\
**Post date:** [March 3, 2022, 4:07pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/5 "2022-03-03T16:07:09Z")

</div>

Ah thanks for pointing that out! I would prefer however if I could launch my simulations without having side effects on the global rng. Is that possible?

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [March 3, 2022, 4:11pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/6 "2022-03-03T16:11:08Z")

</div>

The pattern I’ve used above should be fine. For something more elegant, see: [Sum result of Threads.foreach() - #10 by tkf](https://discourse.julialang.org/t/sum-result-of-threads-foreach/76701/10)

---

<div class="post-metadata">

**Author:** ![jw3126](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jw3126/32/3086_2.png) [@jw3126](https://discourse.julialang.org/u/jw3126)\
**Post date:** [March 3, 2022, 4:12pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/7 "2022-03-03T16:12:53Z")

</div>

Thanks, I guess that solves my MWE. However for my real application, I need something more composable. E.g. there are multiple `Theads.@threads` loops involving `rand` scattered over many functions. What I would like is a thread safe way to pass around rng state.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 3, 2022, 4:14pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/8 "2022-03-03T16:14:46Z")

</div>

> [@lmiq](#):
>
> But what about the overhead associated to the concurrent access to the global random state

There isn’t a global random state. There’s a [per-task random state](https://github.com/JuliaLang/julia/blob/cbf6c1c8965aaba00683bfe4d54be9889368c1f6/stdlib/Random/src/Xoshiro.jl#L95-L114) used [by default](https://github.com/JuliaLang/julia/blob/82dc13005039a8215f1c910f48870ef378403cbc/stdlib/Random/src/RNGs.jl#L358-L359), if I understand correctly.

In general, now that [Julia 1.7 implemented task migration between threads](https://github.com/JuliaLang/julia/pull/40715), it is [problematic to use `threadid` to emulate thread-local state](https://github.com/JuliaLang/julia/pull/39168#issuecomment-757551959), and one should instead [use task-local state](https://docs.julialang.org/en/v1/base/parallel/#Base.task_local_storage-Tuple%7BAny%7D).

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 3, 2022, 4:22pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/9 "2022-03-03T16:22:08Z")

</div>

> [@jw3126](#):
>
> Ah thanks for pointing that out! I would prefer however if I could launch my simulations without having side effects on the global rng. Is that possible?

Yes, because the “global” rng is actually a task-local state.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [March 3, 2022, 4:24pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/10 "2022-03-03T16:24:28Z")

</div>

Uhm… so what explains the best performance of `simulation_perthread` above?

> [@jw3126](#):
>
> What I would like is a thread safe way to pass around rng state.

I think this is a safe way:

```julia
julia> using Random

julia> ntasks = 3 # not necessarily equal to the number of threads
3

julia> const rngs = [MersenneTwister() for i in 1:ntasks]
3-element Vector{MersenneTwister}:
 MersenneTwister(0x38c65f20ec9d5c7eccf63a103295d3c)
 MersenneTwister(0xfd4719ad02ca867d1ba23e9810767b28)
 MersenneTwister(0x8ce6e5b0de492eddd0b904d019c2cd1d)

julia> function simulation_pertask(ntasks)
           s = zeros(ntasks)
           Threads.@threads for it in 1:ntasks
               rng = rngs[it]
               for j in it:ntasks:10^6 
                   s[it] += rand(rng)
               end
           end
           sum(s)
       end
simulation_pertask (generic function with 1 method)

julia> @btime simulation_pertask(3);
  2.467 ms (22 allocations: 2.05 KiB)

```

Note that `threadid()` is never used, and the tasks have a unique index `it`, which is actually independent on the actual number of threads run.

---

<div class="post-metadata">

**Author:** ![jw3126](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jw3126/32/3086_2.png) [@jw3126](https://discourse.julialang.org/u/jw3126)\
**Post date:** [March 3, 2022, 4:25pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/11 "2022-03-03T16:25:16Z")

</div>

```julia
using Random
Random.seed!(0)
@show rand()
Random.seed!(0)
@show rand()
function doit()
    Threads.@threads for i in 1:100
        rand()
    end
end
Random.seed!(0)
doit()
@show rand()

```

gives me

```julia
rand() = 0.4056994708920292
rand() = 0.4056994708920292
rand() = 0.08910091331143466

```

So there seems to be a side effect?

---

<div class="post-metadata">

**Author:** ![jw3126](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jw3126/32/3086_2.png) [@jw3126](https://discourse.julialang.org/u/jw3126)\
**Post date:** [March 3, 2022, 4:26pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/12 "2022-03-03T16:26:59Z")

</div>

That is a cool trick. It may create load balancing headaches however.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 3, 2022, 4:27pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/13 "2022-03-03T16:27:13Z")

</div>

> [@lmiq](#):
>
> `rng = Random.MersenneTwister()`

This is not even the same rng as the default (Xoshiro).

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 3, 2022, 4:29pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/14 "2022-03-03T16:29:49Z")

</div>

> [@jw3126](#):
>
> So there seems to be a side effect?

`Threads.@threads` is still allowed run things in the main task (and it probably does — that task is already running, so why shouldn’t it use it for one of the threads?).

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [March 3, 2022, 4:32pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/15 "2022-03-03T16:32:09Z")

</div>

> [@stevengj](#):
>
> This is not even the same rng as the default (Xoshiro).

Yes, true, that changed as well in 1.7. Thanks (Xoshiro is slower in that minimal test, effectively).

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [March 3, 2022, 4:36pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/16 "2022-03-03T16:36:49Z")

</div>

> [@jw3126](#):
>
> That is a cool trick. It may create load balancing headaches however.

The trick some time ago was doing something like this:

```julia
RNG = [randjump(Random.MersenneTwister(options.seed),big(10)^20)] 
foreach(_ -> push!(RNG, randjump(last(RNG),big(10)^20)), 2:Threads.nthreads())

```

I don’t know if this is still recommended in any case. It is in one of my packages, though, where the parallel use of `rand` was a bottleneck.

edit: Though I think that this was to avoid the overlap of the sequences. I can’t find the thread exactly where I got that from.

---

<div class="post-metadata">

**Author:** ![jw3126](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jw3126/32/3086_2.png) [@jw3126](https://discourse.julialang.org/u/jw3126)\
**Post date:** [March 3, 2022, 4:39pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/17 "2022-03-03T16:39:47Z")

</div>

> [@stevengj](#):
>
> > [@jw3126](#):
> >
> > So there seems to be a side effect?
> 
> `Threads.@threads` is still allowed run things in the main task (and it probably does — that task is already running, so why shouldn’t it use it for one of the threads?).

> [@stevengj](#):
>
> > [@jw3126](#):
> >
> > Ah thanks for pointing that out! I would prefer however if I could launch my simulations without having side effects on the global rng. Is that possible?
> 
> Yes, because the “global” rng is actually a task-local state.

I am confused by these two answers, they seem contradicting to me. So how to do thread safe rng without side effect on the global rng?

---

<div class="post-metadata">

**Author:** ![fredrikekre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikekre/32/1688_2.png) [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)\
**Post date:** [March 3, 2022, 4:45pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/18 "2022-03-03T16:45:13Z")

</div>

The task local RNG is seeded with a random number from the RNG of the current task – ~~scheduling~~ creating TaskB from TaskA therefore advances the RNG of TaskA.

---

<div class="post-metadata">

**Author:** ![jw3126](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jw3126/32/3086_2.png) [@jw3126](https://discourse.julialang.org/u/jw3126)\
**Post date:** [March 3, 2022, 5:03pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/19 "2022-03-03T17:03:47Z")

</div>

Ok I think I get it. If I spawn things in my main task, its rng state changes. If I spawn things in another task, that other tasks rng state changes, but the main task rng state is untouched.

```julia
using Random

println("no spawning")
Random.seed!(0)
@show rand()
println("no spawning")
Random.seed!(0)
@show rand()
println("spawning once from main")
Random.seed!(0)
@sync begin
    Threads.@spawn nothing
end
@show rand()

println("spawning once from main")
Random.seed!(0)
@sync begin
    Threads.@spawn nothing
end
@show rand()
println("spawning twice from main")
Random.seed!(0)
@sync begin
    Threads.@spawn nothing
    Threads.@spawn nothing
end
@show rand()

println("spawning once from child")
task = Task() do 
    Threads.@spawn nothing
    rand()
end
Random.seed!(0)
@sync schedule(task)

@show rand()

```

```julia
no spawning
rand() = 0.4056994708920292
no spawning
rand() = 0.4056994708920292
spawning once from main
rand() = 0.6616126907308237
spawning once from main
rand() = 0.6616126907308237
spawning twice from main
rand() = 0.2895098423219379
spawning once from child
rand() = 0.4056994708920292
0.4056994708920292

```

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [March 3, 2022, 9:27pm UTC](https://discourse.julialang.org/t/random-numbers-and-threads/77364/20 "2022-03-03T21:27:22Z")

</div>

> [@jw3126](#):
>
> That is a cool trick. It may create load balancing headaches however.

Maybe I misinterpreted a point here. You can tune the load balancing in this pattern by setting the number of tasks (or the task size). That is actually a good way to control how the parallel code runs, depending on the problem. And you can (with current Julia) emulate the future behavior of `@threads` using `@spawn`.

For example, here I create a function which may have a bad load balancing, because each iteration sleeps for a fraction of a second:

```julia
julia> function test(lapses,ntasks)
           s = zeros(ntasks)
           Threads.@sync for it in 1:ntasks
               Threads.@spawn for i in it:ntasks:length(lapses)
                   sleep(lapses[i])
                   s[it] += lapses[i]
               end
           end
           sum(s)
       end
test (generic function with 1 method)

julia> lapses = [1e-3*i for j in 1:50 for i in 1:4];

julia> sum(lapses) # total sleep time
0.5000000000000001

julia> @btime test($lapses,1) # one task
  756.380 ms (1027 allocations: 31.58 KiB)
0.5000000000000003

julia> @btime test($lapses,4) # ntasks == nthreads
  258.188 ms (1046 allocations: 33.19 KiB)
0.5000000000000003

julia> @btime test($lapses,20) # ntasks >> nthreads
  50.320 ms (1141 allocations: 42.12 KiB)
0.5

```

Thus, if the tasks are very heterogeneous, you can improve balancing by controlling the number of tasks. Ideally, the code can setup some optimal task size depending on the input given, taking into consideration the problem characteristics, avoiding both excessive or insufficient task spawning, both which can be detrimental for performance.

Anyway, this may be offtopic, but who knows.

[Next page](https://discourse.julialang.org/t/random-numbers-and-threads/77364.md?page=2)