# Why is a, b = rand(2) so bad?

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

``````
2 Likes

Presumably because `rand(2)` has to allocate an array.

6 Likes

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

``````

Only if you try to run it before `using Random` and it gets compiled without full inference information.

2 Likes

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
``````
3 Likes

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:
 j(N::Int64)
@ Main ./REPL:5
 top-level scope
@ REPL:1

julia> using Random

julia> @btime j(10^6)
121.164 ms (5000003 allocations: 106.81 MiB)
0.6667798259191253
``````
11 Likes

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

5 Likes

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?

7 Likes

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).

8 Likes

Thanks, that’s reassuring 2 Likes