minhash.py

1"""2MinHash：用随机排列 / 线性哈希估计 Jaccard 相似度（第18章 §18.3.1）34正文的核心性质：5    h_π(C) = min_{r∈C} π(r)6    P[h_π(A) = h_π(B)] = |A∩B| / |A∪B| = J(A,B)78用 L 个独立哈希得到签名，则9    Ĵ = (1/L) Σ 1[h_i(A) = h_i(B)] → J(A,B)10    std(Ĵ) = sqrt(J(1−J)/L)1112正文给了两种实现：13    1) 随机排列（permutation）：最直观、完全符合定义，但需要保存 L 张14       长度 d 的置换表，内存 O(L·d)；15    2) 线性哈希：h(x) = (a·x + b) mod p，p 取大素数（2^31−1）。16       每个哈希函数只存 (a, b) 两个整数，内存 O(L)，工程上等价。1718本脚本：191) 在小集合上逐位对比两种实现，确认签名一致、碰撞概率约等于 J；202) 扫签名长度 L，画「估计误差 vs sqrt(J(1−J)/L)」，验证正文的方差公式；213) 展示为什么工程上用线性哈希：内存、构建耗时、精度三者对比；224) 接一个「用 MinHash 找相似集合」的小检索，说明它是 LSH 的基础。2324运行：python3 ch03_lsh/minhash.py25"""2627from __future__ import annotations2829import pathlib30import sys3132import numpy as np3334sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent.parent))3536from common import table  # noqa: E4023738P = 2 ** 31 - 1           # 正文用的大素数39SEED = 2026404142# --------------------------------------------------------------------------43# 一、基于随机排列的实现（完全对应定义）44# --------------------------------------------------------------------------45class MinHashPermutation:46    """47    预生成 L 套独立随机排列 π_1..π_L，签名第 i 位取集合内元素在 π_i 下的最小 rank。4849    这里保存的是「元素 → rank」的反向映射而不是排列本身，因为查询只需要 rank。50    注意内存：L × universe_size 个整数，universe 是 10^6 时 L=128 就要 5 亿项。51    """5253    def __init__(self, num_hash: int, universe_size: int, seed: int = 42) -> None:54        self.num_hash = num_hash55        self.universe_size = universe_size56        rng = np.random.default_rng(seed)57        self.ranks = []58        for _ in range(num_hash):59            perm = rng.permutation(universe_size)60            rank = np.empty(universe_size, dtype=np.int32)61            rank[perm] = np.arange(universe_size, dtype=np.int32)62            self.ranks.append(rank)6364    def get_signature(self, s: set) -> np.ndarray:65        items = np.fromiter(s, dtype=np.int64, count=len(s))66        sig = np.empty(self.num_hash, dtype=np.int64)67        for i, rank in enumerate(self.ranks):68            sig[i] = int(rank[items].min())69        return sig7071    def nbytes(self) -> int:72        return self.num_hash * self.universe_size * 4737475# --------------------------------------------------------------------------76# 二、基于线性哈希的实现（工程版）77# --------------------------------------------------------------------------78class MinHashLinear:79    """80    h_i(x) = (a_i·x + b_i) mod p。8182    当 p 足够大（远大于最大元素值）时，(a·x+b) mod p 对 x 而言是近似均匀的83    「随机重排 + 压缩」，取最小值得到的碰撞概率仍等于 J。数学上它不再是严格84    的排列（可能把两个元素映到同一个值），但 p ≫ d 时碰撞概率可忽略。85    """8687    def __init__(self, num_hash: int, seed: int = 42) -> None:88        rng = np.random.default_rng(seed)89        self.num_hash = num_hash90        self.p = P91        self.a = rng.integers(1, P, size=num_hash, dtype=np.int64)92        self.b = rng.integers(0, P, size=num_hash, dtype=np.int64)9394    def get_signature(self, s: set) -> np.ndarray:95        items = np.fromiter(s, dtype=np.int64, count=len(s))96        # [L, |s|]：一次性算完所有哈希，再沿集合维度取最小97        h = (self.a[:, None] * items[None, :] + self.b[:, None]) % self.p98        return h.min(axis=1)99100    def nbytes(self) -> int:101        return self.num_hash * 16      # 两个 int64102103104def jaccard(A: set, B: set) -> float:105    return len(A & B) / len(A | B)106107108def make_pair(target_j: float, size: int = 60,109              universe: int = 4096, seed: int = SEED) -> tuple[set, set]:110    """111    构造一对指定 Jaccard 的集合。112113    A 取 size 个随机元素；B 保留其中 m 个、再补 size−m 个新元素，于是114        J = m / (2·size − m)   ⟹   m = 2·size·J / (1 + J)115    """116    rng = np.random.default_rng(seed)117    pool = rng.choice(universe, size * 2, replace=False)118    A = set(pool[:size].tolist())119    m = int(round(2 * size * target_j / (1 + target_j)))120    B = set(pool[:m].tolist()) | set(pool[size:size + (size - m)].tolist())121    return A, B122123124def make_sets(n_sets: int, universe: int, size: int, n_groups: int,125              keep: float = 0.6, seed: int = SEED) -> tuple[list[set], list[int]]:126    """127    生成有簇结构的集合：同一 group 的集合共享一部分基础元素，其余各自随机。128129    keep 控制共享比例：同组两个集合各自从 base 里独立抽 keep 比例，130    交集约 keep²·size，并集约 size + (1−keep)·size，因此131        同组 J ≈ keep² / (2 − keep)      （keep=0.6 时 ≈ 0.26）。132    完全随机的集合两两 J ≈ size/universe 且几乎一样，看不出 MinHash 的区分能力。133    """134    rng = np.random.default_rng(seed)135    bases = [set(rng.choice(universe, size, replace=False).tolist())136             for _ in range(n_groups)]137    sets, labels = [], []138    for _ in range(n_sets):139        g = int(rng.integers(n_groups))140        base = np.array(sorted(bases[g]))141        n_keep = max(1, int(round(size * keep)))142        shared = set(rng.choice(base, n_keep, replace=False).tolist())143        noise = set(rng.choice(universe, size - n_keep, replace=False).tolist())144        sets.append(shared | noise)145        labels.append(g)146    return sets, labels147148149def main() -> None:150    print("MinHash：用随机排列 / 线性哈希估计 Jaccard 相似度（§18.3.1）\n")151    rng = np.random.default_rng(SEED)152153    # ---------------- 1) 两种实现的一致性 ----------------154    print("① 两种签名实现的一致性（小 universe，便于做完整排列）：")155    universe = 4096156    sets, _ = make_sets(n_sets=20, universe=universe, size=60, n_groups=4)157    L = 64158    mh_perm = MinHashPermutation(L, universe, seed=SEED)159    mh_lin = MinHashLinear(L, seed=SEED)160161    A, B = sets[0], sets[1]162    sig_p, sig_l = mh_perm.get_signature(A), mh_lin.get_signature(A)163    same = float((sig_p == sig_l).mean())164    print(f"  同一集合的两种签名：逐位相同比例 = {same:.3f}")165    print("  （两者只是「等价的估计器」，用的是不同随机源，签名本身不应相同；")166    print("   要比较的是「碰撞比例是否都接近真实 Jaccard」）\n")167168    rows = []169    for tgt in (0.2, 0.4, 0.6, 0.8):170        A2, B2 = make_pair(tgt, seed=SEED + int(tgt * 10))171        J = jaccard(A2, B2)172        sp = float((mh_perm.get_signature(A2)173                    == mh_perm.get_signature(B2)).mean())174        sl = float((mh_lin.get_signature(A2)175                    == mh_lin.get_signature(B2)).mean())176        # 单次估计的方差是 sqrt(J(1−J)/64) ≈ 0.05–0.06，一行的偏差说明不了问题，177        # 所以同时列出 200 个哈希种子下的均值，用来判断「有没有系统性偏差」178        many = []179        for t in range(200):180            mm = MinHashLinear(L, seed=SEED + 3000 + t)181            many.append(float((mm.get_signature(A2) == mm.get_signature(B2)).mean()))182        rows.append([f"{tgt:.1f}", f"{J:.3f}", f"{sp:.3f}", f"{sl:.3f}",183                     f"{np.mean(many):.3f}"])184    table(["目标 J", "真实 J", "排列估计", "线性哈希估计", "线性哈希 200 种子均值"], rows)185    print("\n  逐列看：单次估计在真实 J 上下随机波动（L=64 时 std ≈ 0.05–0.06），")186    print("  而 200 个种子的均值几乎贴着真实 J——两者都是无偏估计器，没有系统性偏差。")187    print("  这正是「碰撞概率 = Jaccard」的直接体现：偏差来自抽样，不来自实现。")188    print("  ② 会把「误差如何随 L 收缩」定量画出来。\n")189190    # ---------------- 2) 误差 vs 理论标准差 ----------------191    print("② 估计误差随签名长度 L 的下降（与 sqrt(J(1−J)/L) 对照）：")192    print("   （同组集合对 J≈0.19–0.25，各用 300 个不同哈希种子重复，取标准差）")193    rows = []194    pairs = [(i, j) for i in range(8) for j in range(i + 1, 8)]195    for L_ in (16, 64, 256, 1024):196        for i, j in pairs:197            J = jaccard(sets[i], sets[j])198            if not (0.15 <= J <= 0.85):      # 太接近 0/1 的对方差近乎 0，比值无意义199                continue200            ests = []201            for t in range(300):202                m = MinHashLinear(L_, seed=SEED + 1000 + t)203                sig_i, sig_j = m.get_signature(sets[i]), m.get_signature(sets[j])204                ests.append(float((sig_i == sig_j).mean()))205            ests = np.array(ests)206            emp_std = float(ests.std())207            theo = float(np.sqrt(J * (1 - J) / L_))208            rows.append([f"({i},{j})", f"{J:.2f}", L_, f"{ests.mean():.3f}",209                         f"{emp_std:.4f}", f"{theo:.4f}",210                         f"{emp_std / theo:.2f}"])211    table(["集合对", "真实 J", "L", "均值估计", "经验 std", "理论 std", "比值"], rows)212    print("\n  经验标准差与 sqrt(J(1−J)/L) 基本吻合（比值接近 1）。")213    print("  这就是正文的方差推导——每一位是伯努利变量，L 位独立，方差 J(1−J)/L。")214    print("  注意：J 越接近 0 或 1，方差越小，估计本来就更容易准；")215    print("  真正难估的是 J≈0.5 的中间区域，这也决定了检索阈值附近的误判率。\n")216217    # ---------------- 3) 工程视角：内存与耗时 ----------------218    print("③ 为什么工业界用线性哈希（universe = 10^6，L = 128）：")219    rows = []220    for name, obj in (("随机排列", MinHashPermutation(128, 10 ** 6, seed=SEED)),221                      ("线性哈希", MinHashLinear(128, seed=SEED))):222        import time223        t0 = time.perf_counter()224        for _ in range(200):225            obj.get_signature(sets[0])226        dt = (time.perf_counter() - t0) / 200227        mb = obj.nbytes() / 1e6228        size_s = f"{mb:.1f} MB" if mb >= 1 else f"{mb * 1000:.0f} KB"229        rows.append([name, size_s, f"{dt * 1e6:.0f} µs / 次"])230    table(["实现", "索引内存", "单集合签名耗时"], rows)231    print("\n  随机排列需要 128 × 10^6 个整数 ≈ 512 MB（还没算 universe 更大时）；")232    print("  线性哈希只需要 128 对 (a, b)，约 2 KB。耗时上前者也慢一个量级，")233    print("  因为它要对每个元素做一次随机内存访问。这就是正文「不适合大规模落地」的原因。\n")234235    # ---------------- 4) 用 MinHash 找相似集合 ----------------236    print("④ 检索演示：给一个查询集合，用 MinHash 签名找最相似的 5 个：")237    sets_big, labels_big = make_sets(n_sets=200, universe=20000, size=80,238                                     n_groups=10)239    mh = MinHashLinear(256, seed=SEED)240    sigs = np.stack([mh.get_signature(s) for s in sets_big])   # [200, 256]241    q_idx = 0242    q_sig = sigs[q_idx]243    sim = (sigs == q_sig[None, :]).mean(axis=1)244    sim[q_idx] = -1.0                                    # 排除查询自身245    order = np.argsort(-sim)[:5]246    rows = []247    for k in order:248        true_j = jaccard(sets_big[q_idx], sets_big[k])249        rows.append([k, labels_big[k], f"{sim[k]:.3f}", f"{true_j:.3f}"])250    table(["候选编号", "所属簇", "签名相似度", "真实 Jaccard"], rows)251    print("\n  Top-5 全部来自查询集合自己的簇，签名相似度与真实 Jaccard 同序——")252    print("  但注意这只是「对 200 个集合全量算签名相似度」，还不是 ANN。")253    print("  真正的子线性检索要靠 §18.3.3 的 LSH Banding 把签名切成 band 建桶。\n")254255    print("⑤ 一句话：MinHash 把集合压缩成 L 个整数，碰撞概率 = Jaccard；")256    print("   签名长度控制精度（std=sqrt(J(1−J)/L)），线性哈希控制成本。")257258259if __name__ == "__main__":260    main()