bugGNU Octave - Bugs: bug #53299, Document more fully rand...

 
 

bug #53299: Document more fully rand state/seed behavior

Submitter:  deego <deego>
Submitted:  Wed 07 Mar 2018 06:57:45 AM UTC
   
 
Category:  Documentation Severity:  2 - Minor
Priority:  5 - Normal Item Group:  Documentation
Status:  In Progress Assigned to:  None
Originator Name:  DAVE GOEL Open/Closed:  * Open
Release:  * 4.2.1 Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sat 10 Mar 2018 05:15:48 PM UTC, comment #20: 

I'm marking this as "In Progress".  I've started re-drafting the documentation, but it is actually quite difficult to succinctly explain the behavioral differences between the random number generators and capture all of the funny corner cases.  This is another indication that we should be deprecating and discouraging the use of the older generators.

Rik <rik5>
Group administrator
Thu 08 Mar 2018 04:47:17 PM UTC, comment #19: 

@Mike: For conv2 bugs, see bug #39314 and bug #34893.  This function has consistently been a problem.  The root issue now seems to be that the inner and outer convolutions are performed in different directions and the summation of floating point numbers can accumulate different small errors because of that.

For some reason, the random numbers that are generated by the "seed" generator are apparently precisely represented in IEEE-854 so that they do not accumulate these errors.  I tested with different random values of the seed itself and the conv2 test always pass.  Thus, it is not the particular value (42) that we use in the BIST tests, but the fact that we use the "seed" generator.

I also generated a random set of numbers, saved them to a file, and then re-started Octave.  I wanted to make sure that it wasn't a strange side effect like initializing the "seed" generator which also set some CPU processor flags.  But reading in the numbers and executing the conv2 test on them works.  Hence, it really is the numbers themselves that have the necessary magic.

Anyways, it would be better to resolve the issue with the conv2 algorithm rather than messing with our random number generation.

Rik <rik5>
Group administrator
Thu 08 Mar 2018 04:28:31 PM UTC, comment #18: 

Rik, rightly, said:But I have have been in design situations where the exact random values in the exact order encountered previously were desired. The only way to guarantee that, it would seem, would be to keep the "seed" generator implementation.

This is important for any work that is, for example, published.
It is essential to be able to reproduce computed results.
If necessary, this would mean capturing the outputs of randn
and saving them for future use. This may be the safest method
for now. But, it would be very good to explain what is going
on in the (new) documentation. It is not really necessary to
keep the "seed" generator. Actually, it may be best to remove
it so that people do not get confused and, in any case, lose the
data that they will need in order to reproduce their results.

Michael Godfrey <godfrey>
Group Member
Thu 08 Mar 2018 01:21:36 PM UTC, comment #17: 

I would say that a first step would be to actually undocument the "seed" generator in the individual randXXX functions.  If the feature is unknown to programmers, it won't get used.

The next step is to collect the information about the "seed" generator and put it in a subheading of the manual, along with the warning that this is all subject to deprecation and shouldn't be used in any new code and that you might get strange results.

I doubt we can implement an API on top of the "state" generators to mimic "seed" generators.  For most people, it is true, they don't care too much about the generator as long as their results are random.  But I have have been in design situations where the exact random values in the exact order encountered previously were desired.  The only way to guarantee that, it would seem, would be to keep the "seed" generator implementation.

Rik <rik5>
Group administrator
Thu 08 Mar 2018 04:10:26 AM UTC, comment #16: 

So can we say that this bug is about documenting the following two peculiarities:

1. rand('state', 'reset') seems to not depend on the system clock, but rand('seed', 'reset') does depend on the integer value of the system clock (same goes for the other rand* functions).

2. rand('seed', ...) switches the internal algorithm for all of the rand functions. This is documented in the section of the manual describing the rand functions, but not in the rand doc strings themselves.

Does that summarize the issues here?

The startup random state bug I found is now bug #53300.

I agree with deego, it might be possible to turn the rand('seed', ...) style calls into wrappers for the new algorithm and just remove the old algorithms completely. There might still be some behavioral compatibility issues to think about. That could be tracked in a separate bug report or project.

