generate an independent seed with a fixed number

Hi to all,

I am simulating a simple problem in Matlab and optimizing it through OptQuest. The function in Matlab is as follows:

function saddleproduct(infile, outfile, replication)
% Read input file
inp = readtable(infile);
disp(infile);
disp(outfile);
disp(replication);
disp(inp);
%Table indices are 1-based
var1 = inp{1,2};
var2 = inp{2,2};
%variance for simulation: Jalali et al. 
W_expected(1) = var1 + var2;
W_expected(2) = 1.5-var1-2*var2-(0.5)*sin(2*pi*(var1^2-2*var2));
W_expected(3) = var1^2+var2^2-1.5;
variances = [(0.45*W_expected(1)+0.3)^2 (0.45*W_expected(2)+1.15)^2                  (0.45*W_expected(3)+0.98)^2]; 
%simulate at current point
w0 = var1 + var2 + normrnd(0, sqrt(variances(1)));
w1 = 1.5 - var1 - 2*var2 -0.5*sin(2*pi*(var1^2 - 2*var2)) + ...
    normrnd(0, sqrt(variances(2)));
w2 = var1^2 + var2^2 -1.5 + normrnd(0, sqrt(variances(3)));
%outputs
product=w0;
sum = w1; 
quotient=w2;
% Output results to file
output = table({'func-product';'func-sum';'func-quotient'}, {product;sum;quotient});
disp(output);
writetable(output, outfile, 'WriteVariableNames', false);
end

The function reads inputs, gets replication number from OptQuest, and writes output to a file. In simulation (i.e., normrnd), Matlab always uses the default seed (rng('default)). I have to change the code in such a way that Matlab uses replication to determine an independent seed for each replication, where replications at the same input has the same number of replications.

How should I use the replication number to create a seed independent for each replication at the same input? Is clock the only option for it? What is your recommendation?

Answer 1

[Peter Perkins]

Ebru, you are doing parallel simulations and presumably want to combine the results under the assumption that those results are (pseudo)independent across replications.

Don't use seeds for that. You can, but there are better ways. Take a look at the parallel generators that support streams and substreams, and use one of those. These strategies are well documented, see

https://www.mathworks.com/help/matlab/math/generate-random-numbers-that-are-repeatable.html

https://www.mathworks.com/help/matlab/math/generate-random-numbers-that-are-different.html

[Ebru Angun]

Thank you Peter Perkins.

Best regards

Ebru

[Ebru Angun]

Hi Peter,

I have further questions about your comment above.

It seems that I should generate a prn stream, say,

s=RandStream('mrg32k3a');

But this has to be done once at the beginning of the run, and should not be reset each time it returns to the function. So, I can save this into a .mat file;

save seed.mat s

and make s a global variable.

I will be using 12 parallel processors.

1- How can I use 12 substreams from this stream?

2- Before Matlab returns to OptQuest, I have to save the current state (which has to be different than the initial state) as the next seed to be used when OptQuest returns to Matlab. How can I do that?

Best regards

Ebru

[Ebru Angun]

Hi Peter,

I detailed a bit more what I have to do:

1- Generate one stream with substreams:

s=RandStream('mrg32k3a');

Then, save the state as the current_state.mat.

2- Each time OptQuest calls Matlab, Matlab has to load the file current_state.mat.

3- Then, Matlab has to assign 12 substreams to 12 parallel processors, which work independently to simulate.

4- At the end of the simulation before returning to OptQuest, the 12 by 1 vector state should be written in current_state.mat, to be used at the next call.

How can I assign 12 substreams to 12 local processors?

Should the state vector be written in current_state.mat before returning to OptQuest? Will this way give independent pseudo-random numbers at each call?

Do I have to assign the current state to a global variable within Matlab function?

Thank you

Ebru

[Peter Perkins]

The whole point of stream numbers and substream numbers is that you don't need to save any states. Each stream or substream is easily referenced by its number. Read those doc sections.

Answer 2

[Jan]

rng('default')

is the default seed at the start of the Matlab session. All subsequent requesrts of random numbers move the seed accordingly. There is no need to seed the generator repeatedly.

If you wat the initial seed to be random also, use

rng('shuffle')

If you want the seed to depend of the loop index:

rng(replication)

assuming that replication is a positive integer.

But I'm confused by the question:

"How should I use the replication number to create a seed independent for each replication at the same input?" - with useing a defined seed, the rng is not "independent", but dependent.

"Is clock the only option for it?" - Why? If the rng should be independent, do not set further seeds.

[Ebru Angun]

Thank you Jan, but this did not work. I tried it already. What happened when I used rng('default') is that each time OptQuest asks a new simulation replication, Matlab always uses the default seed, and of course produces the same output. Hence, using rng does not work. For the same input combination, the replication number is the same, but this number has to produce an independent pseudorandom number each time a new replication is to be made.

To better explain, suppose the current input is (1,1), and the output at that input is w_0 = 1 + 1 + e_0, where e_0 is normal(0, sigma). Now suppose that, for input (1,1), OptQuest decides to have 10 replications. If I use rng('default') to start the seed at its default value, for each replication at (1,1), Matlab generates, say, e_0 = 0.164 for all 10 replications.

In such a case, does it make sense to use

seed_mine = floor(replication * current_seconds_in_clock);

and use it in generating stream

s = RandStream('mt19937ar','Seed',seed_mine,'NormalTransform','Polar');

or maybe in

rng(seed_mine)?

The replications will not be independent since they will be using the same number of replications; yet, hopefully, they will produce different output values. Of course, it may be better to use

seed_mine = floor(current_seconds_in_clock);

So this is the confusing question. I hope that it is explained better now.

Thanks again

Ebru

[Ebru Angun]

I must also add that Matlab and OptQuest (from SimWrapper) are two different programs that communicate. OptQuest passes the information of inputs and replication numbers to Matlab, and Matlab runs the simulation and returns the output values to OptQuest.

[Ebru Angun]

I added one more thing to my questions above that I tried in the mean time. Each time OptQuest returns to Matlab:

c = clock;

seed = floor( c(6) * replication );

rng( seed );

This produces different values for output at the same input.

Is there any way to make it better?

Best regards

Ebru

[Jan]

The 6th element of clock has a low entropy only. The mutliplication by replication is not useful. Remember that only the integer part of the seed is used. Then providing the milliseconds in the 6th element of clock matter, if the replication value is > 1000.

Use rng('shuffle') to get a random shuffle. But using different streams as suggested by Peter is much better.

[Ebru Angun]

Thank you Jan.