Random Selection, in depth: k of n without replacement with sparse Fisher-Yates, Floyd's algorithm, Algorithm S, unbiased integers and hash-based samples

The problem, precisely

There are C(n, k) subsets, and a correct sampler gives each probability 1 / C(n, k). Equivalently, every item is included with probability k/n and the inclusions are jointly right. Some samplers meet the first condition and fail the second. Which method fits depends on three questions: Can you hold or index the n items? Is k small next to n? Does the output need to be sorted, streamed in order, or reproducible?

The sections below take these in turn. One caution applies to all of them: the guarantee is only as good as the integers they draw, which is why bounded integers get their own section.

Partial Fisher-Yates, and a sparse version for huge n

If the items sit in an array you may permute, run Fisher-Yates for only k steps. Swap position i with a uniform position in [i, n), and after k steps the prefix is a uniform random k-subset in random order. This costs O(k) time but needs a mutable array of n, and copying a read-only input costs O(n). The Fisher-Yates article proves the full shuffle. Stopping early is uniform for the same reason: each prefix of the swap sequence is a uniform partial permutation.

When n is huge, say 10¹² user IDs or every row of a table, you cannot allocate the array. You do not need to. A virtual array where position i holds i until something is swapped needs storage only for touched positions, at most k entries:

def sparse_partial_fy(n, k, rng):
    swaps = {}                                # position -> value, only where it differs
    out = []
    for i in range(k):
        j = rng.randrange(i, n)               # uniform in [i, n)
        vi = swaps.get(i, i); vj = swaps.get(j, j)
        swaps[j] = vi                         # position j receives what was at i
        out.append(vj)                        # position i now holds vj; it is never touched again
    return out

With n = 10¹², k = 5 and seed 1, this returned five IDs and stored at most five dictionary entries. The output order is random, which matters if you assign the first item to one bucket and the rest to another.

Floyd's algorithm

Robert Floyd's algorithm, which Jon Bentley published in his Programming Pearls column in 1987, does the same job with only a set of chosen values:

def floyd_sample(n, k, rng):
    chosen = set()
    for j in range(n - k, n):
        t = rng.randrange(j + 1)              # uniform in [0, j]
        chosen.add(j if t in chosen else t)   # collision: take j, which cannot be chosen yet
    return chosen

It makes exactly k random draws. No draw is ever rejected, so the cost does not climb as the set fills, unlike draw-until-new loops. The correctness proof is by induction. Suppose that after the step for j - 1 the set is a uniform random (m - 1)-subset of [0, j). The step for j draws t from j + 1 values. Any particular m-subset S of [0, j] that contains j arises when the previous set was S minus j and t was j or one of the m - 1 already-chosen values. Any S without j arises from m different previous sets, each with exactly one t. Both cases give the same count, so the result is uniform. The output is a set, so if you need random order, shuffle the k results, which costs O(k).

Sequential selection: Algorithm S and skips

Sometimes you must walk the items in order, for example a sorted file you can read only once or a table scan, and emit the sample in that order. Knuth's Algorithm S (TAOCP volume 2, section 3.4.2) selects item i with probability (still needed) / (still remaining):

def selection_sampling(n, k, rng):
    out, needed = [], k
    for i in range(n):
        if needed == 0: break
        if rng.random() * (n - i) < needed:   # P = needed / (n - i)
            out.append(i); needed -= 1
    return out                                # already sorted; exactly k items

It always returns exactly k items: once remaining equals needed, the probability reaches 1. The cost is O(n) draws. Vitter's Algorithm D (1984) brings that to O(k) expected time by drawing the gap to the next selected item from its exact distribution and skipping over the rest. Use it when n is large and reading an item is cheap to skip, such as seeking in a file of fixed-size records.

Do not confuse this with Bernoulli sampling, which keeps each item independently with probability p. Bernoulli sampling is simpler and parallelizes trivially, but the sample size is random: binomial with mean pn. Use it when you need a rate, not an exact count.

Uniform integers without modulo bias

Every sampler above calls "uniform integer in [0, m)". The classic bug is rand32() % m. When m does not divide 2³², the low remainders come up more often. With m = 3·2³⁰ the effect is dramatic. Bucketing 300,000 results into thirds of the range gave 149,723 / 75,419 / 74,858 with modulo, because values below 2³⁰ have two preimages. Lemire's multiply-shift method with rejection gave 99,966 / 100,463 / 99,571:

