Reservoir Sampling

Reservoir Sampling

Reservoir sampling is a family of randomized algorithms designed to select a simple random sample, without replacement, of $k$ items from a population of unknown size $n$ in a single pass.^[17] The algorithm processes the population sequentially, ensuring that at any point during the execution, every item seen so far has an equal probability of being included in the sample. This mathematical guarantee is maintained regardless of how large the population grows or when the stream of data terminates.

The primary problem solved by reservoir sampling is the selection of an unbiased sample from massive data streams or files where the total number of elements cannot be known in advance, or where the dataset is simply too large to fit into main memory. In traditional simple random sampling, the total population size $n$ must be known to assign a uniform probability of $k/n$ to each item.^[18] When $n$ is unknown, or when the data is infinite, traditional methods fail because they either require multiple passes over the data (one to count the items, and another to sample them) or require loading the entire dataset into memory to perform a random shuffle.

By maintaining a fixed-size buffer, known as the "reservoir," these algorithms operate with a strict upper bound on memory consumption. As each new item arrives from the data stream, the algorithm uses a probabilistic mechanism to decide whether the new item should be ignored or whether it should replace an existing item in the reservoir.^[18] This single-pass, memory-bounded characteristic makes reservoir sampling a foundational technique in stream processing, database query approximation, and network telemetry.^[17]

Algorithmic Mechanism and Complexity

The most fundamental version of reservoir sampling is known as Algorithm R. The mechanism operates by maintaining an array of size $k$, which serves as the reservoir. The algorithm processes the incoming data stream in two distinct phases: initialization and iteration.

During the initialization phase, the algorithm unconditionally accepts the first $k$ items from the data stream and places them directly into the reservoir array.^[19] If the data stream contains fewer than $k$ items, the algorithm simply returns all available items, as the entire population fits within the requested sample size.

Once the reservoir is full, the iteration phase begins for all subsequent items. For the $i$-th item in the stream (where $i > k$), the algorithm generates a random integer $j$ between $1$ and $i$, inclusive. If the generated integer $j$ is less than or equal to $k$, the algorithm replaces the item at the $j$-th position in the reservoir with the new $i$-th item. If $j$ is greater than $k$, the $i$-th item is discarded, and the reservoir remains unchanged. This process repeats until the data stream is exhausted.

function ReservoirSample(stream, k):
    reservoir = array of size k
    for i from 1 to k:
        reservoir[i] = stream.next()
    
    i = k + 1
    while stream.hasNext():
        item = stream.next()
        j = randomInteger(1, i)
        if j <= k:
            reservoir[j] = item
        i = i + 1
        
    return reservoir

The space complexity of Algorithm R is strictly $O(k)$, as the algorithm only needs to store the $k$ items in the reservoir and a few integer variables for counting and random number generation. The time complexity is $O(n)$, where $n$ is the total number of items processed from the stream. This linear time complexity arises because the algorithm must generate a new pseudorandom number for every single item in the stream from index $k+1$ to $n$.

Mathematical Proof of Uniformity

The correctness of Algorithm R, specifically the guarantee that every item has an equal probability of $k/n$ of being in the final reservoir, is established through mathematical induction Reservoir sampling - Wikipedia^[1].

For the base case, consider the state of the reservoir after exactly $k$ items have been processed. The algorithm places the first $k$ items directly into the reservoir. Therefore, the probability that any of the first $k$ items is in the reservoir is $1$, which is mathematically equivalent to $k/k$. The base case holds true.

For the inductive step, assume that after processing $i - 1$ items (where $i - 1 \ge k$), every item seen so far has a probability of $k / (i - 1)$ of being present in the reservoir. The algorithm then processes the $i$-th item. According to the rules of Algorithm R, the $i$-th item is selected to enter the reservoir with a probability of $k / i$. This immediately satisfies the uniformity condition for the newly arrived item.

