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, Randi, and Shuffle

Subject: Randperm, Randi, and Shuffle

From: Anthony Hopf

Date: 2 Nov, 2010 22:02:04

Message: 1 of 19

I am trying to replace the following code:

[rix riy riz] = ind2sub(sampbox,randi(myrandistr,prod(sampbox),SPtemp,1));

where sampbox can get large ( [1100 200 200]) and SPtemp can be 125e3. This is not the case all of the time, but it is definitely bogs down the process.

I want to replace my current code because it seems to show a bias when I plot out the locations it is picking within the 3d space.

Originally I wanted to replace this piece of code with randperm, but it was incredibly slow. Luckily I found Jan's addition to FEX called Shuttle. It seemed to really fit the bill because I would specify the number of NOut within his function. Orders of magnitude faster than randperm for sure.

Problem is now Shuffle is fast, 0.32sec for the above extreme, but not as fast as randi, 0.014 for the same case.

This is the code I am using with Shuffle:

        randpts = Shuffle(prod(sampbox),'index',SPtemp);
        [rix riy riz] = ind2sub(sampbox,randpts);

Thank you

Anthony

Subject: Randperm, Randi, and Shuffle

From: Peter Perkins

Date: 3 Nov, 2010 13:30:11

Message: 2 of 19

On 11/2/2010 6:02 PM, Anthony Hopf wrote:
> Problem is now Shuffle is fast, 0.32sec for the above extreme, but not
> as fast as randi, 0.014 for the same case.

