A Proof of Correctness of Durstenfeld's Random Permutation Algorithm

Question: Does anyone have a precise mathematical proof of Durstenfeld's Algorithm?

The first $O(n)$ shuffle or random permutation generator was published by Richard Durstenfeld in 1964. This algorithm came to the notice of programmers because it was included in Knuth 1969, page 125, 2, as Algorithm P. A succinct statement of the algorithm is: $$\text{for }\ k \leftarrow n,n-1,\ldots,2\ \text{ do } \begin{cases} \text{Choose at random a number}\ \ r \in [1,k] \\ \text{Interchange }\pi[r] \text{ and }\pi[k],\ \end{cases}$$ where $\pi[1 \ldots n]$ is the array to be shuffled or randomly permuted.

This is a truly elegant algorithm, computationally and mathematically. Computationally it handles set operations by dividing the array $\pi$ into an unshuffled set $\pi[1,\ldots, k]$, and a shuffled set $\pi[k+1,\ldots, n]$. At each iteration $k$, it chooses a random element $\pi[r]$ from the unshuffled set $\pi[1,\ldots, k]$, and interchanges it with $\pi[k]$. This moves $\pi[r]$ into the new shuffled set $\pi[k,\ldots, n]$, where it remains, unmoved, for all subsequent iterations. This is done in $O(1)$ time and space. Thus, the time and space complexity of Durstenfeld's algorithm is $O(n)$, which is optimal.

Mathematically, the algorithm is elegant because it implicitly uses the well-known lemma that every permutation $\pi$ is a product of transpositions. At each iteration $k$, the algorithm chooses a random number $r_k \in [1,\ldots,k]$ and performs a transposition $(\pi[k], \pi[r_k])$. Thus the random permutation produced can be written as

$$\pi = (\pi[n],\pi[r_n])(\pi[n-1],\pi[r_{n-1}]) \cdots (\pi[k],\pi[r_k])\cdots (\pi[2],\pi[r_2]),$$

where $(\pi[k],\pi[r_k])$ is the transposition or interchange of $\pi[k]$ and $\pi[r_k]$.

Matlab Implementation of Durstenfeld's Algorithm

function p = GRPdurG(p)
% -------------------------------------------------------------
% Randomly permute the elements of p(1:n) using Durstenfeld's
% Random Permutation Algorithm, CACM, Vol 7, No. 7, 1964.
% See Knuth, Section 3.4.2, TAOCP, Vol 2, 3rd Ed.
% Derek O'Connor, 8 Dec 2010.  derekroconnor@eircom.net
%
% USE: n=10^7; p = 1:n; p = GRPdur(p);
% -------------------------------------------------------------
n = length(p);
for k = n:-1:2
r = 1+floor(rand*k);    % random integer between 1 and k
t    = p(k);
p(k) = p(r);            % Swap(p(r),p(k)).
p(r) = t;
end
%  GRPdurG


A Combinatorial Conundrum

A random permutation of the integers $\{1,2, \ldots,10^7\}$ can be generated in a few seconds using the Durstenfeld Shuffle algorithm on a standard PC. The combinatorial space $\mathcal{C}$ in this case is set of all permutations of size $10^7$. The size of this sample space is \begin{equation*} |\mathcal{C}| = n! = (10^7)! \approx 10^{65,657,059} \end{equation*} This is a gigantic number and very far beyond limits of any known random number generator. That is, $|\mathcal{C}| \ggg |\mathcal{G}|$, the state space of the random number generator. For example, the period (size of state space) of the well-known Mersenne Twister RNG (MT) is about $10^{6000}$. Hence random permutation generators using MT must miss nearly all permutations in this space.

None-the-less, the Matlab function above produces random permutations of length $10^7$ in the correct statistical proportions, e.g, $1/e$ are derangements, $1/n$ are cycles, etc., yet it must miss most of the possible permutations. Hence the conundrum. The most important question to me is:

Do random permutation generators sample uniformly from the space of permutations?

assuming they use the best random number generator available. I suspect that this is a difficult mathematical problem.

See this Matlab discussion: What Does RANDPERM Miss?

Historical Note