uint32_t bounded(uint32_t m) {                 /* Lemire, ACM TOMS 2019 */
    uint64_t prod = (uint64_t)next32() * m;
    uint32_t low = (uint32_t)prod;
    if (low < m) {                              /* rare slow path */
        uint32_t threshold = -m % m;            /* (2^32 - m) mod m */
        while (low < threshold) {
            prod = (uint64_t)next32() * m;
            low = (uint32_t)prod;
        }
    }
    return prod >> 32;
}

Real-world m values are usually tiny compared with 2³², so the bias is small. But it is systematic, and it compounds across k draws. Python's randrange and Java's nextInt(bound) already reject, so call them rather than reducing raw bits yourself. Also consider the generator's state size. A 64-bit seed cannot reach most subsets when C(n, k) is astronomically large, which matters for lotteries and cryptographic uses, not for a validation split. Use a CSPRNG such as secrets when someone could profit from predicting the sample.

Worked example: testing uniformity with chi-square

Uniformity is testable on a small case. With n = 6 and k = 3 there are 20 subsets. Draw 200,000 samples, count each subset (sorted), and compute the chi-square statistic against 10,000 expected per subset. With 19 degrees of freedom, a correct sampler stays below about 36.2 on 99 percent of runs. Each sampler was seeded separately:

The last row is a common "fix" for collisions: on a duplicate, step to the next unused value. Every item still has a fair chance of some inclusion, but values right after taken ones are favoured. The subset {2, 3, 4} came up 15,005 times and the rarest subset only 5,418. Spot checks of marginal frequencies can miss this, and the subset-level test catches it at once. Keep a test like this in CI for any sampler you write. The table above came from a single run, and on any one run about 1 percent of correct samplers will exceed the threshold, so retest before you panic.

Reproducible and distributed selection by hashing

Random draws have two operational weaknesses. Re-running with new data reshuffles everything, and shards cannot agree on a sample without coordinating. Hashing fixes both. Keep a key when a salted hash falls below a threshold:

def in_sample(key, salt="val-2026", rate=0.1):
    h = hashlib.sha256(f"{salt}:{key}".encode()).digest()
    return int.from_bytes(h[:8], "big") < rate * 2**64

Over 100,000 user keys this kept 9,873, close to the expected 10,000. Each key's decision never changes, so a user stays in validation when new users arrive, and any shard decides locally. This is Bernoulli-style: the count varies. For an exact k, take the bottom-k: the k keys with the smallest hashes. Each shard keeps its own k smallest, and a merge keeps the k smallest overall, which is the same answer a single machine would give. Change the salt to draw an independent sample, and never reuse the salt that defines your test split for anything else.

Operational guidance

• Use the library when it fits. Python's random.sample(population, k) switches between a partial shuffle and a set of chosen indices depending on n and k. NumPy's Generator.choice(n, k, replace=False) is the vectorized equivalent. Pass a seeded generator, not global state.
• Sample IDs, then fetch. Selecting k row IDs and fetching them beats ORDER BY random() LIMIT k, which assigns a key to every row and sorts.
• Know your SQL sampler. PostgreSQL's TABLESAMPLE SYSTEM samples whole pages, so rows on a page come in clusters. BERNOULLI samples rows. Neither gives an exact count.
• Split on the right unit. Sample users, sessions or documents, not rows, when rows from one entity are correlated, or the validation score will leak.
• Record the seed or salt with every dataset you derive, so the sample can be rebuilt and audited.

Failure modes

• Modulo bias in a hand-rolled bounded integer.
• Draw-until-new loops that slow sharply as k approaches n. Floyd needs no retries.
• Collision repair by probing or re-hashing, which biases the subsets, as the table above shows.
• Sorting by a random key with ties, for example a 16-bit key on a large table, falls back to input order on equal keys.
• Using a set's iteration order as a random order: CPython sets of small integers iterate in nearly sorted order, so shuffle Floyd's output if order matters.
• Global PRNG state shared across threads or workers, so forked data loaders draw identical samples.

What to do next

1. List every place your system samples k of n, and note whether n is known and whether the sample must be reproducible.
2. Replace % m reductions and probe-on-collision code with library calls or the methods above.
3. Add the 20-subset chi-square test for any custom sampler to your test suite.
4. Move train/validation and experiment assignment to salted-hash selection keyed on the entity, and store the salt with the dataset.
5. For huge ID spaces, use Floyd or sparse partial Fisher-Yates instead of materializing the population.