Overflow in `randn` using CUDA.jl's native RNG

[danielwe]

Describe the bug

The native CUDA.jl rng regularly produces Infs when sampling standard normal Float32s. This is not expected: for practical purposes, the standard normal distribution has zero measure near floatmax(Float32), and neither CURAND's rng nor the default CPU rng have similar behavior.

To reproduce

The Minimal Working Example (MWE) for this bug:

using CUDA
using Random

println(maximum(randn(CUDA.default_rng(), Float32, 2^26)))         # Likely Inf
println(maximum(randn(CUDA.CURAND.default_rng(), Float32, 2^26)))  # Likely in (5.3, 6.1)
println(maximum(randn(Random.default_rng(), Float32, 2^26)))       # Likely in (5.3, 6.1)

Manifest.toml

CUDA v3.8.5
GPUArrays v8.3.2
GPUCompiler v0.13.14
LLVM v4.9.1

Expected behavior

No overflow, like with CUDA.CURAND.default_rng() and Random.default_rng().

Version info

Details on Julia:

Julia Version 1.7.2
Commit bf53498635 (2022-02-06 15:21 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Xeon(R) Platinum 8280 CPU @ 2.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, cascadelake)

Details on CUDA:

CUDA toolkit 11.6, artifact installation
NVIDIA driver 465.19.1, for CUDA 11.3
CUDA driver 11.3

Libraries: 
- CUBLAS: 11.8.1
- CURAND: 10.2.9
- CUFFT: 10.7.0
- CUSOLVER: 11.3.2
- CUSPARSE: 11.7.1
- CUPTI: 16.0.0
- NVML: 11.0.0+465.19.1
- CUDNN: 8.30.2 (for CUDA 11.5.0)
- CUTENSOR: 1.4.0 (for CUDA 11.5.0)

Toolchain:
- Julia: 1.7.2
- LLVM: 12.0.1
- PTX ISA support: 3.2, 4.0, 4.1, 4.2, 4.3, 5.0, 6.0, 6.1, 6.3, 6.4, 6.5, 7.0
- Device capability support: sm_35, sm_37, sm_50, sm_52, sm_53, sm_60, sm_61, sm_62, sm_70, sm_72, sm_75, sm_80

2 devices:
  0: NVIDIA A100-PCIE-40GB (sm_80, 35.178 GiB / 39.586 GiB available)
  1: NVIDIA A100-PCIE-40GB (sm_80, 26.063 GiB / 39.586 GiB available)

[simsurace]

I think this comes down to the Random.rand(device_rng, Float32) sometimes producing 0f0, which in the Box-Muller transform will give log(0f0) = -Inf32

minimum(rand(CUDA.default_rng(), 2^26)) # Likely 0f0
minimum(rand(CUDA.CURAND.default_rng(), 2^26)) # Likely > 1.5f-10

Is this a known limitation of the Philox2x32RNG?

[maleadt]

The implementation is based on Random123.jl, maybe @sunoru can comment

[simsurace]

Where does the conversion from Int32 to Float32 happen?

[maleadt]

https://github.com/JuliaLang/julia/blob/701468127a335d130e968a8a180021e4287abf3f/stdlib/Random/src/generation.jl#L25-L27= probably?

[sunoru]

