mvnrnd bug when run in parfor?

We are having trouble running the multivariate normal random number generator (mvnrnd) inside a 'parfor' loop. For example, when we run mvrnd within a 'parfor' loop with 5 iterations (i=5) and mvrnd randomly calculate the variable A. We end up with:

A =    [2, 7, -3, 12]
    [2, 7, -3, 12]
    [2, 7, -3, 12]
    [2, 7, -3, 12]
    [2, 7, -3, 12]

with each iteration of i being identical. Re-running the code will generate a new random set of A, however, like above each iteration i will be identical:

A =    [1, 3, 5, -2]
    [1, 3, 5, -2]
    [1, 3, 5, -2]
    [1, 3, 5, -2]
    [1, 3, 5, -2]

However, if we run the same bit of code within a 'for' loop, A, always changes randomly as we'd expect.

A =    [2, 7, -5, 12]
    [3, -9, 2, 12]
    [17, 1, -1, 1]
    [2, -2, 8, 11]
    [5, 7, -3, 19]

This strikes us as rather odd behavior. Can anyone explain why mvrnd doesn't randomly generate A within a parlor loop like it does in a for loop?

Thanks in advance

Answer 1

[Shashank Prasanna]

This is expected behavior since the random numbers have been seeded to the same initial conditions at the start of each workers. You can change the seed as follows:

parfor i = 1:4
rng(i)
A(i,:) = mvnrnd(zeros(1,4),diag(ones(4,1)));
end

RNG will seed it to a different seed based on the loop index

Answer 2

[Kirk]

Thanks for the response Shashank,

And thank you for pointing out rng(). However, upon closer inspection of our code, I misspoke earlier. The replication of numbers does not occer within each iteration of i, but rather when comparing different runs (i.e. different parfor loops). For example:

%RUN1 parfor i=1:4 A=mvnrnd(mu,sigma) end;

yields 
i=1, A=2
i=2, A=7
i=3, A=-3
i=4, A=12

%Run2 A=mvnrnd(mu,sigma) mvnrnd end;

yields 
i=1, A=2
i=2, A=7
i=3, A=-3
i=4, A=12

Whereas changing parfor() to for() (and running a single worker of course), we get the following:

%RUN1 for i=1:4 A=mvnrnd(mu,sigma) end;

yields 
i=1, A=3
i=2, A=1
i=3, A=-3
i=4, A=11

%Run2 for i=1:4 A=mvnrnd(mu,sigma) end;

yields 
i=1, A=5
i=2, A=-2
i=3, A=-1
i=4, A=9

So... If I am following your point, it appears that parfor starts each individual run with the same seed, whereas for() re-seeds each run? Aslo I can replicate the parfor() behavior by adding (as you suggested) rng(i), to the for loop. However, adding rng('suffle') to the parfor loop does NOT produce the expected "new random numbers" upon each run.

Answer 3

[Edric Ellis]

As Shashank points out, this is expected behaviour. Matlabpool workers have their random number generators set up carefully at the start of the session to get precisely the same results, in an exactly analogous way to what MATLAB does. For example, when I start MATLAB R2013a, the first random number generated by RAND is always the same:

>> rand
ans =
    0.8147

Likewise, immediately after opening a matlabpool, the first random number is always precisely the same:

>> matlabpool('open', 'local', 1); spmd, rand, end
Starting matlabpool using the 'local' profile ... connected to 1 workers.
Lab 1: 
  ans =
      0.3246

If you do not close the matlabpool, then the random number stream will continue to progress and different numbers will be generated. If you want to choose a different starting point, then setting the seed directly on the workers is not the best idea.

There's more detail in this answer which tells you about options for setting up the random number generation on the workers.

[Kirk]

Sorry to beat this to death, and this all makes good sense and is how my machine behaves as well (ans = 0.8147 and ans = 0.3246). But... I am talking about re-running the random number generator (after those first instances), and getting back the same and... over and over. This does not seem correct.

For example, if I create mu as: mu(12,4)

and sigma as: sigma(4,4,12)

Then open 4 workers as: matlabpool(4)

