How can we generate random numbers in parallel in Julia?

Question

I am writing a parallelized julia code for monte carlo simulations. This requires me to generate random numbers in parallel on different cores. In a simple test code on my workstation, I tried to generate random numbers on 4 cores and got the following results:

julia -p 4

julia> @everywhere using Random

julia> @everywhere x = randn(1)

julia> remotecall_fetch(println,1,x[1])
-1.9348951407543997

julia> remotecall_fetch(println,2,x[1])
      From worker 2:    -1.9348951407543997

julia> remotecall_fetch(println,3,x[1])
      From worker 3:    -1.9348951407543997

julia> remotecall_fetch(println,4,x[1])
      From worker 4:    -1.9348951407543997

I do not understand why the numbers fetched from different processes give exactly the same result. I am not sure what the mistake is. My understanding is that using the @everywhere macro lets you run the same piece of code on all the processes in parallel. I am currently julia 1.6.0 on my computer. Thank you.

UPDATES: Thank you for the responses. Basically, what I am looking for is an assignment statement like x = y, where x and y are both local to a worker processes. I tried something like this:

julia -p 4

@sync @distributed for i = 1:2
       x = randn(1)
       println(x)
       end
      From worker 3:    [0.4451131733445428]
      From worker 2:    [-0.4875627629008678]
Task (done) @0x00007f1d92037340

julia> remotecall_fetch(println,2,x)
ERROR: UndefVarError: x not defined
Stacktrace:
 [1] top-level scope
   @ REPL[23]:1

This seems to generate the random numbers independently on each process. However, I do not know how to access the variable x anymore. I tried remotecall_fetch(println, 2,x) but the variable x seems to be not defined on the worker processes. This has been super confusing.

I wish there was nice flowchart or good documentation explaining the scope of variables and expressions in Julia during parallel computations.

Answer 1

remotecall_fetch sends x[1] for evaluation from your local process (having id 1). You can check it by running the following code:

# julia -p 4

julia> @everywhere x = myid() # make sure x holds a worker number

julia> remotecall_fetch(println, 4, x) # x passed from worker 1 (local machine) to println
      From worker 4:    1

julia> @sync @everywhere println(x) # x is evaluated on worker
1
      From worker 3:    3
      From worker 2:    2
      From worker 4:    4
      From worker 5:    5

julia> @sync @everywhere println($x) # x interpolated from local machine
1
      From worker 4:    1
      From worker 5:    1
      From worker 3:    1
      From worker 2:    1

---

Regarding random number generation on remote machines you should make sure to create independent random number streams on each of them. The simplest method that is good enough for most cases is to just use Random.seed! function with different seeds on different workers. If you want to be extra careful then use Future.randjump to ensure that random number generators on worker processes do not have an overlap.

Answer 2

The best way is to have a single random state and slice it to pieces where each worker has a single piece. This can be done as:

using Distributed
addprocs(4)
@everywhere import Future, Random
@everywhere const rng = Future.randjump(Random.MersenneTwister(0), myid()*big(10)^20)

Now you have at each worker a worker local yet in Julia sense global rng variable.

In this example I have used 0 ad the random number seed. The randjump size of big(10)^20 is recommended in the manual since in Julia this step has already precomputed values.

To make use of such rng you could define a function such as:

@everywhere getr(rng=rng) = rand(rng, 5)

which could be called as

fetch(@spawnat 2 getr())

Basically rng since it is global, it should be passed as the outermost parameter to whatever you are calling on the remote worker or defined as const as mentioned in the comments.

Answer 3

You say on top of your question "on different cores". This is the solution that I use for functions that involves random number generations on a thread based parallelism using the Threads.@thread macro. It has the advantage that I can choose, based on the needs when I run the code, the level of stochasticity that I need:

using Random, Test, Statistics

FIXEDRNG = MersenneTwister(123)

println("** Testing generateParallelRngs()...")
x = rand(copy(FIXEDRNG),100)

function innerExpensiveFunction(bootstrappedx; rng=Random.GLOBAL_RNG)
     sum(bootstrappedx .* rand(rng) ./ 0.5)
end
function outerFunction(x;rng = Random.GLOBAL_RNG)
    masterSeed = rand(rng,100:9999999999999) 
    rngs       = [deepcopy(rng) for i in 1:Threads.nthreads()]  # make new copy instances
    results    = Array{Float64,1}(undef,30)
    Threads.@threads for i in 1:30
        tsrng         = rngs[Threads.threadid()]    # Thread safe random number generator: one RNG per thread
        Random.seed!(tsrng,masterSeed+i*10)         # But the seeding depends on the i of the loop not the thread: we get same results indipendently of the number of threads
        toSample      = rand(tsrng, 1:100,100)
        bootstrappedx = x[toSample]
        innerResult   = innerExpensiveFunction(bootstrappedx, rng=tsrng)
        results[i]    = innerResult
    end
    overallResult = mean(results)
    return overallResult
end


# Different sequences..
@test outerFunction(x) != outerFunction(x)

# Different values, but same sequence
mainRng = copy(FIXEDRNG)
a = outerFunction(x, rng=mainRng)
b = outerFunction(x, rng=mainRng)

mainRng = copy(FIXEDRNG)
A = outerFunction(x, rng=mainRng)
B = outerFunction(x, rng=mainRng)
@test a != b && a == A && b == B


# Same value at each call
a = outerFunction(x,rng=copy(FIXEDRNG))
b = outerFunction(x,rng=copy(FIXEDRNG))
@test a == b