One pass, k slots, and two randomness budgets

Take a uniform random sample of a thousand items from a stream that cannot be counted, rewound or held. It arrives one item at a time; there may be a million of them or a hundred; when it stops the sample must hold exactly a thousand, each of which was equally likely to be any item that went past.

Written that way it sounds like it needs the length. Reservoir sampling does it in one pass and $k$ slots, and the algorithm is three lines. What makes it worth an essay here is not the trick but two things measurement adds to it: an off-by-one that produces a perfectly ordinary-looking sample and is invisible to any single run, and a second version with identical output whose randomness cost differs from the first by a factor of 145.

Algorithm R #

Keep the first $k$ items as the reservoir. Then for each later item at index $i$, draw a uniform integer $j$ in $[0, i]$; if $j < k$, overwrite slot $j$ with the new item, and otherwise discard it. That is the whole algorithm: one draw and one conditional store per item, with no state beyond the $k$ slots and the running index.

The correctness argument is one induction. Suppose after processing $i$ items each of them is in the reservoir with probability $k/i$. Item $i$ (zero-indexed, so the $(i{+}1)$-th) enters with probability $k/(i+1)$, directly. An older item survives if either the new item was rejected — probability $1 - k/(i+1)$ — or it was accepted but displaced a different slot, which happens with probability $\frac{k}{i+1}\cdot\frac{k-1}{k}$. Add those and multiply by the $k/i$ it already had:

$$\frac{k}{i}\left(1 - \frac{k}{i+1} + \frac{k-1}{i+1}\right) = \frac{k}{i}\cdot\frac{i}{i+1} = \frac{k}{i+1}$$

which is the same probability the new item got. The invariant holds at every step, so it holds when the stream stops, whenever that is — and the algorithm never needed to know when.

That last clause is the whole point. The uniformity is maintained continuously, so the sample is valid at every prefix. Stop the stream anywhere and the reservoir holds a uniform sample of what has gone past.

What the off-by-one does #

Write uniform integer in [0, i) instead of [0, i] — a mistake that is one character in most languages — and the algorithm still runs, still fills $k$ slots, and still returns a sample that looks exactly like a sample.

The bug is that the arriving item can now displace an older one but can never be rejected outright: when $j$ is drawn below $i$ rather than below $i+1$, the acceptance probability is $k/i$ rather than $k/(i+1)$, so every item is slightly over-favoured relative to the ones before it, and the effect compounds down the stream.

The first position is chosen 9.71% of the time where it is owed 12.5%: it is under-represented by 22%. The last position is chosen 12.97%. Over a whole stream the early items are systematically thinned.

Three things about this failure are worth naming, because it is the archetype for a whole class.

No single run reveals it. A biased sample of four items out of thirty-two is four items out of thirty-two. There is nothing to look at.

No correctness test catches it. The output has the right size, the right type, contains only items from the stream, and contains no duplicates. Every property an assertion would naturally check is satisfied.

Only a frequency measurement over many runs catches it, and only if the number of runs is chosen so that the estimator’s noise is below the bias. At 40,000 runs the noise floor is 1.41% and the bias is 23.5%, which is a comfortable margin. At 400 runs the noise floor is 14% and the two are indistinguishable, so the test would pass on the broken sampler and prove nothing about the correct one.

That last point is the reason the site’s gate brackets its tolerance rather than picking a round number: the threshold has to sit above the measured noise of a correct sampler (3.34%) and below the measured departure of the broken one (23.54%), and 8% is in that gap by a wide margin on both sides. A threshold chosen without measuring both ends is a threshold that has never been shown to be able to fail.

Algorithm L, and what it costs differently #

Algorithm R draws a random number for every item in the stream. For a stream of sixty-five thousand items and a sample of sixteen, that is 65,520 draws to produce a sixteen-element answer.

Almost all of them are rejections. When $i$ is large the acceptance probability $k/(i+1)$ is tiny, so the algorithm’s inner loop is overwhelmingly “draw a number, discard it”. Algorithm L, due to Li, replaces that with drawing how far to skip before the next acceptance — the gap is geometrically distributed, so it can be sampled directly with one uniform and a logarithm.

The output distribution is unchanged. That is the claim, and it is the one worth checking rather than assuming, because a sampler that skips is a sampler with a second opportunity to be subtly wrong.

And then the difference, which needs the bit counter to see at all:

Algorithm R’s consumption is $\Theta(n\log n)$ — one draw per item, each costing about $\log_2 i$ bits — and fits that class with a spread of 1.14 across three orders of magnitude. Algorithm L’s is $\Theta(k\log(n/k))$, which for fixed $k$ is $\Theta(\log n)$, and fits at 1.52.

