# What is wrong with this “naive” shuffling algorithm?

This is a follow-up to a Stackoverflow question about shuffling an array randomly.

There are established algorithms (such as the Knuth-Fisher-Yates Shuffle) that one should use to shuffle an array, rather than relying on "naive" ad-hoc implementations.

I am now interested in proving (or disproving) that my naive algorithm is broken (as in: does not generate all possible permutations with equal probability).

Here is the algorithm:

Loop a couple of times (length of array should do), and in every iteration, get two random array indexes and swap the two elements there.

Obviously, this needs more random numbers than KFY (twice as much), but aside from that does it work properly? And what would be the appropriate number of iterations (is "length of array" enough)?

• I just can't understand why people think this swapping one is 'simpler' or 'more naive' than FY... When I was solving this problem for a first time I have just implemented FY (not knowing it has even a name), just because it seemed the simplest way to do it for me. – user88 Sep 26 '10 at 9:56
• @mbq: personally, I find them equally easy, although I agree that FY seems more "natural" to me. – nico Sep 26 '10 at 10:29
• When I researched shuffling algorithms after writing my own (a practice I have since abandoned), I was all "holy crap, it's been done, and it has a name!!" – J. M. is not a statistician Sep 26 '10 at 13:55
• – DMGregory Sep 11 '18 at 2:01

It is broken, although if you perform enough shuffles it can be an excellent approximation (as the previous answers have indicated).

Just to get a handle on what's going on, consider how often your algorithm will generate shuffles of a $k$ element array in which the first element is fixed, $k \ge 2$. When permutations are generated with equal probability, this should happen $1/k$ of the time. Let $p_n$ be the relative frequency of this occurrence after $n$ shuffles with your algorithm. Let's be generous, too, and suppose you are actually selecting distinct pairs of indexes uniformly at random for your shuffles, so that each pair is selected with probability $1/{k \choose 2}$ = $2/\left( k (k-1) \right)$. (This means there are no "trivial" shuffles wasted. On the other hand, it totally breaks your algorithm for a two-element array, because you alternate between fixing the two elements and swapping them, so if you stop after a predetermined number of steps, there is no randomness to the outcome whatsoever!)

This frequency satisfies a simple recurrence, because the first element is found in its original place after $n+1$ shuffles in two disjoint ways. One is that it was fixed after $n$ shuffles and the next shuffle does not move the first element. The other is that it was moved after $n$ shuffles but the $n+1^{st}$ shuffle moves it back. The chance of not moving the first element equals ${k-1 \choose 2}/{k \choose 2}$ = $(k-2)/k$, whereas the chance of moving the first element back equals $1/{k \choose 2}$ = $2/\left( k (k-1) \right)$. Whence:

$$p_0 =1$$ because the first element starts out in its rightful place;

$$p_{n+1} = \frac{k-2}{k} p_n + \frac{2}{k(k-1)} \left( 1 - p_n \right).$$

The solution is

$$p_n = 1/k + \left( \frac{k-3}{k-1} \right) ^n \frac{k-1}{k}.$$

Subtracting $1/k$, we see that the frequency is wrong by $\left( \frac{k-3}{k-1} \right) ^n \frac{k-1}{k}$. For large $k$ and $n$, a good approximation is $\frac{k-1}{k} \exp(-\frac{2n}{k-1})$. This shows that the error in this particular frequency will decrease exponentially with the number of swaps relative to the size of the array ($n/k$), indicating it will be difficult to detect with large arrays if you have made a relatively large number of swaps--but the error is always there.

It is difficult to provide a comprehensive analysis of the errors in all frequencies. It's likely they will behave like this one, though, which shows that at a minimum you would need $n$ (the number of swaps) to be large enough to make the error acceptably small. An approximate solution is

$$n \gt \frac{1}{2} \left(1 - (k-1) \log(\epsilon) \right)$$

where $\epsilon$ should be very small compared to $1/k$. This implies $n$ should be several times $k$ for even crude approximations (i.e., where $\epsilon$ is on the order of $0.01$ times $1/k$ or so.)

All this begs the question: why would you choose to use an algorithm that is not quite (but only approximately) correct, employs exactly the same techniques as another algorithm that is provably correct, and yet which requires more computation?

# Edit

Thilo's comment is apt (and I was hoping nobody would point this out, so I could be spared this extra work!). Let me explain the logic.

