Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

mvnrnd bug when run in parfor?

Asked by Kirk on 18 Jun 2013

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

0 Comments

Kirk

Products

No products are associated with this question.

4 Answers

Answer by Shashank Prasanna on 18 Jun 2013

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

0 Comments

Shashank Prasanna
Answer by Kirk on 19 Jun 2013

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.

0 Comments

Kirk
Answer by Edric Ellis on 20 Jun 2013
Edited by Edric Ellis on 20 Jun 2013

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.

4 Comments

Edric Ellis on 21 Jun 2013

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 on 24 Jun 2013

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 on 25 Jun 2013

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...

Edric Ellis
Answer by Peter Perkins on 24 Jun 2013

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.

2 Comments

Kirk on 24 Jun 2013

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 on 27 Jun 2013

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".

Peter Perkins

Contact us