# Problems with Octave random number generation, part 1 Classic List Threaded 7 messages Open this post in threaded view
|

## Problems with Octave random number generation, part 1

 Maintainers, I believe there are serious, but subtle, problems with Octave's random number generation.  Resolving them, without introducing further issues, requires careful thought and review beyond just a single coder (myself).  Hence, I'm bringing this issue to the Maintainer's list.  I've also CC'ed Jordi and Philipp Kutin who were active in initially discovering the bugs, and Michael Leitner who was helpful in removing the bias from Octave's randi function. The tip of the iceberg was bug report #41742 (https://savannah.gnu.org/bugs/?41742).  In this report, using rand to generate an integer in the range [1, a] using the straightforward code x = round (rand * a + 0.5) occasionally fails and returns a result of a+1 outside the desired range.  This is only possible if rand should return exactly 1.  But, according to the documentation, rand produces a floating point number in the range (0, 1), i.e., endpoints are excluded.  Unfortunately, it is quite simple to demonstrate that the rand function does include the endpoint.  Test code that I found that always reproduces the problem is octave:1> rand ("state", 1); octave:2> x = max (rand ([20e6, 1], "single")) x = 1 This result is contrary to the documentation, and contrary to the behavior of the rand function in Matlab.  Discussion between Jordi and Philipp ensued at this thread https://octave.1599824.n4.nabble.com/Re-Random-number-generation-quirk-td4662514.html. The Octave code which generates a "single" (32-bit IEEE floating point) random number is found in liboctave/numeric/randmtzig.cc   /* generates a random number on (0,1)-real-interval */   static float randu32 (void)   {     return (static_cast (randi32 ()) + 0.5) * (1.0/4294967296.0);     /* divided by 2^32 */   } The function randi32 (which I have no reason to doubt) generates an unsigned 32-bit integer in the range [0, 2^32 - 1].  Unfortunately, I think there was some confusion by the original programmer between Octave and C implicit conversion rules.  In Octave, arguments are demoted to lower precision before the operand is applied: single + double = single This is exactly reversed in C where arguments are promoted to the more precise class.  In C, single + double = double So this code is actually performed completely with double precision values, and then truncated to single precision float at the end of the function because the return type for the function is float.  In Octave pseudo-code, single ( (double (randi32 ()) + 0.5) * ( 1.0 / 4294967296.0) ) A useful Octave function for mimicking the C code is fs = @(x) single ((floor (x) + 0.5) * (1.0 / 4294967296.0)) Using fs() one can determine the maximum valuable allowable in order for the result not to become exactly equal to 1. octave:28> fs (2^32-129) ans = 0.9999999 octave:29> fs (2^32-128) ans = 1 I assume Philipp did something similar.  He has provided a patch that, after unpacking various macros, does the following static float randu32 (void) {   uint32_t i;   do     {       i = randi32 ();     }   while (i == 0 || i > (2^32 - 129);   return i * (1.0 / 4294967296.0); } My fear, and more mathematical help in this specific would be nice, is that this will introduce a bias.  In Octave, the new code can be modeled as fs2 = @(x) single (floor (x) * (1.0 / 4294967296.0 )) Trying some sample values octave:78> fs2(1) ans = 2.3283064e-10 octave:79> fs2(2) ans = 4.6566129e-10 octave:80> fs2(3) ans = 6.9849193e-10 octave:82> fs2(2^32-128) ans = 1 octave:83> fs2(2^32-129) ans = 0.9999999 octave:84> fs2(2^32-130) ans = 0.9999999 octave:100> fs2(2^32-129) - fs2(2^32-383) ans = 0 What one notices is that the small integer values all get mapped to different floating point values, but at the upper end, large ranges of different integer values get mapped to the same floating point value. Given that IEEE single-precision floating point has only 24 bits for the mantissa, I think something like this would work better   static float randu32 (void)   {     uint32_t i;     do       {         i = randi32 () & static_cast (0xFFFFFF);       }     while (i == 0 || i == 0xFFFFFF);     return i * (1.0f / 16777216.0f);   } In my testing, this does not produce exact values of 1, nor does it map multiple integer values to the same floating point value.  Is anyone aware of problems with this algorithm? Resolving this must happen, but is only part of the issue.  There is a #define in randmtzig.cc #define RANDU randu32() and then RANDU is used in the algorithms for rand_normal and rnd_exponential.  In particular, the use occurs as the argument to the log function.  For example, yy = - std::log (RANDU); The old RANDU = randu32() occasionally put out values of 1 which will be mapped by the log function to 0.  Whereas any new function will always be slightly less than one which means the log function output will be small, but not exactly 0.  Do rand_exponential and rand_normal rely on the old behavior in order to function correctly?  I have no idea. --Rik
