Skip to content

algorithms.md

Latest commit

 

History

History
294 lines (234 loc) · 17.6 KB

algorithms.md

File metadata and controls

294 lines (234 loc) · 17.6 KB
Raw

Algorithms for sampling without replacement

In the standard form of the problem, the answers can be in any order and so the order must also be chosen fairly; there is a variation in which the answers must be in sorted order. Any algorithm can be turned into one with a "sorted" order guarantee with an O(k log k) sort, or a "random" order guarantee with an O(k) Fisher-Yates shuffle.

Algorithm Python C++ Rust Order guarantee Data structures Time
Quadratic rejection sampling Python C++ Rust Random none k^2
Set-based rejection sampling Python C++ Rust Random Set k
Floyd's F2 Python C++ Rust none Set k
Quadratic F2 Python C++ Rust Random none k^2
Iterative choosing Python C++ Rust Sorted none n
Reservoir sampling Python C++ Rust Random none n
SELECT Python C++ Rust Random n-sized list n
HSEL Python C++ Rust Random Dictionary k
RANKSB Python C++ Rust Sorted none k
Vitter's Algorithm D Python C++ Rust Sorted none k
Hidden Shuffle Python C++ Rust Sorted none k
Cardchoose Python C++ Rust Sorted none k log k

Quadratic rejection sampling

Algorithm Python C++ Rust Order guarantee Data structures Time
Quadratic rejection sampling Python C++ Rust Random none k^2

For each output, try choosing an integer in [0, n), but check every integer already in the output to see if it's the same, and reject and re-draw if so. This is quadratic in k, but for small k the small constants can mean this algorithm does very well. It takes O(k^2) expected time for k ≤ n/2, and more as k approaches n and more draws are rejected, as for set-based rejection sampling.

The check is the inner loop, so it matters that the compiler can vectorise it, comparing several entries in one instruction. My C++ checks every entry rather than stopping at a match, since a loop that can stop early doesn't vectorise. My Rust uses the standard library's contains, which gets both: it compares entries 16 at a time with vector instructions, and stops after the first block of 16 that has a match.

Set-based rejection sampling

Algorithm Python C++ Rust Order guarantee Data structures Time
Set-based rejection sampling Python C++ Rust Random Set k

Maintain a set of already-produced values. While the output isn't full, try choosing a random integer in [0, n). If this is in the set of already-produced values, discard it. Otherwise, append it to the output, and add it to the list of already-produced values.

One of the simplest possible algorithms, but the cost of set lookups and insertions is high. It takes O(k) expected time for k ≤ n/2, where each draw is rejected with probability at most a half. As k approaches n, rejections mount: the expected number of draws is n(Hn − Hn_−_k), where H is the harmonic number, so at k = n it takes O(n log n) time.

Floyd's F2

Algorithm Python C++ Rust Order guarantee Data structures Time
Floyd's F2 Python C++ Rust none Set k

I refer you to two excellent explanations:

F2 returns a set, with no order. For random order, the C++ and Rust implementations do not shuffle afterwards: as in Daniel Ting's random-order version, when the chosen value is already present, they move a randomly chosen earlier output to the end and put the new value in its place, which costs an extra draw only on collisions.

The original presentation is as one of Jon Bentley's Programming Pearls, though the algorithm is due to Robert W Floyd: A Sample of Brilliance.

Quadratic F2

Algorithm Python C++ Rust Order guarantee Data structures Time
Quadratic F2 Python C++ Rust Random none k^2

A quadratic variant on Floyd's F2 which needs no external data structure: the output array is the set, scanned to check each draw. When the draw r is already present, it replaces that entry with m and appends r, which keeps the output in uniformly random order without a shuffle. The Rust rand crate used Floyd's F2 with a scan in the same way, but kept the order random either by inserting m before the matching entry, which shifts every entry after it, or by shuffling the whole output afterwards; it has used this version since rand pull request 1277, in 2023. Unlike Ting's random-order Floyd's F2 above, it needs no extra draw on a collision, because the scan has already found where r is. That pull request is the first publication of this way of keeping the order random that I know of. Like quadratic rejection sampling, the small constants mean this does very well for small k.

For random order, each step draws t and replaces the existing t, if there is one, with j. At most one entry matches, so my implementations write every entry, keeping its value unless it matches, rather than stopping at the match; a loop that can stop early doesn't vectorise. For sorted order, only the set matters until the final sort, so where t is already present, a step adds j instead of moving the existing t to j. That gives the same set, and the step only has to check whether t is present, as in quadratic rejection sampling, which is faster than rewriting every entry.

Iterative choosing (Algorithm S)

Algorithm Python C++ Rust Order guarantee Data structures Time
Iterative choosing Python C++ Rust Sorted none n

Iterate through every possible value in [0, n) in order and decide whether to include it in the output, weighting the probability appropriately depending on the number of spaces left in the output and the number of values remaining to test. One of the faster algorithms for producing sorted output where k is a large fraction of n.

  • Donald Knuth
    • The Art of Computer Programming
    • Volume 2: Seminumerical Algorithms.
    • Section 3.4.2
    • Page 137
    • "Algorithm S"

Reservoir sampling (Algorithm R)

Algorithm Python C++ Rust Order guarantee Data structures Time
Reservoir sampling Python C++ Rust Random none n

Initialize the output with the integers [0, k) in random order. Then iterate through the remaining values [k, n) in order and decide whether to replace an existing, randomly chosen output with that integer. Fast for producing random output where k is a large fraction of n.

  • Donald Knuth
    • The Art of Computer Programming
    • Volume 2: Seminumerical Algorithms.
    • Section 3.4.2
    • Page 138
    • "Algorithm R"

