Mining massive datasets with streaming, hashing, and algorithms that fit in one pass

Why this subject matters now

The subject grew up around a specific engineering situation. The web got larger than any machine, and the algorithms that assume random access to the whole input stopped being algorithms anyone could run. What survived that filter is a body of technique built on three moves, which are to hash instead of sort, to sample or sketch instead of store, and to bound the error instead of pretending there is none. Twenty years later the individual systems have churned (MapReduce is effectively retired) but the techniques are load-bearing in more places than ever. Every Redis deployment ships HyperLogLog, and BigQuery's APPROX_COUNT_DISTINCT is HyperLogLog++. Network telemetry runs on Count-Min and its descendants. Every major open LLM pretraining corpus (RefinedWeb, SlimPajama, FineWeb, Dolma) was deduplicated with MinHash-LSH, the exact algorithm derived below, and the retrieval layer of every RAG system is an approximate-nearest-neighbor index whose design is a direct response to the strengths and failures of locality-sensitive hashing. A practitioner today is expected to know not just that these tools exist but what their guarantees say, because the guarantees are the interface. A Bloom filter never lies in one direction, Count-Min never underestimates, a HyperLogLog with \(2^{14}\) registers has a standard error of 0.81 percent, and choosing parameters means doing the small calculation, not guessing. The other reason the subject matters is the cost model it teaches. Modern accelerators changed the arithmetic but not the conclusion. Moving data costs more than computing on it, so the number that decides whether an algorithm is viable is passes over the input and bytes across the network, not FLOPs. That way of counting is the durable content of this page.

Core theory

The setting, one pass in sublinear space, and the cost model

The data-stream model

The formal model is austere. A stream is a sequence \( a_1, a_2, \dots, a_N \) of items drawn from a universe \( U = \{1, \dots, n\} \), arriving in an order the algorithm does not control. The algorithm may keep a state of size \( O(\mathrm{polylog}(n, N)) \), far smaller than either the stream or the universe, and after one pass must answer queries about the whole stream, such as how many distinct items appeared, the frequency of item \(x\), which items exceeded a frequency threshold, or how many ones were in the last \(N\) bits. Two facts shape everything. First, exact answers are impossible in sublinear space. An algorithm that counts distinct elements exactly must distinguish \(2^n\) possible subsets of the universe and therefore needs \(\Omega(n)\) bits of state, an information-theoretic bound made precise by Alon, Matias, and Szegedy. Second, randomization plus approximation dissolves the impossibility. Allow the answer to be within a factor \((1 \pm \varepsilon)\) with probability \(1 - \delta\) and the space drops to something like \( O(\varepsilon^{-2} \log(1/\delta)) \), independent of the stream length. Every sketch on this page is an instance of that trade, a fixed-size randomized summary, linear in a useful sense, whose error is a theorem rather than a hope.

The cost model that decides viability

The second constraint is a hierarchy of costs. The figures below are order-of-magnitude and move every year, but the ratios have been stable for decades, and the ratios are what matter.

Three consequences follow. A disk-resident algorithm must be a sequential scanner, because one seek costs as much as reading a megabyte. This is why every algorithm in this subject is phrased as passes over the data, and why the difference between a two-pass and a three-pass frequent-itemset algorithm is a real engineering distinction rather than a constant factor. A distributed algorithm is priced by what crosses the network, not by what each machine computes, and the communication-cost model below makes that precise. And a sketch that fits in cache is effectively free to update, which is why a 12 KB HyperLogLog inside a database counts billions of rows without appearing in the profile. The recurring discipline, before asking whether an algorithm is fast, is to ask how many times it touches each byte and where that byte lives.

MapReduce, the shuffle, and the communication-cost model

The programming model

MapReduce, published by Dean and Ghemawat at OSDI 2004, is a two-phase template. A map function turns each input element into a list of key-value pairs. The system then performs the shuffle, in which all pairs with the same key, from every mapper on every machine, are routed to the same reducer. A reduce function receives one key with the list of all its values and emits the output. Word count is the canonical instance. Map emits \((w, 1)\) for each word occurrence, and reduce sums the list. The model's contribution was never the expressiveness, which is deliberately minimal, but the contract. The runtime handles partitioning, scheduling, re-execution of failed tasks (possible because map and reduce are pure functions of their inputs), and the routing. The programmer writes two functions and gets fault tolerance across thousands of unreliable machines.

input splits          map tasks                shuffle                reduce tasks
[ split 0 ] ──► map ──► (k,v) pairs ─┐   group by key, all-to-all  ┌─► reduce(k, [v...]) ─► out 0
[ split 1 ] ──► map ──► (k,v) pairs ─┼──► every mapper sends to ──►┼─► reduce(k, [v...]) ─► out 1
[ split 2 ] ──► map ──► (k,v) pairs ─┤    every reducer (disk +    └─► reduce(k, [v...]) ─► out 2
[ split 3 ] ──► map ──► (k,v) pairs ─┘    network materialization)

The shuffle is the real cost. Map output is written to local disk, partitioned by hashed key, pulled across the network by every reducer, and merged. It is an all-to-all exchange materialized through the slowest two resources in the hierarchy. A combiner, an associative pre-reduce run on each mapper's local output, exists purely to shrink this exchange. For word count it collapses a mapper's ten thousand \((\text{the}, 1)\) pairs into one \((\text{the}, 10000)\) before anything crosses the network. Combiners require the reduce operation to be associative and commutative, and the discipline of asking whether an aggregation admits a combiner is the same discipline as asking whether a sketch is mergeable, which is why sketches and MapReduce fit together so well. Every sketch in this page (Bloom, Count-Min, HyperLogLog) merges by cellwise OR, sum, or max, so each mapper can sketch its shard and ship kilobytes instead of the data.

The communication-cost model, with a lower bound

Afrati and Ullman formalized the right way to price a MapReduce round. Let each reducer receive at most \(q\) inputs (the reducer size, bounded by the memory of one machine), and define the replication rate \(r\) as the total number of input copies sent to reducers divided by the input size, the average number of reducers each input must visit. Replication rate is communication, and reducer size is parallelism (smaller \(q\) means more, smaller reducers). The model's content is that for a fixed problem these trade off against each other, and the trade can be lower-bounded exactly.

The worked bound. Consider the all-pairs problem, \(n\) inputs where every one of the \(\binom{n}{2}\) pairs must be examined by some reducer (this is a similarity join, comparing every document with every other). A reducer that receives \(q\) inputs can cover at most \(g(q) = \binom{q}{2} \le q^2/2\) pairs. If there are \(p\) reducers, covering all outputs requires

Total communication is at least \(pq \ge n^2/q\) input copies, so the replication rate obeys

The bound is tight. Partition the \(n\) inputs into \(2n/q\) groups of size \(q/2\), create one reducer per pair of groups, and send each group to the \(2n/q - 1\) reducers that need it. Each input is then replicated \(\Theta(n/q)\) times and each reducer holds \(q\) inputs. The reading is that with a billion documents and reducers that hold a million, every document must cross the network a thousand times, no matter how clever the code is. This is the theorem behind the practical instinct that all-pairs comparison at scale is infeasible and something like LSH, which avoids examining most pairs at the price of missing a bounded fraction, is not an optimization but a necessity.

The same model prices multiway joins. For the triangle join \( R(A,B) \bowtie S(B,C) \bowtie T(C,A) \) on \(k\) reducers, give each attribute a share, so reducers are indexed by \( (h_A(a), h_B(b), h_C(c)) \) with the three hash ranges multiplying to \(k\). A tuple of \(R(A,B)\) knows its \(A\) and \(B\) hashes but not \(C\), so it must be replicated to every value of the \(C\) share. Minimizing total communication by symmetry (or a short Lagrangian argument) gives each share \(k^{1/3}\). With \(k = 64\) reducers, each tuple is sent to \(64^{1/3} = 4\) reducers, and total communication is \(4(|R| + |S| + |T|)\), versus a replication of \(\sqrt{k} = 8\) for the naive cascade of two binary joins on skewed data. The 2010s' shared-nothing SQL engines (Hive, Presto/Trino, Spark SQL, BigQuery) all price plans in exactly this currency.

Spark and what actually survived

MapReduce the system is historical. Google retired it internally (Dataflow/Beam replaced it), and Hadoop MapReduce is legacy. What killed it was the rigidity of the two-phase template. Iterative algorithms, and almost everything on this page is iterative (PageRank, k-means, alternating least squares), had to be expressed as chains of jobs, each writing its entire state to the distributed filesystem and reading it back, so a 50-iteration PageRank paid 50 rounds of replicated disk I/O for data that never changed. Zaharia and colleagues' RDD abstraction (NSDI 2012), the core of Spark, fixed precisely this. A resilient distributed dataset is an immutable partitioned collection that remembers its lineage, the graph of deterministic transformations that produced it, so it can live in memory and be recomputed per-partition after a failure instead of being replicated to disk as insurance. Transformations build a DAG, and the scheduler cuts the DAG into stages at shuffle boundaries and pipelines everything within a stage. The honest summary is that Spark changed where intermediate data lives, not what costs. The stage boundary is still the shuffle, a wide dependency (groupByKey, join, repartition) still materializes an all-to-all exchange, and the communication-cost model above transfers verbatim. MapReduce is dead, but its accounting is not.

Finding similar items with shingles, MinHash, and LSH

Jaccard similarity and shingling

The problem is to find, among \(n\) documents, the pairs that are nearly identical, where \(n\) is large enough that the \(\binom{n}{2}\) comparisons forbidden by the lower bound above are out of the question. The first move is to represent a document as a set, because sets have a similarity measure with exactly the right hashing structure. A \(k\)-shingle is any substring of length \(k\) (characters, or more commonly word \(k\)-grams), and the document becomes the set of its shingles. \(k\) is chosen so that shingles are rare enough to be discriminative, \(k = 5\) words or \(k \approx 9\) characters for web text, small enough that a light rewrite still shares most shingles, large enough that unrelated documents share few. Similarity between shingle sets is Jaccard similarity,

the fraction of either set's combined content that is shared. Two near-duplicate news articles might have \(J \approx 0.9\), and two unrelated pages \(J \approx 0.01\). The sets are still large (a page has thousands of shingles), so the second move compresses each set to a short signature that preserves Jaccard similarity in expectation.

MinHash, with the collision proof

Fix a random permutation \(\pi\) of the universe of shingles. The MinHash of a set \(S\) under \(\pi\) is the element of \(S\) that appears first in the permuted order, \( h_\pi(S) = \argmin_{x \in S} \pi(x) \). The theorem, due to Broder (1997), is that the collision probability of one MinHash equals the Jaccard similarity exactly.

Proof. Partition the elements of \( S \cup T \) into two classes, the \(a = |S \cap T|\) elements in both sets and the \(b = |S \cup T| - |S \cap T|\) elements in exactly one. Elements outside \(S \cup T\) are irrelevant, since they affect neither MinHash. Scan the universe in the order given by \(\pi\), stop at the first element of \(S \cup T\), and call it \(x\). Every element of \(S \cup T\) is equally likely to be \(x\), because \(\pi\) is a uniform permutation. If \( x \in S \cap T \), then \(x\) is simultaneously the first element of \(S\) and the first element of \(T\) in permuted order, so \( h_\pi(S) = h_\pi(T) = x \). If \(x\) is in exactly one set, say \(S\), then \( h_\pi(S) = x \) while \( h_\pi(T) \) is some element appearing strictly later, so the hashes differ. The two MinHashes therefore agree exactly when the first-arriving element of the union lies in the intersection, an event of probability \( a / (a + b) = |S \cap T| / |S \cup T| = J(S,T) \). \(\square\)

