perlmeditation
tlm
I recently came across the following code for generating a "random permutation"<sup><small>1</small></sup> of an array's elements.
<code>
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;
}
</code>
It looks reasonable enough, if perhaps a bit zealous (we know that only <em>N</em> - 1 swaps are required, since fixing a permutation of <em>N</em> elements is a problem with <em>N</em> - 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.</p>
<p>The easiest way to see this is to realize that what this algorithm <em>does</em> sample uniformly is a space of possible <em>execution paths</em>, as determined by the values returned by the expression <code>int rand $n</code> at the various iteration steps. The size of <em>this</em> sample space is <em>N<sup><small>N</small></sup></em>. To each element of this space corresponds a permutation, but the size of the space of all possible permutations is <em>N</em>! , which is not only smaller than <em>N<sup><small>N</small></sup></em> for any <em>N</em> > 1, but more importantly, it is not a <em>divisor</em> of <em>N<sup><small>N</small></sup></em>, which means that it is <em>impossible</em> 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 <em>N</em> = 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.</p>
<p>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.</p>
<p>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.</p>
<readmore>
<p>To beat a dead horse, and to indulge in a free-association to a reference I made in another [id://441020|post], here's a "forking paths" strategy to pin down the algorithm's actual sampling. We modify <code>naive_shuffle</code> (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 <code>int rand $n</code>.</p>
<p>The first modification is to "recursify" the loop:
<code>
# 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;
}
</code>
Here we rely on recursive calls and an additional "state argument" to handle the iteration of $i from 0 to $n-1<sup><small>2</small></sup>.</p>
<p>The next step is more delicate. It involves, essentially, replacing the expression <code>int rand $n</code> with all the possible execution paths implicit in this expression. Here's one way to do it:
<code>
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;
}
}
</code>
A call to
<code>
naive_shuffle_traverse(-1, undef, 1..3);
</code>
causes the program to traverse all possible execution paths of the original <code>naive_shuffle( 1..3 )</code>. The business with the -1 value for <code>$i</code> is because now the first call to the function is to set things up for <code>$j</code>.</p>
<p>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 <code>$tally</code>:
<code>
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;
}
}
</code>
Now, executing
<code>
my %tally = ();
naive_shuffle_count( -1, undef, \%tally, 1..3 );
</code>
results in the following contents for <code>%tally</code>:
<code>
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
</code>
showing that with this algorithm some permutations are "more equal than others."</p>
</readmore>
<p>At any rate, the correct algorithm, if a monk may allow himself to be categorical about it<sup><small>3</small></sup>, is this one:
<code>
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;
}
</code>
Note that the size of the space of possible execution paths for this algorithm is exactly <em>N</em>!.</p>
<p><small>the lowliest monk</small></p>
<hr>
<small>
<p><small><sup>1</sup></small>To put it more pedantic terms, this function is supposed to sample <em>uniformly at random</em> the space of all the permutations of an array's elements.</p>
<p><small><sup>2</sup></small>As an aside, though this modification involves making the sub recursive, <em>the algorithm is still iterative</em>. This is because the subroutine is <em>tail-recursive</em>, which means that its space requirements, in principle at least, do not grow with the size of the input. For more on this see [http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-11.html#%_sec_1.2.1|here], especially [http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-11.html#%_idx_682|here].</p>
<p><small><sup>3</sup></small>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 <code>rand_perm</code> 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.</p>
</small>
<hr>
<p><b>Update:</b> <em>I found a [http://c2.com/cgi/wiki?LinearShuffleSummary|page] that gives a nice discussion of the algorithm used by <code>naive_shuffle</code>, where it is described as "short, elegant, clear, and wrong". Also, much more material [http://c2.com/cgi/wiki?LinearShuffle|here], although I object to this page's description of "Knuth's algorithm" (the one used by <code>random_perm</code>) as being "harder to code". This page does point out something that I failed to mention, namely that <code>naive_shuffle</code> 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.</em></p>