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

Thread Subject:
Randperm is not random! How to REALLY randomise each randperm() call?

Subject: Randperm is not random! How to REALLY randomise each randperm() call?

From: David Painter

Date: 23 Feb, 2010 11:46:04

Message: 1 of 8

Hi,

I've finally realised randperm is subject to the same "limitations" as the rand functions as it contains rand (obviously) and so operates on the default random stream. :'-( This may seem obvious, but it should be documented explicitly in the help file. No warning/hint is given in Matlab R2008B.

Wondering why it is not advisable to seed the random by the block more than once per Matlab session?

s = RandStream.create('mt19937ar','seed',sum(100*clock));
RandStream.setDefaultStream(s);
randperm(10)
s = RandStream.create('mt19937ar','seed',sum(100*clock));
RandStream.setDefaultStream(s);
randperm(10)

I want my randperms to produce different results on each run. How can I do it?

Thanks,

David

Subject: Randperm is not random! How to REALLY randomise each randperm() call?

From: John D'Errico

Date: 23 Feb, 2010 11:55:06

Message: 2 of 8

"David Painter" <david.ross.painter@gmail.com> wrote in message <hm0f5s$ct6$1@fred.mathworks.com>...
> Hi,
>
> I've finally realised randperm is subject to the same "limitations" as the rand functions as it contains rand (obviously) and so operates on the default random stream. :'-( This may seem obvious, but it should be documented explicitly in the help file. No warning/hint is given in Matlab R2008B.
>
> Wondering why it is not advisable to seed the random by the block more than once per Matlab session?
>
> s = RandStream.create('mt19937ar','seed',sum(100*clock));
> RandStream.setDefaultStream(s);
> randperm(10)
> s = RandStream.create('mt19937ar','seed',sum(100*clock));
> RandStream.setDefaultStream(s);
> randperm(10)
>
> I want my randperms to produce different results on each run. How can I do it?

Do NOT reseed the random generator multiple times
within one session. This will make the random
number generator LESS random, not more so.

And randperm DOES return different results on each
call.

randperm(10)
ans =
     6 3 7 8 5 1 2 4 9 10

randperm(10)
ans =
     6 1 7 4 9 5 8 3 10 2

randperm(10)
ans =
     2 10 8 9 1 5 7 6 3 4
 
John

Subject: Randperm is not random! How to REALLY randomise each randperm()

From: Peter Perkins

Date: 23 Feb, 2010 14:13:44

Message: 3 of 8

On 2/23/2010 6:46 AM, David Painter wrote:

> I want my randperms to produce different results on each run. How can I
> do it?

David, that depends on what you mean by "run".

1) If you mean, "each successive time I call it", then as John points out, it does that.

2) If you mean, "every time I start up MATLAB", then the right thing to do depends on what you are really doing.

a) It may be that you are running a Monte-Carlo simulation that depends on RANDPERM (and possibly other random number calls), and you want to run that simulation in several MATLAB sessions and think of the different runs as statistically independent so that you can combine the results, as if you had run one long simulation. Using a seed based on clock (for example) is one way to do that. However, you won't be able to go back and reproduce a run unless you save the seed. That may or may not be important to you, but another option would be to use one of the inherently parallel generator algorithms (e.g., 'mrg32k3a') and use a known stream index in each run.

b) It may be that you are running a MC simulation repeatedly, but not trying to combine the different runs as statistically independent, but are just uncomfortable having the same set of random values used each time. For example, you might be doing MC integration on different functions (you probably aren't doing that, but just as an illustration). You are free to feel uncomfortable, and you can do the same thing as in (a), but I might point out that if your conclusions from a MC simulation depends on the specific sequence of random values actually used, then you probably need to run more MC iterations.

c) You may just want to see a different output from RANDPERM each time you call it, regardless of whether or not you just started up MATLAB. That's fine (though I've never felt that need), go ahead and use clock. But there's no reason I know of to do it more than once per session, and you may just want to put something in your startup file.

Hope this helps.

Subject: Randperm is not random! How to REALLY randomise each randperm() call?

From: Jan Simon

Date: 23 Feb, 2010 15:30:26

Message: 4 of 8

Dear David!

> s = RandStream.create('mt19937ar','seed',sum(100*clock));
> RandStream.setDefaultStream(s);
> randperm(10)

> s = RandStream.create('mt19937ar','seed',sum(100*clock));
> RandStream.setDefaultStream(s);
> randperm(10)

If you believe that mt19937ar is a valid random number generator, there is no need to seed it twice. If you don't believe that mt19937ar produces valid random numbers, use another generator.
I personally do trust Matlab's random number generator and I'm convinced that RANDPERM is not biased.

As far as I understood the only drawback of RANDPERM is the stable SORT: If a random value appears twice, the order of these elements is kept. Fortunately the chance to get two equal values in a RAND vector can be neglected.
Another drawback of RANDPERM is the speed for large arrays, which is limited by the sorting. The Fisher-Yates shuffle is simply faster (60% of RANPERM time for 10 elements, 46% of RANDPERM for 1E6 elements):

function X = RandInplacePermute(X)
for i = 1:numel(X)
  w = ceil(rand * i);
  t = X(w);
  X(w) = X(i);
  X(i) = t;
end