Mike Miller <mtmiller>
Group Member
Thu 08 Mar 2018 02:05:40 AM UTC, comment #15: 


>> foo = rand("seed") => Implements K = round(1e20*rand()), implement rand("state", K), and output K. [or minor variants thereupon.]



Also wanted to mention: K is actually determiniscally determined from the current state. So, not only is the state reproducible going forward, it is also reproducible back in time.

Namely, even if the seed-user tried something really a weird "backwards"-test, he will still get the right answer:


K = rand('seed');

rand();

L = rand('seed');


rand('seed' , K)
rand()

Y = rand('seed')


He expects Y=L, and indeed, should find that under the scheme I proposed above.

deego <deego>
Thu 08 Mar 2018 01:57:26 AM UTC, comment #14: 


>> foo = rand("seed") => Implements K = 1e20*rand(), implement rand("state", K), and output K. [or minor variants thereupon.]


I had meant to say K = round(1e20*rand())



deego <deego>
Thu 08 Mar 2018 01:55:53 AM UTC, comment #13: 

Rik, you wrote:

> It will be hard to switch Octave away from these random number generators. Matlab has been trying for years to get people to abandon them. See http://www.mathworks.com/help/matlab/math/updating-your-random-number-generator-syntax.html.





It is understandable that people want their legacy code to continue to work (and sometimes are actually unable to change it, even if they want.). So, good to maintain "seed".

But, they don't REALLY care what underlying mechanism rand() uses.  The output is random, of course!

How about these (or similar) drop-in replacements. These have the advantage that octave never silently switches behind my back to "seed", while continuing to satisfy code that uses seed:

rand("seed", X) => actually imlements rand("state", X).

foo = rand("seed") =>  Implements K = 1e20*rand(), implement rand("state",  K), and output K.  [or minor variants thereupon.]

rand("seed", "reset") does rand("state", "reset").

[remember that 'state' actually expects a vector, but continues to work if a scalar is supplied. For reproducibility, we simply need to supply the same scalar again.]

For the person that uses 'seed', his code continues to do what he expects: He can reset the state. He can query a seed, and he can supply the same seed back for reproducible behavior. All the while, we are actually using state.

deego <deego>
Thu 08 Mar 2018 01:46:54 AM UTC, comment #12: 

Mike

Thanks for your detailed analysis of "bug1". Also, for noticing the (much more glaring) startup state problem. I have added rand("state", "reset") to my init file out of caution.


deego <deego>
Thu 08 Mar 2018 12:29:12 AM UTC, comment #11: 

It will be hard to switch Octave away from these random number generators.  Matlab has been trying for years to get people to abandon them.  See http://www.mathworks.com/help/matlab/math/updating-your-random-number-generator-syntax.html.

They still support the old generators because some people have code that works, but which can't be modified.

Rik <rik5>
Group administrator
Thu 08 Mar 2018 12:22:42 AM UTC, comment #10: 

@deego - compare these examples



for i = 1:100; rand ('seed', 'reset'); disp (rand ('seed')); endfor

for i = 1:100; rand ('seed', 'reset'); disp (rand ('seed')); pause (.1); endfor

for i = 1:100; rand ('seed', 'reset'); disp (rand ('seed')); pause (1); endfor


I think that should explain why you are seeing what you are, but I have no idea how to convey such a mess in a clear way in the documentation.

Mike Miller <mtmiller>
Group Member
Thu 08 Mar 2018 12:03:34 AM UTC, comment #9: 

@rik - I'll file that as a new bug.

@deego - I do see some inconsistency with respect to your "bug 1".

TL;DR - I think rand('seed', 'reset') depends on the current time stamp. So your test of resetting one generator and comparing the seed of the other probably depends on round-off error with the current system clock. I haven't looked into why, but this is probably a Matlab compatibility thing with the old generator.

First, I have verified that randn('seed', 'reset') does switch the rand generator to the old algorithm. But your observation that one generator resets the other but not the other way around isn't quite correct. It seems to differ with each run.


