random sampling from an (abstract)array is slow

[CarloLucibello]

It is three times slower than this alternative implementation

julia> myrand(v) = (i = ceil(Int,rand()*length(v));  v[i])

julia> @benchmark rand(1:100)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     36.048 ns (0.00% GC)
  median time:      36.996 ns (0.00% GC)
  mean time:        37.104 ns (0.00% GC)
  maximum time:     66.210 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     993
  time tolerance:   5.00%
  memory tolerance: 1.00%

julia> @benchmark myrand(1:100)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     11.824 ns (0.00% GC)
  median time:      12.332 ns (0.00% GC)
  mean time:        12.286 ns (0.00% GC)
  maximum time:     35.627 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     999
  time tolerance:   5.00%
  memory tolerance: 1.00%

[StefanKarpinski]

cc @rfourquet, who made sure this stuff was all very fast at one point. We should add BaseBenchmarks for this so that it can't regress again. Thanks for reporting!

[martinholters]

Probably ceil(Int,rand()*length(v)) is a bad way to get a uniform distribution between 1 and length(v)? It might be faster than what we have, but might be biased?

[CarloLucibello]

Probably ceil(Int,rand()*length(v)) is a bad way to get a uniform distribution between 1 and length(v)? It might be faster than what we have, but might be biased?

Yes, I have to delve deeper into this, but it could be slightly unbiased, with a negligible bias for small length(v) wich gets higher for bigger vectors

[StefanKarpinski]

That version is also slightly buggy since there's a small chance of getting a zero index and then getting a bounds error upon indexing.

[bkamins]

This issue is recurring many times on Discourse (https://discourse.julialang.org/t/rand-1-10-vs-int-round-10-rand/14339/9), so I put some benchmarks.

Here is the code that shows where we have a problem:

# tweak function to make rand(range) return number of loop iterations
function Random.rand(rng::AbstractRNG, sp::Random.LessThan)
    i = 0
    while true
        i += 1
        x = rand(rng, sp.s)
        x <= sp.sup && return i-1
    end
end

And now running it for small ranges gives:

julia> [mean(rand(1:j) for i in 1:10^6) for j in 1:16]
16-element Array{Float64,1}:
 1.0
 1.0
 1.333166
 1.0
 1.600809
 1.333687
 1.142155
 1.0
 1.777679
 1.599978
 1.455112
 1.333272
 1.231023
 1.143564
 1.067084
 1.0

I guess it cannot be helped without using div/mod which is more expensive than generating an additional pseudorandom number (@rfourquet did I get your idea in this design right?).

And now the issue is that for a smal n the formula 1+floor(Int,n*rand()) has a small bias. For example for n=10^6 it is around 10^-10. Now, in order to detect this difference statistically the number of required samples is astronomical, so:

Maybe we want to allow for some small bias in trade off for speed for small n? If not in rand then maybe in something like fastrand? OTOH maybe it is better to do such things in packages. Not sure.

[rfourquet]

@rfourquet did I get your idea in this design right?

Indeed, that was the idea of the change in #27560, as generation with MersenneTwister is fast enough that this approach is usually more performant than using div/mod. That said, the default Sampler for other RNGs is still using div/mod (SamplerRangeInt).

[rfourquet]

I think that rather than fastrand, we could create an object Biased (or Fast), so that a call would look like rand(Biased(1:10)). Like that the rand infrastructure is still available (e.g. we don't have to reimplement array generation for fastrand).

[bkamins]

Nice idea. I think that Fast would be nice, and we could avoid exporting it so user would call it using Random.Fast (like with seed!) as I guess this operation will not be so common.