iv. Sampling and Streaming

Sampling and summarizing data that arrives as a stream: unknown length, one pass, memory far smaller than the data. If you do know the length and have the data, simply shuffle and take the first $k$, or pick $k$ random indices.

Problem	Algorithm	Space	Guarantee
uniform sample of $k$	reservoir sampling (Algorithm R; Algorithm L skips ahead)	$O(k)$	exact uniformity
weighted sample of $k$ (without replacement)	A-Res: keep the $k$ largest keys $u^{1/w}$	$O(k)$	exact
$k$ distinct indices from $[0,n)$	Floyd’s algorithm	$O(k)$	exact
frequent items (heavy hitters)	Misra–Gries ($k$ counters)	$O(k)$	undercounts by $\le \frac{n}{k+1}$
frequency of any item	Count–Min sketch	$O(\frac{1}{\epsilon }\log \frac{1}{\delta })$	overcounts by $\le \epsilon n$ w.p. $1-\delta$
membership	Bloom filter	$m$ bits	no false negatives; false positives $\approx (1-e^{-kn/m}{)}^k$
count distinct elements	HyperLogLog	$m$ small registers ($O(m\log \log n)$ bits)	relative error $\approx 1.04/\sqrt{m}$
approximate quantiles	t-digest, GK summary	$O(\frac{1}{\epsilon }\log n)$	rank error $\epsilon n$
majority element	Boyer–Moore voting	$O(1)$	exact if a strict majority exists: Classic Array Problems

When reservoir sampling applies	Why
input length unknown in advance	cannot pre-allocate or index
data does not fit in memory, single pass	$O(k)$ space, streaming
sampling log lines, telemetry, network packets	the standard production use
random element of a linked list, random node during a DFS	one traversal, $k=1$
random element of a filtered subset	apply the filter while streaming
distributed data	each machine keeps a reservoir and its count; merge by sampling from the reservoirs weighted by their counts

Reservoir sampling and exact samplers

template<class T> vector<T> reservoirR(const vector<T>& stream, int k, mt19937_64& rng) {   // Algorithm R
    vector<T> res;
    for (ll i = 0; i < sz(stream); i++) {
        if (i < k) res.push_back(stream[i]);
        else { ll j = uniform_int_distribution<ll>(0, i)(rng); if (j < k) res[j] = stream[i]; }   // keep w.p. k/(i+1)
    }
    return res;
}
 
template<class T> vector<T> reservoirL(const vector<T>& stream, int k, mt19937_64& rng) {   // Algorithm L (Li)
    uniform_real_distribution<double> U(0, 1); int n = sz(stream);
    vector<T> res(stream.begin(), stream.begin() + min(k, n));
    if (n <= k) return res;
    auto u = [&] { double x; do x = U(rng); while (x == 0); return x; };
    double W = exp(log(u()) / k); ll i = k - 1;
    while (true) {
        i += (ll)floor(log(u()) / log1p(-W)) + 1;    // skip a geometric number of items
        if (i >= n) break;
        res[rng() % k] = stream[i];
        W *= exp(log(u()) / k);
    }
    return res;
}
 
vector<int> weightedReservoir(const vector<double>& w, int k, mt19937_64& rng) {   // A-Res: indices, prob. ~ weight
    uniform_real_distribution<double> U(0, 1);
    priority_queue<pair<double, int>, vector<pair<double, int>>, greater<>> pq;   // min-heap of keys
    for (int i = 0; i < sz(w); i++) {
        double key = pow(U(rng), 1.0 / w[i]);
        if (sz(pq) < k) pq.push({key, i}); else if (key > pq.top().first) pq.pop(), pq.push({key, i});
    }
    vector<int> res; while (!pq.empty()) res.push_back(pq.top().second), pq.pop();
    return res;
}
 
