Locality-sensitive hashing

Locality-sensitive hashing (LSH) is the simplest family of approximate nearest-neighbor methods: design a hash function whose collision probability is monotone in the similarity between two points[1]. Similar points land in the same bucket with high probability; dissimilar points land in different buckets. Lookup is then just a bucket fetch.

The defining property

A family $\mathcal{H}$ of hash functions is $(r_1, r_2, p_1, p_2)$-locality sensitive for a distance $d$ if:

• For $d(x, y) \le r_1$: $\Pr_{h \sim \mathcal{H}}(h(x) = h(y)) \ge p_1$.
• For $d(x, y) \ge r_2$: $\Pr_{h \sim \mathcal{H}}(h(x) = h(y)) \le p_2$.

With $r_1 < r_2$ and $p_1 > p_2$, the gap between $p_1$ and $p_2$ can be amplified by combining several hashes (next section).

MinHash for Jaccard similarity

For set-valued data (documents as sets of $k$-shingles, users as sets of liked items), Jaccard similarity is $J(A, B) = |A \cap B|/|A \cup B|$.

MinHash picks a random permutation $\pi$ of the universe and defines $h_\pi(A) = \min_{a \in A}\pi(a)$. The key fact:

The collision probability equals the Jaccard similarity, by symmetry of permutations over union.

In practice, instead of true permutations, use $k$ universal hash functions and take the minimum hash value of items in $A$ under each. A document is summarized by a length-$k$ MinHash signature; estimated Jaccard is the fraction of matching positions.

SimHash for cosine similarity

For vector-valued data with cosine similarity, SimHash uses random hyperplanes:

A standard result:

where $\theta$ is the angle between $x$ and $y$. Smaller angle (more similar) means higher collision probability. With $k$ random hyperplanes, each point gets a $k$-bit signature; Hamming distance between signatures estimates angular distance.

$p$-stable LSH for $\ell_2$

For Euclidean distance, use random projections plus quantization:

with $w$ from a $p$-stable distribution (Gaussian for $\ell_2$, Cauchy for $\ell_1$), $b$ uniform in $[0, r]$, and bucket width $r$. Points within roughly distance $r$ of each other collide with high probability; distant points do not.

Amplification: AND-OR construction

A single hash has weak collision guarantees. Combine multiple:

• AND: concatenate $k$ hashes into one composite. Two points collide only if all $k$ component hashes agree. The collision probability of similar points becomes $p_1^k$; dissimilar becomes $p_2^k$. The ratio $(p_1/p_2)^k$ grows.
• OR: build $L$ hash tables with independent composites. Two points collide if they agree in any table. The collision probability becomes $1 - (1 - p^k)^L$[2].

Tuning $k$ and $L$ trades recall against precision and against the size of the candidate set retrieved.

Limits

• Memory: $L$ tables of $n$ entries each, plus the original vectors.
• Recall ceiling: high-recall regimes (>95%) require expensive AND-OR combinations; graph methods (HNSW) tend to dominate at high recall.
• Tuning: $k$, $L$, and $r$ all need calibration on the actual data distribution.

LSH is conceptually elegant and memory-light at moderate recall, which made it the workhorse before graph methods matured.

Check it on arrays

import numpy as np

rng = np.random.default_rng(0)

# Cosine LSH (SimHash): k random hyperplanes per signature
n, d = 5000, 64
X = rng.normal(size=(n, d))
X /= np.linalg.norm(X, axis=1, keepdims=True)         # unit norm
queries = rng.normal(size=(50, d)); queries /= np.linalg.norm(queries, axis=1, keepdims=True)

def simhash(X, normals):
    return (X @ normals.T > 0).astype(np.uint8)

def signatures_match(sig_a, sig_b):
    return (sig_a == sig_b).sum(axis=-1)

# Single table, vary k (signature length)
def recall_for_k(k_bits, table_count):
    recall_total = 0.0; frac_total = 0.0; count = 0
    for query in queries:
        true_top10 = np.argsort(-X @ query)[:10]      # by cosine similarity
        # OR over `table_count` tables
        candidate_mask = np.zeros(n, dtype=bool)
        for _ in range(table_count):
            normals = rng.normal(size=(k_bits, d))
            H = simhash(X, normals)
            hq = simhash(query[None, :], normals)[0]
            # Match if signature agrees on ALL bits
            same = (H == hq).all(axis=1)
            candidate_mask |= same
        frac_total += candidate_mask.sum() / n         # average over EVERY query
        if candidate_mask.sum() < 10: continue
        cand_idx = np.where(candidate_mask)[0]
        ann_top10 = cand_idx[np.argsort(-(X[cand_idx] @ query))[:10]]
        recall_total += len(set(ann_top10) & set(true_top10)) / 10
        count += 1
    return recall_total / max(count, 1), frac_total / len(queries)

print(f"{'k_bits':>8}  {'L tables':>8}  {'recall@10':>10}  {'cand. fraction':>15}")
for k_bits in [4, 8, 12]:
    for L in [1, 4, 16]:
        r, frac = recall_for_k(k_bits, L)
        print(f"{k_bits:>8}  {L:>8}  {r:>10.3f}  {frac:>15.4f}")

Longer signatures cut the candidate set; more tables (OR amplification) boost recall.

Exercises

Exercise 1. A SimHash family has collision probability $p_1 = 0.9$ for a similar pair and $p_2 = 0.4$ for a dissimilar pair. Form an AND composite of $k = 5$ bits, then an OR over $L = 10$ such composites. Compute the final collision probability for each pair and explain why the gap widened.

Exercise 2. Two documents have MinHash signatures of length $k = 200$ that agree in $150$ positions. Estimate their Jaccard similarity and the standard error of that estimate. Why does a longer signature give a tighter estimate?

Exercise 3 (code). Verify the SimHash collision identity $\Pr_w(h_w(x) = h_w(y)) = 1 - \theta(x, y)/\pi$. Build pairs at several known angles $\theta$ in a high-dimensional space, draw many random hyperplanes, and check that the empirical collision rate matches $1 - \theta/\pi$.

import numpy as np

rng = np.random.default_rng(0)
d = 64

def unit(v):
    return v / np.linalg.norm(v)

# Build an orthonormal pair (x, y0) so cos(theta) x + sin(theta) y0 has angle theta to x
x = unit(rng.normal(size=d))
y0 = rng.normal(size=d)
y0 = unit(y0 - (y0 @ x) * x)

normals = rng.normal(size=(200_000, d))      # random hyperplanes

for theta in [np.pi/6, np.pi/3, np.pi/2, 2*np.pi/3]:
    y = np.cos(theta) * x + np.sin(theta) * y0
    hx = (normals @ x) > 0
    hy = (normals @ y) > 0
    collide = (hx == hy).mean()
    predicted = 1 - theta / np.pi
    print(f"theta={theta:.3f}  empirical={collide:.3f}  predicted={predicted:.3f}")
    assert abs(collide - predicted) < 0.01

The next lesson covers an inverted-index approach (IVF) that scales to billions of vectors.