• If you make sure to generate actual swaps each time, you're utterly screwed. The problem I pointed out for the case $k=2$ extends to all arrays. Only half of all the possible permutations can be obtained by applying an even number of swaps; the other half is obtained by applying an odd number of swaps. Thus, in this situation, you can never generate anywhere near a uniform distribution of permutations (but there are so many possible ones that a simulation study for any sizable $k$ will be unable to detect the problem). That's really bad.

• Therefore it is wise to generate swaps at random by generating the two positions independently at random. This means there is a $1/k$ chance each time of swapping an element with itself; that is, of doing nothing. This process effectively slows down the algorithm a little bit: after $n$ steps, we expect only about $\frac{k-1}{k} N \lt N$ true swaps to have occurred.

• Notice that the size of the error decreases monotonically with the number of distinct swaps. Therefore, conducting fewer swaps on average also increases the error, on average. But this is a price you should be willing to pay in order to overcome the problem described in the first bullet. Consequently, my error estimate is conservatively low, approximately by a factor of $(k-1)/k$.

I also wanted to point out an interesting apparent exception: a close look at the error formula suggests that there is no error in the case $k=3$. This is not a mistake: it is correct. However, here I have examined only one statistic related to the uniform distribution of permutations. The fact that the algorithm can reproduce this one statistic when $k=3$ (namely, getting the right frequency of permutations that fix any given position) does not guarantee the permutations have indeed been distributed uniformly. Indeed, after $2n$ actual swaps, the only possible permutations that can be generated are $(123)$, $(321)$, and the identity. Only the latter fixes any given position, so indeed exactly one-third of the permutations fix a position. But half the permutations are missing! In the other case, after $2n+1$ actual swaps, the only possible permutations are $(12)$, $(23)$, and $(13)$. Again, exactly one of these will fix any given position, so again we obtain the correct frequency of permutations fixing that position, but again we obtain only half of the possible permutations.

This little example helps reveal the main strands of the argument: by being "generous" we conservatively underestimate the error rate for one particular statistic. Because that error rate is nonzero for all $k \ge 4$, we see that the algorithm is broken. Furthermore, by analyzing the decay in the error rate for this statistic we establish a lower bound on the number of iterations of the algorithm needed to have any hope at all of approximating a uniform distribution of permutations.

• "Let's be generous, too, and suppose you are actually selecting distinct pairs of indexes uniformly at random for your shuffles". I don't understand why that assumption can be made, and how it is generous. It does seem to discard possible permutations, resulting in an even less random distribution. – Thilo Sep 27 '10 at 1:14
• @Thilo: Thank you. Your comment deserves an extended answer, so I placed it in the response itself. Let me point out here that being "generous" does not actually discard any permutations: it just eliminates steps in the algorithm that otherwise would do nothing. – whuber Sep 27 '10 at 15:23
• This problem can be analyzed fully as a Markov chain on the Cayley graph of the permutation group. Numerical calculations for k = 1 through 7 (a 5040 by 5040 matrix!) confirm that the largest eigenvalues in size (after 1 and -1) are exactly $(k-3)/(k-1) = 1 - 2/(k-1)$. This implies that once you have coped with the problem of alternating the sign of the permutation (corresponding to the eigenvalue of -1), the errors in all probabilities decay at the rate $(1 - 2/(k-1))^n$ or faster. I suspect this continues to hold for all larger $k$. – whuber Sep 27 '10 at 17:48
• You can do much better than $5040 \times 5040$ since the probabilities are invariant on conjugacy classes, and there are only $15$ partitions of $7$ so you can analyze a $15 \times 15$ matrix instead. – Douglas Zare Aug 18 '12 at 20:54

I think your simple algorithm will shuffle the cards correctly as the number shuffles tends to infinity.

Suppose you have three cards: {A,B,C}. Assume that your cards begin in the following order: A,B,C. Then after one shuffle you have following combinations:

{A,B,C}, {A,B,C}, {A,B,C} #You get this if choose the same RN twice.
{A,C,B}, {A,C,B}
{C,B,A}, {C,B,A}
{B,A,C}, {B,A,C}


Hence, the probability of card A of being in position {1,2,3} is {5/9, 2/9, 2/9}.

If we shuffle the cards a second time, then:

Pr(A in position 1 after 2 shuffles) = 5/9*Pr(A in position 1 after 1 shuffle)
+ 2/9*Pr(A in position 2 after 1 shuffle)
+ 2/9*Pr(A in position 3 after 1 shuffle)


This gives 0.407.

Using the same idea, we can form a recurrence relationship, i.e:

Pr(A in position 1 after n shuffles) = 5/9*Pr(A in position 1 after (n-1) shuffles)
+ 2/9*Pr(A in position 2 after (n-1) shuffles)
+ 2/9*Pr(A in position 3 after (n-1) shuffles).


Coding this up in R (see code below), gives probability of card A of being in position {1,2,3} as {0.33334, 0.33333, 0.33333} after ten shuffles.

R code

## m is the probability matrix of card position
## Row is position
## Col is card A, B, C
m = matrix(0, nrow=3, ncol=3)
m[1,1] = 1; m[2,2] = 1; m[3,3] = 1

## Transition matrix
m_trans = matrix(2/9, nrow=3, ncol=3)
m_trans[1,1] = 5/9; m_trans[2,2] = 5/9; m_trans[3,3] = 5/9

for(i in 1:10){
old_m = m
m[1,1] = sum(m_trans[,1]*old_m[,1])
m[2,1] = sum(m_trans[,2]*old_m[,1])
m[3,1] = sum(m_trans[,3]*old_m[,1])

m[1,2] = sum(m_trans[,1]*old_m[,2])
m[2,2] = sum(m_trans[,2]*old_m[,2])
m[3,2] = sum(m_trans[,3]*old_m[,2])

m[1,3] = sum(m_trans[,1]*old_m[,3])
m[2,3] = sum(m_trans[,2]*old_m[,3])
m[3,3] = sum(m_trans[,3]*old_m[,3])
}
m

• +1. That demonstrates that the probability for a given card to end up in a given position approximates the expected ratio as the number of shuffles increases. However, the same would also be true of an algorithm that just rotates the array once by a random amount: All cards have an equal probability to end up in all positions, but there is still no randomness at all (the array remains sorted). – Thilo Sep 27 '10 at 1:30
• @Thilo: Sorry I don't follow your comment. An "algorithm rotates by a random amount" but there's still "no randomness"? Could you explain further? – csgillespie Sep 27 '10 at 8:21
• If you "shuffle" an N-element array by rotating it between 0 and N-1 positions (randomly), then every card has exactly the same probability to end up in any of the N positions, but 2 is still always located between 1 and 3. – Thilo Sep 27 '10 at 8:26
• @Thio: Ah, I get your point. Well you can work out the probability (using exactly the same idea as above), for the Pr(A in position 2) and Pr(A in position 3) - dito for cards B and C. You will see that all probabilities tend to 1/3. Note: my answer only gives a particular case, whereas @whuber nice answer gives the general case. – csgillespie Sep 27 '10 at 8:46

One way to see that you won't get a perfectly uniform distribution is by divisibility. In the uniform distribution, the probability of each permutation is $1/n!$. When you generate a sequence of $t$ random transpositions, and then collect sequences by their product, the probabilities you get are of the form $A/n^{2t}$ for some integer $A$. If $1/n! = A/n^{2t}$, then $n^{2t}/n! = A$. By Bertrand's Postulate (a theorem), for $n \ge 3$ there are primes which occur in the denominator and which do not divide $n$, so $n^{2t}/n!$ is not an integer, and there isn't a way to divide the transpositions evenly into $n!$ permutations. For example, if $n=52$, then the denominator of $1/52!$ is divisible by $3, 5, 7, ..., 47$ while the denominator of $1/52^{2t}$ is not, so $A/52^{2t}$ can't reduce to $1/52!$.

How many do you need to approximate a random permutation well? Generating a random permutation by random transpositions was analyzed by Diaconis and Shahshahani using representation theory of the symmetric group in

Diaconis, P., Shahshahani, M. (1981): "Generating a random permutation with random transpositions." Z. Wahrsch. Verw. Geb. 57, 159–179.

One conclusion was that it takes $\frac 12 n \log n$ transpositions in the sense that after $(1-\epsilon) \frac12 n \log n$ the permutations are far from random, but after $(1+\epsilon) \frac 12 n \log n$ the result is close to random, both in the sense of total variation and $L^2$ distance. This type of cutoff phenomenon is common in random walks on groups, and is related to the famous result that you need $7$ riffle shuffles before a deck becomes close to random.

Bear in mind I am not a statistician, but I'll put my 2cents.

I made a little test in R (careful, it's very slow for high numTrials, the code can probably be optimized):

numElements <- 1000
numTrials <- 5000

swapVec <- function()
{
vec.swp <- vec

for (i in 1:numElements)
{
i <- sample(1:numElements)
j <- sample(1:numElements)

tmp <- vec.swp[i]
vec.swp[i] <- vec.swp[j]
vec.swp[j] <- tmp
}

return (vec.swp)
}

# Create a normally distributed array of numElements length
vec <- rnorm(numElements)

# Do several "swapping trials" so we can make some stats on them
swaps <- vec
prog <- txtProgressBar(0, numTrials, style=3)

for (t in 1:numTrials)
{
swaps <- rbind(swaps, swapVec())
setTxtProgressBar(prog, t)
}


This will generate a matrix swaps with numTrials+1 rows (one per trial + the original) and numElements columns (one per each vector element). If the method is correct the distribution of each column (i.e. of the values for each element over the trials) should not be different from the distribution of the original data.

Because our original data was normally distributed we would expect all the columns not to deviate from that.

If we run

par(mfrow= c(2,2))
# Our original data
hist(swaps[1,], 100, col="black", freq=FALSE, main="Original")
# Three "randomly" chosen columns
hist(swaps[,1], 100, col="black", freq=FALSE, main="Trial # 1")
hist(swaps[,257], 100, col="black", freq=FALSE, main="Trial # 257")
hist(swaps[,844], 100, col="black", freq=FALSE, main="Trial # 844")


We get:

which looks very promising. Now, if we want to statistically confirm the distributions do not deviate from the original I think we could use a Kolmogorov-Smirnov test (please can some statistician confirm this is right?) and do, for instance

ks.test(swaps[1, ], swaps[, 234])


Which gives us p=0.9926

If we check all of the columns:

ks.results <- apply(swaps, 2, function(col){ks.test(swaps[1,], col)})
p.values <- unlist(lapply(ks.results, function(x){x$p.value})  And we run hist(p.values, 100, col="black")  we get: So, for the great majority of the elements of the array, your swap method has given a good result, as you can also see looking at the quartiles. 1> quantile(p.values) 0% 25% 50% 75% 100% 0.6819832 0.9963731 0.9999188 0.9999996 1.0000000  Note that, obviously, with a lower number of trials the situation is not as good: 50 trials 1> quantile(p.values) 0% 25% 50% 75% 100% 0.0003399635 0.2920976389 0.5583204486 0.8103852744 0.9999165730  100 trials  0% 25% 50% 75% 100% 0.001434198 0.327553996 0.596603804 0.828037097 0.999999591  500 trials  0% 25% 50% 75% 100% 0.007834701 0.504698404 0.764231550 0.934223503 0.999995887  Here's how I am interpreting your algorithm, in pseudo code: void shuffle(array, length, num_passes) for (pass = 0; pass < num_passes; ++pass) for (n = 0; n < length; ++) i = random_in(0, length-1) j = random_in(0, lenght-1) swap(array[i], array[j]  We can associate a run of this algorithm with a list of$2 \times length \times num\_passes$integers, namely the integers returned by random_in() as the program runs. Each of these integers is in$[0, length-1]$, and so has$length$possible values. Call one of these lists a trace of the program. That means there are$length ^ {2 \times length \times num\_passes}$such traces, and each trace is equally likely. We can also associate with each trace a permutation of the array. Namely, the permutation at the end of the run associated with the trace. There are$length !$possible permutations.$length ! < length ^ {2 \times length \times num\_passes}$so in general a given permutation is associated with more than one trace. Remember, the traces are all equally likely, so for all permutations to be equally likely, any two permutations must be associated with the same number of traces. If that is true, then we must have$length ! \bigm| length ^ {2 \times length \times num\_passes}$. Pick any prime$p$such that$p < length$, but such that$p \nmid length$, which you can do for any$length > 2$. Then$p \bigm| length!$but does not divide$length ^ {2 \times length \times num\_passes}$. It follows that$length ! \nmid length ^ {2 \times length \times num\_passes}$and so all permutations cannot be equally likely if$length > 2$. Does such a prime exist? Yes. If$length$were divisible by all primes$p < length$, then$length-1$must be prime, but then$length-1$would be such a prime that is less than but does not divide$length$. Compare this to Fisher-Yates. In the first iteration, you make a choice among$length$options. The second iteration has$length-1$options, and so on. In other words you have$length !$traces, and$length! \bigm| length!\$. It's not hard to show that each trace results in a different permutation, and from there it is easy to see that Fisher-Yates generates each permutation with equal probability.