Yes. During the conversion, the UInt32 number is firstly truncated to its trailing 23 bits (https://github.com/JuliaLang/julia/blob/701468127a335d130e968a8a180021e4287abf3f/stdlib/Random/src/generation.jl#L124)

But Philox2x32 (as well as other common RNGs) generates uniformly distributed integers, which means a number like the following example would be considered zero when converting to Float32.

rng = Philox2x(UInt32, 0x021242ff, 7)
set_counter!(rng, 0x00a807f2)
println((rand(copy(rng), UInt32), rand(copy(rng), Float32)))
# (0x02000000, 0.0f0)

I'd say this is a known limitation of this conversion approach, but it is very commonly used for its good performance. There is always a possibility around 1/2^23≈1.19e-7 to have 23 trailing zeros.

Random.default_rng() doesn't have this issue because it does not use the Box-Muller transform. I don't know how CURAND does the generation.

[sunoru]

BTW, minimum(rand(Random.default_rng(), Float32, 2^26)) is also likely to be 0.0f0.

[danielwe]

TL;DR:

• The Box-Muller transform should reject and resample if a zero occurs in the "radial" value. Not doing this is the bug.
• Not directly related, but interesting and relates to the smallest nonzero value of rand(rng, ...)
rand(CUDA.default_rng(), ...) only has 23 bits of entropy. For comparison, rand(Random.default_rng(), ...) has 24 bits, and rand(CUDA.CURAND.default_rng(), ...) has even more although it's a little more work to calculate exactly. The "fault" in the case of CUDA.default_rng() lies with the generic transform from random bits to Float32, linked to in the previous comment. See below for an explanation of the transforms used in the other two cases---both of them would be straightforward to implement as the generic case, but I suppose there could be a performance hit.

Full, meandering version

From the previous comment:

https://github.com/JuliaLang/julia/blob/701468127a335d130e968a8a180021e4287abf3f/stdlib/Random/src/generation.jl#L25-L27

If I understand correctly, this transformation generates Float32s by taking 23 random bits as the significand of a Float32 in [1, 2) (an interval that the Float32 representation covers uniformly) and subtracting 1. The smallest nonzero number that can be achieved this way is eps(Float32).

Wikipedia's example implementation of Box-Muller rejects samples below eps for the radial component: https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform#Implementation.

Hence, assuming that floats in Uniform[0, 1) are generated using the transform in generation.jl, simply rejecting exact zeros accomplishes the same as Wikipedia's suggestion to set the threshold at eps. However, if the Uniform[0, 1)-samples are generated in a more sophisticated way that produces nonzero values smaller than eps (see below) I don't see why one would set the threshold there: log never over-/underflows, on the contrary it maps outrageous numbers to sane numbers, so you'd never need to reject anything other than exact zeros (the smallest subnormal representable by Float32 is 2^-149, and log(2f0^-149) ≈ -103).

Now, do the various RNGs use the transform implemented in generation.jl? Not all:

nozero(x) = iszero(x) ? typemax(x) : x
minimum(nozero, rand(CUDA.default_rng(), Float32, 2^25))         # ✅ 1.1920929f-7 == eps(Float32) == 2f0^-23
minimum(nozero, rand(Random.default_rng(), Float32, 2^26))       # ❌ 5.9604645f-8 == eps(0.5f0) == eps(Float32) / 2 == 2f0^-24
minimum(nozero, rand(CUDA.CURAND.default_rng(), Float32, 2^33))  # ❌ 1.1641532f-10 == eps(2f0^-10) == 2f0^-33

Note that the number of samples used in each case is the smallest power of 2 that made the measurement more or less deterministic, so heuristically this exponent should be slightly larger than the number of random bits sampled to generate one Float32. CUDA.default_rng() behaves as surmised; Random.default_rng() seems to use one extra bit and gets twice as close to zero; CUDA.CURAND.default_rng() gets orders of magnitude closer, and based on how frequently the minimum value recurs it seems to make maximal use of 32 random bits.

Another thing to note is that rand(CUDA.CURAND.default_rng(), ...) excludes 0 but includes 1, i.e., it samples in (0, 1], not [0, 1). This is in accordance with the documentation: https://docs.nvidia.com/cuda/curand/host-api-overview.html#generation-functions.

Random.default_rng() is a Xoshiro instance (julia 1.7.2), which implements it's own transform from bits to Float32 here: https://github.com/JuliaLang/julia/blob/701468127a335d130e968a8a180021e4287abf3f/stdlib/Random/src/Xoshiro.jl#L207-L208. This interprets 24 random bits as an integer between 0 and 2^24 - 1, and multiplies this by 2^-24 to produce evenly spaced values from 0 to 1 - 2^-24. This uses one more bit than the generic transform in generation.jl, which makes sense since the density of Float32s in [0, 1) is at least twice as high as in [1, 2) (specifically, it's twice as high in [0.5, 1) and even higher below 0.5). Another way to accomplish the same would be to take the value returned by the transform in generation.jl and flip a coin to decide whether or not to add 2f0^-24.

Meanwhile, CUDA.CURAND.default_rng() makes even more use of the higher resolution at small values. After poking around a little it seems to map 32 random bits one-to-one onto evenly spaced real numbers in the interval (2^-33, 1 - 2^-33). At the lower end these values can be represented exactly in Float32, but for values higher than 2^-8 they must be rounded, so the map becomes many-to-one, and in particular all values larger than 1 - 2^-25 round to 1. I think this can be done using a transform like the one used by Xoshiro, but without discarding bits in advance and instead accepting truncation; something like (Float32(rand(rng, UInt32), RoundUp) + 0.5f0) * 2f0^-32.

[danielwe]

Here's the explanation for why the generic implementation only uses 23 bits (52 for Float64), while Xoshiro uses 24 (53): JuliaLang/julia#33222 (comment)

We still use the 52-bit version for our older RNGs so that they give the same values as before, but for the new xoshiro generator we're using shift-by-11-and-multiply currently.

Is it important that CUDA.default_rng() maintains legacy behavior here?

[simsurace]

The Box-Muller transform should reject and resample if a zero occurs in the "radial" value. Not doing this is the bug.

I believe this sort of conditional resampling may create issues with thread synchronization on the GPU, at least that was something to consider last time I checked, so please correct me if I'm wrong about this. I do use a non-deterministic number of random draws in BinomialGPU.jl, but this did degrade performance when the synchronizing RNG was introduced, even though this degradation was compensated by other improvements. I haven't checked in a while whether this is still an issue.

The following assumes that the thread synchronization still makes the conditional generation of an additional random uniform number bad for performance.

Another way to accomplish the same would be to take the value returned by the transform in generation.jl and flip a coin to decide whether or not to add 2f0^-24.

Because this approach deterministically generates exactly one random number in all threads, it should only slightly degrade performance and not lead to any synchronization issues.

Short of modifying the Int32 to Float32 conversion, this seems the easier route to take, at least in the short term.

[danielwe]

It seems plausible that resampling would be problematic on GPU. My comment referred to what the Box-Muller transform requires in general, I wasn't thinking about GPU limitations.

The conditional addition of 2f0^-24 was just an illustration of what the extra bit of entropy in the 24-bit vs. 23-bit conversion amounts to, it doesn't do anything to solve the problem of infinities other than making them twice as rare. (Sorry for writing such a long and confusing post that this wasn't clear.) However, now that you say it, always adding 2f0^-24 would remove the zeros in rand() while maintaining a uniform distribution in (0, 1). That would solve the problem, so maybe that's the way to go?

[danielwe]

Btw. the corresponding numbers for Float64 and Float16 would be 2.0^-53 and Float16(2)^-11, respectively.

[maleadt]

The Box-Muller transform should reject and resample if a zero occurs in the "radial" value. Not doing this is the bug.

I believe this sort of conditional resampling may create issues with thread synchronization on the GPU, at least that was something to consider last time I checked, so please correct me if I'm wrong about this. I do use a non-deterministic number of random draws in BinomialGPU.jl, but this did degrade performance when the synchronizing RNG was introduced, even though this degradation was compensated by other improvements. I haven't checked in a while whether this is still an issue.

Resampling shouldn't be problematic. I'll create a PR.