comment on

I recently came across the following code for generating a "random permutation"¹ of an array's elements.

sub naive_shuffle {
  my $n = my @p = @_;
  for my $i ( 0..$n-1 ) {
    my $j = int rand $n;
    @p[ $i, $j ] = @p[ $j, $i ]; # swap
  }
  return @p;
}
[download]

It looks reasonable enough, if perhaps a bit zealous (we know that only N - 1 swaps are required, since fixing a permutation of N elements is a problem with N - 1 degrees of freedom; see "correct algorithm" below). But the real problem with this subroutine is that it does not sample the space of permutations uniformly.

The easiest way to see this is to realize that what this algorithm does sample uniformly is a space of possible execution paths, as determined by the values returned by the expression int rand $n at the various iteration steps. The size of this sample space is N^N. To each element of this space corresponds a permutation, but the size of the space of all possible permutations is N! , which is not only smaller than N^N for any N > 1, but more importantly, it is not a divisor of N^N, which means that it is impossible for the algorithm to give equal weight to all the permutations. Some permutations will be generated more often than the others. For example, in the case of N = 3, there are 27 execution paths; of the 6 possible permutations of 3 elements, 3 of them are generated by 5 execution paths each, and the other 3 by 4 paths each. Hence the former are oversampled relative to the latter.

If this argument were not convincing enough, or if we are trying to assess a sampling algorithm for which a rigorous argument as to its correctness is not readily apparent, we can just run the algorithm many times, keeping a tally of the permutations generated, to get an idea of how even-handed the sampling was. The only problem with this approach is that the sampling noise puts a limit on the size of the deviations from uniformity that can be detected, and the cost of reducing this noise grows quadratically with the number of sampling runs. Be that as it may, this is the first thing to try, because it's trivial to implement.

If the lazy approach is not acceptable for whatever reason, an alternative is to modifify the algorithm to systematically follow all the possible execution paths, and keep a count of the permutations obtained in the process of traversing all paths.

To beat a dead horse, and to indulge in a free-association to a reference I made in another post, here's a "forking paths" strategy to pin down the algorithm's actual sampling. We modify naive_shuffle (taking great care not to alter the sampling) so that instead of generating one permutation, it generates all the possible permutations implied by the multiple evaluations of int rand $n.

The first modification is to "recursify" the loop:

# untested
sub naive_shuffle_rec {
  my $i = shift;
  my $n = my @p = @_;
  my $j = int rand $n;
  @p[ $i, $j ] = @p[ $j, $i ]; # swap
  return ++$i < $n ? naive_shuffle_rec( $i, @p ) : @p;
}
[download]

Here we rely on recursive calls and an additional "state argument" to handle the iteration of $i from 0 to $n-1².

The next step is more delicate. It involves, essentially, replacing the expression int rand $n with all the possible execution paths implicit in this expression. Here's one way to do it:

sub naive_shuffle_traverse {
  my $i = shift;
  my $j = shift;
  my $n = my @p = @_;
  @p[ $i, $j ] = @p[ $j, $i ] if defined $j; # swap
  if ( ++$i < $n ) {
    for my $j (0..$n-1) {
      naive_shuffle_traverse( $i, $j, $tally, @p );
    }
  }
  else {
    return;
  }
}
[download]

A call to

naive_shuffle_traverse(-1, undef, 1..3);
[download]

causes the program to traverse all possible execution paths of the original naive_shuffle( 1..3 ). The business with the -1 value for $i is because now the first call to the function is to set things up for $j.

The only remaining modification is an easy one: we need some way to have the routine keep track of the permutations visited so far. For this we use an additional state variable $tally:

sub naive_shuffle_count {
  my $i = shift;
  my $j = shift;
  my $tally = shift;
  my $n = my @p = @_;
  @p[ $i, $j ] = @p[ $j, $i ] if $j > -1;
  if ( ++$i < $n ) {
    for my $j (0..$n-1) {
      naive_shuffle_count( $i, $j, $tally, @p );
    }
  }
  else {
    $tally->{ "@p" }++;
    return;
  }
}
[download]

Now, executing

my %tally = ();
naive_shuffle_count( -1, undef, \%tally, 1..3 );
[download]

results in the following contents for %tally:

  DB<2> x \%tally
0  HASH(0x82d7f94)
   '1 2 3' => 4
   '1 3 2' => 5
   '2 1 3' => 5
   '2 3 1' => 5
   '3 1 2' => 4
   '3 2 1' => 4
[download]

showing that with this algorithm some permutations are "more equal than others."

At any rate, the correct algorithm, if a monk may allow himself to be categorical about it³, is this one:

sub random_perm {
  my $n = my @p = @_;
  for ( my $i = $#p; $i > 0; --$i ) {
    my $j = int rand( $i + 1 );
    @p[ $i, $j ] = @p[ $j, $i ];
  }
  return @p;
}
[download]

Note that the size of the space of possible execution paths for this algorithm is exactly N!.

the lowliest monk

¹To put it more pedantic terms, this function is supposed to sample uniformly at random the space of all the permutations of an array's elements.

²As an aside, though this modification involves making the sub recursive, the algorithm is still iterative. This is because the subroutine is tail-recursive, which means that its space requirements, in principle at least, do not grow with the size of the input. For more on this see here, especially here.

³One of the many enduring coinages of the mathematician Pál Erdös was the phrase "straight from The Book", which was his way of pronouncing a mathematical proof as perfect, and which made playful reference to a book where god keeps all mathematical secrets. Well, if this book exists, it surely has a chapter or two for algorithms, and if so, the algorithm used by rand_perm above, which I learned from Knuth, surely is in there. BTW, the subject of candidates for the chapter on algorithm in The Book would be a fine one for another meditation.

Update: I found a page that gives a nice discussion of the algorithm used by naive_shuffle, where it is described as "short, elegant, clear, and wrong". Also, much more material here, although I object to this page's description of "Knuth's algorithm" (the one used by random_perm) as being "harder to code". This page does point out something that I failed to mention, namely that naive_shuffle gives correct results for the case of N = 2, although it works twice harder to achieve it than needed, like a man who, facing a choice of two things, flips a coin twice, disregards the first flip, and takes the outcome of the second.

In reply to A bad shuffle by tlm

Are you posting in the right place? Check out Where do I post X? to know for sure.
Posts may use any of the Perl Monks Approved HTML tags. Currently these include the following:
<code> <a> <b> <big> <blockquote> <br /> <dd> <dl> <dt> <em> <font> <h1> <h2> <h3> <h4> <h5> <h6> <hr /> <i> <li> <nbsp> <ol> <p> <small> <strike> <strong> <sub> <sup> <table> <td> <th> <tr> <tt> <u> <ul>
Snippets of code should be wrapped in <code> tags not <pre> tags. In fact, <pre> tags should generally be avoided. If they must be used, extreme care should be taken to ensure that their contents do not have long lines (<70 chars), in order to prevent horizontal scrolling (and possible janitor intervention).
Want more info? How to link or How to display code and escape characters are good places to start.


Think about Loose Coupling
	PerlMonks