Bounds of Random.float (exclusive vs inclusive)

Learning
1

I would like to create random floats following a Gaussian (normal) distribution with μ=0 and σ=1.

I’m using the Box-Muller transform, and this is my source code:

let rec gaussian_floats : unit -> float Seq.t =
  let pi_2 = 2.0 *. Float.pi in
  fun () ->
    let rng = Random.split () in
    let rec next () =
      let a = sqrt (-2.0 *. log (Random.State.float rng 1.)) in
      let b = Random.State.float rng pi_2 in
      Seq.Cons (a *. cos b, fun () -> Seq.Cons (a *. sin b, next))
    in
    next

let () =
  gaussian_floats () |> Seq.take 10
  |> Seq.map string_of_float |> Seq.iter print_endline

Following the specifications of Random.float this program could create infinite floats, because:

Random.float bound returns a random floating-point number between 0 and bound (inclusive).

and

# sqrt (-2.0 *. log 0.);;
- : float = infinity

This is something I don’t want to happen (even if the probability is low).

I could do my own if-clause to handle that case, but it’s an extra branch. If the bounds were excluding the upper limit, I could do 1.0 - Random.float 1. Another alternative that even works when both bounds are inclusive is: min_float +. Random.float 1. (because min_float +. 1. = 1.).

But these are extra calculations I would want to avoid (even if they don’t hurt a lot).

Looking into the stdlib, I see that there is already an if-clause to exclude results of zero (assuming the upper bound is not extremely small):

Putting the same if-clause in my own code seems to be redundant, but the specification in the manual forces me to do that.

Also, the upper limit seems to be exclusive, following the current implementation. For any finite positive number x, I would expect Random.float x < x:

# Int64.to_float (Int64.shift_right_logical (-1L) 11)
*. 0x1.p-53 *. max_float < max_float;;
- : bool = true

But this isn’t guaranteed by the specification either.

What should I do? I feel like this is a problem other people may deal with as well.

  • Should the if-clause be removed from the stdlib (because it’s unnecessary)?
  • Could the manual be refined to specify a tighter bound when the upper bound is non-zero and greater than the smallest subnormal number?

Consider:

# 1.0 *. 0x1.p-53 *. (0x1.p-1022);;
- : float = 0.
# 1.0 *. 0x1.p-53 *. (0x1.p-1022 +. 0x1.p-1074);;
- : float = 4.94065645841e-324

Not sure if depending on rounding modes, the limit may be slightly different, but in either case it is so low (lower than 10 to the power of minus 307) that for most application cases, we could safely assume that a random float will never be zero (if the manual would be refined accordingly).

Or am I overlooking some detail here, or is there a good reason why these corner cases should not be specified, e.g. to allow future changes of the implementation? In the latter case, maybe I should implement my own version of Random.float?

2

Apparently the lower and upper bounds became (practically) exclusive in this commit in January 2022 by @xavierleroy. Looking into the corresponding pull request #10742, I found this comment by xavierleroy:

This was also discussed in the RFC: ocaml/RFCs#28 (comment) Eventually, I’d like to export rawfloat (it’s called uniform in PRINGO) with explicit guarantees that it returns a random FP between 0 and 1 excluded.

What is the current state? I think that function has never been exported? How would I currently solve my task to create a random float that is strictly greater than zero?

P.S.: I mean, of course I can add an if x = 0. then …. It’s not really expensive. But I wonder if there is an elegant solution and if that issue is going to be adressed in the stdlib or was forgotten.

3

This does look like a documentation bug in the standard library, I will have a look at it.

Nevertheless, note that a conditional with a branch that is never taken has essentially no cost due to branch prediction. In particular, it seems unlikely to have any performance effect on your gaussian sampling implementation which is:

  • computing a sqrt root, a log, a cos, sin function every two samples
  • using Seq and thus doing three allocations by sample

Also it is unclear to me why you are using Random.split on every sample which are purely sequential.

In other words, in order to optimize meaningfully your code, you should benchmark your code first.
(And probably also compare with a ziggurat implementation). Because yes, languages don’t have a “speed”, it is the combination of a language, a compiler, a processor, an algorithm and a task at hand that has a performance characteristic.

1 Like
4

Not sure if it is: See the comment on exporting rawfloat, which would solve the issue I had. Maybe Random.float has deliberately weaker guarantees. :man_shrugging:

There will be a conversion from int64 to float between the two ifs. Is that still optimized? (Just asking out of curiosity, I’m aware a single if is little cost.)

I did make some (rough) measurements before, and even the floating point operations seem neglible. Also, I learned in the other thread that Seq is expensive. But apparently creating the random data is the bottleneck.

Specifically: It seems to be worth creating two floats out of two floats by remembering the intermediate results.

I don’t think I do? I split it when creating the (ephemeral) sequence, not on every sample, am I wrong?

You mean adding up sample from uniform distributions to obtain a Gaussian one? I tried this with 12 samples (from 6 random int64s), but that was slower, apparently. (edit: I just noticed ziggurat is a different algorithm. I’ll look into that as an alternative, thanks for mentioning it.)

I don’t need an as-fast-as possible algorithm. But if possible, I’d like to omit unnecessary statements if that makes my code better readable (and not slower).

If someone is aware of another suitable algorithm to create a Gaussian distribution (I’ll want integers in the end), I’m thankful for advice.

5

Branch prediction is a (statistical) hardware optimization, what matters is that during the execution only one branch will be taken for a typical run, and thus the speculative execution of the most probable branch will almost always succeed.

Right, I was mislead by the rec keyword on the toplevel definition. Sorry for the noise!

Indeed, Marsaglia’s ziggurat algorithm is a generic rejection sampling method which avoids using any transcendental float functions on the fast path at the cost of a larger memory footprint.