SELECT

Algorithm Python C++ Rust Order guarantee Data structures Time
SELECT Python C++ Rust Random n-sized list n

Initialize an array of the values [0, n). Pick out values at random to include in the output, and delete them by replacing them with a value from one end and shortening the array.

This is SELECT in Goodman and Hedetniemi's Introduction to the Design and Analysis of Algorithms, though in a slight variation I pick off the array at the start, rather than the end.

Each call allocates the array afresh, so its cost depends on the allocator. The C++ and Rust implementations align the array to a 64-byte cache line: malloc only promises 16 bytes, and on x86-64, filling the array took 1.5 times as long in Rust when it started 16 bytes past a 32-byte boundary, which depended on where earlier allocations happened to leave it. On Linux, glibc's malloc serves blocks over 32 MB with a fresh mmap and unmaps them on free, so from n ≈ 8 million every call faults in every page of the array again; this makes filling it about ten times slower than on macOS, where nothing like this shows until the largest n timed.

HSEL

Algorithm Python C++ Rust Order guarantee Data structures Time
HSEL Python C++ Rust Random Dictionary k

This is a variation of the above which uses a hash table to simulate the array of values [0, n), storing only the entries that change. First proposed by Ernvall and Nevalainen in 1982:

Teuhola and Nevalainen give a version that needs 2_k_ storage locations rather than about 5_k_, which draws all k random numbers first and then works out which value each one stands for, from the last back:

It is also presented in Martin Harriman, Sampling Without Replacement: Durstenfeld-Fisher-Yates Permutation for Very Large Permutation Sizes.

My C++ and Rust implementations use a faster variant, which keeps the first k entries of the simulated array in the output array itself, so only swaps that reach beyond k touch the hash table. This variant was described by rici in a Stack Overflow answer in 2015. Thanks to Mathias Weller for bringing it to my attention.

RANKSB

Algorithm Python C++ Rust Order guarantee Data structures Time
RANKSB Python C++ Rust Sorted none k

Divide [0, n) into k bins of nearly equal size. First decide how many of the chosen values fall in each bin, by simulating set-based rejection sampling while keeping only a count for each bin: a draw is rejected with the probability that it would duplicate a value already accepted in its bin, and that depends only on how many values the bin has accepted. Then choose that many values within each bin; bins hold about one value each, so a quadratic method is fine. The output comes out sorted, it takes O(k) expected time for k ≤ n/2, and the book's program needs no storage beyond the output array. As k approaches n it slows down, and for k ≥ n/2 the book pairs it with iterative choosing (which Nijenhuis and Wilf refer to as RKS2); I benchmark RANKSB on its own.

My C++ and Rust implementations follow the program in the book, except that they fill each bin as soon as they know where its values go, rather than first marking out every bin's slots and then scanning for the marks; the Python one keeps the counts in a separate array, which is easier to follow.

  • Albert Nijenhuis and Herbert S. Wilf
    • Combinatorial Algorithms for Computers and Calculators, 2nd edition
    • Academic Press, 1978
    • Chapter 4, "Random k-Subset of an n-Set (RANKSB)"

Vitter's Algorithm D

Algorithm Python C++ Rust Order guarantee Data structures Time
Vitter's Algorithm D Python C++ Rust Sorted none k

A sequential sampling algorithm: like iterative choosing, it produces the output in sorted order, but instead of deciding about each value in [0, n) in turn, it draws the number of values to skip before the next chosen one directly. The skip has a known discrete distribution, and Algorithm D samples it by rejection from a continuous approximation, the smallest of k uniform reals, drawn in one step as in Bentley and Saxe's method for generating sorted random numbers. This takes O(1) expected time per output, so the whole algorithm takes O(k) expected time and needs no extra memory. The cost is floating-point arithmetic: each output needs a few uniform variates and exponentiations. Once n is at most 13 k, it falls back on Vitter's Algorithm A, which draws each skip in O(skip) time, so the total work for the rest of the output is O(n).

My C++ and Rust implementations follow the Pascal-like code in the paper's appendix; the Python one states the method directly, which is easier to follow.

  • Jon Louis Bentley and James B. Saxe

  • Jeffrey Scott Vitter

Hidden Shuffle

Algorithm Python C++ Rust Order guarantee Data structures Time
Hidden Shuffle Python C++ Rust Sorted none k

Picture a Fisher-Yates shuffle of [0, n) that stops after k steps, so the sample is whatever ends up in the first k positions. Call those the low positions and the rest the high positions. Each step reaches a high position with a known probability, independently of the others, and which high position it reaches is uniform and independent of everything else. So first count how many steps reach a high position. Then draw that many high positions with replacement, in sorted order as in Bentley and Saxe's method. The first time a high position is drawn, its value joins the sample. A repeat finds a low value that was swapped there earlier, so the sample gains a low value instead, much as a collision in Floyd's F2 adds a different value. Given how many values are low, they are a uniformly random subset of the low positions by symmetry, and Vitter's Algorithm A chooses them in order. The output comes out sorted, it takes O(k) time in the worst case, and it needs no extra memory. Like Algorithm D, it relies on floating-point arithmetic, including an exponentiation for each high draw.

The code printed in the paper is wrong: it can divide by zero, and the last step has an off-by-one error that biases the output. My implementations apply the corrections from the errata in the authors' repository. Thanks to metasj for bringing this algorithm to my attention.

Cardchoose

Algorithm Python C++ Rust Order guarantee Data structures Time
Cardchoose Python C++ Rust Sorted none k log k

This is my algorithm: a Pólya urn generates a uniform multiset, and stars and bars turns it into a uniform subset. Detailed description, including how it relates to prior art.