make uniformRM (0, 1) for floats equivalent to uniformDouble01M / uniformFloat01M
#166
Comments
I don't really feel that there are good reasons to expect that, especially if the only mismatches are numeric noise.
Which means a major version bump and lots of busy work for everyone to bump their upper bounds. IMHO not worth it. |
|
It would be interesting to see how the change affects statistical quality benchmarks like Diehard (used in splitmix). |
|
I was not originally involved in implementing floating point number generation, it was mostly @Shimuuar and @curiousleo that handled it when we did the rewrite in 2020. So, I had to do a bit of archaeology and tried to figure out why current implementation is the way it is. What I found out was that Ironically, it was pointed out in that conversation that we should record this decision, but I can't find it being recorded anywhere in the documentation. Here is the PR that introduced the change: https://github.com/idontgetoutmuch/random/pull/169/files Looking though changes to documentation still did not clarify it for me which of the problems were solved by changing the formula. @Bodigrim Could you shed some light for my own sanity. What did we get with Here is a draft PR that implements the change: #172 My response to other comments in this ticket:
Diehard tests, from what I know, check floating point generation in the (0, 1) range, while this issue is about rounding error of scaling to a different range once a floating point value in (0,1) has already been generated. So, I can't see how diehard tests could be relevant here.
I'll be releasing major version
I do agree with this argument, however there are other parts of this ticket that I do like. In particular I like the simplification of the scaling formula and avoidance of some corner cases.
I dropped this minor optimization for the sake of clarity of the code. We can bring it back at any time if need be. |
|
I decided to play around a bit with
Posting here my proof, so others can double check that I hopefully did not make any mistakes: sat $ do
let scalingFloat l h x = ite (abs l .< abs h) (l + x * (h - l)) (h + x * (l - h))
(l, h, x) <- (,,) <$> free "l" <*> free "h" <*> free "x"
constrain (x .>= 0)
constrain (x .<= 1)
let y = scalingFloat l h (x :: SFloat)
-- Only care about normal values
constrain (sNot (fpIsInfinite y))
constrain (sNot (fpIsNaN y))
constrain $
ite (abs l .< abs h)
-- l is closer to zero:
(ite (l .< h)
(y .< l) -- l is the true low (find y that is lower than l)
(l .< y) -- l is the high (find y that is higher than l)
)
-- h is closer to zero:
(ite (l .< h)
(y .> h) -- h is the true high (find y that is higher than h)
(y .< h) -- l is the low (find y that is lower than h)
) |
Very interesting finding. I would be very interested in the details too, because this comment argues exactly in the opposite direction of my proposal.
The commit message of 0a36e76 argues: “
Looks good to me.
Actually, I argue that even for the range (0, 1) some bits of randomness might be lost due to rounding.
I could imagine that the “numeric noise” might actually be relevant if a user skews the distribution by applying a function to the drawn numbers. For example applying the n-th root (for large n). Not sure if this is a good example or if there are good examples at all. In general, I would be hesitant to throw away some randomness if there is no good reason.
Fair enough. The optimizer might be clever enough to do this anyways.
Nice! Unfortunately, I do not know |
One difference I can recall is that when |
|
We have solution for problem with loss of precision above. But we still another numeric problem: So code for float would looks like (omitting finiteness checks): uniformRM (l,h) g = do
w <- uniformWord32 g
let x = fromIntegral (clearBit w 31) / 2^32 :: Float
if | testBit w 31 -> return $! l + (h-l)*x
| otherwise -> return $! h + (l-h)*xThis approach should not exceed range on either end. In exchange we may get very slight nonuniformity near x=0.5. But I think that shouldn't be a problem P.S. That's vaguely inspired by https://arxiv.org/abs/1704.07949 |
Indeed, when using infinity for the bounds, the proposed implementation behaves differently: -- scale floating old
sfOld :: Num a => a -> a -> a -> a
sfOld l h x = x * l + (1 - x) * h
-- scale floating as proposed
sfNew :: (Ord a, Num a) => a -> a -> a -> a
sfNew l h x = if abs l < abs h then l + x * (h - l) else h + x * (l - h)
sfOldVsNew
= [(l, h, x, sfOld l h x, sfNew l h x)
| let bounds = [(-1) / 0, 0, 1 / 0]
, l <- bounds
, h <- bounds
, x <- [0, 0.5, 1 :: Float]
]
-- >>> sfOldVsNew
-- [ (-Infinity, -Infinity, 0.0, NaN, NaN)
-- , (-Infinity, -Infinity, 0.5, -Infinity, NaN) -- old better
-- , (-Infinity, -Infinity, 1.0, NaN, NaN)
-- , (-Infinity, 0.0, 0.0, NaN, NaN)
-- , (-Infinity, 0.0, 0.5, -Infinity, -Infinity)
-- , (-Infinity, 0.0, 1.0, -Infinity, -Infinity)
-- , (-Infinity, Infinity, 0.0, NaN, NaN)
-- , (-Infinity, Infinity, 0.5, NaN, NaN)
-- , (-Infinity, Infinity, 1.0, NaN, NaN)
-- , ( 0.0, -Infinity, 0.0, -Infinity, NaN) -- old better(?)
-- , ( 0.0, -Infinity, 0.5, -Infinity, -Infinity)
-- , ( 0.0, -Infinity, 1.0, NaN, -Infinity) -- new better(?)
-- , ( 0.0, 0.0, 0.0, 0.0, 0.0)
-- , ( 0.0, 0.0, 0.5, 0.0, 0.0)
-- , ( 0.0, 0.0, 1.0, 0.0, 0.0)
-- , ( 0.0, Infinity, 0.0, Infinity, NaN) -- old better(?)
-- , ( 0.0, Infinity, 0.5, Infinity, Infinity)
-- , ( 0.0, Infinity, 1.0, NaN, Infinity) -- new better(?)
-- , ( Infinity, -Infinity, 0.0, NaN, NaN)
-- , ( Infinity, -Infinity, 0.5, NaN, NaN)
-- , ( Infinity, -Infinity, 1.0, NaN, NaN)
-- , ( Infinity, 0.0, 0.0, NaN, NaN)
-- , ( Infinity, 0.0, 0.5, Infinity, Infinity)
-- , ( Infinity, 0.0, 1.0, Infinity, Infinity)
-- , ( Infinity, Infinity, 0.0, NaN, NaN)
-- , ( Infinity, Infinity, 0.5, Infinity, NaN) -- old better
-- , ( Infinity, Infinity, 1.0, NaN, NaN)
-- ]However, there is a lot Also, infinite bounds are handled separately anyways in the current and the proposed implementations. If this special handling is kept, no infinite value will reach the scaling calculation.
I do not think that the linked paper is applicable here (or it does not make sense to apply it) because the scaling from [0, 1] to [ However, the paper is a good argument why we should keep as much precision (“numeric noise”) as possible. And there actually are cases where the current as well as the proposed solution loses precision and where the paper does not help: When using negative and positive bounds simultaneously, e.g., Regarding the problem of exceeding the bound that is farther from zero, it might make sense to simply clamp the result to the bound. |
Correct
@Bodigrim I see, your suggestion was made before the explicit check on infinity was added. So, there is no longer a need to keep this computation. Thank you, that clarifies things for me.
@Shimuuar That does look interesting. Question for you: why did you use uniformRM (l,h) g = do
w <- uniformWord32 g
let x = fromIntegral (clearBit w 31) / fromIntegral (maxBound :: Word32) :: Float
return $!
if testBit w 31
then l + (h - l) * x
else h + (l - h) * xI can confirm that it does indeed resolves the floating point caveats. However, for some mysterious reason it is a bit slower than using the calculation that @Flupp suggested, but also I believe it is slightly inferior since we are loosing one bit of randomness, that otherwise would be used to generate more floating point noise. I really like @Flupp's suggestion of just clamping the values, which will effectively eliminate those nasty caveats: scaleFloating l h x =
if abs l < abs h
then let !y = l + x * (h - l)
in if l < h -- `l` is closer to zero
then min y h -- `l` is the true low, ensure `y` is not higher than `h`
else max y h -- `l` is the high, ensure `y` is not lower than `h`
else let !y = h + x * (l - h)
in if l < h -- `h` is closer to zero
then max y l -- `h` is the true high, ensure `y` is not lower than `l`
else min y l -- `h` is the low, ensure `y` is not higher than `l`This clamping has negligent impact on performance, but removes the strange behavior. Win-win in my books. Any final thoughts on this, since I want to make a release with this today. (A Christmas present, so to speak 😄) |
|
While investigating some ideas, I actually noticed an actual problem with the proposed solution: The sub-term |
You have something like this in mind, right: ghci> let maxFloat = 3.402823466e+38 :: Float
ghci> maxFloat - negate maxFloat
InfinityThat is a good point. In fact, I don't think that |
That's why I said "vaguely inspired". But we're doing addition and that's enough to get in trouble. Root of problem is
Just general laziness. '2^32 == 2^32-1
On the contrary it's slightly better wrt to randomness. Highest bit is flag 0 means As for performance. Probably either codegen doesn't do good job at mixing bit operations and FP or CPU doesn't like to mix them. It looks like very low level detail |
@Shimuuar Oh yeah you are right. I do like that! However, it also has the same problem of overflowing: #166 (comment)
Probably. The difference is not that big, so if we could get better quality of the generator from it, the penalty is worth paying. |
|
So we have two algorithms each with its own problem:
Maybe we should pick 2nd but check whether |
|
@Shimuuar are you suggesting something along these lines: uniformFRM :: StatefulGen g m => (Float, Float) -> g -> m Float
uniformFRM (l, h) g
| l == h = return l
| isInfinite l || isInfinite h
-- Optimisation exploiting absorption:
-- (+Infinity) + (-Infinity) = NaN
-- (-Infinity) + (+Infinity) = NaN
-- (+Infinity) + _ = +Infinity
-- (-Infinity) + _ = -Infinity
-- _ + (+Infinity) = +Infinity
-- _ + (-Infinity) = -Infinity
= return $! h + l
| otherwise = do
w <- uniformWord32 g
let x =
fromIntegral (clearBit w 31) / fromIntegral (maxBound :: Word32) :: Float
diff = h - l
return
$! if isInfinite diff
then if testBit w 31
then l * (1 - x) + h * x
else h * (1 - x) + l * x
else if testBit w 31
then l + diff * x
else h + negate diff * x |
|
We could even further clamp the values for the if isInfinite diff
then
let y =
if testBit w 31
then l * (1 - x) + h * x
else h * (1 - x) + l * x
in max (min y (max l h)) (min l h)
... |
|
Yes. Although I'm not sure that trick with And clamping is very good idea. If either |
Good point. This is the final draft that I have: w <- uniformWord32 g
if isInfinite diff
then let !x = fromIntegral w / m
!y = x * l + (1 - x) * h
in max (min y (max l h)) (min l h)
else let !x = fromIntegral (clearBit w 31) / m
in if testBit w 31
then l + diff * x
else h + negate diff * x
where
!diff = h - l
m = fromIntegral (maxBound :: w) :: a |
|
I've adjusted my PR: #172 Please leave feedback and if there will be no objections, I will make a new release with a fix later on today |
uniformRM (l, h)forFloatandDoubleis defined by drawingxfrom [0, 1] usinguniformDouble01M/uniformFloat01Mand then returningx * l + (1 - x) * h.One might expect that
uniformRM (0, 1)may produce the same values asuniformDouble01M/uniformFloat01M. However, this is not the case because of rounding errors when calculating1 - x.For visualization of the rounding problem try:
Note that the
x'values have a lot of repetitions while eachxis different.An easy fix would be changing the calculation to
(1 - x) * l + x * h. However, then the possible results ofuniformRM (-1, 0)are reduced instead.A possible solution could be to first (mathematically) transform the two calculations as follows:
x * l + (1 - x) * h = h + x * (l - h)(1 - x) * l + x * h = l + x * (h - l)Then use calculation 2 when
lis closer to zero thanh, and calculation 2 otherwise. This may be implemented as follows:(This also drops the binding for the second parameter of
uniformRMin order to allow memoization of the bound preprocessing.)This approach also slightly reduces the documented floating point number caveats, because this implementation guarantees that the bound that is closer to 0 is not exceeded.
Note: Obviously, this would break backwards compatibility as it changes the result of
uniformRMfor a given state of the RNG.The text was updated successfully, but these errors were encountered: