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)?
 A: 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

A: 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.
A: 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.
A: 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.
