Optimal-order mixing of the Thorp shuffle

1 Introduction

A single card can become uniformly distributed long before a shuffled deck is close to uniform. The Thorp shuffle makes this distinction particularly clear. Split a deck into two equal halves, pair cards that occupy corresponding positions, and choose the order of each pair by an independent fair coin. Interleaving the pairs produces the new deck. For a deck of size n=2^d, every individual card is uniform after d shuffles. The question is how much longer it takes to randomize the joint order of all n cards.

We prove that a constant multiple of d shuffles suffices. Let q_d denote the law of one complete shuffle as a permutation of the positions, and let U_n be uniform on the symmetric group S_n. Convolution is the law of composition of independent permutations, and

Write t_{\mathrm{mix}}(d) for the least nonnegative integer t at which \|q_d^{*t}-U_n\|_{\mathrm{TV}}\le1/4.

One step in this theorem is one complete physical shuffle. Changing the initial deck translates the permutation law, so its distance from uniform is unchanged. The upper bound is asymptotic in d; the constant 1600 is an explicit choice in the proof, not a proposed sharp constant. Determining a leading constant or a cutoff remains a separate question.

The full-law conclusion controls every statistic of the final deck and, in particular, the joint positions of any deterministic selection of cards, by contraction of total variation under a map. It also shows that this shuffle produces an approximately uniform permutation using O(n\log n) independent fair bits: each of the O(\log n) steps uses n/2 bits. The support bound identifies the same order of randomness as necessary for this construction.

1.1 History and related work

Thorp introduced the model in his 1973 study of the effects of nonrandom shuffling on Faro [19, Section 3.1]. Its simple pairwise rule conceals a difficult question about dependence across the whole deck. For dyadic decks, Morris obtained the first bound polynomial in d, namely O(d^{44}), using evolving sets and a chameleon process that tracks conditional card probabilities [10]. Montenegro and Tetali sharpened this analysis to O(d^{29}) by spectral-profile and evolving-set methods [9, Theorems 6.14–6.15].

A second line of attack used relative entropy. Morris related entropy decay to pairwise card collisions and obtained O((\log n)^4) shuffles for every even n [12]. He then proved O(d^3) for n=2^d by a stronger contraction over a block of d shuffles [13]. Theorem 1.1 attains the logarithmic order discussed in Jonasson’s notes [7, Section 5] for dyadic decks. The corresponding optimal-order question for arbitrary even deck sizes is not addressed by the cube argument here.

Card marginals have also been studied for their own sake and for cryptography. Morris, Rogaway and Stegers analyzed the joint images of ordered input lists in their construction of enciphering schemes on small domains [14, Theorem 1]. Czumaj and Vöcking proved that O((\log n)^2) Thorp shuffles mix any fixed fraction of the cards strictly below one, using non-Markovian couplings and butterfly-network estimates [2]. Our first stage needs a different, averaged guarantee: after 200d shuffles, the joint positions of a uniformly sampled list of 7n/8 labels are close to uniform on average over the list. The passage from these averaged marginals to the full law is a separate part of the proof.

The conditional-flow method has several predecessors. In the original Thorp analysis, Morris used a chameleon process to represent the law of a tagged card conditional on the set-valued trajectory of a collection containing that card; the conditional probabilities average across available pairs and are transported across pairs with only one available position [10, Section 4, Lemma 3]. Hoang, Morris and Rogaway later used conditional squared energy and a sequential total-variation comparison to analyze the swap-or-not shuffle [5, Section 3, Theorem 3 and Lemma 4]. Their uniformly random keys make the next partner uniform over the underlying group. Here the directions are prescribed cyclically. The exact identity retaining only one sweep of noise, together with independence of the two opposite halves during the intervening updates, provides the needed contraction.

For the final passage to the group, we use complete reducibility and Schur orthogonality, together with the tableau description and hook-length formula for symmetric-group degrees [4, 18]. The inverse-degree summability needed here is a special case of a stronger theorem of Liebeck and Shalev [8, Theorem 2.6]; we give an elementary proof of the required case. The final total-variation comparison is the finite-group Fourier method of Diaconis and Shahshahani [3, Section 3, Lemma 14]. The additional step here is to turn bounds on complementary card marginals into operator bounds for matrix Fourier transforms, using a common trimming and the span of one vector’s orbit under a block subgroup.

