Reservoir sampling

Earthquakes of magnitude 5 or more happen about five times a day, and the USGS catalogue for 2015–2025 holds 19,073 of them. Imagine that they come to you as a stream, one at a time, and you must keep a fair random sample of ten: at any moment every event that has arrived so far must be in it with the same probability. How many will arrive in all is unknown, and the whole stream cannot be stored. Press “Play”: the events land on the map, and a reservoir of k of them is updated by one simple rule.

The map draws itself: earthquakes fall along the boundaries of the tectonic plates. Each next event’s chance of entering the reservoir falls (at first every one, then one in hundreds, by the end one in two thousand), but the late events are not short-changed: there are simply many more of them. The “2,000 runs” button checks this directly: the years in two thousand samples are distributed just as in the whole catalogue, including 2021, when there were more earthquakes than usual.

The idea

Put the first k elements into the reservoir. Take element number i > k with probability k/i: choose a random j from 1 to i, and if j ≤ k, put the element in place j, displacing whatever was there. Otherwise skip the element. At any moment the reservoir is a uniform random sample of everything that has arrived.

This is Algorithm R. The step is the same as in the Fisher–Yates shuffle, only the shuffle is never written out in full: if you imagine each new element taking a random place in a permutation of everything that has arrived, the reservoir is its first k places.

Why every element is equally likely

What is more, not only every element but every set of k is equally likely: all C(n, k) samples have the same probability. The listing below checks this on a small example: with n = 6 and k = 2 all 360 equally likely chains of choices are enumerated in full, and each element ends up in the reservoir in exactly 120 of them.

Algorithm L: jumping straight to the next replacement

Algorithm R draws a random number for every element, although nearly all the late ones are skipped. How many times does the reservoir change at all? Element number i replaces something with probability k/i, and the expected number of replacements is the sum of these probabilities:

$$\sum_{i=k+1}^{n} \frac{k}{i} = k\,(H_n - H_k) \approx k \ln \frac{n}{k}.$$

On the earthquake stream with k = 10 that is 75 replacements for 19,073 events. In 1994 Kim-Hung Li proposed Algorithm L, which does not look at skipped elements at all: it computes right away how many elements to skip before the next replacement. The key observation is a threshold. If each element is given a random number from (0, 1) and the reservoir keeps the k elements with the smallest numbers, the sample is also uniform. Let W be the largest of these k numbers. The next element enters the reservoir only if its number is less than W, so each next element is a trial with success probability W. The number of failures before the first success is geometrically distributed and can be drawn with a single random number: ⌊ln u / ln(1 − W)⌋. After a replacement the threshold is recomputed: the new W is the old one multiplied by u^1/k.

The marks are the elements that entered the reservoir. On a logarithmic axis they fall evenly, and that is no coincidence: the stretch from a to b gets k·ln(b/a) replacements on average, the same for each tenfold stretch. Algorithm L goes through a billion elements in the browser in a fraction of a millisecond and draws a few hundred random numbers instead of a billion. If the stream can be fast-forwarded, for example a file of fixed-length records, the skipped elements need not even be read.

Sampling with weights

Sometimes elements are not equal: frequent search queries matter more than rare ones, strong earthquakes more than weak ones. In 2006 Pavlos Efraimidis and Paul Spirakis proposed a simple key for this: give an element with weight w and a random u from (0, 1) the key u^1/w, and keep the k elements with the largest keys in a heap, as in any “k best from a stream” problem. The first place then goes to element number i with probability exactly w_i/Σw, and likewise further on among those remaining. The natural-looking key “u·w” does not give this: it favours heavy elements, but not in proportion to weight. The numbers u^1/w crowd near one for large weights, and they are compared by their logarithms; exactly how is in the second problem, on the energies of real earthquakes.

Code

The stream is the same earthquake catalogue from the sandbox. Algorithms R and L with k = 5 on one generator, and at the end an exact check of equal likelihood by enumerating all chains.

Algorithm R spent 19,068 random numbers, one per event after the first five, and made 44 replacements against an expected 40.7. Algorithm L made 45 replacements and spent only 137 random numbers: three per replacement (the jump length, the place in the reservoir and the new threshold) plus two more, for the first threshold and the last jump, which went past the end of the stream. Their samples differ, but both are fair.

What it costs

Algorithm R is O(n) time and random numbers and O(k) memory: the reservoir and a counter. Algorithm L spends O(k(1 + log(n/k))) random numbers and as much work if the stream allows elements to be skipped without reading them; if not, reading stays O(n), but calls to the generator become thousands of times fewer. Weighted sampling is O(n log k) because of the heap.

The number of replacements in the reservoir, k·(H_n − H_k), deserves a separate look: if the reservoir lives in shared memory or on disk, every replacement is a write. Over a day with a measurement every millisecond, a reservoir of 1028 elements is rewritten about 11.7 thousand times, once per seven thousand-odd measurements. This is the first problem at the end of the article.

Where it works

The command line. GNU shuf -n k since coreutils 8.22 chooses k random lines with a reservoir if the input is large: the file is not read into memory in full. The 2013 proposal on the coreutils mailing list explained the benefit just so: since it is known how many lines are needed on output, there is no point keeping the whole file.

Metrics. The Dropwizard Metrics library keeps not all response-time measurements but a sample. UniformReservoir is Algorithm R on 1028 elements by default: its quantiles describe the program’s whole lifetime. ExponentiallyDecayingReservoir is a sample in which fresh measurements carry more weight; it is built on the “forward decay” method of Cormode and co-authors and roughly reflects the last five minutes.

Distributed samples. If a stream passes through several servers, each keeps its own reservoir and knows how many elements went through it. The overall sample is assembled by merging: how many elements to take from each reservoir is decided at random with probabilities proportional to the stream lengths. The same trick works for sampling from several log files.

History

The problem of sampling from a stream of unknown length was solved as early as 1962 by Fan, Muller and Rezucha. Donald Knuth gave Algorithm R in its present form in the second volume of The Art of Computer Programming, crediting Alan Waterman; McLeod and Bellhouse found it independently in 1983. Jeffrey Vitter, in the 1985 paper “Random sampling with a reservoir”, gave the method its name and worked out how to skip elements: his Algorithm Z spends O(k(1 + log(n/k))) random numbers. Kim-Hung Li obtained the same bound in 1994 with the much simpler Algorithm L, and Efraimidis and Spirakis the weighted variant in 2006.

Problems

In the first problem you must compute how many times the reservoir is rewritten on average for streams of up to 10¹⁸ elements: a table of harmonic numbers and an asymptotic formula. The second is sampling by energy from a real earthquake catalogue, where the keys cannot be compared directly.