Correct bias in random float sampling.

[tpapp]

Fixes #33222.

Also changed some hardcoded test values accordingly.

[mschauer]

I think this should have a news item.

[tpapp]

I added a new sampler OpenOpen01 and made it default. DSFMT turns out to support this natively.

I tried to keep things consistent, but since Random has a ton of clever micro-optimizations, a thorough review would be appreciated. @rfourquet, can you please review?

[andreasnoack]

@nanosoldier runbenchmarks("random", vs = ":master")

[perrutquist]

This can make z equal prevfloat(Float32(2)), so that Float16(z - prevfloat(Float16(1))) rounds to 1.

[mschauer]

Hm, I am surprised, in
#33222 (comment)
I checked the extrema of the new generator. Or we should really use a hard-coded constant Float16(reinterpret(Float32,n(UInt32(i) << 13) | 0x3f800000) - 0.9995117f0) as I proposed.

[perrutquist]

Sorry. I didn't think this fully through. You're right, it can't become 1.

[tpapp]

@perrutquist, thanks for bringing this up anyway, I am always uneasy operating around eps() precision, and @mschauer, thanks for checking this.

[nanosoldier]

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @ararslan

[tpapp]

It would be great if someone could provide guidance about whether to address anything from the benchmarks. I don't have the perspective to say what's just acceptable noise and what I should investigate.

[rfourquet]

a thorough review would be appreciated

Sure I will do, but shouldn't a decision be made first? In the discourse thread it seemed that the preference was for keeping [0,1) by default instead of (0,1).

[tpapp]

That's not how I read that discussion. People (including me) argued that [0,1) was practical and other languages are doing it, so it is "OK". I am not sure anyone is strongly attached to [0,1), I think it was just what other languages happened to be doing, and can be done relatively efficiently.

This issue is kind of orthogonal: it is about fixing a bias. Incidentally it makes random numbers in (0,1), which may be irrelevant or important to people. I got the impression that ruling out 0s is important to some as it simplifies code, but getting those 0s with 2^(-52) probability is not something that is considered crucial.

[andreasnoack]

I think it's fine to provide a (0,1) version of the sampler but I don't see the need for the changing the default. Floats are not reals and the bias is negligible for any practical purposes. At least, I don't think it's worth paying any slowdown for this. (If you start considering Float16 then all kinds of issues will show up. E.g. you'd get many repeated values for even very small sample sizes).

[rfourquet]

IIRC, one mentioned reason to favor [0, 1) is that it's easy to rule out 0 whereas the reverse is not true.
I personally have no opinion on the matter, the current [0,1) is fine for most people and the small bias could be documented. Both options can be provided anyway (but of course only one default). The issue you opened is probably the better place to discuss, but it's a (slightly) breaking change (in particular for Float16 and Float32 where it's easy to get 0) so an explicit decision must be made.

[tpapp]

I am not yet sure that this implies any slowdown (that is not fixable), so maybe it would be better to treat the issue of correctness separately from that.

As I said in the other thread, if there is a bias we can fix, we should just fix it, even if it is "small".

Regarding breaking changes: is there a commitment to the same random stream from a specific seed with the same RNG? I was under the impression that there wasn't.

[rfourquet]

Regarding breaking changes: is there a commitment to the same random stream from a specific seed with the same RNG?

I was refering to the breaking behavior that 0 won't result from a call to rand() or rand(Float16), some applications may rely on the possiblity of getting 0.

[perrutquist]

From the documentation of Random.seed!: " rng will give a reproducible sequence of numbers if and only if a seed is provided"

There's some wiggle room of course, but I would at least call this a "minor change" that should not be in a patch release.

As for code relying on rand(Float16) returning zero, I doubt that there's any in the wild, but if there is, code like

rand(Float16) == 0 # this will happen with probability 0.000977 (approximately)

would have to be rewritten as

rand(Float16) < eps(Float16) # ≈ 0.000977

which gives the same result with or without the bias correction. (And I would argue that this code is also easier to read, as the probability is stated more explicitly.)

[mschauer]

It would be great if someone could provide guidance about whether to address anything from the benchmarks. I don't have the perspective to say what's just acceptable noise and what I should investigate.

It appears that for example

Float16(reinterpret(Float32, (UInt32(i) << 13)  | 0x3f800000) - prevfloat(Float16(1)))

is slower than

