MATLAB Answers

## mvnrnd bug when run in parfor?

Asked by Kirk

### Kirk (view profile)

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

## Products

No products are associated with this question.

## 4 Answers

### Shashank Prasanna (view profile)

Answer by Shashank Prasanna

### Shashank Prasanna (view profile)

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

Answer by Kirk

### Kirk (view profile)

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.

### Edric Ellis (view profile)

Answer by Edric Ellis

### Edric Ellis (view profile)

on 20 Jun 2013
Edited by Edric Ellis

### Edric Ellis (view profile)

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.

Edric Ellis

### Edric Ellis (view profile)

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

### Kirk (view profile)

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

### Edric Ellis (view profile)

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

### Peter Perkins (view profile)

Answer by Peter Perkins

### Peter Perkins (view profile)

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.

Kirk

### Kirk (view profile)

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

### Peter Perkins (view profile)

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

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

### Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

New to MATLAB?