One MinHash is one Bernoulli trial with success probability \(J\). Take \(n_h\) independent permutations (in practice, \(n_h\) independent hash functions standing in for permutations, a substitution justified by Broder, Charikar, Frieze, and Mitzenmacher's min-wise independent families) and stack the results into a signature of \(n_h\) integers per document. The fraction of agreeing positions is an unbiased estimator of \(J\) with standard deviation \( \sqrt{J(1-J)/n_h} \), about \(0.035\) at \(n_h = 200\) and \(J = 0.5\). A verified run of the implementation below on two overlapping integer ranges with true \(J = 0.5\) returned an estimate of \(0.545\) from 200 hashes, one standard deviation off, as expected. The signature matrix (rows are hash functions, columns are documents) has replaced megabyte documents with a few hundred bytes each while preserving the one number pairwise comparison needs.

LSH by banding, and the S-curve

Signatures compress the comparison but not the number of comparisons, since \(\binom{n}{2}\) remains. Locality-sensitive hashing (Indyk and Motwani 1998, with the banding form popularized by Leskovec, Rajaraman, and Ullman's textbook) fixes this by hashing signatures so that only probable matches collide. Split the \(n_h = b \cdot r\) signature rows into \(b\) bands of \(r\) rows. Within each band, hash the \(r\)-row slice of each column into a table. Two documents become a candidate pair if they collide in at least one band, that is, agree on all \(r\) rows of some band. Only candidate pairs are ever verified.

signature matrix (n_h = b·r rows, one column per doc)
┌────────────┐
│ band 1 (r) │──hash──►  table 1 ─┐
│ band 2 (r) │──hash──►  table 2 ─┼─► candidate pair if any collision
│   ...      │                    │
│ band b (r) │──hash──►  table b ─┘
└────────────┘

Derive the candidate probability for a pair with Jaccard similarity \(s\). Each signature row agrees independently with probability \(s\) (the collision theorem), so one band of \(r\) rows agrees entirely with probability \(s^r\). The band fails with probability \(1 - s^r\), all \(b\) bands fail with probability \((1 - s^r)^b\), and hence

This is the S-curve. For small \(s\), \(s^r\) is tiny and the probability is approximately \(b\,s^r \approx 0\), for large \(s\) it saturates at 1, and in between it rises steeply. The threshold, where the rise is steepest, sits approximately where \( b \, s^r = 1 \), i.e.

Concretely, \(b = 20\) bands of \(r = 5\) rows (100 hashes) give threshold \( (1/20)^{1/5} = 0.549 \), and the curve evaluates to the following.

A pair at \(s = 0.8\) is missed with probability \(0.00036\) (a false negative), while a pair at \(s = 0.3\) costs a wasted verification with probability \(0.047\) (a false positive, harmless to correctness but costly in time). To choose \(b\) and \(r\) for a target threshold \(t\) with a hash budget \(n_h\), fix \(n_h = br\) and tabulate the few factorizations of \(n_h\), picking the one whose threshold \((1/b)^{1/r}\) is nearest the target, biasing toward larger \(r\) (a sharper curve, fewer false positives) when verification is expensive and toward smaller \(r\) when false negatives are expensive. The worked problem below does this calculation for a deduplication setting.

SimHash, random hyperplanes for cosine similarity

Jaccard is the right measure for sets. For real-valued vectors (TF-IDF documents, embeddings) the natural measure is the angle. Charikar's SimHash (2002) supplies an LSH family for it. Draw a random vector \(w\) with i.i.d. \(\mathcal N(0,1)\) entries and define the one-bit hash \( h_w(x) = \mathrm{sign}(w \cdot x) \). The claim is that

Proof. Only the component of \(w\) in the plane \( \mathrm{span}(x, y) \) affects the two signs, and because a spherically symmetric Gaussian projected onto any 2-plane is a spherically symmetric 2D Gaussian, the projected direction of \(w\) is uniform on the circle. In that plane, \(w \cdot x\) and \(w \cdot y\) have opposite signs exactly when the line perpendicular to \(w\) separates \(x\) from \(y\), which happens when the direction of \(w\) falls in one of the two arcs of angular width \(\theta\) (one on each side). The probability is \( 2\theta / 2\pi = \theta/\pi \). \(\square\) Numerically, at \( \theta = 60^\circ \) the collision probability is \( 1 - 60/180 = 2/3 \), and a Monte Carlo run with 200,000 random hyperplanes measured 0.6669. Concatenating 64 such bits gives a 64-bit fingerprint whose Hamming distance is a binomial estimate of \( 64\,\theta/\pi \). Manku, Jain, and Das Sarma used exactly this at Google (WWW 2007) to near-duplicate web pages by finding fingerprints within Hamming distance 3.

p-stable projections for Euclidean distance

For \(L_2\) distance, Datar, Immorlica, Indyk, and Mirrokni (2004) built an LSH family from stable distributions. A distribution \(\D\) is \(p\)-stable if for any fixed vector \(a\), the weighted sum \( \sum_i a_i X_i \) of i.i.d. draws \(X_i \sim \D\) is distributed as \( \|a\|_p \, X \) for a single draw \(X \sim \D\). The Gaussian is 2-stable (a linear combination of Gaussians is Gaussian with the right variance), and the Cauchy is 1-stable. Hash by projecting and quantizing,

For two points \(u, v\) at distance \( c = \|u - v\|_2 \), the projection difference \( a \cdot (u - v) \) is distributed as \( c \, Z \) with \(Z\) standard normal, by 2-stability. The two points collide when their projections land in the same width-\(w\) bin, which for the random offset \(\beta\) gives

where \(f\) is the density of \(|Z|\). The first factor is the chance the projected gap equals \(t\), and the second is the chance the random bin boundary does not fall in the gap. \(p(c)\) is monotonically decreasing in \(c\), which is all the LSH machinery needs, and the same AND-OR banding construction then amplifies the gap between near and far. This family is what "LSH" means in most ANN benchmarks, and its weakness, which the modern-practice section returns to, is that the partitions are data-oblivious. The hyperplanes and bins are drawn without looking at the data, so real datasets with low-dimensional structure need many tables to reach high recall.

Frequent itemsets, A-Priori and its descendants

The market-basket model

The data is a file of baskets, each a small set of items from a large universe, such as purchases, words in documents, or side effects per patient. The file is disk-resident and read sequentially. The working assumption is that passes over the file dominate cost, so algorithms are graded by pass count. An itemset is frequent if its support, the number of baskets containing all its items, is at least a threshold \(s\). From frequent itemsets come association rules \( I \to j \) with confidence \( \mathrm{conf} = \mathrm{supp}(I \cup \{j\}) / \mathrm{supp}(I) \) and interest \( \mathrm{conf} - \P[j] \), but the computational problem is the itemsets, and the bottleneck is pairs. With \(10^5\) items there are \(5 \times 10^9\) pairs, and counting them all needs 20 GB of 4-byte counters, which is the memory wall the algorithms below negotiate.

A-Priori and monotonicity

The engine is one observation, monotonicity. If \(I\) is frequent, every subset of \(I\) is frequent, because every basket containing \(I\) contains each subset, and support can only grow when items are removed. The contrapositive is that an itemset with an infrequent subset cannot be frequent and need never be counted. A-Priori (Agrawal and Srikant, VLDB 1994) turns this into a level-wise algorithm. Pass 1 counts single items, yielding the frequent items \(L_1\). Pass 2 counts only pairs with both items in \(L_1\), yielding \(L_2\). Pass \(k\) counts only candidate \(k\)-sets whose every \((k-1)\)-subset is in \(L_{k-1}\). Real basket data prunes heavily. Item frequency is heavy-tailed, most items are infrequent, and the candidate pair count falls from \(\binom{n}{2}\) to \(\binom{|L_1|}{2}\).

PCY, multistage, and multihash

A-Priori's pass 1 uses almost no memory (one counter per item), which Park, Chen, and Yu observed is a waste of the machine. PCY spends the idle memory on a hash table of counters. During pass 1, for each basket, hash every pair it contains into one of \(B\) buckets and increment that bucket. A bucket whose total is below \(s\) cannot contain any frequent pair (the bucket count upper-bounds every member pair's support), so before pass 2 the bucket counts are collapsed into a bitmap, one bit per bucket, and pass 2 counts a pair only if both items are frequent and the pair hashes to a frequent bucket. With counters of the infrequent majority never allocated, the pair-counting table often fits in memory when A-Priori's would not. Multistage adds a second pass with a fresh hash function, rehashing only pairs that survived the first bitmap, shrinking false positives geometrically at the cost of a pass. Multihash runs two hash tables in the same pass with half the buckets each, trading per-table resolution for independence. All three are the same idea as a Bloom filter, met later, hashed counters as a cheap, one-sided pre-filter.

SON and Toivonen, two passes or one and a gamble

SON (Savasere, Omiecinski, Navathe) makes frequent-itemset mining distributed and exactly two-pass. Split the file into chunks that fit in memory, mine each chunk in isolation at the scaled threshold \(p s\) where \(p\) is the chunk's fraction of the file, take the union of all chunk-frequent itemsets as candidates, and count the candidates exactly in a second full pass. No false negatives are possible. An itemset frequent in the whole file, with support at least \(s\) over the file, must reach fraction-scaled support \(ps\) in at least one chunk (if it fell below \(ps\) in every chunk, summing gives total support below \(s\), a contradiction, and this is monotonicity in averaged form). The two phases are one MapReduce round each, and SON is why frequent itemsets parallelize cleanly.

Toivonen's algorithm gambles for one pass plus a usually-small correction. Mine a random sample at a deliberately lowered threshold (say \(0.8\, p s\)) to make false negatives unlikely. Compute the negative border, the itemsets not frequent in the sample all of whose immediate subsets are. Then in one full pass count both the sample-frequent sets and the border. If no border set turns out frequent in the full data, the sample-frequent sets that verified are exactly the frequent itemsets, with a proof. Any itemset frequent in the full data but missed by the sample would have a smallest missed subset, which by construction lies on the border and would have been caught frequent, a contradiction. If some border set is frequent, the pass is inconclusive and must be repeated with a new sample. Lowering the sample threshold trades memory for repeat probability.

Streaming sketches, with their guarantees derived

Reservoir sampling, uniformity proved by induction

The simplest streaming primitive is to maintain a uniform random sample of \(k\) items from a stream of unknown length, using \(O(k)\) memory and one pass. Vitter's Algorithm R stores the first \(k\) items, and on seeing item \(i > k\), draws \( j \sim \mathrm{Unif}\{1, \dots, i\} \) and, if \( j \le k \), replaces slot \(j\) with the new item (otherwise discarding it).

Claim. After \(n \ge k\) items, every item seen so far is in the reservoir with probability exactly \(k/n\).

Proof by induction on \(n\). In the base case \(n = k\), all \(k\) items are stored, probability \(1 = k/k\). For the inductive step, assume each of the first \(n\) items is present with probability \(k/n\). Item \(n+1\) arrives and is admitted with probability \( \P[j \le k] = k/(n+1) \), which is its required retention probability. For any earlier item \(x\) currently in the reservoir, \(x\) is evicted only if item \(n+1\) is admitted (probability \(k/(n+1)\)) and the replaced slot is the one holding \(x\) (probability \(1/k\), since \(j\) is uniform over the \(k\) slots conditioned on admission). So \(x\) survives this step with probability

and its total presence probability is \( \frac{k}{n} \cdot \frac{n}{n+1} = \frac{k}{n+1} \), completing the induction. \(\square\) In an empirical check of the implementation below, sampling \(k = 5\) from a stream of 20 across 200,000 trials, every item's observed frequency fell in \([0.2473,\, 0.2524]\) against the exact \(0.25\). Weighted variants (Efraimidis and Spirakis, who key each item by \(u^{1/w}\) for \(u \sim \mathrm{Unif}(0,1)\) and keep the top \(k\) keys) extend the idea, and reservoir sampling remains the honest baseline against which more elaborate stream summaries are judged. When a question can be answered from a uniform sample, nothing simpler exists.

Bloom filters, with the false-positive rate derived and optimized

A Bloom filter answers set membership with one-sided error, where "no" is always correct and "yes" may be wrong. The structure is a bit array of \(m\) bits and \(k\) hash functions. Inserting \(x\) sets bits \( h_1(x), \dots, h_k(x) \), and querying \(y\) reports present when all \(k\) bits are set. False negatives are impossible because set bits are never cleared. For the false-positive rate, note that after inserting \(n\) items, each of the \(kn\) hash evaluations misses a particular bit with probability \(1 - 1/m\), so a given bit is still zero with probability

using \( (1 - 1/m)^m \to e^{-1} \). A query on an absent item reads \(k\) bits and returns a false positive when all are set, which under the mild and fixable approximation that bits are independent (Mitzenmacher and Upfal treat the dependence carefully) gives

To optimize \(k\) for fixed \(m/n\), write \( p = e^{-kn/m} \), so \( \ln p_{\mathrm{fp}} = k \ln(1 - p) \) with \( k = -\frac{m}{n} \ln p \). Then \( \ln p_{\mathrm{fp}} = -\frac{m}{n} \ln p \ln(1 - p) \), and the product \( \ln p \ln(1-p) \) is symmetric under \( p \leftrightarrow 1 - p \), maximized in magnitude at \( p = 1/2 \). So the optimum keeps each bit set with probability one half (the filter is at maximum entropy, which is the right intuition, since a half-full filter carries the most information per bit), giving

Inverting for sizing gives \( m = n \ln(1/p_{\mathrm{fp}}) / (\ln 2)^2 \), which is \(1.44 \log_2(1/p_{\mathrm{fp}})\) bits per element regardless of element size, 9.6 bits per element for 1 percent and 14.4 for 0.1 percent. In the measured check, \(n = 100{,}000\) keys, \( m = 958{,}506 \) bits (117 KB), and \(k = 7\) predict \( (1 - e^{-0.7303})^7 = 1.00\% \), and querying 100,000 absent keys against the implementation below measured 0.94 percent.

Count-Min sketch, frequency estimation via Markov and a union bound

Cormode and Muthukrishnan's Count-Min sketch (2005) estimates item frequencies from \(d \times w\) counters. Row \(j\) has a hash \( h_j: U \to \{1, \dots, w\} \). On arrival of \(x\), increment \( C[j, h_j(x)] \) for every row, and estimate

Every counter that \(x\) touches contains \(f_x\) plus the counts of colliding items, so \( \hat f_x \ge f_x \) always. The error is one-sided, an overestimate. Bound it. Fix row \(j\) and let the excess be \( X_j = C[j, h_j(x)] - f_x = \sum_{y \ne x} f_y \, \mathbf 1[h_j(y) = h_j(x)] \). With \(h_j\) pairwise independent, \( \P[h_j(y) = h_j(x)] = 1/w \), so

where \( N = \sum_y f_y \) is the stream length. Choose \( w = \lceil e/\varepsilon \rceil \) so \( \E[X_j] \le \varepsilon N / e \). \(X_j\) is nonnegative, so Markov's inequality applies with no variance needed,

The rows use independent hashes, so all \(d\) rows overshoot simultaneously with probability at most \( e^{-d} \), and the minimum over rows fails only if all rows fail. Setting \( d = \lceil \ln(1/\delta) \rceil \) gives

The guarantee is additive in \(N\), which is the sketch's signature. Heavy items are estimated well in relative terms, and light items are buried in the \(\varepsilon N\) noise floor. Measured against theory, \( \varepsilon = 10^{-3}, \delta = 10^{-2} \) gives \( w = 2719, d = 5 \), 13,595 counters (106 KB at 8 bytes). On a Zipf(1.2) stream of \( N = 858{,}625 \) items, the bound permits error up to \( \varepsilon N = 859 \) with 1 percent exceptions. The measured maximum error over 2,000 queried items was 130, mean 12.8, and zero queries exceeded the bound. Theory is loose in the typical direction, since Markov plus a union bound gives a worst-case guarantee, and real skewed streams sit far inside it.

The counting-Bloom comparison is worth making precisely. A counting Bloom filter also keeps hashed counters, but all \(k\) hashes index one shared array, so a query's minimum is over cells that other items' insertions also incremented via any of their \(k\) hashes. Its natural analysis targets membership (is the count nonzero) rather than frequency, and deletions are its selling point. Count-Min's separation into \(d\) independent rows is what makes the clean Markov-per-row, union-over-rows argument work and gives the \((\varepsilon, \delta)\) frequency guarantee. Conservative update (increment only the cells equal to the current minimum) tightens Count-Min further on skewed data at the cost of supporting only increments.

Count sketch and the AMS second-moment estimator

Count-Min's one-sided error comes from every collision adding. The Count sketch (Charikar, Chen, and Farach-Colton 2002) adds a random sign. Row \(j\) carries a second hash \( g_j: U \to \{-1, +1\} \), updates \( C[j, h_j(x)] \mathrel{+}= g_j(x) \), and estimates \( \hat f_x = \mathrm{median}_j g_j(x)\, C[j, h_j(x)] \). Colliding items now cancel in expectation, \( \E[g_j(x) C[j, h_j(x)]] = f_x \) exactly, because each colliding \(y\) contributes \( f_y \E[g_j(x) g_j(y)] = 0 \). The variance of one row's estimate is \( \sum_{y \ne x} f_y^2 / w \le F_2 / w \), so the error scale is \( \varepsilon \sqrt{F_2} \) with \( w = O(1/\varepsilon^2) \) rather than \( \varepsilon F_1 \), which is better on skewed streams (where \( \sqrt{F_2} \ll N \)), at the price of two-sided error and \(1/\varepsilon^2\) width. The median over rows converts the per-row constant-probability guarantee into \(1 - \delta\) via a Chernoff argument. The same signed-sum trick is feature hashing in ML (Weinberger et al. 2009). The Count sketch of a sparse feature vector is a linear, inner-product-preserving projection, which is why "the hashing trick" and this sketch are the same mathematics.

The ancestor is the AMS sketch (Alon, Matias, Szegedy 1996) for the second moment \( F_2 = \sum_x f_x^2 \). Keep \( Z = \sum_x f_x\, s(x) \) for a 4-wise independent sign function \(s\) (updatable in a stream by adding \(s(x)\) per arrival), and estimate \( \hat F_2 = Z^2 \). Unbiasedness follows from

since \( \E[s(x)s(y)] = \mathbf 1[x = y] \) by pairwise independence. For the variance, 4-wise independence makes \( \E[Z^4] = \sum_x f_x^4 + 3 \sum_{x \ne y} f_x^2 f_y^2 \) (only paired-up index patterns survive), so \( \Var[Z^2] = \E[Z^4] - F_2^2 \le 2 F_2^2 \). Averaging \( O(1/\varepsilon^2) \) independent copies and taking a median of \( O(\log(1/\delta)) \) group means gives the standard \( (1 \pm \varepsilon) \) estimate. A small measured run, with 64 sign-counters arranged as a median of 8 means on a 20,000-item skewed stream, estimated \(F_2\) with 20 percent relative error, consistent with the \( \sqrt{2/8} = 50\% \) per-group standard deviation that the median then tightens. Production accuracy needs the hundreds of counters the \(1/\varepsilon^2\) law demands, and \(F_2\) in practice is usually consumed via the Count sketch (whose row variance is exactly an AMS estimate) rather than alone.

HyperLogLog, counting distinct items in kilobytes

The distinct-count problem asks how many different items appear in the stream, in memory that does not grow with the answer. The idea chain has three links. First, hash every item to a uniform bit string. Duplicates hash identically, so the multiset becomes a set of uniform random values and cardinality becomes a property of that random set. Second, there is the pattern observed by Flajolet and Martin, that among \(n\) uniform bit strings the maximum number of leading zeros \(R\) is concentrated near \( \log_2 n \), because a run of \(\rho\) leading zeros has probability \(2^{-\rho}\) and appears at all once \( n \approx 2^\rho \). So \( 2^{R} \) estimates \(n\), but with very large variance, since one lucky hash doubles the estimate. Third, stochastic averaging (Flajolet, Fusy, Gandouet, Meunier 2007). Use the first \(p\) hash bits to route each item to one of \( m = 2^p \) registers, and let register \(j\) keep \( M_j \), the maximum leading-zero count (plus one) seen among its items' remaining bits. Each register estimates its substream's cardinality \( \approx n/m \), and combining them by the harmonic mean suppresses the outlier registers that plague the arithmetic mean,

where \(\alpha_m\) corrects the harmonic mean's multiplicative bias (its derivation via the Mellin transform is the hard part of the 2007 paper and is out of scope here, though the analysis also yields the relative standard error). The accuracy law is

Registers store leading-zero counts, at most \( 64 - p \), so 6 bits suffice, and \( m = 2^{14} \) registers is 12 KB for 0.81 percent standard error at any cardinality up to billions, which is exactly Redis's configuration for PFCOUNT. Two corrections handle the ends of the range. At small \(n\) most registers are zero and the raw estimator biases high, so when \( \hat n \le 2.5m \) the algorithm switches to linear counting, \( \hat n = m \ln(m/V) \) with \(V\) the number of zero registers (the occupancy estimator). At the top of a 32-bit hash range a saturation correction applies, or one uses 64-bit hashes and forgets it, the HyperLogLog++ choice (Heule, Nunkesser, Hall 2013) along with sparse storage of small sketches. In a hand-scale example, \( m = 4 \) (so \( \alpha_4 = 0.673 \)) with registers \( M = [3, 4, 2, 3] \) gives \( \sum 2^{-M_j} = 0.125 + 0.0625 + 0.25 + 0.125 = 0.5625 \), \( \hat n = 0.673 \times 16 / 0.5625 = 19.1 \). In the measured check at realistic scale, \( p = 11 \) (\( m = 2048 \), about 1.2 KB) predicts RSE \( 1.04/\sqrt{2048} = 2.30\% \). One run over \(10^6\) distinct items estimated 991,581 (0.84 percent low), and 20 independent runs at \(10^5\) distinct items had RMS relative error 2.03 percent, inside the predicted 2.30. HyperLogLogs merge by registerwise max, which is why they federate across shards and days for free, and why every serious analytics store exposes them.

Misra-Gries and Space-Saving, heavy hitters deterministically

Frequency estimation has a deterministic corner. The Misra-Gries algorithm (1982), generalizing Boyer-Moore majority, keeps at most \( k - 1 \) candidate counters. On arrival of \(x\), if \(x\) has a counter, increment it. If not and a counter slot is free, start \(x\) at 1. Otherwise decrement every counter, dropping those that reach zero (the arriving item is also discarded). The bound is

Derivation. The estimate never exceeds \(f_x\) because a counter increments only on arrivals of its own item. For the lower bound, count decrements. Each decrement event destroys \(k\) units of count at once, the \(k-1\) stored counters plus the arriving item's own unit, and the total count mass ever created is \(N\), so there are at most \( N/k \) decrement events. Any single item's counter (including its discarded arrivals) loses at most one unit per event, hence at most \(N/k\) in total. \(\square\) Consequently any item with \( f_x > N/k \), any true heavy hitter at the \(1/k\) level, must survive with a positive counter, so there are no false negatives. Measured with \( k = 100 \) counters on the same Zipf stream (\(N = 858{,}625\), bound \( N/k = 8586 \)), the largest underestimate among the twenty most frequent items was 4,331, within the bound, and every item above the \(N/k\) line was retained. Space-Saving (Metwally, Agrawal, El Abbadi 2005) is the practical twin. When full, it evicts the minimum counter and gives its count (plus one) to the newcomer, overestimating by at most the evicted minimum, which it tracks per counter. Its error bound is the same \(N/k\), its estimates are one-sided high rather than low, and its "stream summary" data structure makes updates \(O(1)\). Space-Saving is what production top-K systems (including Redis's TOPK type) actually run.

DGIM, counting ones in a sliding window

Streams often need windowed answers, such as how many ones appeared in the last \(N\) bits, where storing the window itself (\(N\) bits) is too much. Datar, Gionis, Indyk, and Motwani (2002) achieve \( O(\log^2 N) \) bits with a guaranteed at most 50 percent error, improvable to \( 1 + \varepsilon \) by keeping more buckets per size. The structure partitions the ones into buckets, each summarized by its size (a power of two, the number of ones it covers) and the timestamp of its most recent one. The invariants are that sizes are nondecreasing going back in time and that each size has one or two buckets. On a new one, create a size-1 bucket. If that makes three of size 1, merge the two oldest into a size-2 bucket (keeping the newer of their right-end timestamps), and cascade the merge upward as needed. Buckets whose timestamp falls out of the window are dropped. To estimate the ones in the last \(N\) bits, sum the sizes of all buckets with timestamps in the window, but count the oldest such bucket at half size, since the window boundary may fall anywhere inside it. The error is at most half the oldest bucket's size. Because bucket sizes at most double going back and every smaller size has at least one bucket present, the oldest bucket is at most roughly equal to the sum of the newer ones, bounding the relative error by 50 percent.

Worked on a bit stream. Feed \( 1,0,1,1,0,1,1,1,0,0,1,0,1,1,1,0 \) (timestamps 1 to 16). Tracing the merges, after \(t = 4\) (a third size-1 bucket appears) the two oldest merge, giving \( \{(1, t{=}4), (2, t{=}3)\} \). After \(t = 7\) the state is \( \{(1,7), (2,6), (2,3)\} \), and after \(t = 8\) it is \( \{(1,8), (1,7), (2,6), (2,3)\} \). The one at \( t = 11 \) triggers a cascade. Three size-1 buckets merge to \( (2, 8) \), which makes three size-2 buckets, whose two oldest merge to \( (4, 6) \), leaving \( \{(1,11), (2,8), (4,6)\} \). After \(t = 15\) the state is \( \{(1,15), (1,14), (2,13), (2,8), (4,6)\} \), which satisfies both invariants. Now query at \( t = 16 \) for the ones in the last \( N = 10 \) bits (timestamps 7 to 16). The buckets with timestamp \( \ge 7 \) are \( (1,15), (1,14), (2,13), (2,8) \), while the bucket \( (4,6) \) ends before the window and is ignored. The estimate is the full sizes of the newer three, \( 1 + 1 + 2 = 4 \), plus half the oldest contributing bucket, \( 2/2 = 1 \), total 5. The true count of ones in timestamps 7 to 16 is 6, so the estimate errs by 17 percent, comfortably inside the 50 percent guarantee. The same machinery extends beyond bit counting to sums of bounded integers by decomposing into bit planes, and exponentially decayed aggregates, which need only one counter, are the limiting alternative when a sharp window edge is not required.

Link analysis, PageRank as a Markov chain

The stationary-distribution derivation

Model a web surfer as a Markov chain on the graph of pages, following a uniformly random out-link at each step. If page \(j\) has \( d_j \) out-links, the column-stochastic transition matrix has \( M_{ij} = 1/d_j \) when \( j \to i \) and 0 otherwise, and a distribution over pages evolves as \( v' = M v \). PageRank (Page, Brin, Motwani, Winograd 1999) defines a page's importance as its probability mass in the stationary distribution, the \(v\) with

that is, the principal eigenvector of \(M\) with eigenvalue 1. The recursive story it encodes is that a page is important if important pages link to it, with each page dividing its endorsement among its out-links. Perron-Frobenius theory says a unique positive stationary distribution exists, and power iteration \( v \leftarrow M v \) converges to it, provided the chain is irreducible (every page reachable from every other) and aperiodic. The real web violates both, in two characteristic ways. A dead end, a page with no out-links, makes its column of \(M\) all zeros. Each iteration then multiplies total probability by less than one, and mass leaks until \(v \to 0\). A spider trap, a set of pages linking only among themselves, is absorbing. The chain eventually enters and never leaves, so the trap accumulates all the mass and every outside page scores zero, an outcome an adversary can manufacture cheaply.

Teleportation, and power iteration with actual numbers

Both diseases are cured by teleportation. With probability \( \beta \) (typically 0.85) follow a random link, and with probability \( 1 - \beta \) jump to a page chosen uniformly at random. The iteration becomes

equivalently the eigenproblem of the dense "Google matrix" \( A = \beta M + (1-\beta)\frac{1}{n} \mathbf 1 \mathbf 1\T \), which is strictly positive, hence irreducible and aperiodic, with a unique stationary distribution. Moreover its second eigenvalue is at most \(\beta\), so power iteration converges geometrically at rate \(\beta\) (about 50 iterations for \(10^{-4}\) accuracy at \(\beta = 0.85\)). Dead ends are handled either by deleting them recursively and redistributing afterward, or, the cleaner implementation, by renormalizing. Compute \( v' = \beta M v \), then add the lost mass \( (1 - \|v'\|_1)/n \) to every entry, which folds the dead-end leak into the teleport in one line.

Computing it when the graph does not fit

At web scale, \(v\) alone is tens of gigabytes (\(10^{10}\) pages \(\times\) 8 bytes) and \(M\), stored as an edge list or adjacency lists with about \(10^{11}\) edges, is terabytes. The matrix is never materialized densely. Teleportation is a rank-one term applied analytically, and each iteration is one sequential scan of the edge list, where for each source \(j\) with degree \(d_j\) and rank \( v_j \) the scan scatters \( \beta v_j / d_j \) to each destination. If \(v'\) fits in memory, one pass suffices. If not, the block-stripe layout partitions destinations into \(k\) blocks and stores the edges as \(k\) stripes, stripe \(i\) holding, for every source, only its destinations in block \(i\) (with the source's total degree replicated into each stripe, since the scatter needs \(d_j\)). Each stripe is scanned while only block \(i\) of \(v'\) is memory-resident. The cost is reading the edge data \(k\) times per iteration, or once per stripe with the replication overhead, against the alternative of thrashing random writes across a disk-resident \(v'\), which the cost model forbids. The same computation in MapReduce or Spark is the join of the edge list with \(v\) on source, followed by a reduce on destination, and the shuffle is the scatter.

Topic-sensitive PageRank, HITS, and spam

Teleporting uniformly is a choice, and changing it changes the question asked. Topic-sensitive PageRank (Haveliwala 2002) teleports only to a curated set \(S\) of pages on a topic, \( v' = \beta M v + (1-\beta) \frac{1}{|S|} \mathbf 1_S \). The result ranks pages by importance as seen from \(S\), and a handful of topic vectors are precomputed and blended per-query. The same mathematics with \(S\) = one user's trusted pages is personalized PageRank, which reappears throughout graph ML (it is the propagation kernel in APPNP, and random-walk positional encodings descend from it). TrustRank (Gyongyi, Garcia-Molina, Pedersen 2004) is the anti-spam instance. Teleport to human-vetted seeds, so mass can only reach a page through link paths from trusted territory. A page's spam mass, the fraction of its ordinary PageRank not explainable by its TrustRank, flags link farms, which inflate ordinary PageRank by wiring thousands of pages into exactly the spider traps and reciprocal-link structures teleportation was invented to defuse.

Kleinberg's HITS (1999) is the contemporaneous alternative with a two-role model. A page is a good hub if it links to good authorities, and a good authority if linked from good hubs. With adjacency matrix \(E\) (\(E_{ij} = 1\) if \(i \to j\)), iterate \( a \leftarrow E\T h \), \( h \leftarrow E a \), normalizing each. The fixed points are the principal eigenvectors of \( E\T E \) and \( E E\T \), the top singular vectors of \(E\). HITS is query-time (run on the neighborhood of the result set) where PageRank is precomputed and query independent, and HITS's tight coupling to local link structure makes it more spam-sensitive, which is much of why the PageRank lineage won in web search while HITS's hub/authority decomposition survives in citation and e-commerce graph analysis.

Recommendation, from neighborhoods to factor models

Content-based versus collaborative

The utility matrix \(R\) has one row per user, one column per item, and known entries for a tiny fraction of cells (about 1 percent in the Netflix-Prize data). Content-based methods build an item profile from features (genres, text TF-IDF, nowadays an embedding), build a user profile as a weighted average of liked items' profiles, and recommend by profile similarity. They need no other users, work for brand-new items, and cannot surprise, since they only extrapolate a user's own history. Collaborative filtering uses only the matrix. The signal is correlation across users, and its characteristic strength (cross-genre discovery) and weakness (nothing to say about new users or items) both follow.

Neighborhood methods, with the similarity math

User-user CF predicts \( r_{ui} \) from users similar to \(u\), while item-item predicts from items similar to \(i\) that \(u\) already rated. Raw cosine on rating vectors is misleading because users differ in mean level (one user's 3 is another's 5), so ratings are mean-centered per user first. The resulting "adjusted cosine" on items, computed over co-raters, doubles as a Pearson-style correlation. Item-item won in practice (Amazon's 2003 report, Linden, Smith, York) for a statistical reason and a systems reason. An item's rating column is denser and more stable than a user's row, so item similarities are better estimated, and items change more slowly than users, so the similarity matrix can be precomputed offline and served cheaply. In a small worked computation, on the 5-user, 4-item matrix whose rows are \( (5,3,\cdot,1) \), \( (4,\cdot,\cdot,1) \), \( (1,1,\cdot,5) \), \( (1,\cdot,4,4) \), \( (\cdot,1,5,4) \) (dots are missing ratings), the user means are \( (3,\, 2.5,\, 2.33,\, 3,\, 3.33) \). After centering each row by its mean, adjusted cosine gives \( \mathrm{sim}(1, 2) = 0.555 \) over two co-raters and \( \mathrm{sim}(1, 4) = -0.898 \) over four, and the prediction for user \(u\) is the similarity-weighted average \( \hat r_{ui} = \bar r_u + \sum_{j \in N} s_{ij} (r_{uj} - \bar r_u) / \sum_j |s_{ij}| \) over the \(k\) most similar items \(u\) has rated. The co-rater counts matter. A similarity from two co-raters is noise, and production systems shrink it, e.g. multiplying by \( n_{ij}/(n_{ij} + \lambda) \) with \( \lambda \approx 25 \).

Latent factors, and the SGD update derived in full

The factor model (Koren, Bell, Volinsky 2009) posits \( R \approx P Q\T \) with user factors \( p_u \in \R^f \), item factors \( q_i \in \R^f \), plus biases, minimizing over observed entries only (treating blanks as zeros would train the model to predict zero, the classic UV-decomposition mistake),

Here \(\mu\) is the global mean, \(b_u, b_i\) capture "this user is harsh, this film is loved," which alone explain most of the variance, and the factors capture interaction structure. To derive the SGD step, for one observed \( (u, i) \), write the residual \( e_{ui} = r_{ui} - \hat r_{ui} \) with \( \hat r_{ui} = \mu + b_u + b_i + p_u\T q_i \). The per-example loss is \( \ell = e_{ui}^2 + \lambda(\|p_u\|^2 + \|q_i\|^2 + b_u^2 + b_i^2) \). Differentiate with the chain rule, using \( \partial e_{ui} / \partial p_u = -q_i \), to get

and symmetrically for \( b_i \). Descending with learning rate \(\eta\) (absorbing the 2) gives

with the caution that \(p_u\) and \(q_i\) must be updated from each other's old values. For a worked numeric step, take \( \mu = 3.5, b_u = 0.2, b_i = 0.3, p_u = (0.5, -0.3), q_i = (0.8, 0.4) \), true rating \( r = 4 \), \( \eta = 0.05, \lambda = 0.1 \). The prediction is \( 3.5 + 0.2 + 0.3 + (0.5)(0.8) + (-0.3)(0.4) = 4.28 \), so \( e = -0.28 \). The updates are \( b_u \to 0.2 + 0.05(-0.28 - 0.02) = 0.185 \), \( b_i \to 0.3 + 0.05(-0.28 - 0.03) = 0.2845 \), \( p_u \to (0.5, -0.3) + 0.05\,[(-0.28)(0.8, 0.4) - 0.1(0.5, -0.3)] = (0.4863, -0.3041) \), and \( q_i \to (0.8, 0.4) + 0.05\,[(-0.28)(0.5, -0.3) - 0.1(0.8, 0.4)] = (0.7890, 0.4022) \). Every parameter moved to shave the 0.28 overprediction, restrained by the \(\lambda\) pull toward zero.

Implicit feedback, weighted ALS and BPR

Real systems rarely have ratings. They have clicks, plays, and purchases, positive-only signals where the unobserved cells are a mix of dislike and never-seen. Hu, Koren, and Volinsky (2008) model a binary preference \( p_{ui} = \mathbf 1[r_{ui} > 0] \) with a confidence \( c_{ui} = 1 + \alpha r_{ui} \) (more plays, more confidence), and minimize \( \sum_{u,i} c_{ui} (p_{ui} - p_u\T q_i)^2 + \lambda(\cdot) \) over all cells, unobserved included at confidence 1. The sum over all cells kills SGD but suits alternating least squares. Fixing \(Q\), each \(p_u\) solves the ridge system \( (Q\T C^u Q + \lambda I)\, p_u = Q\T C^u p(u) \), and the trick that makes it tractable is \( Q\T C^u Q = Q\T Q + Q\T (C^u - I) Q \), where \( Q\T Q \) is precomputed once per sweep and the correction touches only the user's observed items. Rendle and colleagues' BPR (2009) instead optimizes ranking directly. For sampled triples (user \(u\), consumed item \(i\), unconsumed item \(j\)), maximize \( \ln \sigma(\hat x_{ui} - \hat x_{uj}) \) minus regularization, whose gradient \( (1 - \sigma(\hat x_{uij})) \) times the score gradients pushes consumed items above unconsumed ones per user. BPR's pairwise objective is the direct ancestor of the sampled softmax losses that two-tower models use today. The cold-start problem is structural for all of these. A new item has no interactions, hence no factor, and the standard escapes are content features (hybrid models mapping item features to the factor space), exploration policies that buy information with impressions, and, in current practice, initializing item embeddings from pretrained content encoders.

Evaluation, and why RMSE misleads

The Netflix Prize enshrined RMSE on held-out ratings, and the field spent a decade learning why that was the wrong target. The delivered product is a top-\(k\) list, and RMSE weights errors uniformly over the rating scale and item popularity. A model can win on RMSE by predicting mid-scale ratings of popular items precisely while misordering the handful of items that would appear in anyone's top ten. The observed entries are also missing-not-at-random (people rate what they chose to watch), so held-out-rating accuracy is estimated on a biased slice. Ranking metrics fix the target, recall@\(k\) (fraction of held-out positives appearing in the top \(k\)), MAP, and NDCG, which discounts gain by rank,

Worked, a ranking places items with relevances \( (3, 2, 0, 1) \). With linear gain, \( \mathrm{DCG} = 3/1 + 2/1.585 + 0/2 + 1/2.322 = 4.692 \), the ideal order \( (3, 2, 1, 0) \) gives \( \mathrm{IDCG} = 3 + 1.262 + 0.5 + 0 = 4.762 \), and NDCG \( = 4.692/4.762 = 0.985 \). The residual sin, common to all offline metrics, is that they score against logged behavior generated by the previous recommender, so offline gains and online A/B results routinely disagree. Ranking metrics narrow the gap, interleaving experiments and off-policy estimators narrow it further, and no offline number closes it.

Two-tower retrieval, the production successor

The architecture that actually serves recommendations at scale factorizes the problem into retrieval and ranking. The retrieval model is two towers, a user tower mapping context and history to \( u \in \R^d \) and an item tower mapping item features to \( v_i \in \R^d \), trained so that \( u \cdot v_i \) is high for observed interactions, with a sampled softmax over in-batch negatives corrected for the sampling bias (frequent items appear as negatives more often, so the logQ correction of Yi et al., RecSys 2019, subtracts \( \log q_i \) from the logit). This is matrix factorization with nonlinear towers, and the connection runs deeper. Serving it means "given \(u\), find the top items by inner product over millions of \(v_i\)," which is exactly the approximate nearest-neighbor problem, solved by the vector indexes of the modern-practice section. The heavy ranker, a gradient-boosted or deep model over hundreds of features, sees only the few hundred retrieved candidates. Covington, Adams, and Sargin's YouTube paper (2016) is the canonical description, and every large feed system since is a variation.

Clustering, dimensionality, and two more classics

The curse of dimensionality, concretely

In high dimension, distance loses contrast. For points with i.i.d. coordinates in \( [0, 1]^d \), the pairwise distance concentrates. The mean grows like \( \sqrt{d/6} \) while the standard deviation stays \(O(1)\), so the ratio of farthest to nearest neighbor tends to 1 and "nearest" stops meaning much. Angles degenerate too. Two random unit vectors have \( \E[\cos\theta] = 0 \) with \( \Var \approx 1/d \), so at \( d = 1000 \) essentially every pair is within a few degrees of orthogonal. Clustering algorithms that lean on Euclidean proximity inherit these failures, which is why the scalable classics either assume strong cluster structure (BFR) or work hard to represent shape cheaply (CURE), and why practical pipelines reduce dimension first.

BFR and CURE

BFR (Bradley, Fayyad, Reina) is k-means for a disk-resident dataset read in one pass, under the assumption that clusters are axis-aligned Gaussians. Load a memory-sized chunk. Assign points that are close enough to an existing cluster (Mahalanobis distance below a threshold like \( 2\sigma \) per dimension) to its discard set, retaining not the points but the sufficient statistics \( (N, \mathrm{SUM}, \mathrm{SUMSQ}) \) per dimension, from which mean \( \mathrm{SUM}/N \) and variance \( \mathrm{SUMSQ}/N - (\mathrm{SUM}/N)^2 \) are recoverable. Cluster the leftovers into compression sets (mini clusters, same statistics) and a small retained set of true outliers. The additivity of the statistics is the entire design, since merging two clusters is adding vectors. CURE drops the Gaussian assumption for arbitrary shapes. Each cluster is represented by \( c \approx 10 \) well-scattered sample points shrunk a fraction \(\alpha\) toward the centroid (shrinkage tames outliers), hierarchical agglomeration merges clusters by closest representative pair, and a final pass assigns every point to the cluster of its nearest representative. Representatives near a cluster's boundary let CURE capture elongated and non-convex shapes that centroid methods split.

k-means|| and coresets

k-means++ seeding (sample each next center with probability proportional to \( D^2(x) \), squared distance to the nearest chosen center) gives an \( O(\log k) \) approximation in expectation but is inherently sequential, needing \(k\) passes. Bahmani and colleagues' k-means|| (2012) parallelizes it. Each of \( O(\log \psi) \) rounds samples each point independently with probability \( \ell \, D^2(x) / \sum_y D^2(y) \) with oversampling \( \ell \approx 2k \) (one distributed pass per round, \(\psi\) the initial cost), accumulating \( O(\ell \log \psi) \) candidates, which are then weighted by the number of points they serve and clustered to exactly \(k\) centers by k-means++ on one machine. In practice about five rounds match sequential k-means++ quality. The coreset abstraction generalizes all of this, a weighted subset \(S\) such that for every candidate set of \(k\) centers, cost on \(S\) is within \( 1 \pm \varepsilon \) of cost on the full data. Sizes like \( O(k \log n / \varepsilon^2) \) suffice for k-means (line of work from Har-Peled, Feldman, Langberg, and others), and coresets compose by merge and reduce. Coresets of chunks union into a coreset of the union, then re-reduce, which turns any coreset construction into a streaming and a distributed algorithm simultaneously, the same mergeability that made sketches MapReduce-friendly.

SVD, CUR, and random projection

The SVD \( A = U \Sigma V\T \) underlies the offline dimensionality reductions. Truncating to the top \(k\) singular triplets gives \( A_k = U_k \Sigma_k V_k\T \), and the Eckart-Young theorem states this is the best rank-\(k\) approximation in both spectral and Frobenius norm. The spectral-norm argument fits in one paragraph. For any rank-\(k\) matrix \(B\), its null space has dimension at least \( n - k \), so it intersects the \( (k+1) \)-dimensional span of the top right singular vectors \( v_1, \dots, v_{k+1} \) in some unit vector \(w\). Then \( \|A - B\|_2^2 \ge \|(A - B)w\|^2 = \|Aw\|^2 = \sum_{i \le k+1} \sigma_i^2 (v_i\T w)^2 \ge \sigma_{k+1}^2 \), and \( A_k \) achieves exactly \( \sigma_{k+1} \). (The Frobenius case follows by applying this argument through the singular values one at a time, and the full proof is in Golub and Van Loan.) The practical objections to SVD at scale are that it densifies (singular vectors of a sparse matrix are dense) and that its coordinates are uninterpretable mixtures. CUR decomposition answers both. Sample actual columns \(C\) and rows \(R\) of \(A\) with probability proportional to squared norm (or leverage scores, the row norms of \(V_k\), for the relative-error guarantees of Drineas, Mahoney, Muthukrishnan), set \( U = W^+ \) the pseudoinverse of the intersection \(W\), and \( A \approx C U R \) holds with additive or relative Frobenius error depending on the sampling. The factors are sparse when \(A\) is and are literally rows and columns of the data, readable by a human.

Random projection abandons optimality for obliviousness. The Johnson-Lindenstrauss lemma says that for any \(n\) points in \( \R^d \) and \( \varepsilon \in (0, 1/2) \), a random linear map \( \Pi: \R^d \to \R^k \) with

preserves all \( \binom{n}{2} \) pairwise squared distances to within \( 1 \pm \varepsilon \) with high probability, with \(k\) independent of \(d\). The proof outline runs as follows. For a Gaussian \( \Pi \) scaled by \( 1/\sqrt{k} \), the squared length of a projected unit vector is a \( \chi^2_k / k \) variable, with \( \P\big[\, |\text{length}^2 - 1| \ge \varepsilon \,\big] \le 2\exp(-k(\varepsilon^2/2 - \varepsilon^3/3)/2) \), the Dasgupta-Gupta computation. Applying the concentration bound to each of the \( \binom{n}{2} \) difference vectors and taking a union bound gives the stated \(k\). Achlioptas showed the Gaussian can be replaced by \( \pm 1 \) entries (or a sparse \( \{-1, 0, +1\} \) with two-thirds zeros) with the same guarantee, which is what one implements. In the concrete calculation, \( n = 10^6 \) points at \( \varepsilon = 0.1 \) need \( k \ge 4 \ln(10^6) / (0.005 - 0.000333) = 4 \times 13.816 / 0.004667 = 11{,}842 \) dimensions (the simplified \( 8\ln n/\varepsilon^2 \) gives 11,052), whether \(d\) was two thousand or two billion. At \( \varepsilon = 0.2 \), \( k \approx 2{,}763 \). JL is the preprocessing step that makes the curse-of-dimensionality section's problems finite, and it is also the theory underlying the p-stable LSH projections above. Both project onto random Gaussians, LSH keeping one quantized coordinate per hash, JL keeping enough coordinates to preserve geometry outright.

Submodular maximization and the greedy 1 - 1/e

A set function \( f: 2^V \to \R \) is submodular when marginal gains diminish, meaning for all \( A \subseteq B \) and \( x \notin B \),

Coverage functions (how many elements do these sets cover), facility-location summaries (\( \sum_v \max_{x \in S} \mathrm{sim}(v, x) \), the objective behind representative-subset selection for training data), and influence spread in the independent-cascade model (Kempe, Kleinberg, Tardos 2003, who proved its submodularity) all qualify. Maximizing a monotone submodular \(f\) under \( |S| \le k \) is NP-hard, and the greedy algorithm, add the element with the largest marginal gain \(k\) times, achieves the best possible constant (Nemhauser, Wolsey, Fisher 1978),

Proof. Let \( S^\ast = \{x_1^\ast, \dots, x_k^\ast\} \) be optimal with value \( \mathrm{OPT} \), and \( S_i \) the greedy set after \(i\) steps. By monotonicity and then submodularity (adding the optimal elements to \( S_i \) one at a time, each marginal gain is no larger than that element's gain on \( S_i \) itself),

where \( \Delta_i \) is the best single-element gain available at step \(i\), which greedy takes. Rearranged, \( \Delta_i \ge (\mathrm{OPT} - f(S_i))/k \), so the gap \( g_i = \mathrm{OPT} - f(S_i) \) contracts, \( g_{i+1} \le (1 - 1/k)\, g_i \), giving \( g_k \le (1 - 1/k)^k\, \mathrm{OPT} \le e^{-1}\, \mathrm{OPT} \), i.e. \( f(S_k) \ge (1 - 1/e)\,\mathrm{OPT} \). \(\square\) Two practical notes. Lazy greedy (Minoux) exploits submodularity again. An element's gain can only shrink as \(S\) grows, so stale gains in a priority queue are upper bounds and most re-evaluations are skipped, often a hundredfold speedup with an identical result. And the guarantee is tight, since no polynomial algorithm beats \( 1 - 1/e \) for coverage unless P = NP (Feige 1998), so the greedy baseline is not a heuristic to apologize for but the optimal-in-class answer.

Online matching and AdWords

Search advertising poses matching under irrevocability. Advertisers (with bids and daily budgets) are known, queries arrive one at a time, each must be assigned immediately or dropped, and the future is unknown. Grade algorithms by the competitive ratio, the worst-case ratio of achieved to offline-optimal value. For unweighted bipartite matching, greedy (match any eligible advertiser) is 1/2-competitive, because each greedy match can block at most one optimal match. Every edge of the optimal matching missed by greedy has an endpoint that greedy matched, charging each greedy edge for at most two optimal edges, hence at least half. The bound is tight. Take two advertisers, one bidding on both query types and one on only the first. If greedy gives the first query to the flexible advertiser and only the second query type arrives thereafter, greedy earns 1 of 2. Karp, Vazirani, and Vazirani (1990) proved randomization helps. RANKING, which fixes one random priority order and matches greedily by it, achieves \( 1 - 1/e \approx 0.632 \), optimal for the problem.

AdWords adds budgets and bids (Mehta, Saberi, Vazirani, Vazirani 2007). With equal bids and budgets, the BALANCE rule, assign each query to the eligible advertiser with the most remaining budget, is the right instinct, since spending down budgets evenly keeps options open. Worked with two advertisers \( A_1 \) (bids only on \(x\)) and \( A_2 \) (bids on \(x\) and \(y\)), budgets 4 each, and the stream \( x, x, x, x, y, y, y, y \), the optimum routes all \(x\) to \( A_1 \) and all \(y\) to \( A_2 \), revenue 8. Greedy with unlucky ties sends all four \(x\) to \( A_2 \), exhausting it. The \(y\) queries then find no bidder, revenue 4, ratio 1/2. BALANCE alternates the \(x\) queries (each goes to whichever has more budget, splitting 2 and 2), then serves \(y\) from \( A_2 \)'s remaining 2, for revenue 6, ratio 3/4, which is exactly BALANCE's worst case for two advertisers, and in general BALANCE is \( 1 - 1/e \)-competitive, again optimal. For arbitrary bids, MSVV scale bids by a penalty on spent budget fraction, assigning to maximize \( \mathrm{bid} \times (1 - e^{f - 1}) \) where \(f\) is the fraction spent. The tradeoff function is derived by a factor-revealing LP and recovers \( 1 - 1/e \). The lesson that generalizes is that under adversarial arrival, hedging beats greed by exactly \( 1/2 \to 1 - 1/e \), and the same \( 1 - 1/e \) keeps appearing because both proofs run the same \( (1 - 1/k)^k \) contraction.

Worked problems

Problems 1 through 4 appear in their sections above (LSH parameter choice, A-Priori pass by pass, Bloom filter sizing, PageRank by hand). Three more, in the same spirit.

Implementation

Every sketch below was run against ground truth before the numbers in this page were written down, and the measured errors sit next to the theoretical bounds in the table after the first block. The implementations favor clarity over speed but are honest, using real hash functions (keyed BLAKE2b standing in for the pairwise-independent families the proofs assume), real bit manipulation, and no shortcuts that would change the guarantees.

import hashlib, math, random

def h64(x, seed=0):
    """64-bit keyed hash: stand-in for a random hash function."""
    d = hashlib.blake2b(repr(x).encode(), digest_size=8,
                        key=seed.to_bytes(8, "little")).digest()
    return int.from_bytes(d, "little")

class Bloom:
    def __init__(self, m, k):
        self.m, self.k = m, k
        self.bits = bytearray((m + 7) // 8)
    def _idx(self, x):                      # double hashing: h1 + i*h2
        h1, h2 = h64(x, 1), h64(x, 2) | 1
        return [(h1 + i * h2) % self.m for i in range(self.k)]
    def add(self, x):
        for i in self._idx(x):
            self.bits[i >> 3] |= 1 << (i & 7)
    def __contains__(self, x):
        return all(self.bits[i >> 3] >> (i & 7) & 1 for i in self._idx(x))

class CountMin:
    def __init__(self, eps, delta):
        self.w = math.ceil(math.e / eps)          # width: e / epsilon
        self.d = math.ceil(math.log(1 / delta))   # depth: ln(1/delta)
        self.tab = [[0] * self.w for _ in range(self.d)]
    def add(self, x, c=1):
        for j in range(self.d):
            self.tab[j][h64(x, 10 + j) % self.w] += c
    def query(self, x):                           # >= true count, w.h.p. within eps*N
        return min(self.tab[j][h64(x, 10 + j) % self.w] for j in range(self.d))

class HyperLogLog:
    def __init__(self, p=11):                     # m = 2^p registers
        self.p, self.m = p, 1 << p
        self.M = [0] * self.m
        self.alpha = 0.7213 / (1 + 1.079 / self.m)
    def add(self, x):
        h = h64(x, 99)
        j = h & (self.m - 1)                      # first p bits pick a register
        w = h >> self.p                           # remaining 64-p bits
        rho = (64 - self.p) - w.bit_length() + 1  # leading-zero rank
        self.M[j] = max(self.M[j], rho)
    def estimate(self):
        E = self.alpha * self.m ** 2 / sum(2.0 ** -r for r in self.M)
        if E <= 2.5 * self.m:                     # small-range: linear counting
            V = self.M.count(0)
            if V:
                E = self.m * math.log(self.m / V)
        return E

def reservoir(stream, k, rng=random):
    R = []
    for i, x in enumerate(stream):
        if i < k:
            R.append(x)
        else:
            j = rng.randrange(i + 1)              # uniform over 0..i
            if j < k:
                R[j] = x
    return R

def misra_gries(stream, k):
    """At most k-1 counters; f - N/k <= est <= f."""
    C = {}
    for x in stream:
        if x in C:
            C[x] += 1
        elif len(C) < k - 1:
            C[x] = 1
        else:                                     # decrement-all event
            for y in list(C):
                C[y] -= 1
                if C[y] == 0:
                    del C[y]
    return C

Measured error versus theory, from runs of exactly this code (the Zipf stream is \( \mathrm{Zipf}(1.2) \) truncated to items below \(10^4\), length \( N = 858{,}625 \)).

The pattern worth internalizing is that every measured error is inside its bound, usually far inside, because the bounds are worst-case over adversarial streams and real skewed data is kinder. Engineering practice inverts this into a sizing rule. Trust the bound for provisioning (it will not be exceeded) and expect several times better in operation.

MinHash signatures and LSH banding, the deduplication core. The banding index is a dictionary per band mapping the hashed band slice to the documents containing it. Candidates are verified with exact Jaccard, so false positives cost time, never correctness.

from collections import defaultdict

def shingles(text, k=5):
    words = text.split()
    return {" ".join(words[i:i + k]) for i in range(len(words) - k + 1)}

def signature(S, n_hashes=128):
    # column of the signature matrix: min over the set per hash function
    return [min(h64(x, seed) for x in S) for seed in range(n_hashes)]

def lsh_index(sigs, b=16, r=8):
    """sigs: {doc_id: signature}. Returns candidate pairs (b*r = len(sig))."""
    tables = [defaultdict(list) for _ in range(b)]
    for doc, sig in sigs.items():
        for band in range(b):
            key = tuple(sig[band * r:(band + 1) * r])   # r rows of this band
            tables[band][key].append(doc)
    cands = set()
    for table in tables:
        for bucket in table.values():
            for i in range(len(bucket)):
                for j in range(i + 1, len(bucket)):
                    cands.add((bucket[i], bucket[j]))
    return cands

def jaccard(A, B):
    return len(A & B) / len(A | B)

# verify candidates exactly; threshold matches (1/b)^(1/r) = 0.707
def dedup(docs, thresh=0.7):
    sets = {d: shingles(t) for d, t in docs.items()}
    sigs = {d: signature(s) for d, s in sets.items()}
    return [(a, c) for a, c in lsh_index(sigs)
            if jaccard(sets[a], sets[c]) >= thresh]

PageRank as a matrix-free power iteration over an edge list, the in-memory miniature of the block-stripe computation. Rank flows along edges via gather and scatter-add, and the teleport plus dead-end correction is applied analytically as the mass missing from the sum. Both versions converge to the worked example's \( (0.3797, 0.1989, 0.3839, 0.0375) \) on the 4-node graph. On an accelerator this loop is purely memory-bound. At the 2,992.4 GB/s fp32 copy bandwidth measured on an H100 80GB, an edge scan costing about 12 bytes of traffic per edge (two 4-byte indices plus a float gather) tops out near 250 billion edges per second per iteration, and no amount of compute optimization moves that number. Only edge compression does.

import torch

def pagerank(src, dst, n, beta=0.85, iters=50):
    # src, dst: (E,) int64 edge list
    ones = torch.ones_like(src, dtype=torch.float32)
    deg = torch.zeros(n).scatter_add_(0, src, ones)        # out-degrees (n,)
    v = torch.full((n,), 1.0 / n)
    for _ in range(iters):
        contrib = beta * v[src] / deg[src]                 # (E,) gather
        v = torch.zeros(n).scatter_add_(0, dst, contrib)   # (n,) scatter-add
        v = v + (1.0 - v.sum()) / n                        # teleport + dead ends
    return v

# 4-node worked example: A=0, B=1, C=2, D=3
src = torch.tensor([0, 0, 1, 2, 3, 3])
dst = torch.tensor([1, 2, 2, 0, 0, 2])
v = pagerank(src, dst, 4)
print(v)   # tensor([0.3797, 0.1989, 0.3839, 0.0375])

import jax
import jax.numpy as jnp

def pagerank(src, dst, n, beta=0.85, iters=50):
    ones = jnp.ones_like(src, dtype=jnp.float32)
    deg = jax.ops.segment_sum(ones, src, num_segments=n)   # out-degrees (n,)
    def step(v, _):
        contrib = beta * v[src] / deg[src]                 # (E,)
        v = jax.ops.segment_sum(contrib, dst, num_segments=n)
        v = v + (1.0 - v.sum()) / n                        # teleport + dead ends
        return v, None
    v0 = jnp.full((n,), 1.0 / n)
    v, _ = jax.lax.scan(step, v0, None, length=iters)
    return v

src = jnp.array([0, 0, 1, 2, 3, 3])
dst = jnp.array([1, 2, 2, 0, 0, 2])
print(pagerank(src, dst, 4))   # [0.3797 0.1989 0.3839 0.0375]

Matrix factorization with biases, implementing the SGD derivation exactly. The PyTorch version leans on autograd and expresses the \(\lambda\) penalty as weight decay (note it then also decays \(\mu\), and excluding it via parameter groups is the production nicety), while the JAX version writes the loss explicitly so the gradient matches the derived updates term for term.

import torch
import torch.nn as nn

class MF(nn.Module):
    def __init__(self, n_users, n_items, f=32):
        super().__init__()
        self.P = nn.Embedding(n_users, f)     # user factors (n_users, f)
        self.Q = nn.Embedding(n_items, f)     # item factors (n_items, f)
        self.bu = nn.Embedding(n_users, 1)
        self.bi = nn.Embedding(n_items, 1)
        self.mu = nn.Parameter(torch.tensor(0.0))
        nn.init.normal_(self.P.weight, std=0.1)
        nn.init.normal_(self.Q.weight, std=0.1)
        nn.init.zeros_(self.bu.weight)
        nn.init.zeros_(self.bi.weight)

    def forward(self, u, i):                  # u, i: (B,) int64
        dot = (self.P(u) * self.Q(i)).sum(-1)             # (B,)
        return self.mu + self.bu(u).squeeze(-1) + self.bi(i).squeeze(-1) + dot

model = MF(n_users=1000, n_items=1700)
opt = torch.optim.SGD(model.parameters(), lr=0.05, weight_decay=1e-4)
for u, i, r in batches:                       # observed triples only
    loss = ((model(u, i) - r) ** 2).mean()
    opt.zero_grad(); loss.backward(); opt.step()

import jax
import jax.numpy as jnp

def init(key, n_users, n_items, f=32):
    k1, k2 = jax.random.split(key)
    return dict(P=0.1 * jax.random.normal(k1, (n_users, f)),
                Q=0.1 * jax.random.normal(k2, (n_items, f)),
                bu=jnp.zeros(n_users), bi=jnp.zeros(n_items),
                mu=jnp.array(3.5))

def loss_fn(params, u, i, r, lam=1e-4):
    pred = (params["mu"] + params["bu"][u] + params["bi"][i]
            + (params["P"][u] * params["Q"][i]).sum(-1))  # (B,)
    reg = sum((v ** 2).sum() for k, v in params.items() if k != "mu")
    return ((pred - r) ** 2).mean() + lam * reg

@jax.jit
def step(params, u, i, r, eta=0.05):
    grads = jax.grad(loss_fn)(params, u, i, r)
    # identical to the derived updates: p += eta*(e*q - lam*p), etc.
    return jax.tree.map(lambda p, g: p - eta * g, params, grads)

Random projection, with the JL guarantee checked empirically by projecting and comparing pairwise squared distances before and after. Dense Gaussian projection is a single matmul, the cheapest operation an accelerator has. The same H100 measured 744.6 TFLOPS on bf16 matmuls at \( n = 4096 \), so projecting a million 512-dimensional vectors to \( k = 2763 \) dimensions (the JL size for \( \varepsilon = 0.2 \)) is roughly \( 2 \times 10^6 \times 512 \times 2763 = 2.8 \times 10^{12} \) FLOPs, under ten milliseconds of tensor-core time. The bottleneck is, as always, moving the vectors.

import math
import torch

n, d, eps = 2000, 10_000, 0.2
k = int(8 * math.log(n) / eps ** 2)          # JL dimension: 1520 for n=2000
X = torch.randn(n, d)
Pi = torch.randn(d, k) / math.sqrt(k)        # E[|Pi w|^2] = |w|^2
Y = X @ Pi                                   # (n, k)

pairs = torch.randint(0, n, (5000, 2))
a, b = pairs[:, 0], pairs[:, 1]
a, b = a[a != b], b[a != b]                  # drop degenerate self-pairs
before = ((X[a] - X[b]) ** 2).sum(-1)
after = ((Y[a] - Y[b]) ** 2).sum(-1)
ratio = after / before
print(ratio.min().item(), ratio.max().item())  # inside (1-eps, 1+eps) w.h.p.

import math
import jax
import jax.numpy as jnp

n, d, eps = 2000, 10_000, 0.2
k = int(8 * math.log(n) / eps ** 2)
key = jax.random.PRNGKey(0)
kx, kp, ki = jax.random.split(key, 3)
X = jax.random.normal(kx, (n, d))
Pi = jax.random.normal(kp, (d, k)) / math.sqrt(k)
Y = X @ Pi

pairs = jax.random.randint(ki, (5000, 2), 0, n)
a, b = pairs[:, 0], pairs[:, 1]
a, b = a[a != b], b[a != b]                  # drop degenerate self-pairs
before = ((X[a] - X[b]) ** 2).sum(-1)
after = ((Y[a] - Y[b]) ** 2).sum(-1)
print((after / before).min(), (after / before).max())

Finally, the analytics framing is frequent pairs as SQL. The self-join on basket id is pass 2 of A-Priori, the HAVING on items is the \(L_1\) filter, and a database's hash join plus pre-aggregation is PCY and the combiner wearing different clothes. The Python tab is the same computation as an explicit two-pass program.

-- baskets(basket_id, item): one row per item occurrence, s = 3
WITH freq_items AS (                      -- pass 1: frequent singletons
  SELECT item
  FROM baskets
  GROUP BY item
  HAVING COUNT(*) >= 3
)
SELECT a.item AS item1, b.item AS item2, COUNT(*) AS support
FROM baskets a
JOIN baskets b
  ON a.basket_id = b.basket_id
 AND a.item < b.item                      -- each unordered pair once
WHERE a.item IN (SELECT item FROM freq_items)   -- monotonicity pruning
  AND b.item IN (SELECT item FROM freq_items)
GROUP BY a.item, b.item
HAVING COUNT(*) >= 3                      -- pass 2: frequent pairs
ORDER BY support DESC;

from collections import Counter
from itertools import combinations

def apriori_pairs(baskets, s=3):
    item_counts = Counter(i for b in baskets for i in b)      # pass 1
    L1 = {i for i, c in item_counts.items() if c >= s}
    pair_counts = Counter()
    for b in baskets:                                          # pass 2
        for pair in combinations(sorted(set(b) & L1), 2):
            pair_counts[pair] += 1
    return {p: c for p, c in pair_counts.items() if c >= s}

baskets = [{"m","b","r"}, {"m","b"}, {"b","r","d"}, {"m","d"},
           {"m","b","d","r"}, {"b","d"}, {"m","b","d"}, {"r","d"}]
print(apriori_pairs(baskets))
# {('b','m'): 4, ('b','r'): 3, ('b','d'): 4, ('d','r'): 3, ('d','m'): 3}

How it is done in practice

Vector databases, how similarity search is actually served

The production descendants of LSH are approximate nearest-neighbor indexes, and the two designs that dominate are graphs and quantized inverted files. HNSW (Malkov and Yashunin 2018) builds a hierarchy of proximity graphs. Every vector is a node in the bottom layer with links to \( M \approx 16\text{-}48 \) near neighbors chosen by a diversity heuristic, and each higher layer keeps an exponentially thinning random subset with longer-range links, a structure deliberately analogous to a skip list. A query enters at the sparse top, greedily walks to the nearest node layer by layer, and at the bottom runs a best-first beam search of width ef. Recall is dialed at query time by ef alone. IVF-PQ (Jegou, Douze, Schmid 2011, the core of FAISS) instead clusters the corpus into nlist cells by k-means, searches only the nprobe nearest cells, and compresses residual vectors by product quantization. Split into \(m\) subvectors, quantize each against its own 256-entry codebook, store \(m\) bytes per vector, and score queries against codes through \( m \times 256 \) lookup tables (asymmetric distance) without ever decompressing. The arithmetic that sells it is that \(10^9\) vectors at 768 dimensions of fp32 is 3 TB raw, while PQ at 64 bytes per vector is 64 GB, one machine instead of a cluster.

The reason graphs largely beat LSH in this role is that LSH partitions space without looking at the data, so to reach high recall it must union many independent tables, each probe of which is a hash lookup returning candidates that still need exact scoring. On real corpora with low intrinsic dimension, a navigable graph adapts to the data manifold and reaches the same recall with far fewer distance computations, typically an order of magnitude fewer at recall 0.95 and above on the standard million-to-billion-scale benchmarks. LSH retains three roles the graphs cannot fill. It has worst-case sublinear guarantees independent of data distribution (Indyk-Motwani, and Andoni-Indyk's optimal families), its signatures are mergeable streaming objects, and for set similarity at extreme scale (deduplication) banding examines only colliding pairs, a batch capability no per-query graph index offers. Hybrids own the frontier. DiskANN (Microsoft Research, 2019) pushes graph search to SSD-resident billion-vector corpora with PQ in memory and full vectors on disk, while ScaNN (Google, 2020) sharpened PQ training with an anisotropic loss that weights errors along the query direction, and its successors serve retrieval for production LLM stacks.

Deduplicating LLM pretraining corpora

The highest-stakes deployment of MinHash today is cleaning web-scale text before pretraining. Lee and colleagues' "Deduplicating Training Data Makes Language Models Better" (ACL 2022) established the case. Near-duplicate training text wastes compute, inflates evaluation via train-test leakage, and increases verbatim memorization, and removing it improves perplexity at fixed compute. Their pipeline pairs exact substring deduplication (suffix arrays over the corpus) with MinHash-LSH over n-gram shingles for near-duplicates, and that recipe, with tuned parameters, is now standard. RefinedWeb, SlimPajama (13-gram shingles, Jaccard threshold 0.8), Dolma, and FineWeb (5-gram shingles, 112 hashes in 14 bands of 8 rows, an S-curve threshold of \( (1/14)^{1/8} \approx 0.72 \)) all ship MinHash stages, and the FineWeb report is unusually candid about the second-order effects. Deduplicating too aggressively across snapshots removed high-quality recurring documents and hurt downstream accuracy, so the final pipeline dedups within snapshots but not across them. The numbers involved are the point. FineWeb processed on the order of \(10^{13}\) tokens across tens of billions of documents, a regime where the only algorithms that run at all are the ones in this page, hash, band, bucket, and never compare all pairs.

Streams, logs, and feature stores

The streaming model became infrastructure as the log. Kafka durably stores an append-only partitioned log. Producers write, consumer groups read at their own pace, and the log's replay-ability converts "one pass" from a constraint into a choice. Flink executes the actual streaming computation, with event-time windows and watermarks (the disciplined answer to out-of-order arrival), keyed state per operator, and exactly-once state via distributed snapshots, aligned barriers flowing through the dataflow in the Chandy-Lamport style, checkpointing every operator's state consistently so failure recovery replays from the last snapshot. Inside such pipelines the sketches of this page run as operator state, per-key Count-Min for rate estimation, HyperLogLog for windowed distinct users, and Space-Saving for trending items. The same sketches ship inside storage engines. Redis's PFCOUNT is a 16,384-register, 6-bit-packed HyperLogLog (12 KB, 0.81 percent standard error, exactly \( 1.04/\sqrt{2^{14}} \)), and BigQuery's APPROX_COUNT_DISTINCT is HyperLogLog++. The Apache DataSketches library (from Yahoo) provides mergeable theta and quantile sketches for warehouse rollups. LevelDB and RocksDB consult per-SSTable Bloom filters before touching disk, the purest cost-model play in the list, 10 bits per key in memory to avoid a 100-microsecond read that would return nothing. Feature stores (Feast and its commercial kin) exist because recommendation features are computed twice, batch for training and streaming for serving, and the two paths drift. The store's contract is that both read the same definitions, with the streaming path built from exactly the windowed-aggregation machinery above.

The current research frontier

Three active fronts descend directly from this material. The first is learned augmentation of classical structures. Hsu, Indyk, Katabi, and Vakilian (ICLR 2019, at MIT) attached a learned oracle to Count-Min and Count sketches. A model predicts which items are heavy and gives them dedicated counters, provably shrinking error on power-law streams while the sketch still catches the oracle's mistakes. The same pattern (model proposes, classical structure guarantees) runs through learned Bloom filters and the learned-index line started by Kraska and colleagues at MIT and Google, with follow-up analysis from Mitzenmacher at Harvard establishing when the learned variants actually save space. The second front is billion-scale ANN, where the competition is explicitly cross-institutional. Microsoft's DiskANN line (SSD-resident graphs, with streaming updates in FreshDiskANN) competes against Google's ScaNN quantization line, NVIDIA's GPU-native CAGRA graphs in cuVS, Zilliz's Milvus engine, and academic entrants like NTU's RaBitQ, which put quantization back on a rigorous footing with distance estimates carrying probabilistic error bounds, and the public BigANN benchmarks referee. Filtered and hybrid search (vector similarity constrained by metadata predicates) is where the open problems now concentrate, since neither graphs nor inverted files handle selective filters gracefully. The third front is data curation for pretraining at corpus scale. Beyond MinHash near-dedup, Meta researchers proposed SemDeDup (embedding-space semantic deduplication) and D4, showing duplicate-aware data selection shifts scaling curves. Hugging Face's FineWeb and AI2's Dolma published full ablations of dedup strategies, turning what was folklore into measured pipeline science, and exact-substring dedup via suffix automata over trillions of tokens is its own systems problem. Beneath all three, the sketch literature itself stays lively where streams meet new constraints, with mergeable summaries as a composable algebra (Agarwal, Cormode, and others), sketches under differential privacy (Apple's deployed private Count-Mean Sketch for telemetry), and sketches in programmable network hardware (the UnivMon line from CMU and Johns Hopkins, which put universal frequency-moment estimation into switches).

Open source to read

Each of these repositories embodies one section of this page, and the named file is the right first door.

• apache/spark, the RDD/DAG model in the flesh. Open core/src/main/scala/org/apache/spark/rdd/RDD.scala, where the lineage-and-partitions abstraction, with narrow versus wide dependencies, is the whole NSDI paper in one class.
• facebookresearch/faiss, the reference ANN library. Open faiss/IndexIVFPQ.cpp to see coarse quantization, residual product quantization, and the ADC lookup-table scan in one file.
• nmslib/hnswlib, HNSW in a single header. Open hnswlib/hnswalg.h, where layer assignment, the neighbor-diversity heuristic, and the ef-bounded beam search are all short enough to actually read.
• apache/flink, production stream processing. Open flink-streaming-java/src/main/java/org/apache/flink/streaming/api/datastream/DataStream.java for the operator algebra, then follow windowing and checkpoint barriers outward.
• google/leveldb, the storage-engine angle on this page. Open util/bloom.cc, a production Bloom filter in under a hundred lines, including the double-hashing trick and the bits-per-key arithmetic derived above.
• redis/redis, where src/hyperloglog.c is one of the best-commented implementations anywhere. The header explains the 6-bit register packing, the sparse-to-dense transition, and the bias-corrected estimator.
• ekzhu/datasketch, the Python MinHash/LSH toolkit used by several corpus dedup pipelines. Open datasketch/minhash.py, then lsh.py for the banding index.
• apache/kafka, the log as infrastructure. Open clients/src/main/java/org/apache/kafka/clients/producer/KafkaProducer.java to see partitioning, batching, and acknowledgment, the practical face of "the shuffle is the cost."
• milvus-io/milvus, a full vector database around the ANN cores. Start in internal/proxy/task_search.go and follow a search request through routing, segment fan-out, and result reduction.

Common misconceptions

"Bloom filters can also miss items that were inserted." They cannot, because insertion sets bits and nothing ever clears them, so every inserted item finds all its bits set forever. False negatives appear only in variants that delete (counting filters with counter underflow, or aging schemes), which is precisely why plain Bloom filters do not support deletion.

"The Count-Min sketch has small relative error." Its guarantee is additive, \( \hat f \le f + \varepsilon N \), with \(N\) the whole stream. For an item with \( f \approx \varepsilon N \) the relative error can be 100 percent or more. The sketch is an instrument for heavy hitters. Reading light-item frequencies off it is a category error, and the Count sketch's \( \sqrt{F_2} \) scale is the honest upgrade when the tail matters.

"Vector databases made LSH irrelevant." Graph indexes won the online k-NN query workload, but LSH banding solves a different problem, finding all similar pairs in a batch corpus without a query loop. Every major LLM pretraining corpus was deduplicated with MinHash-LSH in the last three years, and the technique is arguably at its historical peak of deployed importance.

"PageRank is basically in-degree with extra steps." The worked example refutes this in four nodes. Rank depends on who links, recursively, not how many. A single link from a hub outweighs many links from leaf pages, and teleportation gives every page a floor of \( (1-\beta)/n \) that in-degree counting has no analog for. The distinction is exactly what link farms exploit and TrustRank patches.

"MapReduce is obsolete, so studying its cost model is archaeology." The system is obsolete, but the model prices every shuffle in Spark, every distributed join in a SQL warehouse, and every all-to-all in distributed training. The replication-rate lower bound derived above is a statement about the problem, not about Hadoop, and it holds for whatever executes the join next decade.

"A random sample can answer distinct-count questions." Sampling is the wrong tool for cardinality. A 1 percent uniform sample cannot distinguish \(10^6\) items each appearing once from \(10^4\) items each appearing 100 times without seeing item identities across the whole stream. Distinct counting needs hashing precisely because duplicates must collide deterministically, and HyperLogLog is small because it hashes everything and samples nothing.

"Lower RMSE means better recommendations." RMSE is measured on held-out ratings of items users chose to rate, weighted uniformly, while the product is a ranked top-k slate over items they have not seen. A model can improve RMSE by calibrating mid-scale predictions on popular items while degrading the ordering at the top of the list, and missing-not-at-random feedback biases the metric further. Ranking metrics on held-out interactions, and ultimately online experiments, are the ground truth.

"HyperLogLog registers can be combined by averaging the estimates." Merging two HLLs means taking the registerwise maximum, which yields exactly the sketch of the union, as if one sketch had seen both streams. Averaging two estimates double-counts the intersection. The max-merge is why HLLs federate across shards and time windows losslessly, and also why they cannot subtract (there is no registerwise operation for set difference, and intersections come indirectly via inclusion-exclusion, with error amplification).

Self-check

References

1. Leskovec, J., Rajaraman, A., Ullman, J. D. Mining of Massive Datasets, 3rd ed., 2020. mmds.org
2. Muthukrishnan, S. Data Streams: Algorithms and Applications. Foundations and Trends in Theoretical Computer Science, 2005. doi:10.1561/0400000002
3. Dean, J., Ghemawat, S. MapReduce: Simplified Data Processing on Large Clusters. OSDI 2004 and CACM 51(1), 2008. doi:10.1145/1327452.1327492
4. Zaharia, M., Chowdhury, M., Das, T., Dave, A., Ma, J., McCauley, M., Franklin, M. J., Shenker, S., Stoica, I. Resilient Distributed Datasets: A Fault-Tolerant Abstraction for In-Memory Cluster Computing. NSDI 2012. usenix.org
5. Afrati, F. N., Ullman, J. D. Optimizing Joins in a Map-Reduce Environment. EDBT 2010. doi:10.1145/1739041.1739056
6. Broder, A. Z. On the Resemblance and Containment of Documents. SEQUENCES 1997. doi:10.1109/SEQUEN.1997.666900
7. Indyk, P., Motwani, R. Approximate Nearest Neighbors: Towards Removing the Curse of Dimensionality. STOC 1998. doi:10.1145/276698.276876
8. Charikar, M. Similarity Estimation Techniques from Rounding Algorithms. STOC 2002. doi:10.1145/509907.509965
9. Datar, M., Immorlica, N., Indyk, P., Mirrokni, V. Locality-Sensitive Hashing Scheme Based on p-Stable Distributions. SoCG 2004. doi:10.1145/997817.997857
10. Andoni, A., Indyk, P. Near-Optimal Hashing Algorithms for Approximate Nearest Neighbor in High Dimensions. CACM 51(1), 2008. doi:10.1145/1327452.1327494
11. Cormode, G., Muthukrishnan, S. An Improved Data Stream Summary: The Count-Min Sketch and its Applications. Journal of Algorithms 55(1), 2005. doi:10.1016/j.jalgor.2003.12.001
12. Alon, N., Matias, Y., Szegedy, M. The Space Complexity of Approximating the Frequency Moments. JCSS 58(1), 1999 (STOC 1996). doi:10.1006/jcss.1997.1545
13. Flajolet, P., Fusy, E., Gandouet, O., Meunier, F. HyperLogLog: The Analysis of a Near-Optimal Cardinality Estimation Algorithm. AofA 2007. dmtcs.episciences.org/3545
14. Misra, J., Gries, D. Finding Repeated Elements. Science of Computer Programming 2, 1982. doi:10.1016/0167-6423(82)90012-0
15. Datar, M., Gionis, A., Indyk, P., Motwani, R. Maintaining Stream Statistics over Sliding Windows. SIAM Journal on Computing 31(6), 2002. doi:10.1137/S0097539701398363
16. Brin, S., Page, L. The Anatomy of a Large-Scale Hypertextual Web Search Engine. Computer Networks 30, 1998. doi:10.1016/S0169-7552(98)00110-X
17. Kleinberg, J. M. Authoritative Sources in a Hyperlinked Environment. JACM 46(5), 1999. doi:10.1145/324133.324140
18. Koren, Y., Bell, R., Volinsky, C. Matrix Factorization Techniques for Recommender Systems. IEEE Computer 42(8), 2009. doi:10.1109/MC.2009.263
19. Hu, Y., Koren, Y., Volinsky, C. Collaborative Filtering for Implicit Feedback Datasets. ICDM 2008. doi:10.1109/ICDM.2008.22
20. Rendle, S., Freudenthaler, C., Gantner, Z., Schmidt-Thieme, L. BPR: Bayesian Personalized Ranking from Implicit Feedback. UAI 2009. arXiv:1205.2618
21. Nemhauser, G. L., Wolsey, L. A., Fisher, M. L. An Analysis of Approximations for Maximizing Submodular Set Functions I. Mathematical Programming 14, 1978. doi:10.1007/BF01588971
22. Mehta, A., Saberi, A., Vazirani, U., Vazirani, V. AdWords and Generalized Online Matching. JACM 54(5), 2007. doi:10.1145/1284320.1284321
23. Malkov, Y. A., Yashunin, D. A. Efficient and Robust Approximate Nearest Neighbor Search Using Hierarchical Navigable Small World Graphs. IEEE TPAMI 42(4), 2020. arXiv:1603.09320
24. Jegou, H., Douze, M., Schmid, C. Product Quantization for Nearest Neighbor Search. IEEE TPAMI 33(1), 2011. doi:10.1109/TPAMI.2010.57
25. Lee, K., Ippolito, D., Nystrom, A., Zhang, C., Eck, D., Callison-Burch, C., Carlini, N. Deduplicating Training Data Makes Language Models Better. ACL 2022. arXiv:2107.06499

Key takeaway