lsh_banding.py

1"""2LSH Banding：碰撞概率公式、S 型曲线与候选集（第18章 §18.3.3）34正文的推导（签名长度 n = b·r）：56    单个 Band 匹配        P = t^r7    b 个 Band 全不匹配    P = (1 − t^r)^b8    至少一个 Band 匹配    P(candidate) = 1 − (1 − t^r)^b      ← 经典公式910其中 t = P(h(x) = h(y))：MinHash 时是 Jaccard，SimHash 时是 1 − θ/π。11拐点（碰撞概率达到 1/2 的位置）满足12    t_0^r = 1 − 2^(−1/b)   ⟹   t_0 ≈ (ln2 / b)^(1/r) ≈ (1/b)^(1/r)1314本脚本分四段，每一段都对应正文里的一个论断：151) 曲线本身：固定 (r, b) 打印 P(t)，并验证「r 增大曲线右移、b 增大左移」；162) 用随机 bit 实测碰撞率，确认 1−(1−t^r)^b 不是近似而是精确概率；173) 拐点：精确解 1−2^(−1/b) 与工程近似 (1/b)^(1/r) 的差距；184) 真实检索：MinHash 签名 + Banding 建桶，扫 (b, r)，19   统计 Recall@10 / 候选数 / 耗时，展示 b 与 r 的取舍。2021运行：python3 ch03_lsh/lsh_banding.py22"""2324from __future__ import annotations2526import pathlib27import sys28import time2930import numpy as np3132sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent.parent))33sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))3435from common import table  # noqa: E40236from minhash import jaccard, make_sets  # noqa: E4023738SEED = 202639P = 2 ** 31 - 1404142def collide_prob(t: float, r: int, b: int) -> float:43    """正文的碰撞概率公式 1 − (1 − t^r)^b。"""44    return 1.0 - (1.0 - t ** r) ** b454647def hinge_exact(b: int, r: int) -> float:48    """P(t_0)=1/2 的精确解：t_0^r = 1 − 2^(−1/b)。"""49    return float((1.0 - 2.0 ** (-1.0 / b)) ** (1.0 / r))505152def hinge_approx(b: int, r: int) -> float:53    """正文的工程近似 t_0 ≈ (1/b)^(1/r)。"""54    return float((1.0 / b) ** (1.0 / r))555657# --------------------------------------------------------------------------58# MinHash 签名（线性哈希，与 §18.3.1 的实现一致）59# --------------------------------------------------------------------------60class MinHashSigner:61    def __init__(self, num_hash: int, seed: int = SEED) -> None:62        rng = np.random.default_rng(seed)63        self.L = num_hash64        self.a = rng.integers(1, P, size=num_hash, dtype=np.int64)65        self.b = rng.integers(0, P, size=num_hash, dtype=np.int64)6667    def sign(self, s: set) -> np.ndarray:68        items = np.fromiter(s, dtype=np.int64, count=len(s))69        h = (self.a[:, None] * items[None, :] + self.b[:, None]) % P70        return h.min(axis=1)717273def band_keys(sig: np.ndarray, r: int) -> list[tuple]:74    """把签名切成 b 个连续 Band，每个 Band 的 r 位作为倒排表的 Key。"""75    out = []76    for s in range(0, len(sig) - r + 1, r):77        out.append(tuple(int(x) for x in sig[s:s + r]))78    return out798081def main() -> None:82    print("LSH Banding：碰撞概率公式、S 型曲线与候选集（§18.3.3）\n")8384    # ---------------- 1) S 型曲线 ----------------85    print("① 碰撞概率 P(t) = 1 − (1 − t^r)^b：")86    ts = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]87    rows = []88    for t in ts:89        rows.append([f"{t:.1f}"] + [f"{collide_prob(t, r, 16):.3f}"90                                    for r in (2, 4, 8)])91    table(["t", "r=2, b=16", "r=4, b=16", "r=8, b=16"], rows)92    print("\n  固定 b=16，只增大 r：同一条 t 上的碰撞概率单调下降，曲线整体右移。")93    print("  含义是「筛选更严格」——只有更像的对象才能进候选集。")94    rows = []95    for t in ts:96        rows.append([f"{t:.1f}"] + [f"{collide_prob(t, 4, b):.3f}"97                                    for b in (4, 16, 64)])98    table(["t", "r=4, b=4", "r=4, b=16", "r=4, b=64"], rows)99    print("\n  固定 r=4，只增大 b：碰撞概率单调上升，曲线整体左移。")100    print("  b 是「独立抽签次数」，次数越多，只要中一次就进候选，所以召回更高、")101    print("  候选也更多。两者合起来构成一条可控的 S 型判别曲线。\n")102103    # ---------------- 2) 实测碰撞率 ----------------104    print("② 实测：从随机 bit 签名统计碰撞率，验证公式是精确恒等式")105    print("   （每条 2000 对、签名长度 n=256，t 位相同的概率由独立翻转控制）")106    n_bits = 256107    n_pairs = 2000108    rng = np.random.default_rng(SEED)109    base = rng.integers(0, 2, size=(n_pairs, n_bits)).astype(np.int8)110    rows = []111    for t in (0.3, 0.5, 0.7, 0.9):112        flips = rng.random((n_pairs, n_bits)) > t       # 以 1−t 的概率翻转113        other = base ^ flips.astype(np.int8)114        for r, b in ((4, 16), (8, 16)):115            hit = np.zeros(n_pairs, dtype=bool)116            for s in range(0, r * b, r):117                same_band = (base[:, s:s + r] == other[:, s:s + r]).all(axis=1)118                hit |= same_band119            rows.append([f"{t:.1f}", r, b, f"{hit.mean():.4f}",120                         f"{collide_prob(t, r, b):.4f}"])121    table(["t", "r", "b", "实测碰撞率", "公式值"], rows)122    print("\n  逐行相等（误差在 2000 次抽样的 ±0.01 量级）。")123    print("  注意公式里唯一的假设是「各位独立」；MinHash 的 L 个哈希如果用的是")124    print("  同一组随机数，独立性会打折扣，工程上通常直接用 L 个独立种子。\n")125126    # ---------------- 3) 拐点 ----------------127    print("③ 拐点位置：精确解 vs 正文的 (1/b)^(1/r) 近似")128    rows = []129    for b, r in ((16, 4), (16, 8), (32, 4), (64, 4), (64, 8), (256, 4)):130        te = hinge_exact(b, r)131        ta = hinge_approx(b, r)132        rows.append([b, r, f"{te:.4f}", f"{ta:.4f}",133                     f"{(ta - te) / te:+.1%}",134                     f"1−2^(−1/{b}) = {1 - 2 ** (-1 / b):.5f}"])135    table(["b", "r", "精确 t_0", "近似 t_0", "相对误差", "t_0^r 的精确值"], rows)136    print("\n  近似把 t_0^r 从 1−2^(−1/b) 换成了 1/b；b ≥ 16 时两者的差不到 4%，")137    print("  b 越大越准（1−2^(−1/b) ≈ ln2/b）。所以正文说它「较准确」，")138    print("  做参数初筛足够，做严格的召回率承诺时应该用精确解。\n")139140    # ---------------- 4) 真实检索：MinHash + Banding ----------------141    print("④ 真实检索：2000 个集合、40 个簇（同簇 J≈0.47）、Top-10")142    print("   （签名长度 L 固定，只改 Band 的切法；Recall@10 用真值 Jaccard 排序做基准）")143    sets, labels = make_sets(n_sets=2000, universe=20000, size=40,144                             n_groups=40, keep=0.8)145146    # 全量精确 Jaccard，作为真值147    truth = []148    for i in range(len(sets)):149        js = []150        for j in range(len(sets)):151            if i == j:152                js.append(-1.0)153                continue154            inter = len(sets[i] & sets[j])155            js.append(inter / (len(sets[i]) + len(sets[j]) - inter))156        truth.append(np.argsort(-np.asarray(js))[:10])157158    L_FIX = 512159    signer = MinHashSigner(L_FIX, seed=SEED)160    sigs = np.stack([signer.sign(s) for s in sets])      # [2000, 512]161162    rows = []163    for r in (2, 4, 8, 16):164        b = L_FIX // r165        t0 = time.perf_counter()166        tables = [dict() for _ in range(b)]167        for doc in range(len(sets)):168            for bi, key in enumerate(band_keys(sigs[doc], r)):169                tables[bi].setdefault(key, []).append(doc)170        t_build = time.perf_counter() - t0171172        rec = 0.0173        cand_tot = 0174        t0_ = time.perf_counter()175        for q in range(len(sets)):176            cand: set = set()177            for bi, key in enumerate(band_keys(sigs[q], r)):178                cand.update(tables[bi].get(key, ()))179            cand.discard(q)180            cand_arr = np.fromiter(cand, dtype=np.int64, count=len(cand))181            if len(cand_arr) == 0:182                continue183            # 候选内用签名估计相似度排序（真实系统里这一步是「粗排」）184            sim = (sigs[cand_arr] == sigs[q][None, :]).mean(axis=1)185            top = cand_arr[np.argsort(-sim)[:10]]186            rec += len(set(top.tolist()) & set(truth[q].tolist())) / 10187            cand_tot += len(cand_arr)188        dt = (time.perf_counter() - t0_) / len(sets)189        p_theory = collide_prob(0.47, r, b)190        rows.append([r, b, f"{p_theory:.3f}", f"{rec / len(sets):.3f}",191                     f"{cand_tot / len(sets):.0f}",192                     f"{t_build * 1e3:.0f} ms", f"{dt * 1e3:.2f} ms"])193    table(["r", "b", "理论碰撞率(t=0.47)", "Recall@10", "平均候选数",194           "建表耗时", "单查询"], rows)195    print("\n  三列对照着读：r 从 2 增到 16，理论碰撞率按 1−(1−t^r)^b 剧烈下降，")196    print("  实测 Recall 同步下降，候选数也从 61 条收敛到 0 条。")197    print("  r=2、r=4 两档都一样，是因为 61 和 50 个候选里已经包含了几乎全部")198    print("  同簇对象，差别只剩粗排能否排进 Top-10——这正是「候选数到瓶颈后，")199    print("  召回由精排决定」的现象。")200    print("  r=8 时理论碰撞率只剩 0.14，候选数掉到 8 条，召回随之掉到 0.26：")201    print("  Banding 一旦成了瓶颈，后面的粗排再准也救不回来。\n")202203    print("   固定 r=2、只加长签名（同时让 b = L/2），看召回的天花板在哪里：")204    rows = []205    for L in (64, 128, 256, 512, 1024):206        mh = MinHashSigner(L, seed=SEED)207        sg = np.stack([mh.sign(s) for s in sets])208        b2 = L // 2209        tabs = [dict() for _ in range(b2)]210        for doc in range(len(sets)):211            for bi in range(b2):212                tabs[bi].setdefault(tuple(int(x) for x in sg[doc, 2 * bi:2 * bi + 2]), []).append(doc)213        rec = 0.0214        cand_tot = 0215        for q in range(len(sets)):216            cand: set = set()217            for bi in range(b2):218                cand.update(tabs[bi].get(tuple(int(x) for x in sg[q, 2 * bi:2 * bi + 2]), ()))219            cand.discard(q)220            ca = np.fromiter(cand, dtype=np.int64, count=len(cand))221            if len(ca) == 0:222                continue223            sim = (sg[ca] == sg[q][None, :]).mean(axis=1)224            top = ca[np.argsort(-sim, kind="stable")[:10]]225            rec += len(set(top.tolist()) & set(truth[q].tolist())) / 10226            cand_tot += len(ca)227        rows.append([L, b2, f"{rec / len(sets):.3f}", f"{cand_tot / len(sets):.1f}"])228    table(["签名长度 L", "b（r=2）", "Recall@10", "平均候选数"], rows)229    print("\n  候选数几乎不变（50 → 62），召回却从 0.45 一路升到 0.76：")230    print("  瓶颈已经不在 Banding，而在「用签名估计相似度来排序」这一步——")231    print("  L=512 时同簇对象两两估计值的标准差约 0.02，几十个对象挤在同一个")232    print("  分数上，Top-10 里谁进谁不进基本靠运气（L=64 时更糟）。")233    print("  工业系统的做法是在候选内回原始表算精确相似度，而不是继续加长签名。\n")234235    print("   同一份签名只改 b（固定 r=2，签名长度 512 不变）：")236    rows = []237    for b in (8, 16, 32, 64, 128, 256):238        tables = [dict() for _ in range(b)]239        for doc in range(len(sets)):240            for bi in range(b):241                key = tuple(int(x) for x in sigs[doc, 2 * bi:2 * bi + 2])242                tables[bi].setdefault(key, []).append(doc)243        rec = 0.0244        cand_tot = 0245        for q in range(len(sets)):246            cand: set = set()247            for bi in range(b):248                key = tuple(int(x) for x in sigs[q, 2 * bi:2 * bi + 2])249                cand.update(tables[bi].get(key, ()))250            cand.discard(q)251            cand_arr = np.fromiter(cand, dtype=np.int64, count=len(cand))252            if len(cand_arr) == 0:253                continue254            sim = (sigs[cand_arr] == sigs[q][None, :]).mean(axis=1)255            top = cand_arr[np.argsort(-sim)[:10]]256            rec += len(set(top.tolist()) & set(truth[q].tolist())) / 10257            cand_tot += len(cand_arr)258        rows.append([b, f"{collide_prob(0.47, 2, b):.4f}",259                     f"{rec / len(sets):.3f}", f"{cand_tot / len(sets):.0f}"])260    table(["b（r=2）", "理论碰撞率", "Recall@10", "平均候选数"], rows)261    print("\n  b=8 时碰撞率 ≈ 0.86，约有 14% 的近邻根本进不了候选，召回卡在 0.67；")262    print("  b 增大后碰撞率趋于 1，召回升到 0.70 附近就不再涨——剩下的缺口来自")263    print("  签名排序的分辨率，不是候选覆盖。候选数从 44 缓升到 61。")264    print("  所以「加 b」是用候选数换召回，「加 r」是用召回换精度，二者要联调。\n")265266    # 对照：随机采样同样候选数能做到多少267    print("   对照：随机抽相同数量的候选（200 条查询，签名粗排 Top-10）：")268    rows = []269    rng2 = np.random.default_rng(SEED)270    for m in (50, 200, 800, 2000):271        rec = 0.0272        for q in range(200):273            cand = rng2.choice(2000, m, replace=False)274            sim = (sigs[cand] == sigs[q][None, :]).mean(axis=1)275            top = cand[np.argsort(-sim)[:10]]276            rec += len(set(top.tolist()) & set(truth[q].tolist())) / 10277        rows.append([f"{m}", f"{rec / 200:.3f}"])278    table(["随机候选数", "Recall@10"], rows)279    print("\n  同样几百个候选，随机采样只有 10%–38% 的召回，而 Banding 的同桶候选")280    print("  在同等规模下能到 70%——LSH 的价值不在于「少扫」，")281    print("  而在于「扫的那一小撮是对的」。\n")282283    print("⑤ 一句话：P(candidate)=1−(1−t^r)^b 把签名变成判别器，")284    print("   r 管每个 Band 的严格程度，b 管独立机会数，拐点在 t_0≈(1/b)^(1/r)。")285286287if __name__ == "__main__":288    main()