How to change the default RNG for a given script?

General Usage
,
1

I know I can pass RNG to many (all?) stochastic functions, and that StableRNG is slow, but is there a way I can set at the beginning of my model something like:

random_seed = 123          # from a config file
Random.GLOBAL_RNG = StableRNG(123)`

?

This doesn’t work as Random.GLOBAL_RNG is a constant.

1 Like
2

No, unless you want to edit Random.jl, but you can seed!(123), or you just have to pass the rng object down the call chain.

2 Likes
3

And if it’s too tedious to pass down the call chain, you could make your RNG a global const and just access it at the call sites. Although if it’s not then passing it down is a tad nicer.

1 Like
4

And faster:

1 Like
5

One could use a ScopedValue for the best of both worlds. It would be nice if the default rng was a ScopedValue but that’s probably infeasible as the default rng is very special and gets seeded on task creation for reproducible randomness with multi threading.

1 Like
6

I was curious about these statements:

So I tested on Julia v1.11.6

@btime pi_default($N);
 2.505 s (0 allocations: 0 bytes)

@btime pi_global($N);
 2.524 s (0 allocations: 0 bytes)

@btime pi_local($N, $rng);
 2.511 s (0 allocations: 0 bytes)

@btime pi_local($N, $stablerng);
 1.394 s (0 allocations: 0 bytes)
Code 
using Random
using StableRNGs
using BenchmarkTools

const MY_RNG = Xoshiro()

function pi_default(n::Int)
    s = 0
    for _ in 1:n
        s += rand()^2 + rand()^2 < 1
    end
    return 4 * s / n
end

function pi_global(n::Int)
    s = 0
    for _ in 1:n
        s += rand(MY_RNG)^2 + rand(MY_RNG)^2 < 1
    end
    return 4 * s / n
end

function pi_local(n::Int, rng::AbstractRNG)
    s = 0
    for _ in 1:n
        s += rand(rng)^2 + rand(rng)^2 < 1
    end
    return 4 * s / n
end

N = 10^9
rng = Xoshiro()
stablerng = StableRNG(123)

@btime pi_default($N);
@btime pi_global($N);
@btime pi_local($N, $rng);
@btime pi_local($N, $stablerng);
2 Likes
7

Yes, for me StableRNGs are slightly faster in that case.
However, the story is different when generating large-ish arrays of random numbers, as it seems the default rng has a good SIMD implementation for that case.

julia> stab = StableRNG(123);

julia> @btime rand(N) setup=(N=10^5);
  49.200 μs (3 allocations: 781.32 KiB)

julia> @btime rand($stab, N) setup=(N=10^5);
  106.700 μs (3 allocations: 781.32 KiB)

julia> @btime rand(N) setup=(N=10^7);
  13.358 ms (3 allocations: 76.29 MiB)

julia> @btime rand($stab, N) setup=(N=10^7);
  18.210 ms (3 allocations: 76.29 MiB)
2 Likes
8

Separate from the Random vs StableRNGs performance comparison:

One factor making RNGs slow when generating just a few random numbers in a call is that they’re mutable. Reading/writing the RNG state to RAM incurs some overhead for every rand call (not every number that’s generated, which is why only small calls suffer). An immutable RNG state (implemented like an iterator, for example, where the state is carried as a local variable) would not need to do this, so would be faster for small calls (orthogonal to the SIMD question).

It’s been on my list (long and probably never to be completed) to experiment with adding immutable RNGs to Random, but it’s a pretty big refactor to integrate them properly and there would be some debate over the interface.

2 Likes