What, if anything, is wrong with this shuffling algorithm and how can I know?

Question

Just as background, I'm aware of the Fisher-Yates perfect shuffle. It is a great shuffle with its O(n) complexity and its guaranteed uniformity and I'd be a fool not to use it ... in an environment that permits in-place updates of arrays (so in most, if not all, imperative programming environments).

Sadly the functional programming world doesn't give you access to mutable state.

Because of Fisher-Yates, however, there's not a lot of literature I can find on how to design a shuffling algorithm. The few places that address it at all do so briefly before saying, in effect, "so here's Fisher-Yates which is all the shuffling you need to know". I had to, in the end, come up with my own solution.

The solution I came up with works like this to shuffle any list of data:

• If the list is empty, return the empty set.
• If the list has a single item, return that single item.
• If the list is non-empty, partition the list with a random number generator and apply the algorithm recursively to each partition, assembling the results.

In Erlang code it looks something like this:

shuffle([])  -> [];
shuffle([L]) -> [L];
shuffle(L)   ->
  {Left, Right} = lists:partition(fun(_) -> 
                                    random:uniform() < 0.5 
                                  end, L),
  shuffle(Left) ++ shuffle(Right).

(If this looks like a deranged quick sort to you, well, that's what it is, basically.)

So here's my problem: the same situation that makes finding shuffling algorithms that aren't Fisher-Yates difficult makes finding tools to analyse a shuffling algorithm equally difficult. There's lots of literature I can find on analysing PRNGs for uniformity, periodicity, etc. but not a lot of information out there on how to analyse a shuffle. (Indeed some of the information I found on analysing shuffles was just plain wrong -- easily deceived through simple techniques.)

So my question is this: how do I analyse my shuffling algorithm (assuming that the random:uniform() call up there is up to the task of generating apropriate random numbers with good characteristics)? What mathematical tools are there at my disposal to judge whether or not, say, 100,000 runs of the shuffler over a list of integers ranging 1..100 has given me plausibly good shuffling results? I've done a few tests of my own (comparing increments to decrements in the shuffles, for example), but I'd like to know a few more.

And if there's any insight into that shuffle algorithm itself that would be appreciated too.

Answer 1

The answers to this question might help:

http://stackoverflow.com/questions/1685339/verify-knuth-shuffle-algorithm-is-as-unbiased-as-possible

It'd also be worth taking a look at Knuth's analysis of Fisher-Yates (see the wikipedia article you linked for a citation).

Answer 2

Your algorithm is a sort-based shuffle, as discussed in the Wikipedia article.

Generally speaking, the computational complexity of sort-based shuffles is the same as the underlying sort algorithm (e.g. O(n log n) average, O(n²) worst case for a quicksort-based shuffle), and while the distribution is not perfectly uniform, it should approach uniform close enough for most practical purposes.

Oleg Kiselyov provides the following article / discussion:

• Provably perfect random shuffling and its pure functional implementations

which covers the limitations of sort-based shuffles in more detail, and also offers two adaptations of the Fischer–Yates strategy: a naive O(n²) one, and a binary-tree-based O(n log n) one.

Sadly the functional programming world doesn't give you access to mutable state.

This is a complete myth: functional programming avoids side effects, but supports mutable state perfectly well, without invoking side effects.

In this particular case, you can use Haskell's mutable arrays to implement the Fischer–Yates algorithm as described in this tutorial:

• Haskell Shuffling (Brett Hall)

Addendum

The specific foundation of your shuffle sort is actually an infinite-key radix sort: as gasche points out, each partition corresponds to a digit grouping.

The main disadvantage of this is the same as any other infinite-key sorting shuffle: there is no termination guarantee. Although the likelihood of termination increases as the comparison proceeds, there is never an upper bound: the worst-case complexity is O(∞).

Answer 3

General remark

My personal approach about correctness of probability-using algorithms : if you know how to prove it's correct, then it's probably correct; if you don't, it's certainly wrong.

Said differently, it's generally hopeless to try to analyse every algorithm you could come up with : you have to keep looking for an algorithm until you find one that you can prove correct.

Analysing a random algorithm by computing the distribution

I know of one way to "automatically" analyse a shuffle (or more generally a random-using algorithm) that is stronger than the simple "throw lots of tests and check for uniformity". You can mechanically compute the distribution associated to each input of your algorithm.

The general idea is that a random-using algorithm explores a part of a world of possibilities. Each time your algorithm asks for an element randomly choosen in a set ({true, false} when flipping a coin), there are two possible outcomes for your algorithm, and one of them is chosen. You can change your algorithm so that, instead of returning one of the possible outcomes, it explores all solutions in parallel and return all possible outcomes with the associated distributions.

In general, that would require to rewrite your algorithm in depth. If your language supports delimited continuations, you don't have to, you can implementation that "exploration of all possible outcomes" inside the function asking for a random element (the idea is that the random generator, instead of returning a result, capture the continuation associated to your program and run it with all different results). For an example of this approach, see oleg's HANSEI

An intermediary, and probably less arcane, solution is to represent this "world of possible outcomes" as a monad, and use a language such as Haskell with facilities for monadic programming.
Here is an example implementation of a variant¹ of your algorithm, in Haskell, using the probability monad of the probability package :

import Numeric.Probability.Distribution

shuffleM :: (Num prob, Fractional prob) => [a] -> T prob [a]
shuffleM [] = return []
shuffleM [x] = return [x]
shuffleM (pivot:li) = do
        (left, right) <- partition li
        sleft <- shuffleM left
        sright <- shuffleM right
        return (sleft ++ [pivot] ++ sright)
  where partition [] = return ([], [])
        partition (x:xs) = do
                  (left, right) <- partition xs
                  uniform [(x:left, right), (left, x:right)]

You can run it for a given input, and get the output distribution :

*Main> shuffleM [1,2]
fromFreqs [([1,2],0.5),([2,1],0.5)]
*Main> shuffleM [1,2,3]
fromFreqs
  [([2,1,3],0.25),([3,1,2],0.25),([1,2,3],0.125),
   ([1,3,2],0.125),([2,3,1],0.125),([3,2,1],0.125)]

You can see that this algorithm is uniform with inputs of size 2, but non-uniform on inputs of size 3.
The difference with the test-based approach is that we can gain absolute certainty in a finite number of steps : it can be quite big, as it amounts to an exhaustive exploration of the world of possibles (but generally smaller than 2^N, as their are factorisations of similar outcomes), but if it returns a non-uniform distribution we know for sure that the algorithm is wrong. Of course, if it returns an uniform distribution for [1..N] and 1 <= N <= 100, you only know that your algorithm is uniform upto lists of size 100, it may still be wrong.

¹: this algorithm is a variant of your Erlang's implementation, because of the specific pivot handling. If I use no pivot, like in your case, the input size doesn't decrease at each step anymore : the algorithm also consider the case were all inputs are in the left list (or right list), and get lost in an infinite loop. This is a weakness of the probability monad implementation (if an algorithm has a probability 0 of non-termination, the distribution computation may still diverge), that I don't know yet how to fix.

Sort-based shuffles

Here is a simple algorithm that I feel confident I could prove correct :

1. pick a random number (a weight) for each element in your collection
2. if the random weights are not all distinct, restart step 1.
3. "sort" your collection, using the comparison on the random weights associated to each element

You can forget step 2. if you know the probability of "conflict" (two random numbers picked are equal) is low enough, but unless it's exactly 0 your shuffle isn't uniform anymore.

If you pick your weights in [1..N] where N is the length of your collection, you'll have lots of conflicts (Birthday problem). If you pick your weight in [min_int..max_int] or as a random float, the probability of conflict is low in practice (but still subject to the birthday problem, see Piet Delport's comment).

