Fisher-Yates Shuffle, in depth: the proof, three classic biased shuffles, PRNG limits, Sattolo and inside-out variants, and shuffling training data

The algorithm

Think of the array as two regions. The suffix on the right holds elements whose final position is decided. The prefix on the left holds everything still waiting. At each step you pick one element uniformly from the prefix, including the last prefix slot, and swap it into the last prefix slot, which then joins the finished suffix. When the prefix has one element left, that element is the only possible choice and you stop.

import secrets

def fisher_yates(a, randbelow=secrets.randbelow):
    """Shuffle list a in place. randbelow(k) must return a uniform int in [0, k)."""
    for i in range(len(a) - 1, 0, -1):
        j = randbelow(i + 1)        # 0 <= j <= i  -- note the +1
        a[i], a[j] = a[j], a[i]
    return a

Two details carry all the weight. The loop runs from the end down to index 1, because position 0 is decided once everything else is. And the random index is drawn from [0, i] inclusive. Leaving j == i possible means an element is allowed to stay where it is; forbidding it produces a different algorithm, covered below.

Passing the generator in makes the dependency explicit: a seeded one for reproducible experiments, an operating-system source for anything adversarial.

Why it is exactly uniform

The proof is a counting argument. The first iteration has n choices, the second n-1, down to 2 choices at the last iteration, so the algorithm has n x (n-1) x ... x 2 = n! equally likely execution paths. Each path yields a permutation, and no two paths yield the same one: given the final arrangement you can recover which element was placed at position n-1, then n-2, and so on, which uniquely determines every choice. A bijection between n! equally likely paths and n! permutations means every permutation has probability exactly 1/n!.

An equivalent inductive view is useful when you design variants. After the step for index i, the element at position i is uniformly distributed over all elements that were in the prefix, and independent of the order already fixed in the suffix. Any change that breaks either property, uniform over the remaining candidates or independent of earlier steps, breaks the algorithm.

Three ways to get it wrong

Almost every shuffling bug in the wild is one of three patterns, and all of them pass a casual eyeball test.

The naive swap. Loop over every index and swap with a position drawn from the whole array: j = randbelow(n) instead of randbelow(i + 1). That makes nⁿ equally likely paths. For n = 3 that is 27 paths spread over 6 permutations, and 27 is not divisible by 6, so the distribution cannot be uniform. Enumerating it gives three permutations with probability 5/27 and three with 4/27. For any n of 3 or more, nⁿ is not a multiple of n!, so the bias never disappears; it just becomes harder to see.

Sort with a random comparator. Code such as items.sort(key=cmp_to_key(lambda a, b: random.choice([-1, 1]))) or the JavaScript idiom arr.sort(() => Math.random() - 0.5) violates the comparator contract, which requires consistency and transitivity. The output then depends on the sort implementation: insertion-sort-style algorithms leave elements near their starting position, merge sorts favour other patterns. The 2010 European browser-choice screen was a well-known public example of this skew. Sorting by a random key drawn once per element is different and correct when keys never tie, but it costs O(n log n) for no benefit in memory.

Modulo bias. Turning a raw 32-bit random integer into an index with r % (i + 1) favours small indices whenever 2³² is not a multiple of i + 1. The bias is tiny for small ranges and grows with the range; with a 16-bit or 8-bit source it becomes obvious. The fix is rejection sampling: discard draws from the incomplete final block.

def randbelow_from_u32(k, next_u32):
    """Uniform int in [0, k) from a source of uniform 32-bit ints, without modulo bias."""
    limit = (1 << 32) - ((1 << 32) % k)    # largest multiple of k that fits
    while True:
        r = next_u32()
        if r < limit:                      # reject the short final block
            return r % k

The expected number of draws is below 2 for any k, and close to 1 for small k. Daniel Lemire's multiply-and-shift method achieves the same result with fewer divisions and is what several modern standard libraries use; either is fine as long as you do not use a bare modulo.

The random source sets the ceiling

A perfect algorithm cannot produce more distinct outputs than its random source has distinct states. A 52-card deck has 52! orderings, and log₂(52!) is about 225.6, so a generator needs at least 226 bits of state, and a seed with that much entropy, to be able to reach every deal. A generator seeded from a 32-bit value can produce at most 2³², about 4.3 billion, different decks no matter how long its period is. An attacker who sees a few cards can search that seed space.

Python's documentation for random.shuffle makes the general point: even for small lists the number of permutations quickly exceeds the period of most generators, and 2080 is the longest sequence whose permutations fit within the Mersenne Twister's period. For statistics and ML this rarely matters, because you need a representative sample of orders, not reachability of all of them. For games, lotteries, security tokens and anything with an adversary, use a cryptographic source such as secrets in Python, crypto.getRandomValues in browsers or java.security.SecureRandom.

Variants you will need

Partial shuffle for sampling. If you need k random elements without replacement, stop after k iterations. Running the loop from the front makes this natural: at step i swap a[i] with a uniform pick from a[i..n-1], and after k steps the first k slots are a uniform random k-sample in random order. Cost is O(k), not O(n). When n is unknown or the data is a stream, use reservoir sampling instead; it is the streaming relative of the same idea.

Inside-out shuffle. When you are copying from a source into a new array anyway, you can shuffle while copying: for each i, pick j in [0, i], move the element at j to position i, and put the new source element at j. It needs no initial copy and works when the source is an iterator of known or unknown length.

