THE THORP SHUFFLE MIXES IN LINEAR TIME
A reader's version for high-school students
ABSTRACT
How many times must you shuffle a deck before it is truly random?
For the Thorp shuffle the answer is a constant times log2(n), where n is
the number of cards. This note proves that statement from scratch.
The proof has three parts: a counting lower bound, a conditional-flow
upper bound for most cards, and a Fourier/representation argument that
recovers the few missing cards.
Throughout, n = 2^d. Positions are binary strings of length d. One
"shuffle" means one physical Thorp shuffle. One "sweep" means d physical
shuffles, one in each coordinate direction. The main theorem says that
1600d physical shuffles suffice for large d, while at least 2d - O(1)
are necessary. Thus the mixing time is Theta(d) = Theta(log n).
1. THE SHUFFLE AND THE QUESTION
Split a deck of n cards into two equal halves. Pair the cards in
corresponding positions. For each pair flip an independent fair coin and
decide whether to keep or swap the pair. Interleave the pairs to get the
new deck. This is the Thorp shuffle.
Label positions by vectors x = (x_1,...,x_d) in F_2^d. Let R be the
cyclic rotation
R(x_1,...,x_d) = (x_2,...,x_d,x_1).
Let S be the operation that independently fixes or swaps each pair
{x, x+e_1}, where e_1 = (1,0,...,0). One physical shuffle is RS.
Equivalently, if we undo the rotation after each step, then step t uses
the coordinate i_t = 1 + (t mod d). After d steps every coordinate has
been used exactly once. That is one sweep.
For a card, each sweep refreshes each of its d bits once with a fair coin.
So after d physical shuffles one card is exactly uniform. The hard part
is the joint law of all cards.
Let q_d be the law of one physical shuffle as a permutation of the
positions. Let U_n be the uniform law on the symmetric group S_n.
Total variation distance is
||mu - nu||_TV = (1/2) sum_g |mu(g) - nu(g)|.
If this distance is at most 1/4, then no event has its probability
changed by more than 1/4. Let t_mix(d) be the least t with
||q_d^{*t} - U_n||_TV <= 1/4.
THEOREM (main).
As d -> infinity,
|| q_d^{*(1600d)} - U_n ||_TV -> 0.
For every integer t >= 0,
|| q_d^{*t} - U_n ||_TV >= 1 - 2^{tn/2}/n!.
Therefore, for large d,
2d - O(1) <= t_mix(d) <= 1600d,
and in particular t_mix(d) = Theta(d).
The constant 1600 is not claimed to be sharp. The point is that it does
not grow with d.
2. THE LOWER BOUND: COUNT THE COINS
One physical shuffle uses n/2 independent fair coins. After t physical
shuffles there are at most 2^{tn/2} possible coin strings, hence at most
that many possible permutations in the support of q_d^{*t}.
Let A be that support. Then q_d^{*t}(A) = 1, while
U_n(A) <= |A|/n! <= 2^{tn/2}/n!.
So
||q_d^{*t} - U_n||_TV >= q_d^{*t}(A) - U_n(A)
>= 1 - 2^{tn/2}/n!.
If the distance is at most 1/4, then 2^{tn/2} >= 3n!/4. Taking logs,
t >= (2/n) log2(3n!/4).
Use (n/e)^n <= n! <= n^n. Since log2 n = d, this gives
t >= 2d - O(1). That is the lower bound. It is information-theoretic:
you need enough random bits to describe a random permutation.
3. THE COORDINATE MODEL
Write positions as binary strings. Let X_t be the permutation after
undoing the deterministic rotations. Then
X_{t+1} = Q_t X_t,
where Q_t is a matching switch in coordinate i_t = 1 + (t mod d).
The Q_t are independent. A block of d consecutive steps is a sweep:
each coordinate is used once. At the end of a sweep we are back in the
physical frame.
Because uniform measure is invariant under every position permutation,
the distance from uniform does not depend on the initial deck. So we
may start from the identity deck.
For completeness, mixing exists in each fixed dimension. A full sweep
has positive probability of being the identity, and positive probability
of making any specified cube-edge transposition. Edge transpositions on
the connected cube generate S_n. Hence some number of sweeps can reach
every permutation, and after padding with identity sweeps the kernel at
one common length gives positive mass to every group element. It
therefore dominates a positive multiple of U_n, and iteration proves
convergence. The theorem is about the rate, uniform in d.
4. CONDITIONAL FLOW
Fix an ordered list (c_1,...,c_k) of distinct labels, chosen
independently of the shuffle. Expose the paths of the first ell cards.
The next card c_{ell+1} is the tagged card. Let B_s be the set of
positions not occupied by the exposed cards, and let m = |B_s| = n - ell.
Define
p_s(x) = P(X_s(c_{ell+1}) = x | history),
w_s(x) = p_s(x) - (1/m) 1_{B_s}(x).
The vector w_s is supported on B_s, sums to zero, and has squared norm at
most 1.
On a switch with two free positions, the two masses are averaged. On a
switch with one free position, the mass follows the route fixed by the
exposed path. Why? The visited coins are fixed by the exposed paths.
All unvisited coins remain independent and fair. Averaging and transport
also carry the uniform vector on B_s to the uniform vector on B_{s+1}.
Let A_i average a function over every coordinate-i pair. Then
w_{s+1} = A_{i_s} w_s + xi_s,
where xi_s has conditional mean zero. On a one-free-position edge whose
free endpoint is x, the noise is
+-(w_s(x)/2) (delta_x - delta_{x^{(i_s)}}),
and signs on distinct edges are independent.
Define the energy
E_s = E ||w_s||_2^2,
and the recent-noise term
L_s = sum_x w_s(x)^2 1_{x^{(i_s)} not in B_s},
b_s = E L_s.
LEMMA (one-sweep energy identity).
For every t >= d,
E_t = sum_{j=1}^d 2^{-j} b_{t-j}.
Proof. Expand the recurrence over the last d layers. The product of all
d coordinate averages maps every vector to its spatial mean, so it kills
w_{t-d}. Noise from different times has zero expected cross product,
because the later term has conditional mean zero and the intervening
averages are deterministic. At one fixed time, noises on distinct edges
also have zero expected cross product, because their signs are
independent conditional on the past. The noise created at time t-j is
followed by averaging in j-1 other coordinates. Its two masses spread
over disjoint subcubes of size 2^{j-1}, so its squared norm is
2^{-j} w_{t-j}(x)^2. Summing gives the identity.
5. WHY RECENT NOISE CONTRACTS
Fix s >= d-1 and put a = s-d+1. Let H_0,H_1 be the two halves defined
by coordinate i_s. The updates from a to s use every other coordinate
once and never cross these halves. Conditional on the history at time a,
the future switch arrays in the two halves are independent.
A position in one half is free at time s with probability equal to the
fraction of free positions in the opposite half at time a:
P(y in B_s | history at a)
= |B_a cap H_{1-h}| / (n/2), for y in H_h.
This is because each remaining coordinate of that position is refreshed
by a conditionally fair coin during the intervening d-1 updates.
Now choose the initial list uniformly at random, independent of the
shuffle. If m >= n/8, the free set B_a is a uniform m-subset of V.
Its count in either half has mean m/2. A bounded-difference martingale
argument (Hoeffding/Azuma) gives
P(|B_a cap H_h| < m/4) <= e^{-m/32}.
With
beta = 1/16, eta = 2 e^{-n/256},
both half densities are at least beta except on a bad event of
probability at most eta. The bad event may be correlated with the
energy, but the energy is at most 1. Therefore
b_s <= (1 - beta) E_s + eta.
Let gamma = 31/32 = 1 - beta/2. Combining with the energy identity gives
E_t <= gamma^{t-2d} + eta/beta
= (31/32)^{t-2d} + 32 e^{-n/256}.
For t <= 2d the bound follows from E_t <= 1. For t > 2d, use induction
in the energy identity. The geometric part is at most gamma^{t-2d}
because
(1-beta) sum_{j>=1} (2 gamma)^{-j}
= (1-beta)/(2 gamma - 1) = 1.
The constant part is at most eta/beta.
6. FROM ENERGY TO ORDERED CARD MARGINALS
We need a standard sequential comparison. Let U_inj be the uniform law
on tuples of distinct positions.
LEMMA (sequential total variation).
For a random tuple (Z_1,...,Z_k) of distinct positions,
|| L(Z_1,...,Z_k) - U_inj ||_TV
<= sum_{j=1}^k E || L(Z_j | Z_1,...,Z_{j-1})
- U_{V \ {Z_1,...,Z_{j-1}}} ||_TV.
Proof. At step j compare the true joint j-tuple with the law that
samples the first j-1 coordinates from their true law and then appends a
uniform available position. The distance is exactly the jth summand.
Applying the same uniform extension to a uniform prefix gives a uniform
j-tuple, and cannot increase total variation. The triangle inequality
and induction prove the lemma.
PROPOSITION (averaged ordered marginals).
Let T = 200d and 1 <= k <= 7n/8. Choose an ordered list
C = (c_1,...,c_k) uniformly among all lists of distinct labels,
independent of the shuffle. Then
E_C || L((X_T(c_1),...,X_T(c_k)) | C) - U_inj ||_TV
<= eps_n,
where
eps_n = (n^{3/2}/2)
sqrt( (31/32)^{198d} + 32 e^{-n/256} )
= O(n^{-5/2}).
Proof. Fix C and apply the sequential lemma to Z_j = X_T(c_j). In the
jth summand, condition further on the paths through time T of the first
j-1 cards. These paths determine the prefix of endpoints, so the
refinement part of the lemma applies.
Let w_T^{(j)} be the resulting centered flow, with ell = j-1 and
m_j = n-j+1. Its conditional distance is (1/2)||w_T^{(j)}||_1.
Cauchy-Schwarz gives
(1/2) E ||w_T^{(j)}||_1
<= (1/2) sqrt( m_j E ||w_T^{(j)}||_2^2 ).
Each stage leaves at least n/8 free positions, so the energy bound
applies. Since m_j <= n and k <= n, summing over j gives the displayed
bound. Finally, 198 log2(32/31) > 8, so (31/32)^{198d} <= n^{-8},
which gives eps_n = O(n^{-5/2}).
7. EIGHT BLOCKS AND TRIMMING
Partition the labels into eight blocks M_1,...,M_8 of size n/8. Each
complement V \ M_a has size 7n/8. Giving it an independent uniform order
produces exactly the list in the proposition. Therefore some
deterministic partition satisfies
delta_n := sum_{a=1}^8
|| L(mu on V \ M_a) - U_inj ||_TV
<= 8 eps_n = o(1),
where mu = q_d^{*(200d)}.
Fix such a partition. Let H_a be the subgroup of permutations that fix
every label outside M_a. Two permutations agree on the complement
exactly when they lie in the same left coset gH_a. So the complementary
image tuple identifies the coset gH_a, and its uniform law is the coset
pushforward of uniform measure on G = S_n.
For each a, remove every coset C in G/H_a with
mu(C) > 2 |H_a|/|G|.
Let E be what remains after all eight removals. A bad coset satisfies
mu(C) <= 2 ( mu(C) - |H_a|/|G| ).
The sum of positive excesses is the marginal total-variation distance, so
the union bound gives
M := mu(E) >= 1 - 2 delta_n.
For large n, M >= 1/2. Define nu = mu( . | E). Then
||nu - mu||_TV = 1 - M = o(1),
and for every a and every g,
nu(gH_a) <= (2/M) |H_a|/|G| <= 4 |H_a|/|G|.
A coset removed for that subgroup has no retained mass; every other coset
had the required original cap. So all eight coset systems have the same
cap under one common law nu.
8. REPRESENTATION DEGREES
A representation of a finite group G is a homomorphism rho from G to
unitary matrices. It is irreducible if it has no nontrivial invariant
subspace. Its degree D is the matrix size. For a probability measure
nu on G, the Fourier transform is
hat nu(rho) = sum_g nu(g) rho(g).
For the symmetric group S_s, irreducible representations are indexed by
partitions lambda of s. Draw lambda as a Young diagram. A standard
Young tableau fills the boxes with 1,...,s, increasing along rows and
columns. The degree f^lambda is the number of such tableaux. The
hook-length formula is
f^lambda = s! / prod_{z in lambda} h(z),
where h(z) counts the box z plus all boxes to its right and below it.
The row (s) and column (1^s) are the trivial and sign representations,
both of degree 1.
LEMMA (reciprocal representation degrees).
The reciprocal-degree sums satisfy
C_1 := sup_{s>=1} sum_{lambda |- s} (f^lambda)^{-1} < infinity,
and
sum_{lambda |- s, lambda not in {(s),(1^s)}} (f^lambda)^{-1}
-> 0 as s -> infinity.
Proof. Let p(u) be the number of partitions of u. The generating
function gives p(u) <= e^{3 sqrt u}. For a diagram of size s, let L be
the larger of its first row and column lengths. Transpose so the first
row has length L, and put u = s-L.
If 1 <= u <= s/4, then
f^lambda >= C(L,u) >= (L/u)^u >= 3^u.
Indeed, place 1,...,u in the first u top-row boxes; choose any u labels
from the remaining L labels for the lower diagram, filling it in
row-major order; and fill the remaining top row increasingly. Every
lower column lies under one of the first u top boxes, so all inequalities
hold. There are at most 2p(u) shapes with this tail size, including
transposes. Their reciprocal contributions are bounded by the summable
sequence 2p(u)3^{-u}; for each fixed positive u, the binomial lower
bound tends to infinity with s. Dominated convergence makes the whole
contribution tend to zero.
If u > s/4 and L >= s/8, retain the top row and the first
b = floor(s/16) boxes below it, read row by row. The retained boxes form
a Young diagram. Every standard tableau of the retained diagram extends
to the full diagram by filling the remaining boxes in the same row order
with the larger labels. In the retained diagram, fix 1,...,b in the
first b top-row boxes and choose the b lower labels from the remaining L
labels. Since L >= 2b, this gives
f^lambda >= C(L,b) >= 2^b
for all sufficiently large s. There are at most e^{3 sqrt s} diagrams,
so their reciprocal contribution tends to zero.
If L < s/8, every hook has length less than s/4. The hook formula and
s! >= (s/e)^s give f^lambda >= (4/e)^s, again dominating the partition
count. These cases cover all shapes except the row and column. Their
degrees are one, so the full sums tend to two. The finitely many initial
sums are finite, which proves the uniform bound C_1. At s=1 the row and
column coincide and are counted once.
9. THE EIGHT-BLOCK FOURIER BOUND
Let n be a positive multiple of eight and partition the labels of
G = S_n into sets M_1,...,M_8 of size n/8. Let H_a fix every label
outside M_a. Suppose a probability measure nu on G satisfies
nu(gH_a) <= 4 |H_a|/|G| for all g and all a.
PROPOSITION (eight-block Fourier bound).
Every irreducible unitary representation rho of G, of degree D,
satisfies
|| hat nu(rho) ||_op <= A D^{-5/16},
where A = sqrt(32 C_1).
Proof. The subgroups H_1,...,H_8 commute and form
H = H_1 x ... x H_8. Restrict rho to H. The restriction decomposes
orthogonally into spaces
(V_{pi_1} tensor ... tensor V_{pi_8}) tensor A_{pi},
where A_{pi} records arbitrary multiplicity and H acts trivially on it.
Every occurring tensor type has product of factor degrees at most D.
Assign its whole isotypic space to one factor whose degree is at most
D^{1/8}. This gives an orthogonal decomposition
V_rho = E_1 + ... + E_8.
Each E_a is H-invariant, and all H_a-irreducible types occurring in it
have degree at most D^{1/8}.
Fix u_a in E_a and let
W_a = span{ rho(h) u_a : h in H_a }.
Let u_{a,pi} be the projection of u_a onto the full H_a-isotypic
component V_pi tensor A_pi, with b = dim pi. Its orbit lies in the image
of
End(V_pi) -> V_pi tensor A_pi,
Q |-> (Q tensor I) u_{a,pi}.
Its dimension is at most b^2, independently of multiplicity. Grouping
equivalent copies first and counting each irreducible type once gives
dim W_a <= sum_{pi in hat H_a, dim pi <= D^{1/8}} (dim pi)^2
<= D^{3/8} sum_{pi in hat H_a} (dim pi)^{-1}
<= C_1 D^{3/8}.
The uniform bound C_1 applies to every block size.
The space W_a is H_a-invariant, so rho(g)W_a depends only on gH_a. Let
P_W be orthogonal projection. By Schur's lemma and the trace,
E_{g uniform in G} P_{rho(g)W_a} = (dim W_a / D) I.
Indeed this average commutes with rho(G), so it is scalar on each
irreducible; its trace is dim W_a.
Apply the coset cap to the nonnegative coset-constant function
|| P_{rho(g)W_a} z ||^2. For any unit vector z,
| <z, hat nu(rho) u_a> |^2
<= E_{g in nu} | <z, rho(g) u_a> |^2
<= ||u_a||^2 E_{g in nu} || P_{rho(g)W_a} z ||^2
<= 4 (dim W_a / D) ||u_a||^2
<= 4 C_1 D^{-5/8} ||u_a||^2.
For u = sum_a u_a, sum these bounds after taking square roots, and use
sum_a ||u_a|| <= sqrt(8) ||u||. This proves the proposition.
10. FINAL FOURIER COMPARISON
For any probability measure sigma on G,
4 || sigma - U_G ||_TV^2
<= sum_{rho != trivial} D_rho
|| hat sigma(rho) ||_HS^2.
Here ||B||_HS^2 = tr(B^*B) is the ordinary unnormalized Hilbert-Schmidt
norm. To see the inequality, Cauchy-Schwarz bounds the left side by
|G| sum_g |sigma(g) - |G|^{-1}|^2. The regular-character orthogonality
identity gives, for every real mass function a,
sum_rho D_rho || sum_g a(g) rho(g) ||_HS^2
= |G| sum_g a(g)^2.
Apply this to a = sigma - U_G. Its trivial coefficient is zero and
uniform measure has zero transform in every nontrivial irreducible.
Convolution multiplies Fourier matrices. The operator norm is
submultiplicative and ||B||_HS <= sqrt(D) ||B||_op, with no normality
assumption. For sigma = nu^{*8}, the nontrivial non-sign contribution is
at most
sum_{rho != trivial, sign} D_rho
|| hat nu(rho)^8 ||_HS^2
<= A^{16} sum_{rho != trivial, sign} D_rho^{-3}
-> 0,
by the reciprocal-degree lemma, since D^{-3} <= D^{-1}.
The sign representation has degree one and needs a separate argument.
In the final physical shuffle of a positive-time block, condition on
every coin but one. Flipping that last coin changes the resulting
permutation by a transposition, so the original block law mu has zero
expected sign. Since nu = mu 1_E / M,
| hat nu(sign) |
= (1/M) | sum_{g notin E} mu(g) sign(g) |
<= (1-M)/M -> 0.
Its contribution after eight convolution factors also tends to zero.
Thus || nu^{*8} - U_G ||_TV -> 0.
Finally, replacing independent convolution factors one at a time costs
at most their individual total-variation distances:
|| mu^{*8} - U_G ||_TV
<= 8(1-M) + || nu^{*8} - U_G ||_TV
-> 0.
Each convolution factor is 200d physical shuffles, so
mu^{*8} = q_d^{*(1600d)}. Translation invariance gives the same bound
from every deterministic initial deck. Together with the support lower
bound and fixed-dimensional convergence, this proves the theorem.
11. WHAT TO REMEMBER
The lower bound is counting: t shuffles produce at most 2^{tn/2}
outcomes, while a random permutation needs about n log2 n bits.
The upper bound is a contraction argument: most cards become jointly
uniform after 200d physical shuffles, and a representation-theoretic
"eight block" argument recovers the dependencies involving the omitted
1/8 of the cards. Eight such blocks give 1600d. The exact constant is
not important; the order Theta(d) is.
FURTHER READING
E. O. Thorp, Nonrandom shuffling with applications to the game of Faro,
J. Amer. Statist. Assoc. 68 (1973), 842-847.
B. Morris, The mixing time of the Thorp shuffle, SIAM J. Comput. 38
(2008), 484-504.
B. Morris, Improved mixing time bounds for the Thorp shuffle,
Combin. Probab. Comput. 22 (2013), 118-132.
A. Czumaj and B. Vocking, Thorp shuffling, butterflies, and
non-Markovian couplings, ICALP 2014.
P. Diaconis and M. Shahshahani, Generating a random permutation with
random transpositions, Z. Wahrsch. Verw. Gebiete 57 (1981), 159-179.