1.2 Proof overview

Label positions by binary strings of length d, the vertices of a cube. Undoing a deterministic cyclic rotation turns successive shuffles into switch layers along successive coordinate directions: the two positions in one switch differ only in the active bit. A sweep through all d directions takes d physical shuffles. Section 2 gives this change of frame and proves the support lower bound.

For the upper bound, first choose an ordered list of 7n/8 labels uniformly and independently of the shuffle. Expose these cards one at a time. Conditional on the paths already exposed, the next card has a probability distribution on the remaining positions. On a switch with two remaining positions its masses are averaged; on a switch with only one remaining position its mass follows the route fixed by the exposed path. Subtracting the uniform distribution on the remaining positions gives a centered flow.

A full sweep of coordinate averages annihilates every mean-zero vector. The actual conditional update differs from coordinate averaging only by noise with conditional mean zero on switches with one available position. Consequently the expected squared norm at the present time is an exact weighted sum of noise contributions from the preceding sweep, as proved in Section 3. To control those contributions, Section 4 looks back from the next active direction through the d-1 other directions. These updates act independently inside the two opposite halves. A flow value in one half is therefore independent of the availability of its partner in the other half, conditional on the history at the start of this interval. A random list leaves a positive fraction of available positions in both halves with high probability. This yields a scalar energy recurrence and then the averaged marginal estimate of Proposition 4.2 at time 200d.

The second stage recovers the dependencies involving the omitted cards. Partition the labels into eight equal blocks, choosing the partition so that the eight complementary marginals have small total error. Each marginal is a distribution on cosets of the subgroup that permutes its omitted block. Removing a small amount of probability makes all eight coset densities bounded under one common law. Restrict an irreducible representation of S_n to the product of these eight subgroups. Each tensor constituent has at least one factor of small degree. Even with arbitrary multiplicity, the orbit of a single vector under the corresponding block subgroup spans a space of controlled dimension. Averaging projections onto its translates then bounds each Fourier operator. Sections 5 and 6 supply the degree estimate and this transfer; Section 7 combines eight independent time intervals of 200d shuffles and handles parity to complete the proof.

The companion on random coordinate frames gives an independent 32800d proof [17, Theorem 1.1]. The conditional-information companion obtains bounds uniform over specified input lists [15, Theorem 1.1], while the Fourier companion develops general transfers from complementary coset bounds [16, Theorem 1.1]. The present proof uses none of these results: its cyclic-coordinate estimate and eight-block transfer are proved here.

2 The model and the support obstruction

Let d\ge1 be an integer. Identify the positions with V=\mathbb F_2^d, write e_1,\ldots,e_d for its standard basis, and label each card by its initial position. A permutation maps a label to its current position; composition means (gh)(x)=g(h(x)). With

one physical shuffle is RS, where S independently fixes or swaps each pair \{x,x+e_1\}. Thus its action is

with independent fair bits \xi_u.

Let S_1,S_2,\ldots be independent copies of S. If Y_t=RS_tY_{t-1} is the physical permutation, started from Y_0=\mathrm{id}, and X_t=R^{-t}Y_t, then

The Q_t are independent matching switches in the directions i_t=1+(t\bmod d), since R^{-t}e_1=e_{i_t}. A block whose length is a multiple of d ends in the physical frame. We analyze this cyclic-coordinate process, writing \Pi_t=X_t.

During one full sweep each coordinate of a given card is refreshed exactly once by a coin fair conditional on the card’s earlier path. These refreshed bits form a uniform binary string. Restoring the deterministic rotation preserves uniformity, proving the one-card assertion in the introduction.

Uniform measure is invariant under every position permutation. Starting from any fixed deck translates the law started at the identity and leaves its total-variation distance from uniform unchanged. Moreover, t shuffles have at most 2^{tn/2} possible coin strings. Their support A therefore gives

Distance at most 1/4 forces 2^{tn/2}\ge3n!/4, proving the rounded lower bound. The inequalities (n/e)^n\le n!\le n^n give its asymptotic form.