Open this post in threaded view
|

## Re: Problems with Octave random number generation, part 1

 Just for the sake of comparison here is the code used in numpy, the convention (and seems to be widespread is to sample [0,1) ) https://github.com/numpy/numpy/blob/master/numpy/core/include/numpy/random/distributions.hhttps://github.com/numpy/numpy/blob/master/numpy/random/src/distributions/distributions.chttps://github.com/numpy/numpy/blob/master/numpy/core/include/numpy/random/bitgen.hThere we see static NPY_INLINE float next_float(bitgen_t *bitgen_state) { return (next_uint32(bitgen_state) >> 9) * (1.0f / 8388608.0f); } and in the mt library we see https://github.com/numpy/numpy/blob/master/numpy/random/src/mt19937/mt19937.hstatic inline double mt19937_next_double(mt19937_state *state) { int32_t a = mt19937_next(state) >> 5, b = mt19937_next(state) >> 6; return (a * 67108864.0 + b) / 9007199254740992.0; } On Tue, Dec 17, 2019 at 8:27 PM Rik <[hidden email]> wrote: > > Maintainers, > > I believe there are serious, but subtle, problems with Octave's random > number generation.  Resolving them, without introducing further issues, > requires careful thought and review beyond just a single coder (myself). > Hence, I'm bringing this issue to the Maintainer's list.  I've also CC'ed > Jordi and Philipp Kutin who were active in initially discovering the bugs, > and Michael Leitner who was helpful in removing the bias from Octave's > randi function. > > The tip of the iceberg was bug report #41742 > (https://savannah.gnu.org/bugs/?41742).  In this report, using rand to > generate an integer in the range [1, a] using the straightforward code > > x = round (rand * a + 0.5) > > occasionally fails and returns a result of a+1 outside the desired range. > This is only possible if rand should return exactly 1.  But, according to > the documentation, rand produces a floating point number in the range (0, > 1), i.e., endpoints are excluded.  Unfortunately, it is quite simple to > demonstrate that the rand function does include the endpoint.  Test code > that I found that always reproduces the problem is > > octave:1> rand ("state", 1); > octave:2> x = max (rand ([20e6, 1], "single")) > x = 1 > > This result is contrary to the documentation, and contrary to the behavior > of the rand function in Matlab.  Discussion between Jordi and Philipp > ensued at this thread > https://octave.1599824.n4.nabble.com/Re-Random-number-generation-quirk-td4662514.html. > > The Octave code which generates a "single" (32-bit IEEE floating point) > random number is found in liboctave/numeric/randmtzig.cc > >   /* generates a random number on (0,1)-real-interval */ >   static float randu32 (void) >   { >     return (static_cast (randi32 ()) + 0.5) * (1.0/4294967296.0); >     /* divided by 2^32 */ >   } > > The function randi32 (which I have no reason to doubt) generates an > unsigned 32-bit integer in the range [0, 2^32 - 1].  Unfortunately, I think > there was some confusion by the original programmer between Octave and C > implicit conversion rules.  In Octave, arguments are demoted to lower > precision before the operand is applied: > > single + double = single > > This is exactly reversed in C where arguments are promoted to the more > precise class.  In C, > > single + double = double > > So this code is actually performed completely with double precision values, > and then truncated to single precision float at the end of the function > because the return type for the function is float.  In Octave pseudo-code, > > single ( (double (randi32 ()) + 0.5) * ( 1.0 / 4294967296.0) ) > > A useful Octave function for mimicking the C code is > > fs = @(x) single ((floor (x) + 0.5) * (1.0 / 4294967296.0)) > > Using fs() one can determine the maximum valuable allowable in order for > the result not to become exactly equal to 1. > > octave:28> fs (2^32-129) > ans = 0.9999999 > octave:29> fs (2^32-128) > ans = 1 > > I assume Philipp did something similar.  He has provided a patch that, > after unpacking various macros, does the following > > static float randu32 (void) > { >   uint32_t i; > >   do >     { >       i = randi32 (); >     } >   while (i == 0 || i > (2^32 - 129); > >   return i * (1.0 / 4294967296.0); > } > > My fear, and more mathematical help in this specific would be nice, is that > this will introduce a bias.  In Octave, the new code can be modeled as > > fs2 = @(x) single (floor (x) * (1.0 / 4294967296.0 )) > > Trying some sample values > > octave:78> fs2(1) > ans = 2.3283064e-10 > octave:79> fs2(2) > ans = 4.6566129e-10 > octave:80> fs2(3) > ans = 6.9849193e-10 > octave:82> fs2(2^32-128) > ans = 1 > octave:83> fs2(2^32-129) > ans = 0.9999999 > octave:84> fs2(2^32-130) > ans = 0.9999999 > octave:100> fs2(2^32-129) - fs2(2^32-383) > ans = 0 > > What one notices is that the small integer values all get mapped to > different floating point values, but at the upper end, large ranges of > different integer values get mapped to the same floating point value. > > Given that IEEE single-precision floating point has only 24 bits for the > mantissa, I think something like this would work better > >   static float randu32 (void) >   { >     uint32_t i; > >     do >       { >         i = randi32 () & static_cast (0xFFFFFF); >       } >     while (i == 0 || i == 0xFFFFFF); > >     return i * (1.0f / 16777216.0f); >   } > > In my testing, this does not produce exact values of 1, nor does it map > multiple integer values to the same floating point value.  Is anyone aware > of problems with this algorithm? > > Resolving this must happen, but is only part of the issue.  There is a > #define in randmtzig.cc > > #define RANDU randu32() > > and then RANDU is used in the algorithms for rand_normal and > rnd_exponential.  In particular, the use occurs as the argument to the log > function.  For example, > > yy = - std::log (RANDU); > > The old RANDU = randu32() occasionally put out values of 1 which will be > mapped by the log function to 0.  Whereas any new function will always be > slightly less than one which means the log function output will be small, > but not exactly 0.  Do rand_exponential and rand_normal rely on the old > behavior in order to function correctly?  I have no idea. > > --Rik > > > > >