vector<ll> floydSample(ll n, int k, mt19937_64& rng) {   // k distinct uniform values of [0, n), O(k) expected
    unordered_set<ll> S;
    for (ll j = n - k; j < n; j++) { ll t = uniform_int_distribution<ll>(0, j)(rng); if (!S.insert(t).second) S.insert(j); }
    return vector<ll>(all(S));
}

Proof: Algorithm R, the k = 1 case, A-Res and Floyd are exact

Algorithm R. Claim: after $n\ge k$ items each is in the reservoir with probability $\frac{k}{n}$. Induction: true for $n=k$. Item $n$ enters with probability $\frac{k}{n}$. An earlier item, present with probability $\frac{k}{n-1}$, survives if the new item is rejected ($1-\frac{k}{n}$) or accepted into another slot ($\frac{k}{n}\cdot \frac{k-1}{k}$): $\frac{k}{n-1}(1-\frac{k}{n}+\frac{k-1}{n})=\frac{k}{n-1}\cdot \frac{n-1}{n}=\frac{k}{n}$. The same induction shows every $k$-subset is equally likely. For $k=1$: keep the $i$-th item (0-based) with probability $\frac{1}{i+1}$.

Algorithm L produces the same distribution: in Algorithm R the reservoir is the $k$ items with the smallest of $n$ independent uniform labels; $W$ tracks the $k$-th smallest label so far (distributed as the max of $k$ uniforms, hence the update by $u^{1/k}$), and the number of items until a label falls below $W$ is geometric with parameter $W$, sampled directly. It draws $O(k\log \frac{n}{k})$ random numbers instead of $n$.

A-Res (Efraimidis and Spirakis). $\mathrm{Pr}[U^{1/w}\le x]=x^w$. For $k=1$, item $i$ has the largest key with probability ${\int }_0^1{\prod }_{j\ne i}{x}^{w_j}d(x^{w_i})=w_i{\int }_0^1{x}^{W-1}dx=\frac{w_i}{W}$; conditioning on the winner and repeating gives weighted sampling without replacement for $k>1$.

Floyd. Claim: after the iteration for $j$, $S$ is a uniform random subset of $[0,j]$ of size $m+1$, where $m=j-(n-k)$. Assume $S$ is a uniform $m$-subset of $[0,j-1]$ and $t$ is uniform in $[0,j]$. A given $(m+1)$-subset $S^{\prime }$ of $[0,j]$ containing $j$ arises from $S=S^{\prime }\setminus \{j\}$ with $t\in S\cup \{j\}$: $m+1$ values of $t$. A subset $S^{\prime }$ not containing $j$ arises from $S=S^{\prime }\setminus \{x\}$ and $t=x$, for each of the $m+1$ elements $x\in S^{\prime }$. Either way the probability is $\frac{m+1}{(j+1)\left(\genfrac{}{}{0ex}{}{j}{m}\right)}$, the same for every $S^{\prime }$. Checked by frequency tests: all four samplers match the exact inclusion probabilities within statistical error.

Sketches

ull splitmix(ull x) { x += 0x9e3779b97f4a7c15ULL; x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL; x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL; return x ^ (x >> 31); }
 
struct MisraGries {                                      // heavy hitters with k counters
    int k; map<ll, ll> c;
    MisraGries(int k) : k(k) {}
    void add(ll x) {
        if (c.count(x) || sz(c) < k) { c[x]++; return; }
        for (auto it = c.begin(); it != c.end();) { if (--it->second == 0) it = c.erase(it); else ++it; }   // decrement all
    }
    ll estimate(ll x) const { auto it = c.find(x); return it == c.end() ? 0 : it->second; }   // f - n/(k+1) <= est <= f
};
 
struct CountMin {                                        // width w ~ e / eps, depth d ~ ln(1 / delta)
    int w, d; vector<vector<ll>> t; vector<ull> seed;
    CountMin(int w, int d, ull s) : w(w), d(d), t(d, vector<ll>(w, 0)), seed(d) { for (int i = 0; i < d; i++) seed[i] = splitmix(s + i); }
    void add(ll x, ll c = 1) { for (int i = 0; i < d; i++) t[i][splitmix(x ^ seed[i]) % w] += c; }
    ll estimate(ll x) const { ll r = LLONG_MAX; for (int i = 0; i < d; i++) r = min(r, t[i][splitmix(x ^ seed[i]) % w]); return r; }
};
 