For completeness, mixing exists in each fixed dimension. A full coordinate sweep has positive probability of being the identity, and of making any specified single cube-edge transposition. Edge transpositions on the connected cube generate S_n. Every permutation can therefore be obtained in some number of sweeps; identity sweeps pad these representations to one common length. The kernel at that length assigns positive mass to every group element and hence dominates a positive multiple of U_n. Iteration proves convergence. This finite-dimensional observation does not assert the numerical bound 1600d for every d.

3 A conditional flow with one sweep of memory

Choose an ordered list (c_1,\ldots,c_k) of distinct labels independently of the shuffle. We will ultimately average over uniform lists, but first fix the list. For 0\le\ell<k\le n, expose the trajectories of its first \ell cards, and call the next card the tagged card. Let \mathcal F_s be the sigma-field generated by the list and the exposed paths through time s. Let B_s be the positions not occupied by those \ell cards, so m=|B_s|=n-\ell is constant. Define

The vector w_s is supported on B_s, has total sum zero, and has squared norm at most one. Averaging and transport also underlie Morris’s chameleon representation [10, Section 4, Lemma 3]. The next lemma identifies the unexposed coins directly for labeled trajectories.

Let A_i average a function across every coordinate-i pair, and write x^{(i)}=x+e_i. Put

Here the expectation includes the initial list when it is random. Before the next exposed positions are revealed, the transport choices on edges with one free position are independent fair coins. Hence

On a one-free-position edge whose free endpoint is x, its noise is

where \delta_x is unit mass at x, and the signs are independent on distinct such edges. Other edges contribute no noise.

4 Why the recent noise contracts

We bound the recent-noise term b_s in terms of the energy E_s, then use the resulting decay to control ordered card marginals.

Fix s\ge d-1, put a=s-d+1, and 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 do not cross these halves. Conditional on \mathcal F_a, the initial free sets B_a\cap H_h and centered flows w_a|_{H_h} are fixed, and the future switch arrays in the two halves are independent. The pair (B_s\cap H_h,w_s|_{H_h}) is determined by this initial state and the switch array inside H_h.

Every position at time a in a half is marginally uniform in that half after the intervening d-1 updates: each of its remaining coordinates is refreshed by a conditionally fair coin. Summing this endpoint probability over the free positions at time a gives, for every y\in H_{1-h},

For x\in H_h, the random variables w_s(x)^2 and \mathbf1_{\{x^{(i_s)}\in B_s\}} depend on opposite future arrays. Their conditional independence therefore yields

The values of the flow and the free set within one half need not be independent.

Now take a uniform random initial list, independent of the shuffle, and suppose m\ge n/8. For every fixed realization of the shuffle, the image of the uniformly chosen free set is uniform among the m-subsets of V. Thus B_a is unconditionally a uniform m-subset. Its count in either half has mean m/2, and

For completeness we give the bounded-difference exponential-moment argument of Hoeffding and Azuma [6, 1] in this setting. Expose the uniform subset in random order and use the Doob martingale for its final half-count. Changing the next draw can be coupled with the remaining completion by exchanging the two possible drawn points; the increment has absolute value at most one. A centered random variable in an interval of length two has exponential moment at most e^{\theta^2/2}: for \psi(\theta)=\log\mathbb E e^{\theta Z}, the second derivative is a tilted variance at most one, and integration from zero proves the bound. Conditional multiplication over the m martingale increments gives \mathbb E e^{\theta(M_m-M_0)}\le e^{m\theta^2/2}. Markov’s inequality and minimization in \theta give the lower-tail bound e^{-u^2/(2m)}; use u=m/4 to obtain (4.2).

With

both half densities in (4.1) are at least \beta, except on an \mathcal F_a-measurable event of probability at most \eta. The exceptional event may be correlated with the energy. Since the energy is at most one, (4.1), summed over the halves, nonetheless gives

This uses an unconditional exceptional probability and does not claim concentration conditional on the exposed paths.

Let \gamma=31/32=1-\beta/2. Equations (3.2) and (4.3) imply

For t\le2d the bound follows from E_t\le1. For t>2d, use induction in (3.2); all indices t-j are at least d. The geometric contribution is at most \gamma^{t-2d}, because

The constant contribution is at most (1-\beta)\eta/\beta+\eta=\eta/\beta. This proves (4.4).