Open this post in threaded view
|

## Re: Problems with Octave random number generation, part 1

Open this post in threaded view
|

## Re: Problems with Octave random number generation, part 1

 In reply to this post by Juan Pablo Carbajal-2 That's in essence Rik's approach, but not rejecting 0. Thus, it would be incompatible with Matlab, which specifies (0,1). -- Sent from: https://octave.1599824.n4.nabble.com/Octave-Maintainers-f1638794.html
Open this post in threaded view
|

## Re: Problems with Octave random number generation, part 1

 In reply to this post by Rik-4 On 12/17/19 2:27 PM, Rik wrote: > Maintainers, > > I believe there are serious, but subtle, problems with Octave's random > number generation.  Resolving them, without introducing further issues, > requires careful thought and review beyond just a single coder (myself). > Hence, I'm bringing this issue to the Maintainer's list.  I've also CC'ed > Jordi and Philipp Kutin who were active in initially discovering the bugs, > and Michael Leitner who was helpful in removing the bias from Octave's > randi function. > > The tip of the iceberg was bug report #41742 > (https://savannah.gnu.org/bugs/?41742).  In this report, using rand to > generate an integer in the range [1, a] using the straightforward code > > x = round (rand * a + 0.5) > > occasionally fails and returns a result of a+1 outside the desired range. > This is only possible if rand should return exactly 1.  But, according to > the documentation, rand produces a floating point number in the range (0, > 1), i.e., endpoints are excluded.  Unfortunately, it is quite simple to > demonstrate that the rand function does include the endpoint.  Test code > that I found that always reproduces the problem is > > octave:1> rand ("state", 1); > octave:2> x = max (rand ([20e6, 1], "single")) > x = 1 > > This result is contrary to the documentation, and contrary to the behavior > of the rand function in Matlab.  Discussion between Jordi and Philipp > ensued at this thread > https://octave.1599824.n4.nabble.com/Re-Random-number-generation-quirk-td4662514.html. > > The Octave code which generates a "single" (32-bit IEEE floating point) > random number is found in liboctave/numeric/randmtzig.cc > >    /* generates a random number on (0,1)-real-interval */ >    static float randu32 (void) >    { >      return (static_cast (randi32 ()) + 0.5) * (1.0/4294967296.0); >      /* divided by 2^32 */ >    } > > The function randi32 (which I have no reason to doubt) generates an > unsigned 32-bit integer in the range [0, 2^32 - 1].  Unfortunately, I think > there was some confusion by the original programmer between Octave and C > implicit conversion rules.  In Octave, arguments are demoted to lower > precision before the operand is applied: > > single + double = single > > This is exactly reversed in C where arguments are promoted to the more > precise class.  In C, > > single + double = double > > So this code is actually performed completely with double precision values, > and then truncated to single precision float at the end of the function > because the return type for the function is float.  In Octave pseudo-code, > > single ( (double (randi32 ()) + 0.5) * ( 1.0 / 4294967296.0) ) > > A useful Octave function for mimicking the C code is > > fs = @(x) single ((floor (x) + 0.5) * (1.0 / 4294967296.0)) > > Using fs() one can determine the maximum valuable allowable in order for > the result not to become exactly equal to 1. > > octave:28> fs (2^32-129) > ans = 0.9999999 > octave:29> fs (2^32-128) > ans = 1 > > I assume Philipp did something similar.  He has provided a patch that, > after unpacking various macros, does the following > > static float randu32 (void) > { >    uint32_t i; > >    do >      { >        i = randi32 (); >      } >    while (i == 0 || i > (2^32 - 129); > is >    return i * (1.0 / 4294967296.0); > } > > My fear, and more mathematical help in this specific would be nice, is that > this will introduce a bias.(1.0/4294967296.0) In the discussion between Jordi, Mike and Philipp is the link to the code: http://hg.savannah.gnu.org/hgweb/octave/file/c579bd4e12c9/liboctave/numeric/randmtzig.c#l392Yes, I don't like the notion of tossing out a range of numbers using the while loop.  RNG's have been extensively analyzed for randomness, and the conditional statement above has the *potential* to effect that. Then again, one analysis might be that good pseudo-RNG should have subintervals that look properly random. Small technicality is that    return ((float)randi32 () + 0.5) * (1.0/4294967296.0);    /* divided by 2^32 */ is actually multiplying by the inverse, not dividing.  Fine for real numbers and a CPU-saving technique in a general sense, but there will always be a bit of a guess as to exactly what a hardware floating point processor is doing.  It probably doesn't make a difference in this case, though. But along that line, is there a slightly smarter way of doing the arithmetic closer to integer math?  For example, if there is a way of avoiding the rounding problems by generating a range that is as delta + [0,M-1] * alpha it might avoid the rare value going out of range (which is kind of what you did in your proposed alternative later). > In Octave, the new code can be modeled as > > fs2 = @(x) single (floor (x) * (1.0 / 4294967296.0 )) > > Trying some sample values > > octave:78> fs2(1) > ans = 2.3283064e-10 > octave:79> fs2(2) > ans = 4.6566129e-10 > octave:80> fs2(3) > ans = 6.9849193e-10 > octave:82> fs2(2^32-128) > ans = 1 > octave:83> fs2(2^32-129) > ans = 0.9999999  return ((float)randi32 () + 0.5) * (1.0/4294967296.0);    /* divided by 2^32 */ > octave:84> fs2(2^32-130) > ans = 0.9999999 > octave:100> fs2(2^32-129) - fs2(2^32-383) > ans = 0 > > What one notices is that the small integer values all get mapped to > different floating point values, but at the upper end, large ranges of > different integer values get mapped to the same floating point value. > > Given that IEEE single-precision floating point has only 24 bits for the > mantissa, I think something like this would work better > >    static float randu32 (void) >    { >      uint32_t i; > >      do >        { >          i = randi32 () & static_cast (0xFFFFFF); >        } >      while (i == 0 || i == 0xFFFFFF); > >      return i * (1.0f / 16777216.0f); >    } > > In my testing, this does not produce exact values of 1, nor does it map > multiple integer values to the same floating point value.  Is anyone aware > of problems with this algorithm? Creative, but mapping out that whole upper level of pseudo-random numbers doesn't feel good.  It may be as good as anything though because, yeah, we could change the name of the function to randu24() and it would have the same meaning because of the imposed range of (0,1). In some way, that range must implicitly toss away the exponent and the floating point bit.  It's effectively now fractional point, not floating point. Are we doing a "horse-is-out-of-the-barn" type of thing by generating this (0:1) floating point range and then expecting to scale it?  Say I take a vector of numbers from the RNG with range [1/2^24 : 2^24-1/2^24] (or think of mantissa [1:0xFFFFFE] if you like) and multiply by floating point value greater than 1?  How about less than 1?  Probably some strange dithering effect on the floating point numbers in either case, unless the scaling is a power of 2. > Resolving this must happen, but is only part of the issue.  There is a > #define in randmtzig.cc > > #define RANDU randu32() > > and then RANDU is used in the algorithms for rand_normal and > rnd_exponential.  In particular, the use occurs as the argument to the log > function.  For example, > > yy = - std::log (RANDU); > > The old RANDU = randu32() occasionally put out values of 1 which will be > mapped by the log function to 0.  Whereas any new function will always be > slightly less than one which means the log function output will be small, > but not exactly 0.  Do rand_exponential and rand_normal rely on the old > behavior in order to function correctly?  I have no idea. This is probably where the significance of [0:1) versus (0:1) makes a difference.  Several math functions throw an error or require special treatment with an input of 0.  It would be a mistake to think of floating point in terms of real number (0:1).  Instead it means [delta:1-delta] where delta is the resolution of the floating point number. Dan