1 Like
6

@octachron already gave detailed replies, but just to summarize:

  • The current implementation of Random.float 1.0 guarantees that the result is never 0.0 nor 1.0, which is nice if you’re going to take logarithms afterwards, while not introducing any significant bias. This is not documented, however.
  • In the general case, Random.float x is x *. Random.float 1.0, so FP rounding can produce the values 0.0 and x, albeit with very low probability (if x is not 0.0).
  • If that’s a problem, you can do rejection sampling:
  let rec random_float x =
    let r = Random.float x in
    if r > 0 && r < x then x else random_float x
  • Perhaps the Random.float standard library function should be reimplemented and documented to perform this rejection sampling, instead of the current “let’s draw a non-zero 52-bit integer and scale it to a float in (0,1) then to a float in [0,x]” algorithm.
2 Likes
7

Alright, thanks for your responses.

I’ll then just use rejection sampling for now (i.e. repeat drawing a random value if I get a value that hits the boundary) to be sure I don’t rely on implementation details that are not specified by the documentation.

If the documentation gets updated and/or a version of rawfloat gets exported, I’ll simplify my routine in the future.

8

I finally decided to get rid of Seq and use a dispenser instead, as I don’t need these two features of Seq:

  • Marking an end of the sequence using Seq.Nil (my “sequence” is infinite),
  • Support for persistence (I don’t need to store the current state of iteration, as random values should always be treated as ephemeral).

So I guess this should save me two out of three allocations (Seq.Cons and the closure for obtaining the next element).

I did see there was a recommendation to not use dispensers in the manual though (here):

When possible, we recommend using sequences rather than dispensers (functions of type unit -> 'a option that produce elements upon demand). Whereas sequences can be persistent or ephemeral, dispensers are always ephemeral, and are typically more difficult to work with than sequences.

But I guess in my case, it should be better to use a dispenser. At least I don’t understand why I shouldn’t use one. (If there is one, I’m happy to learn.) This is my code currently:

let random_gaussian_float_dispenser : unit -> unit -> float =
  let pi_2 = 2.0 *. Float.pi in
  fun () ->
    let rng = Random.split () in
    let rec rand_nonzero () =
      let x = Random.State.float rng 1. in
      if x = 0. then rand_nonzero () else x
    in
    let state = ref None in
    fun () ->
      match !state with
      | None -> begin
          let a = sqrt (-2.0 *. log (rand_nonzero ())) in
          let b = Random.State.float rng pi_2 in
          state := Some (a, b);
          a *. cos b
        end
      | Some (a, b) -> begin
          state := None;
          a *. sin b
        end

Unfortunately, it’s a bit less funtional and uses references.

I guess I could also make it to fill an array, but as floats are boxed when taken out (right?), I wonder if accessing values from that array would then create an allocation anyway.

My intention isn’t to write extremely efficient code here, I just want it to look neat and not be unnecessarily slow. Using dispensers, I can do something like:

let next_random_float = random_gaussian_float_dispenser () in
fun () -> Iarray.map (fun x -> x *. next_random_float ()) stddevs

Side note: I noticed that the documentation of Stdlib.Seq and Seq differs. I guess that’s unintended?

9

You might be interested in GitHub - c-cube/gen: Simple, efficient iterators for OCaml · GitHub , as an aside. It’s a library for dispensers (although I named then “generators” back then) :slightly_smiling_face:

10

Unless x is a subnormal number, you cannot have x * y = x with y < 1 and rounding to nearest. (Not sure I can give a good intuition. Let us say that, relatively, multiplying by x amplify the distance between y and 1, and no amount of rounding of x * y will counterbalance it, unless you are rounding toward infinity.)

11

Yeah I see. Thanks for sharing, I’ll keep that in mind. (I don’t think I’ll need it though in this particular scenario, as my needs are quite limited, and I’d like to keep dependencies minimal.)

Also note that my last approach isn’t really a dispenser in terms of Gen.t or Seq.to_dispenser, because I haven’t forseen a termination (i.e. I use unit -> float instead of unit -> float option).

12

Good to know ! I wasn’t aware of that. But I’m sure our demanding users want to do Random.float x with x being denormal :slight_smile:

13

Not sure what the default rounding mode is (how do I set it?), but I get:

let y = Float.pred 1. in
  List.map
    (fun x -> classify_float x, x *. y = x)
    [Float.pred 0x1p-1022; 0x1p-1022; Float.succ 0x1p-1022];;

- : (fpclass * bool) list =
[(FP_subnormal, true); (FP_normal, true); (FP_normal, false)]

So with x being 0x1p-1022, which is the smallest non-subnormal (i.e. normal) number, we can have x *. y = x. Only if x is slightly bigger, the equation is violated.

But not sure if my machine performs rounding to nearest?

14

No, you are right. I tend to have an extensive definition of the subnormal range (i.e., any number that can be reached by rounding a subnormal result, so this includes the smallest normal number), but this makes for a confusing statement and I should have been clearer.

15

I’d like to give an updated summary then, to conclude: (based on post #6)

  • The current implementation of Random.float 1.0 guarantees that the result is never 0.0 nor 1.0.
  • With rounding to nearest, the current implementation also guarantees that the result of Random.float x is never 0.0 for any x that is greater than the smallest normal number (i.e. greater than 0x1p-1022).
  • None of these observations is guaranteed by the documentation as of now, so I should not rely on that, and I can do rejection sampling as in:
  • The stdlib function Random.float could possibly be reimplemented and documented to include this rejection sampling in the future (as it it isn’t a huge performance cost).
  • If Random.float is not changed, the current behavior (giving the above guarantees) could or could not be documented (depending on whether that should remain an implementation detail or not).
1 Like