>> aa = rand ('seed'); randn ('seed', 'reset'); bb = rand ('seed');
>> isequal (aa, bb)
ans = 0      ## But sometimes this returns true, try it multiple times
>> aa = randn ('seed'); rand ('seed', 'reset'); bb = randn ('seed');
>> isequal (aa, bb)
ans = 0      ## But sometimes this returns true, try it multiple times


Now let's test it in loops with different timings:


>> function r = t () aa = rand ('seed'); randn ('seed', 'reset'); bb = rand ('seed'); r = isequal (aa, bb); endfunction
>> sum = 0; for i = 1:1000; sum += t (); endfor; sum
sum =  999      ## That's bad!
>> sum = 0; for i = 1:20; pause (1.1); sum += t (); endfor; sum
sum = 0


So I think this "bug 1" only occurs when the previous rand('seed', 'reset') was done within the same second as the next one.

I don't think this is documented anywhere, but again, this would be great to eliminate by phasing out the "old" algorithm completely.

Mike Miller <mtmiller>
Group Member
Wed 07 Mar 2018 11:45:13 PM UTC, comment #8: 

@Mike: That's a potentially nasty bug.  I could see a test suite starting Octave, querying a random seed and placing it in a file along with the results of tests.  This should allow anyone doing debugging to pull up the file, restore Octave to the original state, and then reproduce a failing test. 

If Octave is lying about the true value of the "state" vector then that isn't going to work.

Rik <rik5>
Group administrator
Wed 07 Mar 2018 11:42:42 PM UTC, comment #7: 

In some glorious future, it would be even better to switch the random number generators to depend on the C++ runtime library.  As of C++11 the random number functionality that is part of the language is extremely well thought out.

Rik <rik5>
Group administrator
Wed 07 Mar 2018 11:37:33 PM UTC, comment #6: 

I'm not sure this is worth pursuing, but I have noticed this odd behavior when testing the rand functions:


## just started Octave (with --norc), rand generator is in some indeterminate state
>> init_state = rand ('state');
>> init_seed = rand ('seed');
>> x = rand (100); ## generate a few numbers to shift the state
>> isequal (init_seed, rand ('seed'))  ## NOT using 'seed'
ans = 1
>> isequal (init_state, rand ('state'))
## now let's try to go back to the initial state
>> rand ('state', init_state);
>> y = rand (100);
>> isequal (x, y)
ans = 0
## uh oh, they should be equal, let's try again
>> rand ('state', init_state);
>> z = rand (100);
>> isequal (x, z)
ans = 0
>> isequal (y, z)
ans = 1


As far as I can tell, setting the rand state does reset the rand generator to a known state and the same sequence is produced. However, when Octave first starts, the vector that is returned from rand('state') is not the actual initial state.

Mike Miller <mtmiller>
Group Member
Wed 07 Mar 2018 11:27:54 PM UTC, comment #5: 

Since I spent a lot of time looking at rand* recently with respect to the test suite, I agree with Rik's assessment that, while confusing, the functions do seem to be behaving as I expect.

There are the "old" vs "new" algorithms, which switch silently under the hood when a user sets "state" or "seed", but not when simply asking for the current value of "state" or "seed". And there is no way to query which algorithm is currently in use.

The only new bit that Rik states that I didn't know is under the old algorithms, the rande/randg/randp functions depend on the others.

In summary, never use "seed".

Aside: I would really like to know why this test passes only when setting "seed" to get a random array from the old algorithm:

https://hg.savannah.gnu.org/hgweb/octave/file/2d68dc548561/libinterp/corefcn/conv2.cc#l274

> It would be nicer to simply remove the old random number generator functionality entirely.


Absolutely, as well as move the state out of global program state and into a class so that multiple independent random number generators can be instantiated.

Mike Miller <mtmiller>
Group Member
Wed 07 Mar 2018 11:25:38 PM UTC, comment #4: 

Also, I do agree: One way to resolve this might be to remove the seed-related code entirely.

