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

 
 

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

Submitted by:  Dave Goel <deego>
Submitted on:  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

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project Administrator
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.

Dave Goel <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())

Dave Goel <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.

Dave Goel <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.

Dave Goel <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>
Project Administrator
Thu 08 Mar 2018 12:22:42 AM UTC, comment #10:

@deego - compare these examples

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>
Project Administrator
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.

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

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>
Project Administrator
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>
Project 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>
Project 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:

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>
Project Administrator
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>
Project Administrator
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.)

Dave Goel <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.)

Dave Goel <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

Examining the test code

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>
Project 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.

Rik <rik5>
Project 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.

Dave Goel <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
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project 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.3