# Should mvnrnd Always Advance the State of the Global Stream

3 views (last 30 days)
Paul on 26 May 2020
Edited: Paul on 27 May 2020
Consider the following:
>> mu=[1 1]; Sigma=eye(2);
rng('default')
preu1 = rand(1,3);
n1 = mvnrnd(mu,Sigma);
u1 = rand(1,3);
rng('default')
preu2 = rand(1,3);
n2 = mvnrnd(mu,Sigma);
u2 = rand(1,3);
rng('default')
n3 = mvnrnd(mu,0*Sigma);
u3 = rand(1,3);
>> [isequal(u1,u2) isequal(u2,u3) isequal(u3,preu2)]
ans =
1×3 logical array
1 0 1
Apparently mvnrnd doesn't actually call randn if it detects that the input covariance is zero. That may be good for efficiency, but is it the best behavior for repeatability? This behavior seems to be a contradiction to the general direction of how to manage the Global Stream to reproduce results. doc mvnrnd is silent on how it handles this special case.

Steven Lord on 26 May 2020
There's no requirement that these two lines of code that call the same function with different inputs need follow identical code paths.
myfun(A, B)
myfun(A, C)
The lines of code in your example look similar to the human eye, but I could rewrite them to be less similar but equivalent to what you wrote.
x1 = mvnrnd(mu, eye(2))
x2 = mvnrnd(mu, zeros(2))

Paul on 26 May 2020
I realize there is not a general requirement that a function with two inputs always follow identical code paths. And I realize that doc mvnrnd doesn't specifically guarantee that it will call randn for all inputs. But my question is if it would be better that mvnrnd always does advance the state of the Global Stream in all use cases both for usability and to meet user expectations. From the documentation:
"By controlling the default random number stream and its state, you can control how the RNGs in Statistics and Machine Learning Toolbox software generate random values. For example, to reproduce the same sequence of values from an RNG, you can save and restore the default stream's state, or reset the default stream. For details on managing the default random number stream, see Managing the Global Stream(MATLAB)."
These statements, along with the content in Managing the Global Stream, certainly suggest a reasonable interpretation that mvnrnd always advances the state of the Global Stream regardless of the inputs.
But apparently it does not, which is fine from a computational perspective. And I'm willing to accept for now that it's not a bug. But my question is: should it advance the Global Stream, from a usability perspective?
As it stands, consider an application with a sequence of calls to various RNGs, including mvnrnd, and assume that anything that controls the dimensions of the outputs of the RNGs is fixed. However, the values of the input parameters that are used to call the RNG functions are provided by the user. Would it not be reasonable to expect that executing the sequence twice, and executing rng('default') prior each execution, result in the same outputs for those RNGs that used the same parameters in both runs? As shown in the example, if on the second run the user sets Sigma to 0 for a call to mvnrnd, then all rand* calls downcode from that call to mvnrnd that use the Global Stream will be affected as well. As it stands, it appears that the only mitigation is for the programmer to know of and catch the special cases to advance the Global Stream manually. Is that what we want? Is there another option?
Steven Lord on 26 May 2020
But my question is: should it advance the Global Stream, from a usability perspective?
I would argue no.
Let's play two hands of poker by calling a hypothetical poker() function with the number of players and the name of the game to play. [Inside poker() we start off with a deck of cards shuffled with randperm.]
Our first hand is
poker(2, 'Texas hold em')
Since there are two of us playing (the first input to the poker() function), we need a total of nine cards: my two hole cards, your two hole cards, and the five common cards for the flop, turn, and river. I'm assuming neither of us chooses to fold at any point.
Now a friend stops by and decides to join. I reset the global generator used by randperm and call
poker(3, 'Texas hold em')
That means we deal eleven cards: my two hole cards, your two hole cards, our friend's two hole cards, and the five common cards. [You and I are cheating, since we know all but the final two common cards from our previous hand, but since this is a hypothetical example on MATLAB Answers that's fine.]
Do you expect that after the first call to poker() we should have discarded two cards (advanced the global stream) expecting / predicting that our friend will stop by and decide to play, so the state of the deck is the same after each hand?
If you want to be certain the stream is in a certain state after the mvnrnd call, set it to that state after the call.
Paul on 26 May 2020
Steve,
I've reread your example several times and I'm still not seeing how it's on point. Among other things, adding that third player changed the requirements on the sequence that is generated for your simulated poker game.
Continuing with the poker theme, suppose you're simulating the last five tables at the WSOP before the next cut and so each table plays a hand simultaneously. So you might have a function that includes code that looks like this to simulate what happens.
for n = 1:ntrials
for table = 1:5
shuffle_and_deal(mu(table),Sigma(:,:,table))
end
end
Assume that the suffle_and_deal function calls mvnrnd with those input arguments. So the result of the simulation depends on the state of the Global Stream at the start of the simulation. Run the simulation once. Reset the Global Stream to whatever its state was at the start of the first run. Now set Sigma(:,:,1) = zeros(2). Rerun the simulation. Would you expect the output of the simulation be any different for tables 2-5? Why should changing only the parameter of the distribution for table 1 have any effect on the other tables? After all, nothing changed for those tables.
If you want to be certain the stream is in a certain state after the mvnrnd call, set it to that state after the call
I'm surprised you said this. It seems counter to the concept of the random number architecture. I'm pretty sure that I've even seen somehwere in the doc that advises to not do things like this, because it shouldn't be needed. It also seems intractable for code that has a sequence of, for example, 1000 calls to mvnrnd. Set the state of the Global Stream after each call? And do so in such a way that you get the same answer for the entire sequence if the state of the Global Stream as the start of the calling squence is reset? I'd rather just a) wrap mvnrnd in my own function and catch the special case, or b) not use mvnrnd at all and just roll my own.
Also, normrnd from the same toolbox has a different behavior when called with standard deviation of zero:
>> s1=RandStream.getGlobalStream.State;x=mvnrnd(1,0);s2=RandStream.getGlobalStream.State;isequal(s1,s2)
ans =
logical
1
>> s1=RandStream.getGlobalStream.State;x=normrnd(1,0);s2=RandStream.getGlobalStream.State;isequal(s1,s2)
ans =
logical
0
Based on the doc for the toolbox, and the doc for Managing the Global Stream, and how normrnd works, and how the other non-toolbox rng functions like randn and rand work, what would be your expectation as to the state of the Global Stream after a call to mvnrnd with Sigma = 0?
Paul

R2019a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!