Random number iterator?

General Usage Performance
,
1

Is it possible to get a random number iterator that doesn’t allocate? I was hoping I could do something like this:

julia> using Random

julia> sum(Iterators.take(MersenneTwister(0), 4))
ERROR: MethodError: no method matching iterate(::MersenneTwister)

Of course I could write a generator,

sum(rand() for _ in 1:4)

but n calls to rand() is usually slower than rand(n):

julia> @btime sum(rand() for _ in 1:100_000);
  146.216 μs (0 allocations: 0 bytes)

julia> @btime sum(rand(100_000));
  105.032 μs (2 allocations: 781.30 KiB)
2

That’s just how it works right?.. Trading space for speed. The iterator will necessarily call rand() to minimize allocation… You can’t not allocate but also want batch rand(1000)…

1 Like
3

Maybe I misunderstand, but why not

julia> using BenchmarkTools, Random

julia> a = zeros(1_000);

julia> @btime rand!($a);
  1.667 μs (0 allocations: 0 bytes)
4

The extent of this would depend on which version you are in. In Julia 1.6, the dicrepancy would mostly be fixed by passing an explicit rng to rand in the generator version. But you your timings make me guess you are on 1.7. There I don’t think you can easily beat the array version, because it uses simd, unlike the scalar version. An idea I would want to explore in a package is wrapping Xoshiro in a way that even the scalar version of rand uses simd, by having an internal cache (like what MersenneTwister does).

EDIT: but to your question about having a rand iterator, you can check out Rand() from the RandomExtensions package. It won’t help performance in this instance, but as it hanldes the Sampler thing in the same way as in arrays, it can occasionally help speed-up things where a “pre-computation” can be shared between multiple calls of rand.

3 Likes
5

I think we could make a really fast operator version that generates 32 at a time and feeds them to you.

7

Thanks! Rand() from RandomExtensions is indeed the kind of thing I was looking for. Here it is in action:

julia> using Random, RandomExtensions, Base.Iterators, BenchmarkTools

julia> rng = MersenneTwister(0);

julia> itr = Rand(rng)
Rand{MersenneTwister, Random.SamplerTrivial{Random.CloseOpen01{Float64}, Float64}}(MersenneTwister(0), Random.SamplerTrivial{Random.CloseOpen01{Float64}, Float64}(Random.CloseOpen01{Float64}()))

julia> @btime sum(take($itr, 100_000));
  169.361 μs (0 allocations: 0 bytes)

julia> @btime sum(rand() for _ in 1:100_000);
  146.219 μs (0 allocations: 0 bytes)

julia> @btime sum(rand(100_000));
  105.492 μs (2 allocations: 781.30 KiB)

I see that RandomExtensions is experimental, but it looks cool. I don’t know anything about pseudo-random number generation algorithms, so I wasn’t sure if there was some conceptual reason why one couldn’t make a random number iterator. Whether that iterator can take advantage of SIMD is a separate question, I suppose.

1 Like