Changes in random number generation performance in Julia

By: Blog by Bogumił Kamiński

Re-posted from: https://bkamins.github.io/julialang/2024/02/16/rng.html

Introduction

Today I want to present a small benchmark of random number generation
performance improvements between current Julia release 1.10.1 and
current LTS version 1.6.7.

The idea for the benchmark follows a discussion with a friend who needed
to run some compute intensive Julia code on its LTS version.

The post was written under Julia 1.10.1 and Julia 1.6.7.

The benchmark

Let us start with presenting the benchmark functions

function test_rand1()
    s = 0
    for i in 1:10^9
        s += rand(1:1_000_000)
    end
    return s
end

function test_rand2()
    s = 0.0
    for i in 1:10^9
        s += rand()
    end
    return s
end

They are relatively simple. I wanted to compare the performance of:
(1) integer generation from some range and (2) generation of floating point numbers from [0, 1) interval
as these are two most common scenarios in practice.

Let us see the results. First comes Julia 1.6.7:

julia> @time test_rand1()
  4.949335 seconds (13 allocations: 35.406 KiB)
499993991047124

julia> @time test_rand1()
  4.663646 seconds
499998112691460

julia> @time test_rand2()
  2.175424 seconds
5.000141761909688e8

julia> @time test_rand2()
  2.238839 seconds
4.9999424544883996e8

And now we have Julia 1.10.1:

julia> @time test_rand1()
  2.355028 seconds
500001818410630

julia> @time test_rand1()
  2.287840 seconds
499998082399284

julia> @time test_rand2()
  1.123886 seconds
5.000026226340503e8

julia> @time test_rand2()
  1.117811 seconds
4.9999201274214923e8

So we see that things run roughly two times faster.

Some additional remarks

What is the reason for this difference?
The major point is that between Julia 1.6.7 and Julia 1.10.1
a default random number generator was changed. Let us see
(below I use copy to ensure explicit instantiation of the random number generator object under Julia 1.10.1).

Again, first we test Julia 1.6.7:

julia> using Random

julia> copy(Random.default_rng())
MersenneTwister(0x2fe644ceb724000ca5e5b4409dc3c6ea, (0, 4502994048, 4502993046, 986, 2502992778, 986))

and next we check Julia 1.10.1:

julia> using Random

julia> copy(Random.default_rng())
Xoshiro(0x1273707731737276, 0x187b3d2e82fb1d48, 0x13f9fd1a82642acb, 0xa7dcba727da742e6, 0x3ed2b4d410aa4b31)

So indeed, we see that MersenneTwister was replaced by Xoshiro generator (to be exact Xoshiro256++).

This has one important consequence, apart from random number generation speed that is related to seeding
of the generator. Let us check, Julia 1.6.7:

julia> Random.seed!(1)
MersenneTwister(1)

julia> rand()
0.23603334566204692

vs Julia 1.10.1:

julia> Random.seed!(1)
TaskLocalRNG()

julia> rand()
0.07336635446929285

This means that when you use the default random number generator you should not expect reproducibility of results between these two Julia versions.
This lack of stability is documented as not ensured across Julia versions.

If you need to ensure such reproducibility you can use e.g. the StableRNGs.jl package.

Conclusions

The topic of changes in random number generation in Julia is probably well known to people doing compute intensive simulations.
However, I thought it is worth to present these results for new users, who might be using different versions of Julia to execute the same code
and wonder why the performance or the results themselves are different across them.