If you pick real numbers (as random infinite lazy lists of booleans) as weights, it is 0 in theory, and you don't need to check for distinctness.

Here is a shuffle implementation in OCaml, using lazy real numbers :

type 'a stream = Cons of 'a * 'a stream lazy_t

let rec real_number () =
  Cons (Random.bool (), lazy (real_number ()))

let rec compare_real a b = match a, b with
| Cons (true, _), Cons (false, _) -> 1
| Cons (false, _), Cons (true, _) -> -1
| Cons (_, lazy a'), Cons (_, lazy b') ->
    compare_real a' b'

let shuffle list =
  List.map snd
    (List.sort (fun (ra, _) (rb, _) -> compare_real ra rb)
       (List.map (fun x -> real_number (), x) list))

There are other approaches to "pure shuffling". A nice one is apfelmus's mergesort-based solution.

Algorithmic considerations : the complexity of the previous algorithm depends on the probability that all picked weights are distinct. If you pick them in [min_int..max_int], you have a 2^{-N} probability of conflict, which means an average number of weight choices of 1, and a O(n log n) algorithm (the only complexity comes from the sorting), assuming picking a random number is O(1).

If you pick random boolean streams, you never have to restart picking, but the complexity is then related to "how many elements of the streams are evaluated on average". I conjecture it is O(log n) in average (hence still a O(n log n) in total), but have no proof.