I believe that Durstenfeld's Algorithm was incorrectly attributed to Fisher and Yates by Knuth, in (2) below. See here O'Connor Scribd Note on Fisher-Yates. The original Durstenfeld Algorithm is here: CACM Algorithm

References

1. Durstenfeld, R. 1964: ACM Algorithm 235: Random Permutation, Communications of the ACM, Vol 7, No 7.
2. Knuth, D. 1969, 1998: Seminumerical Algorithms 1st & 3rd Eds. The Art of Computer Programming Series, Volume 2.
-
What's the question? Also, if you are really convinced you know better than Knuth, who is known for his meticulous historical investigations, then you can write him about it and get some (moral or financial) credit. Fisher and Yates wrote before the era of electronic computing, so usual analysis of execution time does not really apply to their case. –  Marc van Leeuwen Dec 17 '11 at 17:32
@Marc van Leeuwen. Have you read my scribd note above? If you have, what part of it do you find wrong? –  Derek O'Connor Dec 17 '11 at 18:13
I think your attack against Fisher-Yates is a strawman argument, in that you provide a particular computer implementation that is not implied by the original description, and then show that this procedure is not efficient. If you implement the shuffle by taking a deck of cards and repeatedly removing a card from the remaining deck uniformly at random, then you have $O(n)$ time. Transposing to an electronic computer requires some adaptation, which is what Durstenfeld did, but I'm sure Knuth considered it fundamentally the same idea. It's not like bubble versus quicksort. –  Marc van Leeuwen Dec 17 '11 at 18:54
Frank Yates, in Statistical Tables, 6th Edition, Oliver \& Boyd, London, 1963, Example 12, says that the new method by C.R Rao is better, i.e., faster, than the original F&Y method. Yet the Rao method is not remotely similar to Durstenfeld. Your "electronic era" comment maybe a straw man, but certainly it is a red herring: Durstenfeld's Algorithm can be carried out easily with paper and pencil. –  Derek O'Connor Dec 17 '11 at 19:26
The attribution to Fisher and Yates is for what they wrote in 1938, not for what Yates wrote in 1963. One can attribute a method to a fundamental idea even if the person having had that idea did not give the method, or for historical reasons would even be incapable of understanding the method. We talk of Euclidean division of polynomials even though Euclid did not know polynomials. And even for integers, our Euclidean algorithm computes remainders, whereas Euclid did not do division, just repeated subtraction; this is fundamentally less efficient. Should we attribute these differently? –  Marc van Leeuwen Dec 18 '11 at 15:49

I think the Wikipedia article contains the proof, but in any case it is easy. The proof of the correctness of the Durstenfeld algorithm is the same as for the original procedure by Fisher and Yates, by induction on the number $n$ of elements to permute (the cases $n=0$ is trivial). One starts with the items in any order; since they are going to be randomly permuted, this order is of no importance. To get a random permutation of $n$ distinct items, the first item must be chosen uniformly at random. Once this element is chosen, the remaining elements, which we can arrange in any order (and the only difference between Fisher-Yates and Durstenfeld is which order we start out with for the remainder) must be permuted uniformly at random, which is done by applying the algorithm recursively to the list of those remaining items.
Alternatively, one can show (also by induction), that calling the random number generator $n$ times to generate an element of the Cartesian product $[n]\times[n-1]\times\cdots\times[1]$, there is a bijection between the sequences produced by the generator and the resulting permutations of the sequence, so supposing the Cartesian product is sampled uniformly at random, so is the permutation produced.
The difficulty is not the probabilistic argument, but the requirements imposed on the random generator. When producing a random permutation of a large number of items (say $n=10000$), it is very likely that the Cartesian product to be sampled is simply larger than the state space of the random generator, which means that a given generator cannot possibly generate all $n!$ permutations, let alone generate them uniformly. This has nothing to do with the algorithm used for producing the permutation. Whether this is a problem depends on what the permutation is going to be used for; if it is just to assure some mixing in a statistical procedure I think it doesn't matter if there are certain permutations that will never show up, but if it is to produce a shuffle for an on-line game where large amounts of money are involved, it would be embarrassing if a player could prove that the system used cannot possibly produce his winning combination.