simhash.py

1"""2SimHash：随机超平面、谷歌 64 位指纹与分块倒排检索（第18章 §18.3.2）34正文的三层内容，本脚本逐层验证：561) 随机超平面 LSH：h_r(x) = sign(xᵀr)，r ~ N(0, I)，夹角 θ 的两条向量满足7       P[h(x) = h(y)] = 1 − θ/π8   长度 L 的签名里，汉明距离 H 满足 H/L ≈ θ/π，于是 θ ≈ π·H/L。92) 谷歌 SimHash：每个特征先做普通哈希得到 L 位编码 b_ij，转成10       r_ij = 2b_ij − 1 ∈ {−1, +1}11   再按权重累加 s = Σ w_i r_i，取 sign(s)。正文证明了它与随机超平面12   LSH 完全等价——第 j 位就是 r_jᵀx，只是 r_j 由哈希"隐式"生成，13   不需要存储 L×d 的投影矩阵。143) 谷歌的检索方案：64 位指纹切成 m = k+1 块建倒排索引。若两条指纹的15   汉明距离 ≤ k，鸽笼原理保证至少有一块完全相同，所以只查 m 条倒排链16   就能召回全部近重复，不用扫描全库。1718本脚本：191) 实测 P(same) 与 1−θ/π、H/L 与 θ/π 两条等式（随机超平面）；202) 实现谷歌 SimHash（splitmix64 + 带权累加），与显式随机投影逐位比对，21   并验证"汉明距离反推夹角"的精度；223) 20 万条 64 位指纹上对比暴力扫描与 4×16 bit 分块倒排：候选数、23   召回率、耗时，并演示 m < k+1 时鸽笼原理失效、必然漏召回；244) 说明 Faiss IndexLSH 为什么只是"压缩"而不是子线性检索。2526运行：python3 ch03_lsh/simhash.py27"""2829from __future__ import annotations3031import pathlib32import sys33import time3435import numpy as np3637sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent.parent))3839from common import table  # noqa: E4024041SEED = 202642BITS = 64          # 谷歌 SimHash 的指纹长度43N_DB = 200_000     # 检索演示的库规模444546# --------------------------------------------------------------------------47# 一、64 位哈希（splitmix64）：模拟工业界的"普通哈希函数"48# --------------------------------------------------------------------------49GOLD = np.uint64(0x9E3779B97F4A7C15)50MUL1 = np.uint64(0xBF58476D1CE4E5B9)51MUL2 = np.uint64(0x94D049BB133111EB)525354def splitmix64(x: np.ndarray) -> np.ndarray:55    """56    splitmix64：把整数打散成 64 位伪随机数。5758    选它的原因：一行位运算就能在 numpy 里向量化，且雪崩性质好，59    可以当作"普通哈希函数"的替身。正文强调谷歌方案里每个特征的编码60    来自普通哈希而非随机数生成器，这里正是按这个思路实现的。61    """62    with np.errstate(over="ignore"):          # uint64 溢出是预期行为63        z = x + GOLD64        z = (z ^ (z >> np.uint64(30))) * MUL165        z = (z ^ (z >> np.uint64(27))) * MUL266        return z ^ (z >> np.uint64(31))676869def hash_bits(fids: np.ndarray, bits: int = BITS) -> np.ndarray:70    """特征编号 → [m, bits] 的 0/1 编码 b_ij。"""71    h = splitmix64(fids.astype(np.uint64))72    shifts = np.arange(bits, dtype=np.uint64)73    return ((h[:, None] >> shifts[None, :]) & np.uint64(1)).astype(np.int8)747576def google_simhash(fids: np.ndarray, weights: np.ndarray,77                   bits: int = BITS) -> tuple[int, np.ndarray]:78    """79    正文的谷歌 SimHash：8081        b_ij ∈ {0,1}  →  r_ij = 2b_ij − 182        s_j = Σ_i w_i·r_ij83        指纹第 j 位 = 1 if s_j ≥ 0 else 08485    返回 (指纹整数, 64 维累计权重 s)。86    """87    r = 2 * hash_bits(fids, bits).astype(np.float64) - 1.0     # [m, L]88    s = weights @ r                                            # [L]89    bits_arr = (s >= 0).astype(np.uint8)90    fp = int(sum(int(b) << j for j, b in enumerate(bits_arr)))91    return fp, s929394def popcount(x: np.ndarray) -> np.ndarray:95    """对 uint64 数组做 popcount（numpy 没有内置位计数）。"""96    u8 = x.view(np.uint8).reshape(-1, 8)97    return np.unpackbits(u8, axis=1).sum(axis=1)9899100def hamming(fp: np.ndarray, q: int) -> np.ndarray:101    return popcount(fp ^ np.uint64(q))102103104def build_block_tables(db_fp: np.ndarray, m_blocks: int) -> list[dict]:105    """106    谷歌方案的倒排表：64 位指纹切成 m 块，每块作为 Key 建一张表。107108    真实系统里这四张表是离线建好的；这里抽成函数，方便对比不同 m。109    """110    width = BITS // m_blocks111    mask = (1 << width) - 1112    tables: list[dict] = [{} for _ in range(m_blocks)]113    for doc, fpv in enumerate(db_fp):114        v = int(fpv)115        for b in range(m_blocks):116            tables[b].setdefault((v >> (b * width)) & mask, []).append(doc)117    return tables118119120def query_block_tables(q_int: int, tables: list[dict],121                       m_blocks: int) -> set:122    """把查询指纹按同样的切法查 m 条倒排链，取并集。"""123    width = BITS // m_blocks124    mask = (1 << width) - 1125    cand: set = set()126    for b in range(m_blocks):127        cand.update(tables[b].get((q_int >> (b * width)) & mask, ()))128    return cand129130131# --------------------------------------------------------------------------132def batch_simhash(fids: np.ndarray, weights: np.ndarray,133                  bits: int = BITS, chunk: int = 2000) -> np.ndarray:134    """135    批量版谷歌 SimHash：fids / weights 形状 [n, m]，返回 [n] 的 uint64 指纹。136137    逐条调用会做 n 次 Python 循环，20 万条太慢；这里分块向量化：138    [chunk, m, L] 的符号矩阵与权重做 einsum，再按位打包成整数。139    分块是为了控制中间张量的内存（2000×80×64 的 float32 约 41 MB）。140    """141    n = len(fids)142    out = np.empty(n, dtype=np.uint64)143    powers = (np.uint64(1) << np.arange(bits, dtype=np.uint64))144    for s in range(0, n, chunk):145        e = min(s + chunk, n)146        B = hash_bits(fids[s:e].ravel(), bits).reshape(e - s, fids.shape[1], bits)147        R = 2.0 * B.astype(np.float32) - 1.0148        S = np.einsum("nml,nm->nl", R, weights[s:e].astype(np.float32))149        bits_mat = (S >= 0).astype(np.uint64)150        out[s:e] = (bits_mat * powers[None, :]).sum(axis=1)151    return out152153154def rotate(x: np.ndarray, theta: float, u: np.ndarray) -> np.ndarray:155    """把单位向量 x 在 span{x,u} 平面内旋转 θ 弧度。"""156    return np.cos(theta) * x + np.sin(theta) * u157158159def main() -> None:160    print("SimHash：随机超平面、谷歌 64 位指纹与分块倒排检索（§18.3.2）\n")161    rng = np.random.default_rng(SEED)162163    # ---------------- 1) 随机超平面的两条等式 ----------------164    print("① 随机超平面：P(同侧) = 1 − θ/π，H/L ≈ θ/π（D=64，2 万个随机方向）")165    D = 64166    x = rng.normal(size=D)167    x /= np.linalg.norm(x)168    u = rng.normal(size=D)169    u -= (u @ x) * x                 # 与 x 正交170    u /= np.linalg.norm(u)171172    n_proj = 20_000173    R = rng.normal(size=(n_proj, D))174    rows = []175    for frac in (0.0, 1 / 6, 1 / 3, 0.5, 2 / 3, 5 / 6, 1.0):176        theta = np.pi * frac177        y = rotate(x, theta, u)178        hx = (R @ x) >= 0179        hy = (R @ y) >= 0180        emp = float((hx == hy).mean())181        rows.append([f"{frac:.2f}π",182                     f"{emp:.4f}", f"{1 - frac:.4f}",183                     f"{1 - emp:.4f}", f"{frac:.4f}"])184    table(["夹角 θ", "实测 P(同侧)", "1−θ/π", "实测 H/L", "θ/π"], rows)185    print("\n  两列比值几乎逐行相等：每一位就是一个 Bernoulli(1−θ/π)，")186    print("  L 位独立，汉明距离 H 的期望是 L·θ/π，std 是 sqrt(L·θ/π·(1−θ/π))。")187    print("  注意 θ=0 与 θ=π 两行是退化情形：两个向量共线，符号永远相同或永远相反。\n")188189    # ---------------- 2) 谷歌 SimHash 与显式随机投影 ----------------190    print("② 谷歌 SimHash 与显式随机超平面是否逐位一致：")191    m_feat = 200192    fids = np.arange(m_feat)193    weights = rng.uniform(0.2, 3.0, size=m_feat)194195    fp, s = google_simhash(fids, weights)196    bits_google = np.array([(fp >> j) & 1 for j in range(BITS)], dtype=np.uint8)197198    # 显式版本：哈希位矩阵转置后，第 j 列就是随机投影向量 r_j，199    # 那么 s_j = r_jᵀx（x 为特征权重向量）——与谷歌版应逐位相同200    Bm = hash_bits(fids, BITS)                       # [m, L]201    Rj = (2.0 * Bm - 1.0).T                          # [L, m]202    s_explicit = Rj @ weights203    bits_explicit = (s_explicit >= 0).astype(np.uint8)204205    agree = float((bits_google == bits_explicit).mean())206    print(f"  特征数 m = {m_feat}，指纹长度 L = {BITS}，随机权重 w_i ~ U(0.2, 3.0)")207    print(f"  两种实现的逐位一致率 = {agree:.4f}")208    print(f"  最大绝对差 max|s_google − s_explicit| = {np.abs(s - s_explicit).max():.2e}")209    print("  正文的等价性推导在这里是恒等式而非近似：谷歌版算的就是 r_jᵀx，")210    print("  只是 r_j 的第 i 个分量 r_ij 由「特征 i 的哈希第 j 位」给出。")211    print("  工程收益是存储：不必保存 L×m 的投影矩阵，只需一个 64 位哈希函数。\n")212213    # 夹角估计精度214    print("  用汉明距离反推夹角（两条 200 维文档向量，人为旋转到不同夹角）：")215    rows = []216    for frac in (0.1, 0.2, 0.3, 0.4, 0.5):217        theta = np.pi * frac218        # 构造 y：在特征权重空间里把 x 旋转 θ，再用同一套哈希指纹比较219        v = rng.normal(size=m_feat)220        v -= (v @ weights) / (weights @ weights) * weights221        v *= np.linalg.norm(weights) / np.linalg.norm(v)222        wy = np.cos(theta) * weights + np.sin(theta) * v223        fp_x, _ = google_simhash(fids, weights)224        fp_y, _ = google_simhash(fids, wy)225        h = float(popcount(np.array([np.uint64(fp_x ^ fp_y)]))[0])226        rows.append([f"{frac:.1f}π", f"{h:.0f} / {BITS}",227                     f"{np.pi * h / BITS:.3f}", f"{theta:.3f}"])228    table(["真实 θ", "汉明距离 H", "π·H/L 估计", "真实 θ"], rows)229    print("\n  L=64 时每位的标准差约 0.5/√64 ≈ 0.06（弧度约 0.2），")230    print("  所以角度估计有可见误差——这也是工业上 SimHash 只用来做粗筛，")231    print("  候选还要回到原始向量或原始文档做精确比对的原因。\n")232233    # ---------------- 3) 分块倒排检索 ----------------234    print("③ 谷歌的分块倒排检索：暴力扫指纹 vs m = k+1 块查表")235    print(f"   （{N_DB:,} 条 64 位指纹，用随机文档生成）")236237    # 库：随机稀疏文档的 SimHash 指纹；再构造"近重复查询"238    n_docs = N_DB239    vocab = 5000240    kmax = 120241    n_feat = rng.integers(40, kmax + 1, size=n_docs)242    # 用矩形 [n, kmax] 存稀疏特征：不足 kmax 的位置填 fid=0、w=0，243    # 对结果没有影响（权重为 0），但能把整批计算放进 numpy244    feat_ids = np.zeros((n_docs, kmax), dtype=np.int64)245    feat_w = np.zeros((n_docs, kmax), dtype=np.float32)246    for i in range(n_docs):247        k = int(n_feat[i])248        feat_ids[i, :k] = rng.choice(vocab, k, replace=False)249        feat_w[i, :k] = rng.uniform(0.2, 3.0, size=k)250    db_fp = batch_simhash(feat_ids, feat_w)251252    # 查询：以第 0 条为基准，翻转 k 个随机 bit，保证真实汉明距离正好是 k253    q_fp = int(db_fp[0])254    n_q = 20255256    def make_query(k_flip: int, t: int, m_blocks: int | None = None) -> int:257        """258        构造与库中第 0 条指纹汉明距离恰为 k_flip 的查询。259260        m_blocks 给定时，若 m_blocks ≤ k_flip，就强制「每个块至少被翻一位」——261        这正是鸽笼原理失效的最坏情形。不做这个强制的话，随机翻位有262        ~1 − m!/m^k 的概率把两位翻进同一块，恰好留下一个完整块，263        于是 m = k 也会偶然召回，看不出问题（第一版就踩了这个坑）。264        """265        rngq = np.random.default_rng(SEED + 100 * k_flip + t)266        bits_q = [(q_fp >> j) & 1 for j in range(BITS)]267        if k_flip:268            if m_blocks is not None and m_blocks <= k_flip:269                width = BITS // m_blocks270                chosen = [int(b * width + rngq.integers(width))271                          for b in range(m_blocks)]272                rest = k_flip - m_blocks273                if rest > 0:274                    pool = [j for j in range(BITS) if j not in chosen]275                    chosen += rngq.choice(pool, rest, replace=False).tolist()276            else:277                chosen = rngq.choice(BITS, k_flip, replace=False).tolist()278            for j in chosen:279                bits_q[j] ^= 1280        return sum(b << j for j, b in enumerate(bits_q))281282    rows = []283    for k_flip in (0, 1, 2, 3, 4):284        m_blocks = k_flip + 1285        tables = build_block_tables(db_fp, m_blocks)      # 离线建表，不计入查询耗时286        hit = 0287        cand_tot = 0288        t_brute = 0.0289        t_ivf = 0.0290        for t in range(n_q):291            q_int = make_query(k_flip, t)292            t0 = time.perf_counter()293            h_all = hamming(db_fp, q_int)294            truth = np.where(h_all <= k_flip)[0]295            t_brute += time.perf_counter() - t0296297            t0 = time.perf_counter()298            cand = query_block_tables(q_int, tables, m_blocks)299            # 候选上再算精确汉明距离（正文的最后一步）300            cand_arr = np.fromiter(cand, dtype=np.int64, count=len(cand))301            _ = hamming(db_fp[cand_arr], q_int)302            t_ivf += time.perf_counter() - t0303304            hit += len(set(truth.tolist()) & cand) / max(len(truth), 1)305            cand_tot += len(cand)306        rows.append([k_flip, m_blocks, f"{hit / n_q:.4f}",307                     f"{cand_tot / n_q:,.0f}",308                     f"{t_brute / n_q * 1e3:.2f} ms",309                     f"{t_ivf / n_q * 1e3:.2f} ms"])310    table(["k", "块数 m", "近重复召回率", "平均候选数",311           "暴力扫描", "查表+重排"], rows)312    print("\n  建表是离线的，这里只计查询：m 次字典查找 + 候选上的精确汉明距离。")313    print("  候选数远小于 20 万——16 bit 的块有 6.5 万个可能取值，随机指纹几乎")314    print("  不会撞桶，所以查询常只带回 1 条真实近重复；k 增大时块变短，桶才变挤。")315    print("  这就是正文说的「只访问 4 条倒排链」，计算量降低数万倍。")316    print("  召回率始终保持 1.000，因为 m = k+1 时鸽笼原理保证至少一块相同。\n")317318    print("  反例：块数与 k 不匹配时，鸽笼原理失效：")319    rows = []320    for k_flip in (1, 2, 3, 4):321        for m_blocks in (1, 2, 4, 5):322            tables = build_block_tables(db_fp, m_blocks)323            hit = 0324            for t in range(n_q):325                q_int = make_query(k_flip, t, m_blocks)326                truth = set(np.where(hamming(db_fp, q_int) <= k_flip)[0].tolist())327                cand = query_block_tables(q_int, tables, m_blocks)328                hit += len(truth & cand) / max(len(truth), 1)329            rows.append([k_flip, m_blocks, "✓" if m_blocks > k_flip else "✗",330                         f"{hit / n_q:.3f}"])331    table(["k", "块数 m", "m > k?", "近重复召回率"], rows)332    print("\n  只有 m ≥ k+1 的行召回是 1.000；m ≤ k 时每个块都可能被改到，")333    print("  必然漏掉近重复（这里对 m ≤ k 强制每块各翻一位，即最坏情形）。")334    print("  代价是块越短桶越挤：m 增大 → 候选数上升，这正是 §18.3.3 Banding")335    print("  要用 b、r 两个参数平衡的东西。\n")336337    # ---------------- 4) IndexLSH：压缩，不是子线性 ----------------338    print("④ Faiss IndexLSH 的本质：全量汉明比较的压缩检索")339    print("   （2000 条 32 维簇状向量，组内夹角小；随机高维向量两两近正交，")340    print("    SimHash 在这种数据上本来就分不开，不适合演示压缩的代价）")341    N = 2000342    dim = 32343    centers = rng.normal(size=(20, dim)) * 3.0344    labels = rng.integers(0, 20, size=N)345    Xv = (centers[labels] + rng.normal(0, 0.35, (N, dim))).astype(np.float32)346    Xv /= np.linalg.norm(Xv, axis=1, keepdims=True)347    pick = rng.choice(N, 100, replace=False)348    Qv = Xv[pick] + rng.normal(0, 0.01, (100, dim)).astype(np.float32)349    Qv /= np.linalg.norm(Qv, axis=1, keepdims=True)350351    k = 10352    rows = []353    for L in (32, 64, 128):354        proj = rng.normal(size=(L, dim)).astype(np.float32)355        code = (Xv @ proj.T >= 0)356        qcode = (Qv @ proj.T >= 0)357        rec = 0.0358        for i in range(len(Qv)):359            # IndexLSH 的默认查询：对全部编码算汉明距离，取最近的 k 个候选360            hd = (code != qcode[i]).sum(axis=1)361            cand = np.argsort(hd)[:k]362            true = np.argsort(-(Xv @ Qv[i]))[:k]363            rec += len(set(cand.tolist()) & set(true.tolist())) / k364        rows.append([L, f"{rec / len(Qv):.3f}",365                     f"{N * L / 8 / 1024:.1f} KB", "全量 2000 次汉明比较"])366    table(["签名长度 L", "Recall@10", "索引大小", "查询方式"], rows)367    print("\n  召回随 L 上升，但查询始终是 O(N)——索引只是把 32 维 float32")368    print("  （128 B）压成 L 位（L/8 B），省的是内存和带宽，不是扫描量。")369    print("  真正子线性的是 IndexBinaryMultiHash / 谷歌分块倒排：多张哈希表，")370    print("  只访问命中的桶。二者的区别正是正文「压缩编码 vs 子线性检索」那句话。\n")371372    print("⑤ 一句话：SimHash 用 sign(xᵀr) 把夹角编码成 bit，H/L ≈ θ/π；")373    print("   谷歌版用普通哈希隐式生成 r，再用 k+1 块倒排把「找汉明距离 ≤ k」")374    print("   变成几次字典查找。前者是编码，后者才是索引。")375376377if __name__ == "__main__":378    main()