Then run mvnrnd as:

parfor i = 1:64
  A = mvnrnd(mu,sigma);
end

Then examine A,and re-run for parfor loop:

parfor i = 1:64
  A = mvnrnd(mu,sigma);
end

A does not change. Is that really what you mean? If I am just a dull witted numbskull, feel free to point it out I am beginning to seriously question it myself!

[Edric Ellis]

Sorry, I'm not that experienced with MVNRND - could you post a completely standalone script that reproduces the problem. Note that in the example above, your PARFOR loops are not actually going to modify 'A' - but perhaps you just missed something. You need something like

parfor ...
  A(i, :) = mvrnd(...);
end

otherwise PARFOR treats 'A' as a loop temporary.

[Kirk]

Thanks Edric. Sorry for the delay, but I've having a bugger of time getting the three dimensional covariance matrix (sigma) to build in a script. As of now I am still getting an error:

"Error using mvnrnd (line 153) Each page of SIGMA must be a symmetric positive semi-definite matrix."

Even though as far as I can tell, the sigma I am building in the script below is exactly the same as the sigma I am using (which runs with the semi-definite" error. I have even tried to script a simple random (4,4,14) marix, but I continue to get the semi-definite matrix error.

All I can think to do at this point is to paste in the script and hopefully there is enough information so that you can see what I'm trying to do?

%define number of iterations
j = 500;

%create 4x4x12 covarianve matrix sigma
s(:,:,1) = 1.0e+03 * [0.0113    0.0157    0.0052   -0.0009;
    0.0157    0.0241   -0.0199   -0.0008;
    0.0052   -0.0199    2.6108    0.0017;
   -0.0009   -0.0008    0.0017    0.0015];

s(:,:,2) = 1.0e+03 * [0.0126    0.0163   -0.1059    0.0002;
    0.0163    0.0232   -0.1781    0.0006;
   -0.1059   -0.1781    2.8378   -0.0166;
    0.0002    0.0006   -0.0166    0.0003];

s(:,:,3) = 1.0e+03 * [0.0107    0.0091    0.0578    0.0032;
    0.0091    0.0116   -0.0240    0.0054;
    0.0578   -0.0240    5.2269   -0.0509;
    0.0032    0.0054   -0.0509    0.0048];

s(:,:,4) = 1.0e+04 * [0.0011    0.0007    0.0105    0.0002;
    0.0007    0.0005    0.0003    0.0003;
    0.0105    0.0003    1.4131   -0.0126;
    0.0002    0.0003   -0.0126    0.0007];

s(:,:,5) = 1.0e+04 * [0.0005    0.0003   -0.0021    0.0002;
    0.0003    0.0004   -0.0072    0.0004;
   -0.0021   -0.0072    1.5599   -0.0089;
    0.0002    0.0004   -0.0089    0.0022];

s(:,:,6) = 1.0e+04 * [0.0003    0.0002   -0.0003   -0.0002;
    0.0002    0.0003   -0.0024    0.0001;
   -0.0003   -0.0024    1.8137   -0.0327;
   -0.0002    0.0001   -0.0327    0.0026];

s(:,:,7) = 1.0e+04 * [0.0005    0.0004    0.0007    0.0000;
    0.0004    0.0003    0.0022    0.0001;
    0.0007    0.0022    1.1774   -0.0076;
    0.0000    0.0001   -0.0076    0.0008];

s(:,:,8) = 1.0e+03 * [0.0047    0.0035    0.0255   -0.0041;
    0.0035    0.0034   -0.0019   -0.0025;
    0.0255   -0.0019    6.8238   -0.1159;
   -0.0041   -0.0025   -0.1159    0.0181];

s(:,:,9) = 1.0e+03 * [0.0032    0.0023    0.0712    0.0012;
    0.0023    0.0025    0.0181    0.0037;
    0.0712    0.0181    7.6450   -0.0224;
    0.0012    0.0037   -0.0224    0.0323];

s(:,:,10) = 1.0e+03 * [0.0074    0.0036    0.1689    0.0014;
    0.0036    0.0031    0.0419    0.0039;
    0.1689    0.0419    8.0279   -0.1315;
    0.0014    0.0039   -0.1315    0.0175];

s(:,:,11) = 1.0e+03 * [0.0202    0.0150    0.2215    0.0003;
    0.0150    0.0130    0.1175    0.0024;
    0.2215    0.1175    9.6524   -0.1475;
    0.0003    0.0024   -0.1475    0.0095];

s(:,:,12) = 1.0e+03 * [0.0147    0.0139    0.0892   -0.0011;
    0.0139    0.0146    0.0597   -0.0007;
    0.0892    0.0597    4.1552   -0.0130;
   -0.0011   -0.0007   -0.0130    0.0012];

%create climate matix mu
m = [-6.20430107526882,-16.9964157706093,260.027432408801,1.53162000000000;-3.16044061302682,-14.7408319649699,332.993786420239,1.45796000000000;3.47491039426523,-8.00896057347670,473.478101409595,3.05816000000000;13.3595679012346,-0.0648148148148152,659.211580205352,4.04283333333333;20.0909842845327,6.14971050454921,727.234638870377,10.0623076923077;25.3233618233618,12.0754985754986,794.972939707542,11.9321384615385;27.8031660692951,14.5086618876941,892.013383265893,9.12918333333333;26.7810862972153,13.4008822718500,888.608283147339,8.04593846153846;22.4351851851852,9.20524691358025,808.803439710225,9.29216666666667;13.4587813620072,2.47963506028022,576.051870644398,8.44434545454545;2.98611111111111,-6.20833333333333,427.142796993370,3.28506666666667;-4.26360377973281,-13.5777126099707,284.962097288256,1.39007272727273;];

tic
parfor i = 1:j
A = mvnrnd(m,s);
end
toc

[Edric Ellis]

Your PARFOR loop still assigns to the whole of "A", thus it is considered to be a "temporary" variable - you need to assign to slices of "A" to get your results...

Answer 4

[Peter Perkins]

Your cov matrices that are pretty badly badly conditioned, but th reason why M<VNRND refuses to continue is that one of them is not even positive semi-definite:

>> eig(s(:,:,7))
ans =
     -0.23708
       7.0797
       8.6212
        11775

You can't expect MVNRND to accept a matrix that is not a valid cov matrix, and that one isn't.

[Kirk]

Good catch. Thank you Peter. The difference must be in the number of digits displayed in the command window. When I check the variable 'sigma' (which is what I built the variable 's' from, the eigen function reports back that is is a positive semi-definite. The variable 's', on the other hand (as you correctly pointed out) is not.

>> sigma(:,:,7)

ans =

     1.0e+04 *

      0.0005    0.0004    0.0007    0.0000
      0.0004    0.0003    0.0022    0.0001
      0.0007    0.0022    1.1774   -0.0076
      0.0000    0.0001   -0.0076    0.0008

>> eig(sigma(:,:,7))

ans =

     1.0e+04 *

      0.0000
      0.0007
      0.0009
      1.1775

>> s(:,:,7)

ans =

     1.0e+04 *

      0.0005    0.0004    0.0007         0
      0.0004    0.0003    0.0022    0.0001
      0.0007    0.0022    1.1774   -0.0076
           0    0.0001   -0.0076    0.0008

>> eig(s(:,:,7))

ans =

     1.0e+04 *

     -0.0000
      0.0007
      0.0009
      1.1775

This component of my problem seems to be as issue of copy and pasting the display from the Variable window into a script. i.e. not getting all of the significant digits that matlab stored.

Checking with class(), they are both (sigma and s) doubles.

>> class(sigma(:,:,7))

ans =

double

>> class(s(:,:,7))

ans =

double

Ideas?

[Peter Perkins]

I'm not sure what you're asking. These are pretty badly conditioned matrices, and sigma(:,:,7) may be so badly conditioned that mvnrnd will still refuse to work on it even though eig believes that all the eigen values are non-negative. If you're asking how to fix that, that's a hrd question worthy of it's own thread. If you're just asking how to copy and paste and not lose precision, you probably just want to use "format long" or "format long g".