We next turn the energy estimate into a bound on the positions of an ordered list. We use the sequential comparison of Morris [11, Section 5.3, Lemma 12], also given by Morris, Rogaway and Stegers [14, Appendix A, Lemma 2]. It applies to any law on distinct-position tuples; U_{\mathrm{inj}} denotes the uniform law on such tuples.

The bound averages over the input list. We will use that average to choose a single partition into eight blocks for which all eight complementary marginals have small total error.

The remaining question concerns the full permutation, including dependencies involving the omitted cards. We now pass to representations of S_n. Restricting a representation to the subgroups that permute labels within the eight blocks separates it into tensor factors. At least one factor is small enough that the span generated by a single vector remains controlled, even when that factor occurs with multiplicity. This converts the marginal estimates into bounds on the Fourier matrices of a slightly trimmed law. First we prove the dimension estimate needed to sum those bounds; then we carry out the trimming and orbit-span argument.

5 The representation degrees needed for the lift

The transfer needs two properties of symmetric-group representation degrees. A uniform bound on their reciprocal sum will control orbit dimensions for the label-block subgroups in Proposition 6.1. Decay of the sum away from the trivial and sign representations will control the full-group Fourier error in Section 7. We prove both properties here.

Let \lambda\vdash s be a partition of s, represented by its Young diagram, whose row lengths are the parts of \lambda. Write f^\lambda for the degree of the corresponding irreducible representation of S_s. A standard Young tableau fills the diagram with 1,\ldots,s, increasing along each row and column. We use the classical facts that f^\lambda counts these tableaux, is unchanged by transposing the diagram, and satisfies the hook-length formula

where h(z) counts the box z together with the boxes to its right and below it. The row (s) and column (1^s) give the trivial and sign representations, respectively [18, 4].

Stronger asymptotics for every fixed positive inverse-degree exponent were proved by Liebeck and Shalev [8, Theorem 2.6]. The elementary proof below gives all the summability used here.

6 From eight large marginals to the whole group

Assume d\ge3, put G=S_n, and let \mu=q_d^{*(200d)}. Choose a uniform ordered partition (M_1,\ldots,M_8) of the labels into blocks of size n/8. Each complement is a uniform subset of size 7n/8. Giving it an independent uniform order yields exactly the list in (4.6), and changing that order does not change the total-variation distance. Therefore some deterministic partition has

Fix that partition. Let H_a be the subgroup of all permutations that fix every label outside M_a. Two maps g,g' agree on the complement exactly when g'=gh for some h\in H_a. Thus the complementary image tuple identifies the left coset gH_a, and its uniform law is the coset pushforward of uniform measure on G.

For each a, remove every coset C\in G/H_a with \mu(C)>2|H_a|/|G|, and let E be what remains after all eight removals. A bad coset satisfies \mu(C)\le2(\mu(C)-|H_a|/|G|). The sum of positive excesses is its marginal total-variation distance, so the union bound gives

For all sufficiently large n, we have M\ge1/2 and may define the conditional probability law \nu=\mu(\cdot\mid E). Then

A coset removed for that subgroup has no retained mass; every other coset had the required original cap. This proves the same cap for all eight systems under one common law. Consequently, for every nonnegative function that is constant on the cosets of any one H_a, its expectation under \nu is at most four times its uniform expectation. The next proposition applies this comparison to projections onto orbit spans.

We use finite-dimensional complex Hilbert spaces. The operator norm is \|B\|_{\mathrm{op}}, and the Hilbert–Schmidt norm is \|B\|_{\mathrm{HS}}^2=\operatorname{tr}(B^*B) with ordinary, unnormalized trace. Write \widehat H_a for the set of equivalence classes of irreducible representations of H_a. For a probability measure \nu on G and a unitary representation \rho, define \widehat\nu(\rho)=\sum_g\nu(g)\rho(g).

7 Eight convolution factors and the final comparison

The final comparison uses the finite-group Fourier method of Diaconis and Shahshahani [3, Section 3]. The shuffle law is not constant on conjugacy classes, so we retain matrix norms throughout. By (6.2) and Proposition 6.1, the nearby law \nu has the required Fourier operator bound. We first show that \nu^{*8} is close to uniform, then compare it with the law of eight independent intervals of 200d shuffles, \mu^{*8}=q_d^{*(1600d)}.