struct Bloom {                                           // m bits, k hash functions (double hashing)
    int m, k; vector<ull> bits; ull s;
    Bloom(int m, int k, ull s) : m(m), k(k), bits((m + 63) / 64), s(s) {}
    void add(ll x) { ull h1 = splitmix(x ^ s), h2 = splitmix(h1) | 1; for (int i = 0; i < k; i++) { ull b = (h1 + i * h2) % m; bits[b >> 6] |= 1ULL << (b & 63); } }
    bool maybe(ll x) const { ull h1 = splitmix(x ^ s), h2 = splitmix(h1) | 1; for (int i = 0; i < k; i++) { ull b = (h1 + i * h2) % m; if (!(bits[b >> 6] >> (b & 63) & 1)) return false; } return true; }
};
 
struct HyperLogLog {                                     // m = 2^b registers
    int b, m; vector<uint8_t> reg; ull s;
    HyperLogLog(int b, ull s) : b(b), m(1 << b), reg(1 << b, 0), s(s) {}
    void add(ll x) { ull h = splitmix(x ^ s); int idx = h >> (64 - b); ull rest = h << b; int r = rest ? __builtin_clzll(rest) + 1 : 64 - b + 1; reg[idx] = max<int>(reg[idx], r); }
    double estimate() const {
        double sum = 0; int zeros = 0;
        for (int r : reg) sum += ldexp(1.0, -r), zeros += r == 0;
        double alpha = 0.7213 / (1 + 1.079 / m), E = alpha * m * m / sum;
        if (E <= 2.5 * m && zeros) E = m * log((double)m / zeros);   // small-range correction (linear counting)
        return E;
    }
};

Proof: Misra–Gries, Count–Min and the Bloom filter bound

Misra–Gries. A counter is only incremented for its own item, so estimates never exceed true counts. Each “decrement all” step discards $k+1$ occurrences (the new item and one from each of the $k$ counters), so there are at most $\frac{n}{k+1}$ such steps, and an item loses at most one count per step: $f_x-\frac{n}{k+1}\le {\hat{f}}_x\le f_x$. Hence every item with $f_x>\frac{n}{k+1}$ survives (and a second pass can verify candidates).

Count–Min. Every row counts $x$ plus the colliding items, so each row overestimates. In one row, the expected collision mass is at most $\frac{n}{w}$ (pairwise independent hashing), so by Markov $\mathrm{Pr}[\text{row error}>\epsilon n]\le \frac{1}{w\epsilon }=\frac{1}{e}$ for $w=\frac{e}{\epsilon }$; the minimum over $d=\ln \frac{1}{\delta }$ independent rows fails with probability $\le e^{-d}=\delta$.

Bloom filter. An inserted element sets all its $k$ bits, so there are no false negatives. Assuming independent uniform hash positions, a fixed bit stays 0 after $n$ insertions with probability $(1-\frac{1}{m}{)}^{kn}\approx e^{-kn/m}$, so a non-member’s $k$ bits are all 1 with probability about $(1-e^{-kn/m}{)}^k$, minimized at $k=\frac{m}{n}\ln 2$, where it is $2^{-k}\approx 0.6185^{m/n}$ (about 1% at 9.6 bits per element). HyperLogLog’s estimator and its $1.04/\sqrt{m}$ standard error are the analysis of Flajolet, Fusy, Gandouet and Meunier (2007); here it was checked empirically (errors within a few standard errors on random streams). Count–Min and Misra–Gries were checked for their deterministic bounds, the Bloom filter for zero false negatives and a false-positive rate near the formula.

Related notes

• Randomized Algorithms
• Randomized Hashing and Verification
• Classic Array Problems (Boyer–Moore voting)