To prove that the previously processed items maintain their uniform probability, one must calculate the probability that an item already in the reservoir survives the $i$-th step. An item in the reservoir is replaced if and only if two independent events occur: the $i$-th item is selected to enter the reservoir (probability $k / i$), and the specific reservoir slot holding the existing item is chosen as the replacement target (probability $1 / k$). The probability that a specific item in the reservoir is replaced is the product of these two probabilities, which is $(k / i) \times (1 / k) = 1 / i$.

Consequently, the probability that a specific item in the reservoir is not replaced (i.e., it survives the $i$-th step) is $1 - (1 / i) = (i - 1) / i$. To find the total probability that an item from the first $i - 1$ elements is in the reservoir after the $i$-th step, its prior probability of being in the reservoir is multiplied by its probability of surviving the current step. Using the inductive assumption, this is $(k / (i - 1)) \times ((i - 1) / i) = k / i$.

Because both the new $i$-th item and all previously retained items have exactly a $k / i$ probability of being in the reservoir after $i$ steps, the inductive step is proven. By the time the stream terminates at $n$ items, every item has an inclusion probability of exactly $k / n$.

Discovery and Early History

Fan, Muller, and Rezucha's 1962 paper presented sequential sampling algorithms for when the population size is known, not unknown. Their research, published in the Journal of the American Statistical Association under the title "Development of Sampling Plans by Using Sequential (Item by Item) Selection Techniques and Digital Computers," addressed the physical constraints of early digital computing ^[2].^[20]

During the 1960s, large datasets were typically stored on magnetic tape drives. Magnetic tapes are inherently sequential storage mediums; reading data requires physically winding the tape over a read head. Rewinding a tape to perform multiple passes over a dataset was prohibitively slow and mechanically intensive. Furthermore, the exact number of records on a given tape was often unrecorded in any central directory, meaning statistical programs had to handle an unknown population size $n$. Fan, Muller, and Rezucha explored methods to draw random samples under these strict sequential and memory-limited conditions A note on sampling a tape-file^[3].

The specific algorithmic formulation known today as Algorithm R was discovered independently by Alan G. Waterman. Donald Knuth formalized the algorithm, named it Algorithm R, and attributed it to Waterman in his seminal 1969 publication, The Art of Computer Programming, Volume 2: Seminumerical Algorithms The Art of Computer Programming^[4]. Knuth's detailed analysis and publication brought the algorithm into mainstream computer science, establishing it as the standard method for single-pass uniform sampling. The algorithm was later analyzed and proven further by A. I. McLeod and D. R. Bellhouse in their 1983 paper "A Convenient Algorithm for Drawing a Simple Random Sample" A Convenient Algorithm for Drawing a Simple Random Sample^[5].

Later Refinements and Optimizations

While Algorithm R is mathematically sound, it exhibits a significant performance bottleneck when applied to massive datasets: it requires the generation of a pseudorandom number for every single item in the data stream. As the number of processed items $i$ grows much larger than the sample size $k$, the probability of accepting a new item ($k / i$) approaches zero. Consequently, the algorithm expends vast amounts of CPU time generating random numbers only to immediately reject the incoming items.

In 1985, Jeffrey S. Vitter at Brown University addressed this inefficiency by publishing "Random Sampling with a Reservoir" in ACM Transactions on Mathematical Software Random Sampling with a Reservoir^[6]. Vitter recognized that instead of deciding whether to accept or reject each item individually, an algorithm could randomly generate the number of items to skip before the next item is accepted into the reservoir. Vitter introduced Algorithms X, Y, and Z to compute this skip distance. The skip distance follows a geometric distribution, which changes dynamically as the stream progresses. Algorithm Z uses rejection sampling to generate these skip distances efficiently, reducing the time complexity from $O(n)$ to $O(k(1 + \log(n/k)))$. This optimization allows the algorithm to process billions of records rapidly by skipping large blocks of data without inspecting them or generating random numbers for them Faster methods for random sampling^[7].