... and I think your algorithm works

After more reflexion, I think (like douplep), that your implementation is correct. Here is an informal explanation.

Each element in your list is tested by several random:uniform() < 0.5 tests. To an element, you can associate the list of outcomes of those tests, as a list of booleans or {0, 1}. At the beginning of the algorithm, you don't know the list associated to any of those number. After the first partition call, you know the first element of each list, etc. When your algorithm returns, the list of tests are completely known and the elements are sorted according to those lists (sorted in lexicographic order, or considered as binary representations of real numbers).

In substance, your algorithm is equivalent to my "sort over real numbers, as lazy lists of booleans". The action of partitioning the list, reminiscent of quicksort's partition over a pivot element, is actually a way of separating, for a given position in the binary decomposition of those numbers, the elements with valuation 0 from the elements with valuation 1.

The sort is uniform because the real number picked are all different. Indeed, two elements with real numbers equal upto the n-th decimal are on the same side of a partition occuring during a recursive shuffle call of depth n. The algorithm only terminates when all the list resulting from partitions are empty or singletons : all elements have been separated by at least one test, and therefore have one distinct binary decimal.

Probabilistic termination

A subtle point about your algorith (or my equivalent sort-based method) is that the termination condition is probabilistic. Fisher-Yates always terminate after a known number of steps (the number of elements in the array). With your algorithm, the termination depends on the output of the random number generator.

There are possible outputs that would make your algorithm diverge, not terminate. For example, if the random number generate always output 0, each partition call will return the input list unchanged, on which you recursively call the shuffle : you will loop indefinitely.

However, this is not an issue if you're confident that your random number generator is fair : it does not cheat and always return independent uniformly distributed results. In that case, the probability that the test random:uniform() < 0.5 always returns true (or false) is exactly 0 :

• the probability that the first N calls return true is 2^{-N}
• the probability that all calls return true is the probability of the infinite intersection, for all N, of the event that the first N calls return 0; it is the infimum limit¹ of the 2^{-N}, which is 0

¹: for the mathematical details, see http://en.wikipedia.org/wiki/Measure_(mathematics)#Measures_of_infinite_intersections_of_measurable_sets

More generally, the algorithm does not terminate if and only if some of the elements get associated to the same boolean stream. This means that at least two elements have the same boolean stream. But the probability that two random boolean streams are equal is again 0 : the probability that the digits at position K are equal is 1/2, so the probability that the N first digits are equal is 2^{-N}, and the same analysis applies.

Therefore, you know that your algorithm terminates with probability 1. This is a slightly weaker guarantee that the Fisher-Yates algorith, which always terminate. In particular, you're vulnerable to an attack of an evil adversary that would control your random number generator.

With more probability theory, you could also compute the distribution of running times of your algorithm for a given input length. This is beyond my technical abilities, but I assume it's good : I suppose that you only need to look at O(log N) first digits on average to check that all N lazy streams are different, and that the probability of much higher running times decrease exponentially.