The main difference appears for large arrays, because RANDPERM needs 4 temporary vectors of the same size (see source of RANDPERM):
1. RAND(1, n)
2. SORT creates one data copy for the ignored sorted output and one for the indices.
3. X = X(randperm(numel(X)) creates a further copy. I think this could help to re-use the memory:
  X(:) = X(randperm(numel(X))

Anyhow, I think RandInplacePermute is too trivial to be published on the FEX - different opinions?

Kind regards, Jan

Subject: Randperm is not random! How to REALLY randomise each randperm() call?

From: David Painter

Date: 23 Feb, 2010 23:19:07

Message: 5 of 8

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <hm0fmq$ft2$1@fred.mathworks.com>...
> "David Painter" <david.ross.painter@gmail.com> wrote in message <hm0f5s$ct6$1@fred.mathworks.com>...
> > Hi,
> >
> > I've finally realised randperm is subject to the same "limitations" as the rand functions as it contains rand (obviously) and so operates on the default random stream. :'-( This may seem obvious, but it should be documented explicitly in the help file. No warning/hint is given in Matlab R2008B.
> >
> > Wondering why it is not advisable to seed the random by the block more than once per Matlab session?
> >
> > s = RandStream.create('mt19937ar','seed',sum(100*clock));
> > RandStream.setDefaultStream(s);
> > randperm(10)
> > s = RandStream.create('mt19937ar','seed',sum(100*clock));
> > RandStream.setDefaultStream(s);
> > randperm(10)
> >
> > I want my randperms to produce different results on each run. How can I do it?
>
> Do NOT reseed the random generator multiple times
> within one session. This will make the random
> number generator LESS random, not more so.
>
> And randperm DOES return different results on each
> call.
>
> randperm(10)
> ans =
> 6 3 7 8 5 1 2 4 9 10
>
> randperm(10)
> ans =
> 6 1 7 4 9 5 8 3 10 2
>
> randperm(10)
> ans =
> 2 10 8 9 1 5 7 6 3 4
>
> John


Thanks to everyone for their helpful replies. I am guilty of posting while tired!

John, your're right. Randperm produces different results on each call within a session (open Matlab, exit Matlab). The issue is how to have different sequences for different sessions. One call to seed the stream based on the clock is sufficient because clock is likely to be different on each call to:

s = RandStream.create('mt19937ar','seed',sum(100*clock))

Is that right?

Subject: Randperm is not random! How to REALLY randomise each randperm()

From: Peter Perkins

Date: 24 Feb, 2010 16:07:01

Message: 6 of 8

On 2/23/2010 6:19 PM, David Painter wrote:

> One call to seed the stream
> based on the clock is sufficient because clock is likely to be different
> on each call to:
>
> s = RandStream.create('mt19937ar','seed',sum(100*clock))
>
> Is that right?

It is, although

    s = RandStream('mt19937ar','seed',sum(100*clock))

s marginally shorter to type (since you're creating only one random stream), and don't forget that you need to make your new stream the default random number stream,

   RandStream.setDefaultStream(s)

You may find that simply reseeding the existing default stream

    reset(RandStream.getDefaultStream,sum(100*clock))

at the beginning of each MATLAB session (e.g., in your startup file) is even easier.

Subject: Randperm is not random! How to REALLY randomise each randperm()

From: Greg Heath

Date: 25 Feb, 2010 09:20:26

Message: 7 of 8

On Feb 24, 11:07 am, Peter Perkins
<Peter.Perk...@MathRemoveThisWorks.com> wrote:
> On 2/23/2010 6:19 PM, David Painter wrote:
>
> > One call to seed the stream
> > based on the clock is sufficient because clock is likely to be different
> > on each call to:
>
> > s = RandStream.create('mt19937ar','seed',sum(100*clock))
>
> > Is that right?
>
> It is, although
>
> s = RandStream('mt19937ar','seed',sum(100*clock))
>
> s marginally shorter to type (since you're creating only one random stream), and don't forget that you need to make your new stream the default random number stream,
>
> RandStream.setDefaultStream(s)
>
> You may find that simply reseeding the existing default stream
>
> reset(RandStream.getDefaultStream,sum(100*clock))
>
> at the beginning of each MATLAB session (e.g., in your startup file) is even easier.

For today (Feb 25, 2010)

min(sum(clock)) = 2010 + 2 + 25 + 0 + 0 + 0 = 2037
max(sum(clock)) = 2010 + 2 + 25 + 23 + 59 + 59 = 2178

Therefore, Jseed = 100*sum(clock) only allows

100*(23+59+59) = 14,100 integer seeds per day

whereas the total number of available integer seeds is
4,294,967,296 over the interval

[ 0, 2^32-1] = [0 4,294,967,295]

Perhaps

Jseed = 3e7*sum(clock-clock0)

where clock0 = 2010+2+25 (depends on the data)

( or one of many other variations) might yield a better
diversity of random streams.

Hope this helps.

Greg

Subject: Randperm is not random! How to REALLY randomise each randperm()

From: Peter Perkins

Date: 25 Feb, 2010 13:09:25

Message: 8 of 8

On 2/25/2010 4:20 AM, Greg Heath wrote:

> Therefore, Jseed = 100*sum(clock) only allows
>
> 100*(23+59+59) = 14,100 integer seeds per day
>
> whereas the total number of available integer seeds is
> 4,294,967,296 over the interval
>
> [ 0, 2^32-1] = [0 4,294,967,295]
>
> Perhaps
>
> Jseed = 3e7*sum(clock-clock0)
>
> where clock0 = 2010+2+25 (depends on the data)
>
> ( or one of many other variations) might yield a better
> diversity of random streams.

Greg, this might be OK, but if you really need that different streams of random numbers, I would instead recommend using one of the generator types that are specifically intended to support parallel random number generation or substreams. Using clock is probably more like something you'd want to do at the beginning of each MATLAB session just to get something different than last time. If you need thousands of seeds, I'd guess you're running hugely parallel simulations, and you should really be using a tool made for that job. So for example,

s = RandStream('mrg32k3a');
s.Substream = [[pick some number between 1 and 2^51]]
RandStream.setDefaultStream(s);

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us