In 1994, Kim-Hung Li published "Reservoir-Sampling Algorithms of Time Complexity O(n(1+log(N/n)))" (where $N$ denotes the population and $n$ denotes the sample size), introducing Algorithm L Reservoir-Sampling Algorithms of Time Complexity O(n(1+log(N/n)))^[8]. Algorithm L achieves the same asymptotic time complexity as Vitter's Algorithm Z but simplifies the mathematical implementation by using the inverse transform method instead of rejection sampling. Algorithm L maintains a floating-point variable $W$, which represents the largest random key in a simulated random sort. The skip distance is calculated directly using the logarithm of a uniform random variable divided by the logarithm of $1 - W$. After skipping the calculated number of items, the next item replaces a randomly chosen element in the reservoir, and $W$ is updated by multiplying it by $\exp(\log(u) / k)$, where $u$ is a new uniform random number in the range $(0, 1)$. Algorithm L is highly efficient and is frequently the variant implemented in modern standard libraries.

Important Variants

Standard reservoir sampling assumes that every item in the stream has an equal probability of selection. However, many statistical, machine learning, and network routing problems require weighted random sampling, where each item arrives with an associated weight $w_i > 0$, and the probability of the item being included in the final sample must be proportional to its weight.

The theoretical basis for unequal probability sampling was explored by M. T. Chao in his 1982 paper "A general purpose unequal probability sampling plan" A general purpose unequal probability sampling plan^[9]. However, adapting this to a single-pass streaming model with strict memory bounds required further innovation. In 2006, Pavlos S. Efraimidis and Paul G. Spirakis published "Weighted random sampling with a reservoir" in Information Processing Letters, introducing Algorithm A-Res Weighted random sampling with a reservoir^[10].

Algorithm A-Res operates using a key-based priority sampling method. For each incoming item $i$ with weight $w_i$, the algorithm generates a random number $u_i$ uniformly distributed between $0$ and $1$. The algorithm then calculates a key $k_i = u_i^{1/w_i}$. The reservoir is maintained as a priority queue (typically a min-heap) that stores the $k$ items with the largest keys. Because $u_i$ is a fraction, taking the $w_i$-th root pushes the value closer to $1$ for larger weights. Mathematically, this is equivalent to drawing from an exponential distribution, ensuring that heavier items have a proportionally higher probability of generating a large key and surviving in the reservoir Weighted Random Sampling over Data Streams^[11]. Efraimidis and Spirakis also introduced Algorithm A-ExpJ, which optimizes A-Res by generating skip distances for weighted streams, analogous to how Algorithm L optimizes Algorithm R. Recent research continues to expand on these methods, such as the 2026 work by Meligrana and Fazzone on weighted sampling with replacement Weighted Reservoir Sampling with Replacement from Data Streams^[12].

Distributed systems require another variant known as distributed reservoir sampling. In frameworks like MapReduce, data is partitioned across hundreds of independent nodes. Each node processes its local stream and maintains a local reservoir of size $k$. To form a global sample of size $k$, the local reservoirs cannot simply be concatenated. Instead, the system must merge them while accounting for the number of items processed by each node. A mathematically elegant method for this is random sort sampling (or bottom-k min-hash). During the initial pass, each node assigns a uniform random number $r_i$ to each item and keeps the $k$ items with the smallest $r_i$ values. When merging, the central coordinator collects all local reservoirs and simply selects the $k$ items with the smallest $r_i$ values overall. The described random sort sampling (bottom-k min-hash) method produces a sample without replacement, but the cited paper is about sampling with replacement Reservoir-based Random Sampling with Replacement from Data Stream^[13].

Sliding window reservoir sampling is a variant used when the population is an infinite stream, but the sample must only represent the most recent $W$ items. This requires a mechanism to expire items from the reservoir once they fall outside the active window, which complicates the probability calculations. Algorithms to maintain uniform samples over sliding windows rely on generating replacement candidates proactively, ensuring the reservoir remains populated even when older items expire.

