matlabrandomparallel-processingparfor

matlab random number generation in parfor loop


I want to generate same normal random numbers for loop and parfor loop. Following MATLAB documentation I tried three different method:

Method 1: Using rng(seed, 'twister')

N = 1;

for ind = 1:10
    rng(ind, 'twister')
    samps(ind) = normrnd(0,1,N,1);
end

plot(samps); hold on

parfor ind = 1:10
    rng(ind, 'twister')
    samps(ind) = normrnd(0,1,N,1);
end

scatter(1:length(samps), samps)
legend({'using for loop', 'using parfor loop'})

By using rng twister method

Method 2: Using RandStream Method

N = 1;
 
sc = parallel.pool.Constant(RandStream('Threefry'));
for ind = 1:10
    stream = sc.Value;
    stream.Substream = ind;
    samps(ind) = normrnd(0,1,N,1);
end

plot(samps); hold on

sc = parallel.pool.Constant(RandStream('Threefry'));
parfor ind = 1:10
       stream = sc.Value;
       stream.Substream = ind;
       samps(ind) = normrnd(0,1,N,1);
end

scatter(1:length(samps), samps)
legend({'using for loop', 'using parfor loop'})

Using Threefry stream-substream method

Method 3: Using another RandStream method

N = 1;
 
stream = RandStream('mrg32k3a');
for ind = 1:10
    stream.Substream = ind;
    samps(ind) = normrnd(0,1,N,1);
end

plot(samps); hold on

stream = RandStream('mrg32k3a');
parfor ind = 1:10
       set(stream,'Substream',ind);
       samps(ind) = normrnd(0,1,N,1);
end

scatter(1:length(samps), samps)
legend({'using for loop', 'using parfor loop'})

Another randStream method

My question is why the last two methods are not working. If I understood MATLAB documentation correctly, they should have worked. Also, is there anything wrong with using rng(seed, 'twister') method with different seeds to produce statistically independent samples?

Edit: Here is the link to the MATLAB documentation that I used.


Solution

  • As @Cris points out, normrnd uses the current global stream. Here's how to make your 3rd case match between client and workers - you need to set the global stream explicitly. (Note that I'm reverting the global stream immediately after using it too).

    N = 1;
    
    stream = RandStream('mrg32k3a');
    for ind = 1:10
        stream.Substream = ind;
        prev = RandStream.setGlobalStream(stream);
        samps(ind) = normrnd(0,1,N,1);
        RandStream.setGlobalStream(prev);
    end
    
    stream = RandStream('mrg32k3a');
    parfor ind = 1:10
        set(stream,'Substream',ind);
        prev = RandStream.setGlobalStream(stream);
        pf_samps(ind) = normrnd(0,1,N,1);
        RandStream.setGlobalStream(prev);
    end
    

    With this code, samps and pf_samps are identical.

    Also, you should not use "twister" in parallel, since the random samples you obtain by setting the seed in this way will not have as good qualities as if you use a "parallel" generator like mrg32k3a. (I'm no expert, but my understanding is that the twister "seed" values are not able to select statistically-independent sequences of samples in the way that streams/substreams can)

    In recent versions of MATLAB, faster generators that support multiple streams are available. See the doc https://www.mathworks.com/help/matlab/math/creating-and-controlling-a-random-number-stream.html#brvku_2 . Philox and Threefry in particular are good parallel generators.