Sattolo's algorithm. Change the draw to randbelow(i), so j is strictly less than i. Now no element may stay in place at its step, and the result is a uniformly random single cycle of length n, one of (n-1)! cyclic permutations. That suits a pointer-chasing memory benchmark that must visit every node before repeating, and is the wrong answer for a plain shuffle.

def sample_k(a, k, randbelow):
    """First k slots become a uniform random k-subset in random order. O(k)."""
    n = len(a)
    for i in range(k):
        j = i + randbelow(n - i)
        a[i], a[j] = a[j], a[i]
    return a[:k]

def sattolo(a, randbelow):
    """Uniform random cyclic permutation: following a[i] visits every index once."""
    for i in range(len(a) - 1, 0, -1):
        j = randbelow(i)            # strictly less than i
        a[i], a[j] = a[j], a[i]
    return a

Testing a shuffle

Shuffles fail silently, so test them statistically. For a small n you can enumerate all n! permutations, run the shuffle many times, and apply a chi-square goodness-of-fit test against the uniform distribution. With n = 4 there are 24 permutations and 23 degrees of freedom; the critical value at significance 0.001 is about 49.7. A correct shuffle should exceed it roughly once in a thousand runs, while the naive swap blows through it with a few hundred thousand trials.

import itertools, random
from collections import Counter
from math import factorial

def chi_square(shuffle, n=4, trials=240_000, seed=1):
    rng = random.Random(seed)
    counts = Counter()
    for _ in range(trials):
        a = list(range(n))
        shuffle(a, rng.randrange)
        counts[tuple(a)] += 1
    expected = trials / factorial(n)
    return sum((counts.get(p, 0) - expected) ** 2 / expected
               for p in itertools.permutations(range(n)))

def naive(a, randbelow):
    for i in range(len(a)):
        j = randbelow(len(a))
        a[i], a[j] = a[j], a[i]

print("fisher-yates", round(chi_square(fisher_yates), 1))   # typically 10-40
print("naive       ", round(chi_square(naive), 1))          # far above 49.7

For larger n, use a position-frequency matrix instead: count how often element e lands at position p; every cell should be near trials/n. The broader toolkit for reasoning about such tests, including why a single failed run is not proof of a bug, is covered in probability in algorithms and expected value in algorithms.

Shuffling training data

In training pipelines the shuffle is usually hidden inside a library, which moves the risks rather than removing them.

PyTorch. DataLoader(shuffle=True) uses a RandomSampler that draws a fresh permutation of indices each epoch from the loader's generator. In distributed training, DistributedSampler shuffles with a seed derived from the epoch, so you must call sampler.set_epoch(epoch) at the start of each epoch. Forget it and every epoch sees the same order, a bug that shows up only as slightly worse convergence.

Buffer shuffles. Streaming APIs such as tf.data.Dataset.shuffle(buffer_size) and many WebDataset-style pipelines keep a buffer of B examples and emit a random one as each new example arrives. That is not a uniform permutation unless B is at least the dataset size: an example near the start of the stream can only travel so far. If your shards are sorted by source, date or length, a small buffer leaves batches correlated. The standard mitigation is two-level: shuffle the shard list each epoch, interleave several shards, then apply the buffer.

Data too big for memory. A uniform external shuffle is possible in two passes. Pass one writes each record to one of m bucket files chosen uniformly at random. Pass two loads each bucket, Fisher-Yates shuffles it in memory, and concatenates the buckets. Because the bucket assignment is uniform and each bucket is shuffled uniformly, the concatenated result is a uniform permutation. Pick m so each bucket fits in RAM. How much mixing a training run actually needs, and when a deliberate curriculum beats a uniform order, is covered in dataset ordering and shuffling.

# Two-pass external shuffle: uniform even when the data does not fit in memory.
import random, pickle

def external_shuffle(records, m, rng, open_bucket):
    buckets = [open_bucket(b, "wb") for b in range(m)]
    for r in records:                                # pass 1: random scatter
        pickle.dump(r, buckets[rng.randrange(m)])
    for f in buckets:
        f.close()
    for b in range(m):                               # pass 2: shuffle each bucket
        items = []
        with open_bucket(b, "rb") as f:
            while True:
                try:
                    items.append(pickle.load(f))
                except EOFError:
                    break
        fisher_yates(items, rng.randrange)
        yield from items

Reproducibility. Log the seed and the generator algorithm with every run, and derive per-epoch and per-worker seeds deterministically from the base seed. Resuming from a checkpoint should restore the sampler state, not just the model, or the resumed run silently replays or skips data.

Trade-offs

Sort-by-random-key earns its place in distributed engines, where a 64-bit or 128-bit random key and an existing sort give a uniform shuffle with negligible tie probability. On one machine Fisher-Yates wins. For other structures that lean on randomness for balance, see the treap; for generating permutations in order rather than at random, see permutations by backtracking.

What to do next

1. Search your codebase for % n-style index generation, sort( with a random comparator, and hand-written shuffle loops; check each uses an inclusive [0, i] draw.
2. Replace hand-written loops with the standard library shuffle where possible, and pass an explicit generator object rather than relying on global state.
3. Add the chi-square test above to CI for any custom shuffle or sampler, with a fixed seed and a threshold at a strict significance level.
4. Use a cryptographic source for any shuffle an adversary could profit from predicting.
5. In distributed training, confirm set_epoch is called and that resumed runs restore sampler state; log seeds with every run.
6. If you use a streaming buffer shuffle, measure how far examples move, and add shard-level shuffling before the buffer if your shards are sorted.