Real-World Applications

Reservoir sampling is heavily utilized in big data processing frameworks. Apache Spark, a distributed computing engine, provides a sample transformation for Resilient Distributed Datasets (RDDs). When a user requests a random sample from an RDD without replacement, Spark often processes partitions of unknown sizes. Spark's sample transformation uses Bernoulli sampling based on a fraction, not reservoir sampling pyspark.RDD.sample^[14].

In the domain of search engines and distributed databases, Elasticsearch employs reservoir sampling to optimize expensive aggregation queries. When querying billions of documents, calculating exact metrics (such as percentiles or unique counts) consumes excessive CPU and memory resources. Elasticsearch's random_sampler aggregation uses uniform random sampling with a fixed probability (via geometrically distributed skip steps), not reservoir sampling, to collect a statistically representative subset of the documents matching the query. The engine then performs the requested aggregations on this small sample and scales the results by the inverse of the sampling probability, delivering highly accurate approximate answers in a fraction of the time Elasticsearch random sampler aggregation.

In-memory data stores also rely on these algorithms. Redis, a widely used key-value store, provides the SRANDMEMBER command, which allows users to retrieve $k$ random elements from a set. When the requested count is positive, Redis must return distinct elements. For very large sets, allocating a temporary array to shuffle the entire dataset would block the single-threaded Redis server. Instead, Redis SRANDMEMBER does not use reservoir sampling; it samples from in-memory hash tables or intsets since the population size is known and fully in memory Redis SRANDMEMBER^[15].

Network telemetry and monitoring protocols, such as sFlow and NetFlow, use reservoir sampling to monitor high-speed network traffic. Core internet routers process hundreds of millions of packets per second. It is physically impossible to store, analyze, or log every single packet. By implementing reservoir sampling in the router's hardware or firmware, the device can maintain a fixed-size sample of packets that accurately represents the overall traffic distribution, allowing network administrators to detect DDoS attacks or analyze bandwidth usage without degrading routing performance. The technique is also applied in computer graphics for Monte Carlo integration, where spatiotemporal reservoir resampling is used to render dynamic direct lighting in real-time ray tracing Reservoir Sampling - PBR Book.^[21]

Significance and Tradeoffs

The primary significance of reservoir sampling is its ability to provide a mathematically rigorous guarantee of uniform probability while strictly bounding memory usage to $O(k)$. In streaming environments, where data is either infinite or vastly exceeds available random-access memory, this bounded space complexity is an absolute requirement. Furthermore, the skip-based variants (Algorithm Z and Algorithm L) reduce the time complexity to a level where the overhead of sampling becomes negligible, allowing the algorithm to run inline with high-throughput data ingestion pipelines.

Despite its mathematical elegance, reservoir sampling involves specific tradeoffs. The algorithm is inherently designed for sequential access. If a dataset is small enough to fit entirely in main memory, and its exact size $n$ is known in advance, reservoir sampling is computationally inferior to algorithms designed for random access memory. In such cases, generating $k$ unique random indices, or using the Fisher-Yates shuffle to randomize the first $k$ elements of the array, is significantly faster and requires fewer CPU cycles Algorithm 235: Random permutation^[16].

Another tradeoff involves the computational overhead of random number generation. While Algorithm L minimizes the number of random variables generated, it requires floating-point arithmetic, specifically logarithms and exponentiation. In highly constrained embedded systems or low-level hardware implementations, calculating floating-point logarithms can introduce unacceptable latency.

Finally, weighted reservoir sampling variants, such as Algorithm A-Res, impose strict limitations on the input data. The algorithms require all weights to be strictly positive ($w_i > 0$). They cannot process items with zero weight or negative weights. If a dataset contains negative weights, the data must be transformed or shifted before reservoir sampling can be applied, which can distort the statistical properties of the sample if not handled with rigorous mathematical care.