For one, it is twice as slow.

Two, the switching behavior is inconsistent (these bug reports), and less than properly documented.

Three, there's no easy way to query for current behavior in use. And, calls to weird codes can render the behavior secretly switched.
(For example, the old randint from forge switches the behavior to seed, and the user is none the wiser.)

deego <deego>
Wed 07 Mar 2018 11:22:29 PM UTC, comment #3: 


Rik,

Thanks for your comments.
Bug 2:  understood your point. In this case, the behavior should  at least become part of the documentation of rand.

Bug 1: Ok, so, resetting one seed resets the other. This should be documented in rand(n)'s helpstring. But, wait, I don't think even that solves it. It's still weird: (Like I reported), Resetting rand's seed sets that of randn but not the other way.  This behavior is still not even self-consistent. (And, if for some reason, /that/ is the behavior we want, it should definitely be made part of the docstring.)

deego <deego>
Wed 07 Mar 2018 11:09:53 PM UTC, comment #2: 

Bug 1 is also probably explained by the first sentence where it says it is not possible to mix the two generators.  From the normal rand documentation


Older versions of Octave used a different random number generator.
The new generator is used by default as it is significantly faster
than the old generator, and produces random numbers with a
significantly longer cycle time.  However, in some circumstances it
might be desirable to obtain the same random sequences as produced
by the old generators.  To do this the keyword "seed" is used to
specify that the old generators should be used, as in

     rand ("seed", val)

which sets the seed of the generator to VAL.  The seed of the
generator can be queried with

     s = rand ("seed")

However, it should be noted that querying the seed will not cause
'rand' to use the old generators, only setting the seed will.  To
cause 'rand' to once again use the new generators, the keyword
"state" should be used to reset the state of the 'rand'.


Examining the test code


aa = sum(rand('seed')); randn('seed', 'reset'); bb = sum(rand('seed')); assert(aa!=bb);


The first fragment queries the rand "seed" which does not switch generators.  The second fragment executes a reset of the randn "seed" and also switches the random number generators from the new to the old.  When you switch generators entirely there is no guarantee on the behavior of the seeds.  I guess that could be documented.  It would be nicer to simply remove the old random number generator functionality entirely.

Rik <rik5>
Group administrator
Wed 07 Mar 2018 10:55:07 PM UTC, comment #1: 

At least for bug 2, see https://www.gnu.org/software/octave/doc/interpreter/Special-Utility-Matrices.html#Special-Utility-Matrices.


The generators operate in the new or old style together, it is not possible to mix the two. Initializing any generator with "state" or "seed" causes the others to switch to the same style for future calls.

The state of each generator is independent and calls to different generators can be interleaved without affecting the final result. For example,

rand ("state", [11, 22, 33]);
randn ("state", [44, 55, 66]);
u = rand (100, 1);
n = randn (100, 1);

and

rand ("state", [11, 22, 33]);
randn ("state", [44, 55, 66]);
u = zeros (100, 1);
n = zeros (100, 1);
for i = 1:100
  u(i) = rand ();
  n(i) = randn ();
end

produce equivalent results. When the generators are initialized in the old style with "seed" only rand and randn are independent, because the old rande, randg and randp generators make calls to rand and randn.

The generators are initialized with random states at start-up, so that the sequences of random numbers are not the same each time you run Octave.7 If you really do need to reproduce a sequence of numbers exactly, you can set the state or seed to a specific value.

If invoked without arguments, rand and randn return a single element of a random sequence.

The original rand and randn functions use Fortran code from RANLIB, a library of Fortran routines for random number generation, compiled by Barry W. Brown and James Lovato of the Department of Biomathematics at The University of Texas, M.D. Anderson Cancer Center, Houston, TX 77030.



Rik <rik5>
Group administrator
Wed 07 Mar 2018 06:57:45 AM UTC, original submission:  

version ## =>  4.2.1

## BUG 1: rand and randn have independent seeds, yet resetting rand's seed resets that of randn. If not a bug, this is at least undocumented behavior.
## Curiously it doesn't work the other way.