Float16(reinterpret(Float32, (UInt32(i) << 13)  | 0x3f800000) -  0.9995117f

As I proposed hardcoding the constants two times before, let me just propose this once again.

[tpapp]

Thanks for pointing this out again, I will experiment and modify the PR accordingly.

[mschauer]

Let's hope this is actually the problem.

[perrutquist]

With an explicit conversion Float32(prevfloat(Float16(1)))) is as fast as hard-coded 0.9995117f0. (Not sure why this should be necessary.)

[perrutquist]

Another thing that might be a (very slight) performance issue is that dSFMT seems to be producing one more bit of entropy for OpenOpen01 than for CloseOpen01. (Haven't verified this, but it looks that way from the fact that only half of the rands changed when filling an array.) Mixing in this bit probably costs a couple of clock cycles per SIMD vector.

[tpapp]

@perrutquist: are you sure about this? I thought that was just rounding.

For the "native" scalar samplers, here is a script that benchmarks everything, comparing the proposed samplers for OpenOpen01, the one with constant folding, and the CloseOpen01 sampler (which is unchanged) for comparison:

### NOTE run on https://github.com/JuliaLang/julia/pull/33251

using Random: AbstractRNG, SamplerTrivial, Sampler, GLOBAL_RNG, OpenOpen01, UInt10, UInt23,
    CloseOpen12, OpenOpen01_64, CloseOpen01
using BenchmarkTools

RNG = GLOBAL_RNG
S16 = Sampler(RNG, Float16, Val(1))
S32 = Sampler(RNG, Float32, Val(1))
S64 = Sampler(RNG, Float64, Val(1))
C16 = SamplerTrivial(CloseOpen01{Float16}())
C32 = SamplerTrivial(CloseOpen01{Float32}())
C64 = SamplerTrivial(CloseOpen01{Float64}())

###
### Float16
###

function rand_const(r::AbstractRNG, ::SamplerTrivial{OpenOpen01{Float16}})
    z = reinterpret(Float32, (rand(r, UInt10(UInt32)) << 13) | 0x3f800000)
    Float16(z - 0.9995117f0)
end

function rand_conv(r::AbstractRNG, ::SamplerTrivial{OpenOpen01{Float16}})
    z = reinterpret(Float32, (rand(r, UInt10(UInt32)) << 13) | 0x3f800000)
    Float16(z - Float32(prevfloat(Float16(1))))
end

###
### Float32
###

function rand_const(r::AbstractRNG, ::SamplerTrivial{OpenOpen01{Float32}})
    reinterpret(Float32, rand(r, UInt23()) | 0x3f800000) - 0.99999994f0
end

###
### Float64
###

function rand_const(r::AbstractRNG, ::SamplerTrivial{OpenOpen01_64})
    rand(r, CloseOpen12()) - 0.9999999999999999
end

@btime rand($RNG, $S16);
@btime rand_const($RNG, $S16);
@btime rand_conv($RNG, $S16);
@btime rand($RNG, $C16);

@btime rand($RNG, $S32);
@btime rand_const($RNG, $S32);
@btime rand($RNG, $C32);

@btime rand($RNG, $S64);
@btime rand_const($RNG, $S64);
@btime rand($RNG, $C64);

With these results:

julia> @btime rand($RNG, $S16);
  12.395 ns (0 allocations: 0 bytes)

julia> @btime rand_const($RNG, $S16);
  9.083 ns (0 allocations: 0 bytes)

julia> @btime rand_conv($RNG, $S16);
  9.065 ns (0 allocations: 0 bytes)

julia> @btime rand($RNG, $C16);
  9.066 ns (0 allocations: 0 bytes)

julia> @btime rand($RNG, $S32);
  7.159 ns (0 allocations: 0 bytes)

julia> @btime rand_const($RNG, $S32);
  7.155 ns (0 allocations: 0 bytes)

julia> @btime rand($RNG, $C32);
  7.145 ns (0 allocations: 0 bytes)

julia> @btime rand($RNG, $S64);
  6.419 ns (0 allocations: 0 bytes)

julia> @btime rand_const($RNG, $S64);
  6.372 ns (0 allocations: 0 bytes)

julia> @btime rand($RNG, $C64);
  6.478 ns (0 allocations: 0 bytes)

The only regression is for the Float16 (the rest is random fluctuation, I confirmed this with multiple runs). Cf

Body::Float16
1 ─ %1  = Random.UInt10(Random.UInt32)::Core.Compiler.Const(UInt10{UInt32}(), false)
│   %2  = Random.rand(r, %1)::UInt32
│   %3  = (%2 << 13)::UInt32
│   %4  = (%3 | 0x3f800000)::UInt32
│         (z = Random.reinterpret(Random.Float32, %4))
│   %6  = z::Float32
│   %7  = Random.Float16(1)::Core.Compiler.Const(Float16(1.0), false)
│   %8  = Random.prevfloat(%7)::Core.Compiler.Const(Float16(0.9995), false)
│   %9  = (%6 - %8)::Float32
│   %10 = Random.Float16(%9)::Float16
└──       return %10

julia> @code_warntype rand(RNG, S32)
Variables
  #self#::Core.Compiler.Const(rand, false)
  r::Core.Compiler.Const(Random._GLOBAL_RNG(), false)
  #unused#::Core.Compiler.Const(SamplerTrivial{OpenOpen01{Float32},Float32}(OpenOpen01{Float32}()), false)

Body::Float32
1 ─ %1 = Random.UInt23()::Core.Compiler.Const(UInt23{UInt32}(), false)
│   %2 = Random.rand(r, %1)::UInt32
│   %3 = (%2 | 0x3f800000)::UInt32
│   %4 = Random.reinterpret(Random.Float32, %3)::Float32
│   %5 = Random.prevfloat(1.0f0)::Core.Compiler.Const(0.99999994f0, false)
│   %6 = (%4 - %5)::Float32
└──      return %6

which shows that constant folding is not happening for Float16 for some reason. I am not sure whether to open an issue about this.

In any case, I will fix the Float16 case. Thanks for the help.

[perrutquist]

No. I'm not sure. Looking at the dSFMT code, it seems it sets the least significant bit to 1 before subtracting 1.0, rather than directly subtracting prevfloat(1.0) as we do. This would waste one bit of entropy and one clock cycle.

However, looking at the random number vector in the test code, there seems to be randomness in the eps()/2 bit, so maybe I'm looking at the wrong piece of code.

I think this is the relevant part of the dSFMT code:

#if defined(HAVE_SSE2)
...
inline static void convert_o0o1(w128_t *w) {
    w->si = _mm_or_si128(w->si, sse2_int_one);
    w->sd = _mm_add_pd(w->sd, sse2_double_m_one);
}
...
#else
...
inline static void convert_o0o1(w128_t *w) {
    w->u[0] |= 1;
    w->u[1] |= 1;
    w->d[0] -= 1.0;
    w->d[1] -= 1.0;
}
...
#endif

[mschauer]

@nanosoldier runbenchmarks("random", vs = ":master")

[ararslan]

Triggering Nanosoldier requires commit access to Julia.

@nanosoldier runbenchmarks("random", vs=":master")

[nanosoldier]

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @ararslan

[tpapp]

I am wondering if the benchmark improvements/regressions are just noise: in some cases, nothing was changed, in other cases I compared directly and there should be no performance hit.

@rfourquet raised an important point above: it is still not decided whether a change like this would be accepted. It would be great if triage could discuss this.

[stevengj]

I'm skeptical of this. Random floats are supposed to produce samples uniformly in [0,1), but of course that is only up to roundoff error, and any practical application will surely have many other accumulated roundoff errors that swamp a bias of eps(T)/2 And even in Float32 arithmetic you would need an impractical number of samples (~10¹⁴) to conceivably statistically detect the bias.

It doesn't seem reasonable to even slightly slow down all random-number generation—or just change the statistical distribution, for that matter—to correct an effectively invisible bias. Especially since virtually any nontrivial calculation using rand (even just rand() < 0.3) will have a similar "bias" that cannot be eliminated.

[tpapp]

FWIW I don't think there is a slowdown, those numbers are just benchmarking noise. This of course should be investigated if triage is open to the idea, but there should be literally no cost to fixing this, as we are just subtracting a different constant (and dsfmt also supports this natively).

Also, as @mschauer demonstrated in the discussion of #33222, the bias is actually detectable for Float16. I agree that for most calculations with Float64 this should not matter, but why no do it right if we can?

Incidentally, a possible benefit is that algorithms that need to branch for 0 as a special case will be able to just omit that. Cf https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505, https://discourse.julialang.org/t/random-number-in-0-1/16009, and similar topics.

[JeffBezanson]

This makes triage nervous. It's ok to add the open-open sampler and make it available, but we don't want to change the default without a very compelling argument. We also discussed that closed-closed is very defensible, since it's what you'd get if you sampled ideal real numbers and then rounded (some numbers would round to exactly 0 or 1). Anyway, we should not rush to change this.

[tpapp]

I understand. I will revert the changes that make the new sampler default, and then this PR will just add OpenOpen01.

[ViralBShah]

Closing based on this comment: #33222 (comment)