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

See The need for rand speed | Blog by Bogumił Kamiński

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:
 [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
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 :slight_smile:

2 Likes