Reservoir sampling is a family of randomized online algorithms for choosing a simple random sample, without replacement, of k items from a population of unknown size n in a single pass over the items. The size of the population n is not known to the algorithm and is typically too large for all n items to fit into main memory. The population is revealed to the algorithm over time, and the algorithm cannot look back at previous items. At any point, the current state of the algorithm must permit extraction of a simple random sample without replacement of size k over the part of the population seen so far.
Consider the case when we want to maintain a sample of one item. A simple algorithm is to store the first item with probability 1 {\displaystyle 1} and the k {\displaystyle k} -th item with probability 1 k {\displaystyle {\frac {1}{k}}} . We want to show that the sample is uniformly distributed over all of the items that we have seen at each step since.
The proof is by induction on the number of items n {\displaystyle n} . If there is 1 item we always choose it and we are done. If there are already n {\displaystyle n} items, then the probability of one of the items already in the store remaining the store becomes n n + 1 × 1 n = 1 n + 1 {\displaystyle {\frac {n}{n+1}}\times {\frac {1}{n}}={\frac {1}{n+1}}} so the probability of any of the n + 1 {\displaystyle n+1} items remaining in the sample is now 1 n + 1 . {\displaystyle {\frac {1}{n+1}}.} The proof of correctness for the general algorithm where we have an arbitrary number of slots is an extension of this argument.
Reservoir sampling was introduced as early as 1962 by Fan, Muller, and Rezucha, who described a sampling algorithm which did not need to know the population size n in advance. A simpler algorithm, known as Algorithm R, was discovered independently by Alan G. Waterman and McLeod & Bellhouse; it can also be viewed as a variant of the Fisher–Yates shuffle. Algorithm R works as follows.
Initialize an array R {\displaystyle R} indexed from 1 {\displaystyle 1} to k {\displaystyle k} , containing the first k items of the input x 1 , . . . , x k {\displaystyle x_{1},...,x_{k}} . This is the reservoir.
For each new input x i {\displaystyle x_{i}} , generate a random number j uniformly in { 1 , . . . , i } {\displaystyle \{1,...,i\}} . If j ∈ { 1 , . . . , k } {\displaystyle j\in \{1,...,k\}} , then set R [ j ] := x i . {\displaystyle R[j]:=x_{i}.} Otherwise, discard x i {\displaystyle x_{i}} .
Return R {\displaystyle R} after all inputs are processed.
The correctness proof shows by induction on i that after the i'th step of the algorithm, the reservoir contains a uniformly random sample of the inputs x 1 , … , x j {\displaystyle x_{1},\ldots ,x_{j}} .
When i = k {\displaystyle i=k} , Algorithm R returns all inputs, so this is the base case of the induction.
For the inductive step, it is given that just before the ( i + 1 ) {\displaystyle (i+1)} -th input is processed, the reservoir forms a uniformly random sample of x 1 , … , x i {\displaystyle x_{1},\ldots ,x_{i}} . We must show that just after the ( i + 1 ) {\displaystyle (i+1)} -th input is processed, the reservoir forms a uniformly random sample of x 1 , … , x i + 1 {\displaystyle x_{1},\dots ,x_{i+1}} .
Apply Algorithm R to the ( i + 1 ) {\displaystyle (i+1)} -th input. Input x i + 1 {\displaystyle x_{i+1}} is included with probability k i + 1 {\displaystyle {\frac {k}{i+1}}} by definition of the algorithm. For any other input x r ∈ { x 1 , . . . , x i } {\displaystyle x_{r}\in \{x_{1},...,x_{i}\}} the probability that it is included in the reservoir just before the ( i + 1 ) {\displaystyle (i+1)} -th input is processed is k i {\displaystyle {\frac {k}{i}}} according to the inductive hypothesis. The probability that it is still included in the reservoir after x i + 1 {\displaystyle x_{i+1}} is processed (in other words, that x r {\displaystyle x_{r}} is not replaced by x i + 1 {\displaystyle x_{i+1}} ) is k i × ( 1 − 1 i + 1 ) = k i + 1 {\displaystyle {\frac {k}{i}}\times \left(1-{\frac {1}{i+1}}\right)={\frac {k}{i+1}}} as desired. Thus, all inputs have an equal probability k i + 1 {\displaystyle {\frac {k}{i+1}}} of being included in the reservoir. By induction, Algorithm R produces a uniform random sample of the inputs.
The Algorithm R described in the previous section has linear runtime, since it generates a random number for every element in the input stream to determine whether it should be included in the reservoir. It was observed by Jeffrey Vitter that it is possible to obtain faster algorithms with sublinear runtime when the algorithm is allowed to skip past sections of the input stream without examining them. Vitter gave a reservoir sampling algorithm which runs in time O ( k ( 1 + log ⁡ ( n / k ) ) {\displaystyle O(k(1+\log(n/k))} , and proved that this is optimal. Kim-Hung Li described a simple algorithm, Algorithm L, achieving the same asymptotic runtime, and which is described in this section.
The starting point for Algorithm L is the observation that if n {\displaystyle n} random numbers u 1 , . . . , u n ∼ U [ 0 , 1 ] {\displaystyle u_{1},...,u_{n}\sim U[0,1]} are generated uniformly and independently, then the indices of the smallest k {\displaystyle k} of them form a uniform sample of the k {\displaystyle k} -subsets of { 1 , . . . , n } {\displaystyle \{1,...,n\}} .
The process can be done without knowing n {\displaystyle n} :
Keep the smallest k {\displaystyle k} of u 1 , . . . , u i {\displaystyle u_{1},...,u_{i}} that has been seen so far, as well as w i {\displaystyle w_{i}} , the index of the largest among them. For each new u i + 1 {\displaystyle u_{i+1}} , compare it with u w i {\displaystyle u_{w_{i}}} . If u i + 1 < u w i {\displaystyle u_{i+1}<u_{w_{i}}} , then discard u w i {\displaystyle u_{w_{i}}} , store u i + 1 {\displaystyle u_{i+1}} , and set w i + 1 {\displaystyle w_{i+1}} to be the index of the largest among them. Otherwise, discard u i + 1 {\displaystyle u_{i+1}} , and set w i + 1 = w i {\displaystyle w_{i+1}=w_{i}} .
Now couple this with the stream of inputs x 1 , . . . , x n {\displaystyle x_{1},...,x_{n}} . Every time some u i {\displaystyle u_{i}} is accepted, store the corresponding x i {\displaystyle x_{i}} . Every time some u i {\displaystyle u_{i}} is discarded, discard the corresponding x i {\displaystyle x_{i}} .
This algorithm still needs O ( n ) {\displaystyle O(n)} random numbers, thus taking O ( n ) {\displaystyle O(n)} time. However, it can be improved by observing that after each step, the number of following items which will be skipped by the algorithm follows a geometric distribution with parameter u w i {\displaystyle u_{w_{i}}} . Thus, the sampling algorithm can simply generate a sample from this distribution to determine how many items to skip over, rather than examining them one at a time. A further simplification is to store only the maximum value u w i {\displaystyle u_{w_{i}}} at each point, and avoid storing u 1 , … , u i {\displaystyle u_{1},\ldots ,u_{i}} . This results in Algorithm L, which is implemented as follows:
This algorithm computes three random numbers for each item that becomes part of the reservoir, and does not spend any time on items that do not. Its expected running time is O ( k ( 1 + log ⁡ ( n / k ) ) ) {\displaystyle O(k(1+\log(n/k)))} , which is optimal.
Reservoir sampling can also be performed with replacement, meaning that the same input item may appear more than once in the final sample. In sampling without replacement, each input item can appear at most once. For a sample of size one, however, the resulting sample distribution is the same.
Park, Ostrouchov, Samatova, and Geist (2004) derived a skip-based algorithm for reservoir sampling with replacement.
For a reservoir of size one, suppose the currently selected item is at 1-based position n in the stream. The position of the next replacement can be sampled directly from a single uniform random value q:
This gives a particularly simple algorithm:
Here random() denotes a uniform random value strictly between 0 and 1. An implementation whose random-number generator can return an endpoint must exclude or resample such values. The variable n is the absolute 1-based position in the stream, rather than a count of items to skip. The algorithm therefore generates replacement positions directly instead of making a random replacement decision for every input item.
For a reservoir of size one, each new replacement position requires only one random variate. By comparison, Algorithm L normally uses three. When Algorithm L is specialized to k = 1, choosing the reservoir slot requires no randomness, leaving two nontrivial random variates per replacement.
A weighted generalization was later given by Meligrana and Fazzone.
For a single-item weighted reservoir with current cumulative weight W, the cumulative-weight threshold for the next replacement is
where q is uniform on (0,1). The stream is scanned until its cumulative weight reaches that threshold.
If we associate with each item of the input a uniformly generated random number, the k items with the largest (or, equivalently, smallest) associated values form a simple random sample. A simple reservoir-sampling thus maintains the k items with the currently largest associated values in a priority queue.
The expected running time of this algorithm is O ( n + k log ⁡ k log ⁡ ( n / k ) ) {\displaystyle O(n+k\log k\log(n/k))} and it is relevant mainly because it can easily be extended to items with weights.
The methods presented in the previous sections do not allow to obtain a priori fixed inclusion probabilities. Some applications require items' sampling probabilities to be according to weights associated with each item. For example, it might be required to sample queries in a search engine with weight as number of times they were performed so that the sample can be analyzed for overall impact on user experience. Let the weight of item i be w i {\displaystyle w_{i}} , and the sum of all weights be W. There are two ways to interpret weights assigned to each item in the set:
In each round, the probability of every unselected item to be selected in that round is proportional to its weight relative to the weights of all unselected items. If X is the current sample, then the probability of an item i ∉ X {\displaystyle i\notin X} to be selected in the current round is w i / ( W − ∑ j ∈ X w j ) {\displaystyle \textstyle w_{i}/(W-\sum _{j\in X}{w_{j}})} .
The probability of each item to be included in the random sample is proportional to its relative weight, i.e., w i / W {\displaystyle w_{i}/W} . It is not achievable when k = n {\displaystyle k=n} , since each item is then included with probability 1.
The following algorithm was given by Efraimidis and Spirakis that uses interpretation 1:
This algorithm is identical to the algorithm given in Reservoir Sampling with Random Sort except for the generation of the items' keys. The algorithm is equivalent to assigning each item a key r 1 / w i {\displaystyle r^{1/w_{i}}} where r is the random number and then selecting the k items with the largest keys. Equivalently, a more numerically stable formulation of this algorithm computes the keys as − ln ⁡ ( r ) / w i {\displaystyle -\ln(r)/w_{i}} and select the k items with the smallest keys.[failed verification]
The following algorithm is a more efficient version of A-Res, also given by Efraimidis and Spirakis:
This algorithm follows the same mathematical properties that are used in A-Res, but instead of calculating the key for each item and checking whether that item should be inserted or not, it calculates an exponential jump to the next item which will be inserted. This avoids having to create random variates for each item, which may be expensive. The number of random variates required is reduced from O ( n ) {\displaystyle O(n)} to O ( k log ⁡ ( n / k ) ) {\displaystyle O(k\log(n/k))} in expectation, where k {\displaystyle k} is the reservoir size, and n {\displaystyle n} is the number of items in the stream.
Warning: the following description is wrong, see Chao's original paper and the discussion here.
Following algorithm was given by M. T. Chao uses interpretation 2: and Tillé (2006).
For each item, its relative weight is calculated and used to randomly decide if the item will be added into the reservoir. If the item is selected, then one of the existing items of the reservoir is uniformly selected and replaced with the new item. The trick here is that, if the probabilities of all items in the reservoir are already proportional to their weights, then by selecting uniformly which item to replace, the probabilities of all items remain proportional to their weight after the replacement.
Note that Chao doesn't specify how to sample the first k elements. He simple assumes we have some other way of picking them in proportion to their weight. Chao: "Assume that we have a sampling plan of fixed size with respect to S_k at time A; such that its first-order inclusion probability of X_t is π(k; i)".
Similar to the other algorithms, it is possible to compute a random weight j and subtract items' probability mass values, skipping them while j > 0, reducing the number of random numbers that have to be generated.
In machine learning applications, fairness is a critical consideration, especially in scenarios where data streams exhibit class imbalance. To address this, Nikoloutsopoulos, Titsias, and Koutsopoulos proposed the Kullback-Leibler Reservoir Sampling (KLRS) algorithm as a solution to the challenges of Continual Learning, where models must learn incrementally from a continuous data stream.
The KLRS algorithm was designed to create a flexible policy that matches class percentages in the buffer to a target distribution while employing Reservoir Sampling techniques. This is achieved by minimizing the Kullback-Leibler (KL) divergence between the current buffer distribution and the desired target distribution.
KLRS generalizes earlier methods like Reservoir Sampling and Class-Balancing Reservoir Sampling, as verified by experiments using confidence intervals, demonstrating its broader applicability and improved performance.
The KLRS algorithm operates by maintaining a buffer of size and updating its contents as new data points arrive in a stream. Below is the pseudocode for the KLRS algorithm:
Reservoir sampling makes the assumption that the desired sample fits into main memory, often implying that k is a constant independent of n. In applications where we would like to select a large subset of the input list (say a third, i.e. k = n / 3 {\displaystyle k=n/3} ), other methods need to be adopted. Distributed implementations for this problem have been proposed.