At the top of the table the ratio is 145. Two algorithms, identical output distribution, identical space, and one spends 1.36 million random bits where the other spends nine thousand.

Where Algorithm L stops winning #

The table above sweeps the stream length at a fixed sample size, and the ratio grows without limit. Sweeping the other way is more interesting, because it shows where the second algorithm is not worth having.

Algorithm L’s cost is $\Theta(k\log(n/k))$ draws, each costing 32 bits for a uniform plus a few for the slot choice. Algorithm R’s is $n$ draws of about $\log_2 n$ bits. So L wins when $k\log(n/k) \ll n$, which is when the sample is small relative to the stream, and loses when $k$ approaches $n$.

At $n = 32$ and $k = 4$ — the size the frequency figures use — Algorithm R spends 166 bits per run and Algorithm L spends 585. L is three and a half times worse at that size, because it pays a full 32-bit uniform per skip and there are barely any items to skip over. The crossover for these parameters is around $n = 400$.

That is worth stating plainly because it is the ordinary shape of an optimisation and it is easy to lose: Algorithm L is an improvement for large streams and a pessimisation for small ones, the boundary is computable, and neither algorithm is simply better. A library that always uses L is wrong for short streams in exactly the way a library that always uses R is wrong for long ones — and since the whole premise is that the stream length is unknown, choosing between them is genuinely awkward. The practical answer is to start with R and switch to L once the index passes a threshold, which costs nothing because the two maintain the same invariant on the same reservoir.

Both columns, from two closed forms #

The table has four rows and both of its algorithms are recoverable from it in closed form, which is worth doing because the forms say where the crossover between them lives and the table does not.

Divide bits by draws. Algorithm R spends 12.7, 15.3, 17.9 and 20.7 bits a draw across the four rows — that is $\log_2 n$ (10, 12, 14, 16) with about a 29% rejection overhead on top, which is the convention the section below makes explicit. Algorithm L spends 22.8, 22.7, 22.7 and 22.7. R’s cost per draw grows with the stream and L’s does not, so R is paying twice over: more draws, each wider.

The draw counts fit as cleanly. Algorithm R takes $n - k$ of them, exactly. Algorithm L’s expected number of replacements is $k\ln(n/k)$, and at three draws per skip — one for the exponential, one for the gap, one for the slot — that predicts 200, 266, 333 and 399 against the measured 185, 272, 353 and 413. Within eight per cent at every row, with no constant fitted beyond the three.

So the two bills are

$$\text{R} \approx 1.29\,n\log_2 n, \qquad \text{L} \approx 68\,k\ln(n/k)$$

and the second contains $k$ where the first does not. That is the fact the crossover turns on, and it is the one the sweep at fixed $k$ cannot show.

Raising the sample size raises L’s cost and leaves R’s alone. Algorithm R draws once per item whatever $k$ is; Algorithm L draws once per replacement, and there are more replacements when there are more slots to replace. So the stream length at which L becomes worth using is not a property of the stream — it moves with the sample, later for a large one and earlier for a small one.

That has a practical edge, because the switch a library would implement is a threshold on the index. The section above proposes exactly that — start with R, switch to L once the index passes a threshold — and the arithmetic says the threshold cannot be a constant. A sampler taking sixteen items and one taking sixteen thousand cross at very different points, and a library hard-coding one number is running the wrong algorithm for every sample size but one. The threshold has to be computed from $k$, which is the one parameter the caller does supply.

It also explains why the small-stream case runs the way it does. At $n = 32$ and $k = 4$ the ratio $n/k$ is eight, so $\ln(n/k)$ is barely two and L’s per-skip cost of about seventy bits is being paid against a stream where R spends twenty bits a draw over twenty-eight draws. L is paying a fixed toll per replacement in a regime where replacements are nearly every item.

Which is where a crossing moved to’s standing point in a third field: a crossover quoted as a number is a crossover with its other parameters held fixed and unstated, and the useful form is the condition rather than the value. Counting the coin flips is what makes either side of this one measurable at all.

Both samplers are checked again at twice the stream and twice the reservoir, because a uniformity measurement at one shape is a measurement of one shape.

Why the bits matter, when they do #

It would be easy to dismiss this as bookkeeping. Random bits are cheap; a pseudo-random generator produces them at gigabytes per second; a million of them is nothing.

That is true in the common case and it is worth being specific about the cases where it is not.

When the randomness is real. Entropy from an operating system’s pool, or from a hardware source, is rate-limited in a way a PRNG is not. An algorithm needing 1.36 million bits per pass over a stream is not going to get them from getrandom().

