Consistently generating same random sequence with for and parfor loop

Dear all,

I encountered a problem, when I was running monte-carlo simulation. In order to speed up the simulation, I decided to change 'for' loop to 'parfor'. In order to debug and to later check the details of trajectory, I fixed the seed for random number generator within the loop. To my surprise, though the seeds for random number generator in 'for' and 'parfor' loop are same, the random number sequence are completely different. This creates a problem, since I cannot check the trajectories in 'parfor'.

% sample code with 'for' loop
rng(5)
seedForSimulation = randi(100,1,10);
randomNrs = NaN(1,10);
Seeds = NaN(1,10);
for i = 1:10
    rng(SeedForSimulation(1,i));
    Seeds(i) = SeedForSimulation(1,i);
    randomNrs(i) = randi(10);
end   

display(Seeds)
display(randomNrs)

The result that I get is :

Seeds = 23 88 21 92 49 62 77 52 30 19
randomNrs = 6 7 1 9 4 1 10 9 7 1

Whereas the same code with 'parfor'

% code
rng(5)
seedForSimulation = randi(100,1,10);
randomNrs = NaN(1,10);
Seeds = NaN(1,10);
parfor i = 1:10
    rng(SeedForSimulation(1,i));
    Seeds(i) = SeedForSimulation(1,i);
    randomNrs(i) = randi(10);
end   

display(Seeds)
display(randomNrs)

The result that I get :

Seeds = 23 88 21 92 49 62 77 52 30 19
randomNrs = 5 4 9 5 4 4 10 7 8 4

Though the seeds are same within 'for' and 'parfor', I get different random numbers. Can someone explain me what is Matlab doing ? How can I circumvent the problem.

Answer 1

[Robert]

You can produce the same random numbers in a parfor loop if you explicitly set the random number generator type to be the same in the for and parfor loops.

You can do this when you set the seed or before the loop.

parfor i = 1:10
  rng(seedForSimulation(1, i), 'twister');
  Seeds(1, i) = seedForSimulation(1, i);
  randomNrs(1, i) = randi(10);
end

This will return the exact same result if you replace parfor with for. You could also use:

spmd % run once on each parallel worker
  rng(0, 'twister')
end
parfor i = 1:10
  rng(seedForSimulation(1, i)); % no need to specify again
  Seeds(1, i) = seedForSimulation(1, i);
  randomNrs(1, i) = randi(10);
end

This will also return the same result as would the for loop.

[Robert]

My test code:

% code
rng(5, 'combRecursive')
seedForSimulation = randi(10,1,10);
randomNrs = NaN(2,10);
Seeds = NaN(2,10);
for i = 1:10
    rng(seedForSimulation(1,i));
    Seeds(1, i) = seedForSimulation(1,i);
    randomNrs(1, i) = randi(10);
end
parfor i = 1:10
    rng(seedForSimulation(1,i));
    Seeds(2, i) = seedForSimulation(1,i);
    randomNrs(2, i) = randi(10);
end
display(Seeds)
display(randomNrs)

[sulav duwal]

Dear Robert, thank you for the attempt ! It did not work for me. I found out that one has to explicitly mention the generator for example :

parfor i = 1:10
  rng(seedForSimulation(1,i),'twister')
  %Some codes
end

and

for i = 1:10
  rng(seedForSimulation(1,i),'twister')
  %Some codes
end

It seems that the default random number generator in 'parfor' is 'combRecursive', whereas the default random number generator otherwise is 'twister'.

[Robert]

Good catch. I guess I had set my rng source types while playing around with this code and the side effects persisted even when I took 'combRecursive' out of the for and parfor loops.

I will edit my answer to reflect our conclusion: the parfor and for loops need to be explicitly set to use the same random source (no matter which you choose).

[Robert]

My new test code

rng default
seedForSimulation = randi(10,1,10);

rng_source = 'twister'; % or 'combRecursive' or any other
% not needed if using 'twister' since that is default in base workspace
rng(0, rng_source)
% not needed if 'combRecursive' since that is default for parallel workers
spmd
  rng(0, rng_source)
end

randomNrs = NaN(2,10);
Seeds = NaN(2,10);
for i = 1:10
  rng(seedForSimulation(1,i));
  Seeds(1, i) = seedForSimulation(1,i);
  randomNrs(1, i) = randi(10);
end
parfor i = 1:10
  rng(seedForSimulation(1,i));
  Seeds(2, i) = seedForSimulation(1,i);
  randomNrs(2, i) = randi(10);
end

display(Seeds)
display(randomNrs)

whose output is

Seeds =
   9    10     2    10     7     1     3     6    10    10
   9    10     2    10     7     1     3     6    10    10
randomNrs =
   1     8     5     8     1     5     6     9     8     8
   1     8     5     8     1     5     6     9     8     8

[sulav duwal]

Thanks Robert !!