I am trying to find the Expected Value of maximum of two U[0, 1] RVs. Why is a call to rand(2) so much worse than two to rand(), rand()?
julia> function h(N)
s = 0.
for _ in 1:N
a, b = rand(2)
s += (a>b) ? a : b
end
s/N
end
h (generic function with 1 method)
julia> function g(N)
s = 0.
for _ in 1:N
a, b = rand(), rand()
s += (a>b) ? a : b
end
s/N
end
g (generic function with 1 method)
julia> using BenchmarkTools
julia> @btime h(10^6)
47.601 ms (1000000 allocations: 91.55 MiB)
0.6665339890735786
julia> @btime g(10^6)
8.913 ms (0 allocations: 0 bytes)
0.6667565533140116
Interestingly, the non-allocating rand!() function somehow still leads to a lot of allocations in this context
julia> using Random
julia> function j(N)
s = 0.
buf = zeros(2)
for _ in 1:N
a, b = rand!(buf)
s += (a>b) ? a : b
end
s/N
end
j (generic function with 1 method)
julia> using BenchmarkTools
julia> @btime j(10^6)
149.425 ms (5000003 allocations: 106.81 MiB)
0.6666123720042385
That’s because rand!(buf) has to access the global default RNG:
julia> function j(N, rng)
s = 0.
buf = zeros(2)
for _ in 1:N
a, b = rand!(rng, buf)
s += (a>b) ? a : b
end
s/N
end
j (generic function with 2 methods)
julia> @btime j(10^6, $(Random.default_rng()))
16.276 ms (1 allocation: 96 bytes)
0.6669343726647938
To be more explicit, compare the two fresh sessions
julia> using BenchmarkTools, Random
julia> function j(N)
s = 0.
buf = zeros(2)
for _ in 1:N
a, b = rand!(buf)
s += (a>b) ? a : b
end
s/N
end
j (generic function with 1 method)
julia> @btime j(10^6)
16.949 ms (1 allocation: 96 bytes)
0.6665871100597152
and
julia> using BenchmarkTools
julia> function j(N)
s = 0.
buf = zeros(2)
for _ in 1:N
a, b = rand!(buf)
s += (a>b) ? a : b
end
s/N
end
j (generic function with 1 method)
julia> j(1)
ERROR: UndefVarError: rand! not defined
Stacktrace:
[1] j(N::Int64)
@ Main ./REPL[2]:5
[2] top-level scope
@ REPL[3]:1
julia> using Random
julia> @btime j(10^6)
121.164 ms (5000003 allocations: 106.81 MiB)
0.6667798259191253
You’re absolutely correct, that’s what I did. I have never heard of that issue (and I’ve been using julia since at least 0.6)… very interesting, thank you
I’m a bit confused… Is this expected? Shouldn’t the function j get recompiled automatically? Do I have to worry about the order of imports in my code to get good performance?
You do have to worry about order of imports in that running j before using Random results in an error.
If this happen, then a version of j has been compiled that treats rand! as an unknown/dynamic global variable. Because it wasn’t defined, it throws an error.
Once you define it/bring it into the global scope, the function doesn’t recompile. It now just works.
To recompile, j would need to have a backedge pointing to rand!. But it can’t have a backedge to rand! if the problem in the first place was that there was no rand!.
Anyway, as long as you bring rand into scope before actually compiling and running j, you’ll be fine.
julia> using BenchmarkTools
julia> function j(N)
s = 0.
buf = zeros(2)
for _ in 1:N
a, b = rand!(buf)
s += (a>b) ? a : b
end
s/N
end
j (generic function with 1 method)
julia> using Random
julia> @btime j(10^6)
6.222 ms (1 allocation: 80 bytes)
0.6665835627414475
In other words, as long as you aren’t running your code in an order that causes you to get errors – which you should notice and hopefully already be avoiding anyway – no issue an no need to change what you’re doing.
PS. My benchmark was on Julia master. The new RNG seems to help this benchmark a lot vs 1.6 and earlier (which I assume others were using).