When the randomness must be secret. As the previous essay argues, a randomised algorithm’s guarantee holds only while its coins are unknown. Coins that come from a cryptographic generator cost more than a xorshift, and the count then translates directly into cycles.

When the algorithm is distributed. Shared randomness across machines has to be transmitted or derived from a common seed, and the volume is a network cost.

When the same stream is sampled many times. Drawing a hundred independent samples from one stream multiplies the cost by a hundred, and the difference between nine thousand bits and 1.36 million becomes the difference between a megabyte of entropy and a hundred and seventy.

And when the draws themselves are the expensive part. Algorithm R’s 65,520 draws are not free even from a fast generator: that is 65,520 iterations of a loop that mostly does nothing. Algorithm L’s 413 draws include a logarithm each, which is more expensive per draw and 158 times fewer of them. The bit count is the honest way to compare the two, and it comes out of the same counter for both.

The general position this site takes is the same one it takes about comparisons and cache misses: a resource that is not counted is a resource that cannot be traded. Nobody chooses Algorithm L over Algorithm R for its bit consumption if nothing reports the bit consumption.

The rejections are charged, and that is a choice #

One detail of how the bits are counted is a cost model rather than a fact, and it is stated here because every figure that reports a bit count prints it.

Drawing a uniform integer below $k$ when $k$ is not a power of two is done by rejection: take $\lceil \log_2 k\rceil$ bits, and if the value is out of range, start again. The discarded attempts are charged. Two hundred draws below 3 cost 2.53 bits each rather than the 1.58 bits of information a draw below 3 actually carries, because the procedure spends more than the information it extracts.

Charging only the successful attempt would report the information-theoretic lower bound on the draw rather than what the algorithm spent, and those are different quantities with different uses. The lower bound answers “how much entropy does this algorithm fundamentally need”; the charged count answers “how much will it pull from the generator”. This site reports the second, because it is the one that a rate-limited entropy source or a cryptographic generator actually delivers.

The difference is small — about 60% overhead on the worst-case draw and much less on average — and it is the sort of thing that would silently change a reported constant if the convention drifted. So the convention is asserted: assertCoinsAreCounted requires two hundred draws below 3 to cost more than 400 bits, which fails if rejections ever stop being charged.

What the algorithm costs in the resources that were already counted #

For completeness, since two of the three axes are unremarkable and it is worth saying so.

Space is $k$ slots plus a constant, for both algorithms, at every stream length. That is the headline property and it does not vary: a sample of a thousand from a stream of a billion holds a thousand items. There is no version of this problem where the space grows.

Time is one pass, $\Theta(n)$ for Algorithm R — every item is examined even when no random number would have been needed. Algorithm L is $\Theta(k\log(n/k))$ in operations only if the stream can be skipped over without touching the skipped items, which is true for an array or a seekable file and false for a socket. Over a socket both algorithms are $\Theta(n)$ in reads and only the randomness differs, which is the case the table above is measured in.

Randomness is where they diverge, by the factor of 145 above.

Three resources, two algorithms, one difference. That is exactly the shape the space phase found between two merge sorts with identical comparison counts and different allocation totals, and it is the argument for having the third counter at all.

The parallel is worth making explicit because it is the same discovery twice. Two merge sorts, one allocating a single buffer and one allocating per merge, have byte-identical comparison, read, write and swap counts and identical peak space, and their allocation totals differ by a factor of $\log_2 n$. Nothing in the time counters sees it and nothing in the peak sees it; only a counter for the quantity in question does. Here, two reservoir samplers have identical output distributions, identical space and identical passes, and their randomness differs by a factor of 145. A resource with no counter is a resource where two very different algorithms look like the same algorithm.

What is being sampled, and what is not #

A closing precision, because “uniform sample” is doing specific work.

What reservoir sampling gives is a uniform sample without replacement of size $k$: every $k$-subset of the stream is equally likely. That implies each item is present with probability $k/n$, which is the quantity the figures measure, and the figures measure the implication rather than the full statement — checking that all $\binom{n}{k}$ subsets are equally likely would need vastly more runs and is not what has been done here. The essay claims what was measured.

What it does not give is a weighted sample, and weighting is the most common thing people actually want — sample proportional to size, to importance, to recency. The generalisation exists (the A-Res and A-ExpJ algorithms, which keep a priority queue keyed by $u^{1/w}$ for a uniform $u$ and weight $w$) and it is a different algorithm with a different proof, not a parameter of this one. Reaching for reservoir sampling when the requirement is weighted is a substitution that produces a plausible sample of the wrong distribution, which is the same failure mode as the off-by-one and just as quiet.