You are generating a random sample drawn uniformly from the integers
1:n. You need to decide if you want that sample drawn with replacement
(randi), or without replacement (the shuffle function on the FEX.

Hope this helps.

Subject: Randperm, Randi, and Shuffle

From: Bruno Luong

Date: 3 Nov, 2010 13:47:04

Message: 3 of 19

May be you could try this:

a = randi(myrandistr,prod(sampbox),round(SPtemp*1.1),1);
a = unique(a);
a(SPtemp+1:end)=[];

Bruno

Subject: Randperm, Randi, and Shuffle

From: Jan Simon

Date: 3 Nov, 2010 15:02:04

Message: 4 of 19

Dear Anthony,

> Problem is now Shuffle is fast, 0.32sec for the above extreme, but not as fast as randi, 0.014 for the same case.

As said before: The RANDI values are not necessarily unique, while the values from SHUFFLE are. Using UNIQUE(RANDI) will definitely help, but sorting 125.000 values will cost some time also. And unfortuanetly UNQIUE(RANDI) will include the same long sorting as RANDPERM does.

Kind regards, Jan

Subject: Randperm, Randi, and Shuffle

From: Bruno Luong

Date: 3 Nov, 2010 18:01:05

Message: 5 of 19

"Jan Simon" <matlab.THIS_YEAR@nMINUSsimon.de> wrote in message <iarthc$dbb$1@fred.mathworks.com>...
> Dear Anthony,
>
> > Problem is now Shuffle is fast, 0.32sec for the above extreme, but not as fast as randi, 0.014 for the same case.
>
> As said before: The RANDI values are not necessarily unique, while the values from SHUFFLE are. Using UNIQUE(RANDI) will definitely help, but sorting 125.000 values will cost some time also. And unfortuanetly UNQIUE(RANDI) will include the same long sorting as RANDPERM does.

According to my test, UNIQUE will take 5-8 times longer than RAND. OP seems to indicates a ratio of 20 between SHUFFTE and RAND (not tested by me).

So UNIQUE might some amount, and certainly not spectacular.

Bruno

Subject: Randperm, Randi, and Shuffle

From: Jan Simon

Date: 4 Nov, 2010 02:05:04

Message: 6 of 19

Dear Bruno,

You are right: The sorting in UNIQUE belongs to the smaller data set of 125e3 values, while RANDPERM would have to operate on [1:prod(sampbox)].

Therefore Bruno's code is faster than SHUFFLE inspite of the sorting:

> a = randi(myrandistr, prod(sampbox), round(SPtemp*1.1), 1);
> a = unique(a);
> a(SPtemp +1:end) = [];

But it is biased also, because the reply of UNIQUE is sorted. But lukily SHUFFLE works efficient here:
  a = randi(myrandistr, prod(sampbox), round(SPtemp*1.1), 1);
  a = Shuffle(unique(a));
  a(SPtemp +1:end) = [];

How safe is 1.1*SPtemp? At least you should include a test of UNIQUE(a) has enough elements.

I have not implemented an equivalent method in SHUFFLE, because the computing time grows extremly, if the number of chosen values is gets near to the number of all elements. With other words: It is very hard to draw N unique random numbers from [1:(N+tiny)] using the RANDI approach.
But for [1:N*10] the number of collision would be small. Then this approach might be fast:
% Pseudocode, to be implemented in C:
% Choose N values from [1:M]:
  x = zeros(1, N);
  r = zeros(1, N);
  for i = 1:N
     ready = false;
     while ~ready
        a = randi(M); % Draw
        if isempty(find(x(1:i) == a))
          r(i) = a; % Unsorted list as output
          x(1:i) = sort([x(1:i-1), a]); % Sorted list for FIND
          ready = true;
       end
  end

This could use a binary search to find [a] in the list [x] of already drawn numbers, and reuse the found index for sorting.
This competes with the simple method of creating Y=ONES(1,M,'uint8') at first, draw a number X, if Y(X) is 1, accept X and set Y(X) to zero.
Finally there are three methods to draw N out of [1:M]:
1. Fisher-Yates Shuffle: efficient for M=N+(tiny value)
2. The dull draw and check method: efficient for M=N+(fair value)
3. Draw and sort (the pseudocode above): efficient for M=N+(giantic value)
Implementing these 3 methods is easy, but finding a stable method to choose the best method seems to be cruel: It depends on M, N, the available RAM and the size or the processor cache.
I will not start to work on this before 2011, so use Bruno's suggestion.

BTW: Did you say, that RANDI is biased?
This would be worth a new thread!

Going to bed now, expecting N*M*tiny+giantic nightmares, Jan

Subject: Randperm, Randi, and Shuffle

From: Bruno Luong

Date: 4 Nov, 2010 06:20:04

Message: 7 of 19

>
> But it is biased also, because the reply of UNIQUE is sorted. But lukily SHUFFLE works efficient here:

This is good remark. To get unbiased selection, here is a more complicated code that select random points after calling unique. (it's needed my file MAXK on FEX)

n = 1100*200*200;
m = 125e3;

a = ceil(n*rand(1,round(1.1*m)));
a = unique(a);
l = length(a);
i = rand(1,l);
[~, loc] = maxk(i, l-m);
a(loc)=[];

This add a marginal additional runtime.

When m << n, the probability of the event (l < m) is tiny (I'm too lazy to calculate that). My guess is it never happens in practice.

Bruno

Subject: Randperm, Randi, and Shuffle

From: Jan Simon

Date: 4 Nov, 2010 11:05:05

Message: 8 of 19

Dear Bruno,

> When m << n, the probability of the event (l < m) is tiny (I'm too lazy to calculate that). My guess is it never happens in practice.

I'm deeply surprised. I do not know any personal detail of you, but I have the strong impression, that you are definitely not lazy.

Jan

Subject: Randperm, Randi, and Shuffle

From: Derek O'Connor

Date: 9 Dec, 2010 01:49:04

Message: 9 of 19

This is slightly off-topic but may be of interest:

function p = GenPerm(n,m);
%
% Generate a random cyclic permutation using Sattolo's
% modification of Knuth's Algorithm P, Section 3.4.2
% TAOCP, Vol 2, 2nd Ed.
%
% m = 0 is the Fisher-Yates-Knuth Shuffle Algorithm P
% m = 1 is Sattolo's Algorithm
%
% Prodhinger : On the analysis of an algorithm to generate
% a random cyclic permutation.
% http://math.sun.ac.za/~hproding/abstract/abs_161.htm
%
% Wilson : Random and exhaustive generation of permutations
% and cycles
% http://arxiv.org/abs/math/0702753
% Note: Wilson uses zero-based arrays.
%
% Derek O'Connor, 8 Dec 2010. derekroconnor@eircom.net

p = 1:n; % Identity permutation

for k = n:-1:2
   r = 1 + floor(rand*(k-m)); % random integer between 1 and k-m
   t = p(r);
   p(r) = p(k); % Swap(p(r),p(k))
   p(k) = t;
end;
return; % GenCycPerm

Subject: Randperm, Randi, and Shuffle

From: Jos (10584)

Date: 9 Dec, 2010 14:56:06

Message: 10 of 19

"Derek O'Connor" <derekroconnor@eircom.net> wrote in message <idpcig$7a8$1@fred.mathworks.com>...
> This is slightly off-topic but may be of interest:
>
> function p = GenPerm(n,m);
> %
> % Generate a random cyclic permutation using Sattolo's
> % modification of Knuth's Algorithm P, Section 3.4.2
> % TAOCP, Vol 2, 2nd Ed.
> %
> % m = 0 is the Fisher-Yates-Knuth Shuffle Algorithm P
> % m = 1 is Sattolo's Algorithm
> %
> % Prodhinger : On the analysis of an algorithm to generate
> % a random cyclic permutation.
> % http://math.sun.ac.za/~hproding/abstract/abs_161.htm
> %
> % Wilson : Random and exhaustive generation of permutations
> % and cycles
> % http://arxiv.org/abs/math/0702753
> % Note: Wilson uses zero-based arrays.
> %
> % Derek O'Connor, 8 Dec 2010. derekroconnor@eircom.net
>
> p = 1:n; % Identity permutation
>
> for k = n:-1:2
> r = 1 + floor(rand*(k-m)); % random integer between 1 and k-m
> t = p(r);
> p(r) = p(k); % Swap(p(r),p(k))
> p(k) = t;
> end;
> return; % GenCycPerm

You might also like to take a look at (the code) of RANDPERMFULL
http://www.mathworks.com/matlabcentral/fileexchange/23257-randpermfull

Jos

Subject: Randperm, Randi, and Shuffle

From: Derek O'Connor

Date: 9 Dec, 2010 18:25:23

Message: 11 of 19

"Jos (10584) " <#10584@fileexchange.com> wrote in message <idqqm6$rdv$1@fred.mathworks.com>...
>
> You might also like to take a look at (the code) of RANDPERMFULL
> http://www.mathworks.com/matlabcentral/fileexchange/23257-randpermfull
>
> Jos

Dear Jos,

I ran and profiled both GenPerm and randpermfull with n = 10^7.

GenPerm : 8 secs.

  0.08 1 20 p = 1:n;
                    1 22 for k = n:-1:2
  2.05 9999999 23 r = 1 + floor(rand*(k-m));
  2.25 9999999 24 t = p(r);
  1.07 9999999 25 p(r) = p(k);
  1.26 9999999 26 p(k) = t;
  1.25 9999999 27 end;


randpermfull : 24 secs

 0.06 1 48 R = 1:N ;
 0.40 1 61 K = R + ceil((N - R) .* rand(1,N)) ;
                     1 67 for j=1:N-1,
22.20 9999999 69 R([j K(j)]) = R([K(j) j]) ;
  1.37 9999999 73 end

You can see that most of the time is spent in the swap operation.
Also randpermfull uses twice the memory (R & K) that GenPerm uses (p).

System: Dell Precision 690, Intel 2xQuad-Core E5345 @ 2.33GHz
16GB RAM, Windows7 64-bit Professional.
MATLAB Version 7.6.0.324 (R2008a)

Derek O'Connor

Subject: Randperm, Randi, and Shuffle

From: Jos (10584)

Date: 10 Dec, 2010 15:47:11

Message: 12 of 19

"Derek O'Connor" <derekroconnor@eircom.net> wrote in message <idr6uj$n9r$1@fred.mathworks.com>...
> "Jos (10584) " <#10584@fileexchange.com> wrote in message <idqqm6$rdv$1@fred.mathworks.com>...
> >
> > You might also like to take a look at (the code) of RANDPERMFULL
> > http://www.mathworks.com/matlabcentral/fileexchange/23257-randpermfull
> >
> > Jos
>
> Dear Jos,
>
> I ran and profiled both GenPerm and randpermfull with n = 10^7.
>
...

Thanks Derek. Your swapping method is indeed more efficient. I decided to modify the code of RANDPERMFULL, the new version should be up in a few days.

Jos

Subject: Randperm, Randi, and Shuffle

From: Derek O'Connor

Date: 13 Dec, 2010 01:10:28

Message: 13 of 19

"Jos (10584) " <#10584@fileexchange.com> wrote in message <idti1v$55l$1@fred.mathworks.com>...
> "Derek O'Connor" <derekroconnor@eircom.net> wrote in message <idr6uj$n9r$1@fred.mathworks.com>...
> > "Jos (10584) " <#10584@fileexchange.com> wrote in message <idqqm6$rdv$1@fred.mathworks.com>...
> > >
> > > You might also like to take a look at (the code) of RANDPERMFULL
> > > http://www.mathworks.com/matlabcentral/fileexchange/23257-randpermfull
> > >
> > > Jos
> >
> > Dear Jos,
> >
> > I ran and profiled both GenPerm and randpermfull with n = 10^7.
> >
> ...
>
> Thanks Derek. Your swapping method is indeed more efficient. I decided to modify the code of RANDPERMFULL, the new version should be up in a few days.
>
> Jos

Jos,
The profiler can add a lot to the run times so I ran GenPerm and randfullperm
without the profiler to see the actual performance:

GenPerm(10^7) = 2 secs , GenPerm(10^8) = 24 secs

randfullperm(10^7) = 20 secs , randfullperm(10^8) = 211 secs


Derek O'Connor

Subject: Randperm, Randi, and Shuffle

From: Jan Simon

Date: 13 Dec, 2010 02:38:04

Message: 14 of 19

Dear Derek, dear Jos,

The new randpermfull is much faster. I get factor 11 for a 1e6 vector under Matlab 2009a. Nice!

If speed really matters (does it ever?), it might be worth to mention that the C-Mex SHUFFLE is three times faster than GenPerm and randpermfull. Because it uses the Knuth-Shuffle, it is easy to modify it for the Sattolo-algorithm.

Kind regards, Jan

Subject: Randperm, Randi, and Shuffle

From: Derek O'Connor

Date: 14 Dec, 2010 20:27:05

Message: 15 of 19

"Jan Simon" <matlab.THIS_YEAR@nMINUSsimon.de> wrote in message <ie40uc$hcp$1@fred.mathworks.com>...
> Dear Derek, dear Jos,
>
> The new randpermfull is much faster. I get factor 11 for a 1e6 vector under Matlab 2009a. Nice!
>
> If speed really matters (does it ever?), it might be worth to mention that the C-Mex SHUFFLE is three times faster than GenPerm and randpermfull. Because it uses the Knuth-Shuffle, it is easy to modify it for the Sattolo-algorithm.
>
> Kind regards, Jan
--------------------------------
Dear Jan,
I'm still fiddling around with the Fisher-Yates Shuffle Algorithm
and tried the swap operation that Jos and others have used.

As you can see below, the results are disastrous -- although the function looks nicer.

I have two questions:

1. What in Matlab is causing this disastrous slow-down? Perhaps some
Matlab aficionado can enlighten us.

2. Why do most, if not all, languages lack an interchange operation X <=> Y?

Most good textbooks use it at the pseudo-code level. I wonder how many bugs
have been caused by the 3-statement version.
I know that most students get it wrong.


Derek O'Connor



     function p = GenRandPerm(n)
       p = 1:n;
       for k = n:-1:2
           r = 1 + floor(k*rand);
1 t = p(r); % used this statement
2 p(r) = p(k); % p([k r]) = p([r k]);
3 p(k) = t; % instead of 1,2,3
       end;
       return;

     n [ 1e4 1e5 1e6 1e7 ]
     R [ 37 23 24 14 ]

     R = 1 - statement time /3 - statement time
     
     
System: Dell Precision 690, Intel 2xQuad-Core E5345 @ 2.33GHz
        16GB RAM, Windows7 64-bit Professional.

MATLAB Version 7.6.0.324 (R2008a)

Subject: Randperm, Randi, and Shuffle

From: Derek O'Connor

Date: 15 Dec, 2010 07:20:21

Message: 16 of 19


>
> function p = GenRandPerm(n)
> p = 1:n;
> for k = n:-1:2
> r = 1 + floor(k*rand);
> 1 t = p(r); % used this statement
> 2 p(r) = p(k); % p([k r]) = p([r k]);
> 3 p(k) = t; % instead of 1,2,3
> end;
> return;
>
> n [ 1e4 1e5 1e6 1e7 ]
> R [ 37 23 24 14 ]
>
> R = 1 - statement time /3 - statement time
>
>
> System: Dell Precision 690, Intel 2xQuad-Core E5345 @ 2.33GHz
> 16GB RAM, Windows7 64-bit Professional.
> Matlab R2008a

Why does this site mangle the format of what I see in this reply window?

Matlab! you need to provide a better way to format code, tables, and simple mathematics on this site. Also, provide a preview. ''If we can send a man to the moon, surely we can ..., etc."


Derek O'Connor

Subject: Randperm, Randi, and Shuffle

From: Jan Simon

Date: 16 Dec, 2010 13:16:05

Message: 17 of 19

Dear Derek,

Let's examine where the time is lost in the swap operation:

  function p = GenRandPerm1(n)
  p = 1:n;
  for k = n:-1:2
      r = 1 + floor(k*rand);
      t = p(r);
      p(r) = p(k);
      p(k) = t;
  end

GenRandPerm1 is fast and Matlab's JIT can squeeze this loop effciently.

  function p = GenRandPerm2(n)
  p = 1:n;
  for k = n:-1:2
      r = 1 + floor(k*rand);
      p([k r]) = p([r k]);
  end

GenRandPerm2 is about 20 times slower. But the two vektors [k, r] and [r, k] needs to be constructed in each iteration. Although they are tiny, pre-allocation helps:

  function p = GenRandPerm3(n)
  v1 = zeros(1, 2);
  v2 = zeros(1, 2);
  p = 1:n;
  for k = n:-1:2
      r = 1 + floor(k*rand);
      v1(1) = k;
      v1(2) = r;
      v2(1) = r;
      v2(2) = k;
      p(v1) = p(v2);
  end

25% faster than GenRandPerm2! But still p is indexed by a vector and "p(v2)" creates a vector dynamically, which implies a memory allocation. Omit this:

  function p = GenRandPerm4(n)
  v1 = zeros(1, 2);
  v2 = zeros(1, 2);
  p = 1:n;
  for k = n:-1:2
      r = 1 + floor(k*rand);
      v1(1) = k;
      v1(2) = r;
      v2(1) = p(r);
      v2(2) = p(k);
      p(v1) = v2;
  end

Now we get it 60% faster than GenRandPerm2. But GenRandPerm1 is still 7 times faster! It seems like "p(v1)" includes a lot of work. The vector indexing is not implemented really efficiently, see also:
http://www.mathworks.com/matlabcentral/newsreader/view_thread/295653#793525

I do think, that Matlab got its name, because it is very efficient for operating on "Mat"rices, otherwise it would be calle "FOR"lab or "Sca"lab. The fact that a C-Mex can calculate the Fisher-Yates in 30% of the time of the bast GenRandPerm1 implies, that it was a very good idea to implement a C-Mex interface. It is from a scientif point of view an efficient strategy to use a language, which matchs the problem.
Of course I could formulate this posting much more eloquently in my native language, but it would be less efficient, because just a small part of the readers would understand it. (Brr, but at least it is an allegory - or is it a metapher?!).

This is the cause why I flood the FEX with C-Mex files. And I'm aware that the majority of FEX users prefer M-versions, even if they are less efficient.

Kind regards, Jan

Subject: Randperm, Randi, and Shuffle

From: Derek O'Connor

Date: 18 Dec, 2010 09:35:20

Message: 18 of 19

"Jan Simon" <matlab.THIS_YEAR@nMINUSsimon.de> wrote in message <ied3ek$qtt$1@fred.mathworks.com>...
> Dear Derek,
>
> Let's examine where the time is lost in the swap operation:
>
> function p = GenRandPerm1(n)
> p = 1:n;
> for k = n:-1:2
> r = 1 + floor(k*rand);
> t = p(r);
> p(r) = p(k);
> p(k) = t;
> end
> ---- snip ----------

Dear Jan,
Thanks for your explanation.

While I like (somewhat) Matlab's vectorization and index manipulation
facilities, I don't always understand their application. Fancy index
manipulation can be very powerful but it tends to obscure the code,
for me at least.

My old Fortran 66/77 habits still persist and I tend to write code in loop or component-form (as the numerical linear algebra people call it) first.
After I test the component form I try to vectorize, and test again.

As an example of this, here is yet another version of the Shuffle Algorithm
(can there be many more?). In this version the random number generation
is taken from the loop and vectorized:

function p = GenRandPermVecd(n)
p = 1:n;
r = 1 + floor( rand(n,1).*(1:n)' );
for k = n:-1:2
      t = p(r(k));
      p(r(k)) = p(k);
      p(k) = t;
 end

Notice that working memory has doubled from p to p+r. Testing this
for n=10^7 and n=10^8 gave times that are about 5% to 10% longer
than the component form. So this is an example where vectorization
gives the worst of both worlds.

Best wishes,

Derek

Subject: Randperm, Randi, and Shuffle

From: Jan Simon

Date: 18 Dec, 2010 21:44:05

Message: 19 of 19

Dear Derek,

> function p = GenRandPermVecd(n)
> p = 1:n;
> r = 1 + floor( rand(n,1).*(1:n)' );

You got me. I tried exactly the same and got the same slowdown of 10%.
A drawback is the memory consumption:
  1. p = 1:n; % n doubles
  2. a = rand(n,1) % n doubles
  3. b = 1:n % n doubles (should reuse p here)
  4. c = (a * b) % n doubles
  5. d = floor(c) % n doubles
  6. r = 1 + d % n doubles

Sme of the temporarily used memory will be used twice, but without the JIT acceleration this will so much time, that the actually faster vectorized call of RAND does not accelerate the complete program.

I'd be interested in the timing of a corresponding Fortran77 Mex function. I've been very impressed each time I see the speed of naively implemented matrix matrix multiplications in Fortran compared to damn tricky unrolled loop pointer based C versions.

Kind regards, Jan

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