I’m now fully convinced: there does not exist a principled open-open (0,1) floating point rand.
Concretely, rand() will never be able to pass out random variables from any real number interval, despite how badly we wish it would, no matter how perfect our implementation, no matter how high the precision of the floating point type. It gives you a rounded random variable. You can either round it up (leading to a right-closed interval) or down (left-closed).
Now, one could choose to compromise on the principle in favor of pragmatism. ← This is where our discussion goes in circles. Obviously opinions will vary, both on the appetite for sullying the principled approach in the first place and on one’s favored balance point.
I think we are going in circles because the discussion is about goals that the article linked on Hacker News specifies, which sound nice, but are in practice irrelevant. Specifically:
With that exposition done, I can now present my criteria for assessing the quality of functions that claim to generate uniformly-distributed IEEE-754 double-precision floating-point values over (0, 1] or [0, 1] or [0, 1):
What is the probability of seeing exactly zero? The rounding basin of 0 is only 2^-1075 wide, so the probability should be as close as possible to this. Note that probability zero (i.e. impossible) is extremely close to this.
What is the smallest non-zero value that can be seen? Smaller is better.
What is the largest double-precision value less than 1 that cannot be seen? Smaller is better (and negative is best).
How many distinct double-precision values can be seen? Larger is better.
The probability of 0 being low is irrelevant, and having it as 0 ideally is nonsensical. All actual values from a continuous probability distribution have zero probability, so this could apply to all values.
Why do we care about seeing some small value with extremely small probability? It does not make an iota of difference in practice.
Ditto for values close to 1.
Why? I can just use Float128 etc if I really care about that (but I don’t see why I should).
Does the author actually works with random numbers in practice?
As of now, Float32 produces on a grid of size 2^-24 (or is it 23?) Anyway, it’s a crazy bad coarse grid. I think this should be at a minimum 2^-32 resulting in the smallest nonzero number being around 1e-10 this seems of very real world relevance because on my desktop machine it takes less than a minute to make the plots I made above and those kinds of plots are highly relevant to lots of real world probability questions… People want to sample 10M numbers and get a sense of the tail of a distribution, and since the only place that has resolution is near 0 they will map the left tail of [0,1) to the tail of interest if they have knowledge of Float32
The people we are harming are the ones who know enough about Float32 to understand how to use it who have a reasonable expectation it should be able to resolve stuff down into the left tail but find that the standard library is quite inadequate.
Yes, you can get around the problem by rescaling [0,1) into [0,3e-6) or something but why not have the RNG do that?
It especially describes how one should shift singularities around to make use of the higher precision close to zero in floating point numbers.
For a nice intro text that imo gets this right, I like occil. One can also take a look at a C++ implementation that may be extra relevant because it refers to actual simulations failing due to this issue. It also refers to this nice blog post that explains how the C++11 standard has issues.
That being said: “a paper is nice” is unfortunately not relevant unless one one has the context of which authors correspond to which feuding schools, what kind of thing is unwritten, implied and well-known. Literature can’t be approached by reading single articles or books, unfortunately. One needs an apprenticeship to grasp the unwritten subtext, and this is not my field.
PS. I especially like how Monahan is acutely aware of how these issues relate to the number of samples one wishes to draw, and how he will not re-run any simulations because
which is much too large when I once generated 5,000,000 exponentials (using H2) on an IBM machine in a single experiment.
Why is the left tail special here? A lot of distributions are in the positive part of the real line (gamma, exponential, Pareto, etc).
The author’s claim is that using 52 bits to generate random Float64s “is not respecting the pigeonhole principle”, and this is a bug. I fail to understand why.
Personally, I’m approaching this from a documentation standpoint. Our current docs say:
concrete floating point types sample from [0, 1)
This should say more, even with the status quo implementations. Perhaps something like:
The types Float16, Float32, and Float64 return values in [0, 1) as though a random variable were drawn from a theoretically ideal uniform distribution and then rounded down to the previous multiple of eps(T)/2 (this value may decrease in the future).
I’d want to be able to say:
The types Float16, Float32, and Float64 return values in [0, 1) as though a random variable were drawn from a theoretically ideal uniform distribution and then rounded down to the previous floating point number.
Many applications don’t care about the left tail of the distribution. Seldom do mine. Either I’m sampling Bernoulli events with probabilities >1e-5 (usually just to produce test data for some function I actually care about) or I’m drawing random values on (some approximation of) [a,b]. I’ve never tried to sample some vanishing tail of a distribution.
A user whose code breaks for a valid value (e.g., 0) has made a mistake so long as the documentation accurately asserts the possible range, whether the probability that value is 2^{-23} or 2^{-1074}.
It’s better than not if we can fill out the tail without a significant runtime hit. But many uses really don’t need those extra bits and some of those uses are performance-sensitive. That said, usually someone does something with random numbers such that RNG may not be the most expensive part, so we might have a little latitude to regress speed without breaking any real use too badly. But slowing down vectorized generation by more than 2x or 3x is probably a tough sell.
While my benchmark above showed a 7x slowdown on Float32, that was before vectorizing the RNG so there’s some hope we can get it within that range. The comparison for Float64 was actually closer than Float32 (3x-4x slowdown). The precise implementation definitely has room for improvement in any case.
That Rademacher library linked above could make an interesting template for an implementation. With all the effort we’re considering, perhaps we should permit a custom range. This would place the functionality under a new function (e.g. randf), at which point there would be less reason to change the current rand except maybe to extend it to the “native” 16/32/64 bits of entropy.
The left tail is special only on floating point numbers in Uniform(0,1). Only a caveman would sample pareto by mapping u -> (1-u)^-1; obviously one uses u -> u^-1 because of cancellations. Only that doesn’t work because our numbers come pre-cancelled!
Likewise, only a caveman would produce U(-1,1) via u -> 2*u -1, the right way is u -> u * rand(Bool) ? +1 : -1 or u -> u * sign(reinterpret(UInt64, u) & 0x01). Just that the latter is also wrong because the lsb of the mantissa is zero with 75% in our distribution!
Or any person who does not need literally every floating point value. Most real world data does not have a precision of 2^{-24} (72dB SNR) so this wouldn’t even be “noticeable” in most applications even in the significand-only Float32 implementation that people are conceding might stand to be better. It’s all about use cases.
Julia’s gone over 5 years since v1.0, and Python and MATLAB much longer (as far as I can tell), without anyone raising a stink prominently enough to have addressed this already (outside of specialized libraries like Rademacher). Otherwise we’d just look at what was done and implement it without all this debate.
Because the left tail of [0,1) is by design different from the right tail in floating point numbers. We know we have resolving power issues so all the numerical analysis texts will tell you if you need to get close to but not equal to some number, make that number be 0.0f0 by transforming your problem.
Indeed. I am still missing the real world application here. In any random simulation, the noise will overwhelm all possible errors by orders of magnitudes. In the best case scenario (a central limit theorem applies), sampling errors will decrease with \sqrt{} of the sample size. It only gets worse than that (and far worse with extreme values and similar).
You are considering users who
know about floating point cancellations, and still use Float32,
work with random numbers,
yet at the same time remain totally ignorant about how these are generated in the majority of languages today. Note that I am not talking about theoretically “ideal” implementations, but how things are done,
care a lot about a subtle numerical error that will be overshadowed by sampling error in most cases.
Consider eg
julia> using Statistics, Distributions
julia> function simulate_tail_percentiles(::Type{T}, f::F, N) where {T,F}
quantile([f(rand(T)) for _ in 1:N], [0.99, 0.999])
end
simulate_tail_percentiles (generic function with 3 methods)
julia> vcat([simulate_tail_percentiles(Float32, u -> 1/u, 2^20) for _ in 1:10])
# what you propose
10-element Vector{Vector{Float64}}:
[98.7810230255127, 974.203323364292]
[99.67911529541016, 983.9872406006666]
[100.55886840820312, 987.4672012329156]
[101.73590087890625, 1056.1007080078125]
[100.40810775756836, 986.8589416504312]
[101.4587345123291, 994.2953857422232]
[100.84794998168945, 1012.302479553254]
[101.39082717895508, 939.8244155885089]
[101.44968223571777, 1034.891146850589]
[100.05346870422363, 935.6764144897485]
julia> vcat([simulate_tail_percentiles(Float32, u -> 1/(1-u), 2^20) for _ in 1:10])
# "caveman"
10-element Vector{Vector{Float64}}:
[101.67270088195801, 1022.1072372439019]
[100.59202575683594, 952.9247070312676]
[101.27407836914062, 1043.7471038818571]
[101.08370399475098, 989.7363830566842]
[99.90362548828125, 981.0722305298173]
[101.8829345703125, 1023.1262237550194]
[101.62635803222656, 943.1151260376026]
[100.20376396179199, 983.232608032256]
[99.4205379486084, 999.5982208252868]
[101.0535659790039, 1010.6546463013206]
Sampling error overwhelms everything, even with much higher N.
I feel very strongly that whatever decision is reached, it should optimize for ergonomics / principle / performance / etc. etc. etc. for the uniform distribution. so affine transformations like the caveman’s 2*r -1 should work just fine
I care much less about log(r) and whatnot, in those instances it is more obvious to me and presumably most other users that care may be needed around singularities for specialized distributions
Let’s consider a simple textbook StatMech 101 question. Consider the question of computing monte Carlo estimates of material properties with the Lennard Jones potential E(r) = 1/r^{12} - 1/r^6 in some dimensionless measure of radius. This is like a computational statmech 101 question. You might like to know at a certain temperature how close together will two molecules get? So you place one molecule at 0 and randomly sample other molecules in [0,1) and accept a proposed sample with probability exp(-E/kT). Your goal is to get a distribution of radii.
Then later in the semester you come up with better methods, but there’s nothing about this method that’s terrible… Except that the random numbers are on a grid whose resolution is about 1e-8. Admittedly you’re unlikely to have two molecules get that close, but that’s not the issue, suppose at your temperature a typical radius is .0001, then using Float32 randoms your next closer value you could possibly get to is around .00010001 even though eps(.0001) is eps(.0001f0)=7.275958f-12 so the next representable number is around .0001000000001 or something (hard to count those zeros… but you get the point).
Now you probably don’t need resolution down to this 1e-12 range, but you might very well want resolution down to 1e-10 ish… so around 32 bits
EDIT: I accidentally used a double in the above calcs, fixed
And, yes, there are all kinds of ways to do better as the programmer… you could perturb your current sample by a normal random number with scale around 1e-4 and then you’d get all sorts of numbers near 0… and this is probably a reasonable way to do your computation, but I personally find it shocking that the error in the CDF of Float32 looks like the red line below at the left tail when for 11% time cost it could look like the green line
At a fine enough scale, every floating point calculation is wrong. If Float32 is inadequate, Float64 is 536,870,912 times more precise (sometimes 2^{29} just doesn’t get the point across). “rand(Float32) isn’t accurate enough” is a strawman when Float64 is the default.
That said, there has been extremely little pushback to extending Float32 from 24 to 32 bits and Float64 from 53 to 64. We have a number of proposed and mostly-performant ways of producing such “full-width” values correctly. What’s missing for those is optimization and a PR.
The only reason I can imagine this thread is still going is because we’re debating whether we should go further. We’d all like to have a bazillion bits of precision, but the actual solution must be considered with the performance tradeoffs. Some of us (even those of us more comfortable with the status quo) have been trying to produce more bits in a way that does not blow up the run time. This is somewhat more challenging since we must compete with vectorized generation.
The time for philosophical discussion of “what’s right” ended a few-hundred posts ago. The time for debating whether rand should produce some other range than [0,1) ended five years ago when breaking changes were outlawed.
Any discussion of a solution without runtime considerations will only continue this interminable debate, one that has done little to move people’s opinions. The task at hand isn’t to convince people that more bits are better, it’s to convince people that more can be cheap. What we need are implementations.
Do we have then reached consensus on the following:
There will be two samplers (names not representative): FastUniform and IdealUniform. FastUniform will make no tradeoffs that sacrifice speed. IdealUniform will make no tradeoffs that sacrifice accuracy / precision. Behavior will be documented, including recommendations which to use for which applications. Exact zeros for IdealUniform are a separate interminable discussion.
Yes, that will need to happen.
Yes, but I have minor nitpicks.
Disagree
0voters
Where possible without compromising performance, FastUniform will have slightly improved accuracy compared to the current rand(). We can draw inspiration from e.g. amd rocrand for Float32. We will ensure that this works all well on gpu, and we will throw no CPU arch under the performance bus (e.g. avx2 or various arm). We may throw mersenneTwister under the bus, though (e.g. by using all 64 bit for Float64, even though dsfmt natively produces 54 bit).
Yes, that sounds right.
Yes, but I have minor nitpicks.
Disagree
0voters
We will make a reasonable effort to make IdealUniform fast. If it turns out that we can make it almost as fast as FastUniform, then we will, after a grace period, change the default random to this. GPU packages will override this default. This would have the unfortunate effect that unqualified rand() calls on GPU have a subtly different distribution than on CPU.
Yes, that sounds right.
Yes, but I have minor nitpicks.
Disagree
0voters
I’d especially like input by @maleadt on whether the last one would feel like throwing gpu under the bus. You sometimes lose comparability to CPU; you probably gain a rand distribution that is closer to curand / rocrand (compared to the current rand distribution).
IMO it would be best to see the implementations in a package before making a decision about these things, and in the meantime just document the current status quo.
Don’t get me wrong, I would not push back against any improvement, but talking about decisions like
seems a bit premature without anything concrete.
If it turns out that “ideal” is not fast enough, we can still benefit from this project, eg by extending the RNG API to orthogonalize the production of bit strings (32, 64, 128, etc) from transforming these into floating point values. That way, developments in the actual underlying RNG would not affect the choice of the user about what is expected from a random floating point number.