aa = sum(rand('state')); randn('state', 'reset'); bb = sum(rand('state')); assert(aa==bb);

aa = sum(randn('state')); rand('state', 'reset'); bb = sum(randn('state')); assert(aa==bb);


aa = sum(rand('seed')); randn('seed', 'reset'); bb = sum(rand('seed')); assert(aa!=bb);
## ^^^ Bug 1a: Unexpected seed change for rand upon resetting randn.

aa = sum(randn('seed')); rand('seed', 'reset'); bb = sum(randn('seed')); assert(aa==bb);
## ^^^ But, resetting rand does NOT reset randn's seed!

## Bug 1 is not even self-consistent. Resetting randn's seed resets that of rand, but not the other way. Furthermore, resetting state of one does not affect that of another.






## BUG 2: Switching rand to use seed switches randn's behavior as well, and vice versa!
## Again, if a bug, it is undocumented behavior.

## Now, we can't ask rand or randn to report if they are using seed or state, so we will verify this indirectly:


## First ensure we are using state.
randn('state', 'reset');
rand('state', 'reset');

## Sanity checks and warmup:

## These variables keep track of the state.
nt = sum(randn('state')) ## n for randn and t for state =>nt
rt = sum(rand('state'))
nd = randn("seed")    ## n for randn and d for seed => nd
rd = rand("seed")

## We see that all 4 are different above.

## Now, let's verify that drawing a randn changes the randn state, but NOT the other 3.
randn(1);
assert(nt != sum(randn('state')));
assert(rt == sum(rand('state')));
assert(nd == sum(randn('seed')));
assert(rd == sum(rand('seed')));


nt = sum(randn('state')); ## Set nt back to the current randn state.

## Similarly, drawing a rand changes the rand state but not the other 3

rand(1);
assert(nt == sum(randn('state')));
assert(rt != sum(rand('state')));
assert(nd == sum(randn('seed')));
assert(rd == sum(rand('seed')));

rt = sum(rand("state")); ## update since we drew a rand.

## Now, switch randn to seed, and draw a randn.

randn('seed','reset'); 
nd = randn('seed');  ## update since we had reset above.
rd = rand('seed'); ## This update is needed because of Bug 1 above. rand seed has also changed.
randn(1); ##

## Verify that nd has changed, but not others.
assert(nt == sum(randn('state')));
assert(rt == sum(rand('state')));
assert(nd != sum(randn('seed')));
assert(rd == sum(rand('seed')));

nd = sum(randn('seed')); ## Update nd since state changed.

## We'd expect rand to still use state, but it seems that rand has switched to using seed as well!
rand(1);  ## draw a rand. You'd expect rt to change, but it is rd that changes! (((bug 2a)))

## Verify that rd has changed, but not the others
assert(nt == sum(randn('state')));
assert(rt == sum(rand('state')));
assert(nd == sum(randn('seed')));
assert(rd != sum(rand('seed')));

rd = sum(rand('seed'));

## Now, switch rand back to state, and update the corresponding variable
rand('state' ,'reset'); rt = sum(rand('state'));

## For things to be somewhat consistent with bug 2a, we expect randn to switch back to state as well. Verify:
randn(1); ## draw a randn.


## Verify (((bug 2b))) that nt has changed, but not
assert(nt != sum(randn('state')));
assert(rt == sum(rand('state')));
assert(nd == sum(randn('seed')));
assert(rd == sum(rand('seed')));

## We successfully reach here.
 


deego <deego>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by deego (Submitted the item)
  • -email is unavailable- added by deego
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-03-27 mtmiller CategoryOctave Function Documentation
    2018-03-10 rik5 StatusConfirmed In Progress
        SummarySeveral rand/randn bugs with state/seed behavior. Document more fully rand state/seed behavior
    2018-03-08 mtmiller Severity3 - Normal 2 - Minor
        Item GroupIncorrect Result Documentation
        StatusNone Confirmed
    2018-03-07 deego Carbon-Copy- Added -email is unavailable-

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code