Problems Generating Scrambled Quasi-Monte Carlo Numbers in Parfor Loops

I am having an issue generating scrambled quasi-monte carlo numbers in parfor loops.

The problem is that when I generate multiple sets of these numbers within a parfor loop, the numbers in each set end up being identical. I am including a very simple example below.

D = 3;
M = 1000;
numbers = cell(1,4);

mystream = qrandstream(scramble(sobolset(D),'MatousekAffineOwen'));
myfun = @(x) qrand(mystream,x);

parfor i = 1:4
    numbers{i} = myfun(M);
end

To demonstrate the issue, after running this code, the numbers in numbers{1}, numbers{2}, numbers{3} and numbers{4} are identical as:

>>numbers{1}(1:3,:)
ans =
          0.76          0.05          0.77
          0.33          0.96          0.23
          0.60          0.72          0.52

>> numbers{2}(1:3,:)
ans =
          0.76          0.05          0.77
          0.33          0.96          0.23
          0.60          0.72          0.52

I'm wondering whether anyone can think of a fix for this issue. I figure there must be something that I can do, as the problem does not occur when I use a normal random stream of numbers.

I should mention that it will not be possible for me to exploit something like the 'Skip' or 'Leap' properties of Quasi-Random Number streams. The reason is that I use the snippet of code above in a larger MATLAB program that I run in parallel...

Answer 1

[Peter Perkins]

Berk, presumably you have executed something like a matlabpool command prior to running this code.

Your code creates a randomly-scrambled quasi-random stream. The parfor loop then, in effect, pushes copies of that one stream out to each worker. Same thing on all four workers, but separate copies. Then each worker draws from its own stream, separately. Obviously, because each worker is drawing from identical copies of the same stream, they get the same numbers.

I can't say what the right solution is, because I don't know what you're doing. But somehow you need to give the different workers different streams. Perhaps by using different scrambles (though I confess that I don't know the statistical implications of that), perhaps by using skips or leaps.

I don't understand your comment about skips and leaps, because those are tailor made for what you seem to need to do. By way of analogy, when you're using pseudorandom numbers in parallel, you want to either use statistically independent streams, or use some leapfrog/blocking/partitioning scheme on one stream. Skips and leaps are equivalent to the latter, and really are the right way to do what you want.

Two other things:

1) The anonymous function is not necessary in this simple code, though you may have other reasons for using it in your real code. In this example,

     numbers{i} = qrand(mystream,M);

would have worked just as well.

2) Your statement about pseudorandom numbers is not correct. This code

     % mystream = qrandstream(scramble(sobolset(D),'MatousekAffineOwen'));
     mystream = RandStream('twister','seed',0);

     parfor i = 1:4
         numbers{i} = rand(mystream,M,D);
     end
     numbers{:}

demonstrates that, and also demonstrates how in many cases you can write the same code for quasi- or pseudo-random streams. Presumably you were doing this:

     parfor i = 1:4
         numbers{i} = rand(M,D);
     end

which uses the global streams on each worker, which MATLAB sets up differently for each worker, to avoid the whole problem of identical values.

Hope this helps.

Answer 2

[Oleg Komarov]

Give a look at:

http://blogs.mathworks.com/loren/2011/07/07/simpler-control-of-random-number-generation-in-matlab/

Especially comments 8, and 9.

Also http://blogs.mathworks.com/loren/2008/11/13/new-ways-with-random-numbers-part-ii/