dpp_greedy_map.py

1"""对应正文第 7 章 §7.5.2「DPP 算法」(Determinantal Point Process,行列式点过程)。23本实现做了什么4----------------51. 按正文公式构造核矩阵:B_i = r_i f_i,||f_i||_2 = 1,6       S_ij = (1 + <f_i, f_j>) / 2        (把余弦从 [-1,1] 压到 [0,1] 保半正定)7       L  = Diag(r) · S · Diag(r)8       L' = Diag(exp(alpha r)) · S · Diag(exp(alpha r)),alpha = theta / (2(1-theta))9   并数值验证 L 的半正定性(最小特征值 >= 0)与正文公式 log_det_L 的分解恒等式;102. 实现两种贪心 MAP 推断,二者选出的集合必须完全一致:11   (a) 朴素版:每轮对每个候选显式构造子矩阵 L_{Y∪{i}} 并做一次 O(k^3) 的12       Cholesky 分解求 log det,复杂度 O(K^2 N k^3);13   (b) 快速版:复刻正文 code_dpp 的增量递推(Chen et al., NeurIPS 2018),14       维护 c_i 与 d_i^2,每轮只做 O(N k) 的更新:15           e_i    = (L_ji - <c_j, c_i>) / d_j16           c_i'   = [c_i, e_i]17           d_i'^2 = d_i^2 - e_i^218           j      = argmax_i log(d_i^2)19   打印两者耗时与加速比;203. 几何解释(正文图 7.6):det(L_Y) = Vol^2({B_i}_{i in Y}),21   打印选中集合的 log det(L_Y)、log det(S_Y)(纯多样性项)与超平行体体积;224. theta 扫描:theta=1 退化为纯相关性(alpha 置 0 时 L' 的相关性权重消失,23   正文代码对 theta==1 令 alpha=0),theta 越小多样性权重越大。24   打印平均相关性 / 相关性损失 / log det(S_Y) / ILD / 类目熵。2526关键公式(§7.5.2)27-----------------28    P(Y) ∝ det(L_Y) = Vol^2({B_i}_{i in Y})29    log det L = sum_i log(r_i^2) + log det(S)30    log det L' = lambda sum_i log(r_i^2) + (1-lambda) log det(S)31    j = argmax_{i in Z\\Y} [log det(L_{Y∪{i}}) - log det(L_Y)] = argmax_i log(d_i^2)3233复杂度(正文结论):sum_{k=1..K} O(k(N-k)) = O(K^2 N),选满 N 个时为 O(N^3)。34"""3536import math37import random38import time3940import numpy as np4142CATEGORIES = ["影视", "美食", "萌宠", "游戏", "音乐", "知识"]434445# ---------------------------------------------------------------------------46# [MOCK] 合成候选集:同类目物品 embedding 围绕同一类目中心,因此天然高相似47# ---------------------------------------------------------------------------48def build_candidates(n_items, dim, rng):49    centers = {}50    for c in CATEGORIES:51        v = np.array([rng.gauss(0, 1) for _ in range(dim)])52        centers[c] = v / np.linalg.norm(v)53    cats, embs, rels = [], [], []54    for i in range(n_items):55        cat = CATEGORIES[i % len(CATEGORIES)]56        v = centers[cat] + 0.30 * np.array([rng.gauss(0, 1) for _ in range(dim)])57        embs.append(v / np.linalg.norm(v))58        cats.append(cat)59        bonus = 0.25 if cat in ("影视", "美食") else 0.060        rels.append(min(1.0, max(0.05, 0.5 + 0.3 * rng.gauss(0, 1) + bonus)))61    return cats, np.array(embs), np.array(rels)626364# ---------------------------------------------------------------------------65# 核矩阵构造66# ---------------------------------------------------------------------------67def similarity_matrix(embeddings):68    """S_ij = (1 + <f_i, f_j>) / 2,f 已 L2 归一化。"""69    f = embeddings / np.linalg.norm(embeddings, axis=1, keepdims=True)70    return (1.0 + f @ f.T) / 2.0717273def kernel_matrix(embeddings, r, theta):74    """L' = Diag(exp(alpha r)) · S · Diag(exp(alpha r)),alpha = theta/(2(1-theta))。"""75    S = similarity_matrix(embeddings)76    alpha = 0.0 if theta >= 1.0 else theta / (2.0 * (1.0 - theta))77    w = np.exp(alpha * r)78    return w[:, None] * S * w[None, :], S, alpha798081# ---------------------------------------------------------------------------82# (a) 朴素贪心:每个候选显式算 log det(L_{Y∪{i}})83# ---------------------------------------------------------------------------84def greedy_map_naive(L, K):85    n = L.shape[0]86    selected = []87    for _ in range(K):88        best_j, best_gain = None, -float("inf")89        base = logdet_cholesky(L[np.ix_(selected, selected)]) if selected else 0.090        for i in range(n):91            if i in selected:92                continue93            idx = selected + [i]94            gain = logdet_cholesky(L[np.ix_(idx, idx)]) - base95            if gain > best_gain:96                best_j, best_gain = i, gain97        selected.append(best_j)98    return selected99100101def logdet_cholesky(sub):102    """用 Cholesky 分解算 log det,比直接求行列式数值更稳。"""103    if sub.shape[0] == 0:104        return 0.0105    jitter = 1e-10 * np.eye(sub.shape[0])106    low = np.linalg.cholesky(sub + jitter)107    return 2.0 * float(np.sum(np.log(np.diag(low))))108109110# ---------------------------------------------------------------------------111# (b) 快速贪心:正文 code_dpp 的 Cholesky 增量递推112# ---------------------------------------------------------------------------113def fast_greedy_map_inference(item_embeddings, r_u, theta, K):114    """复刻正文 code_dpp。返回选中物品下标列表(长度 K)。"""115    n = item_embeddings.shape[0]116    if K <= 0 or K > n:117        raise ValueError(f"K 必须为正整数且不超过物品总数 {n}")118119    L, _, _ = kernel_matrix(item_embeddings, r_u, theta)120121    c = [[] for _ in range(n)]          # c[i]:物品 i 的辅助向量,长度 = 已选个数122    d_squared = L.diagonal().copy()     # d_i^2 初始为 L 的对角元123    selected = []124125    while len(selected) < K:126        mask = np.ones(n, dtype=bool)127        mask[selected] = False128        valid = np.where(mask)[0]129        if valid.size == 0:130            break131        # argmax log(d_i^2) 等价于 argmax d_i^2(log 单调),这里按正文写法保留 log132        safe = np.maximum(d_squared[valid], 1e-300)133        j = int(valid[int(np.argmax(np.log(safe)))])134        selected.append(j)135        if len(selected) == K:136            break137138        d_j = math.sqrt(max(d_squared[j], 1e-300))139        cj = np.array(c[j]) if c[j] else None140        for i in range(n):141            if i in selected:142                continue143            dot = float(cj @ np.array(c[i])) if cj is not None and c[i] else 0.0144            e_i = (L[j, i] - dot) / d_j145            c[i].append(e_i)146            d_squared[i] -= e_i ** 2147    return selected148149150# ---------------------------------------------------------------------------151# 指标152# ---------------------------------------------------------------------------153def set_metrics(cats, embeddings, r, S, picked):154    m = len(picked)155    sel_emb = embeddings[picked]156    sims = sel_emb @ sel_emb.T157    ild = 0.0158    for x in range(m):159        for y in range(x + 1, m):160            ild += 1.0 - sims[x, y]161    ild /= m * (m - 1) / 2162    sel_cats = [cats[i] for i in picked]163    counts = {c: sel_cats.count(c) for c in set(sel_cats)}164    ent = -sum((v / m) * math.log(v / m) for v in counts.values())165    return {166        "rel": float(np.mean(r[picked])),167        "ild": ild,168        "n_cat": len(counts),169        "entropy": ent,170        "logdet_S": logdet_cholesky(S[np.ix_(picked, picked)]),171    }172173174def main():175    random.seed(42)176    np.random.seed(42)177    rng = random.Random(42)178179    n_items, dim, top_k = 60, 16, 10180    cats, embs, rels = build_candidates(n_items, dim, rng)181182    # ---------------- (1) 核矩阵性质验证 ----------------183    print("== (1) 核矩阵构造与性质验证 ==")184    S = similarity_matrix(embs)185    L_plain = rels[:, None] * S * rels[None, :]          # L = Diag(r) S Diag(r)186    eig_S = np.linalg.eigvalsh(S)187    eig_L = np.linalg.eigvalsh(L_plain)188    print(f"候选集 N={n_items},emb 维度 d={dim},目标长度 K={top_k}")189    print(f"S 取值范围           : [{S.min():.4f}, {S.max():.4f}]  (已压到 [0,1])")190    print(f"S 最小特征值         : {eig_S.min():+.3e}   (>= 0 → 半正定)")191    print(f"L 最小特征值         : {eig_L.min():+.3e}   (>= 0 → 半正定)")192    sub = list(range(8))193    lhs = logdet_cholesky(L_plain[np.ix_(sub, sub)])194    rhs = float(np.sum(np.log(rels[sub] ** 2))) + logdet_cholesky(S[np.ix_(sub, sub)])195    print("正文公式 log det L = sum_i log(r_i^2) + log det S 验证(取前 8 个物品子阵):")196    print(f"  左边 log det(L_Y) = {lhs:+.6f}   右边 = {rhs:+.6f}   "197          f"|差| = {abs(lhs - rhs):.3e}")198199    # ---------------- (2) 两种贪心的一致性与加速比 ----------------200    print("\n== (2) 贪心 MAP 推断:朴素实现 vs Cholesky 增量实现 ==")201    theta = 0.7202    L, _, alpha = kernel_matrix(embs, rels, theta)203    print(f"theta={theta}  →  alpha = theta/(2(1-theta)) = {alpha:.4f}")204205    t0 = time.perf_counter()206    naive = greedy_map_naive(L, top_k)207    t_naive = time.perf_counter() - t0208209    t0 = time.perf_counter()210    fast = fast_greedy_map_inference(embs, rels, theta, top_k)211    t_fast = time.perf_counter() - t0212213    print(f"朴素贪心选中: {naive}")214    print(f"快速贪心选中: {fast}")215    print(f"两者是否一致: {naive == fast}")216    print(f"朴素耗时 {t_naive * 1000:.2f} ms,快速耗时 {t_fast * 1000:.2f} ms,"217          f"加速比 {t_naive / t_fast:.2f}x")218    print("复杂度:朴素每轮要对 N-k 个候选各做一次 O(k^3) Cholesky;"219          "快速版每轮只做 O(N k) 的 c/d 递推 → 总 O(K^2 N)")220221    print("\n不同规模下的耗时对比(单位 ms):")222    print(f"{'N':>6} {'K':>4} {'朴素':>10} {'增量':>10} {'加速比':>8} {'集合一致':>8}")223    for n_probe, k_probe in ((40, 8), (80, 10), (120, 12), (200, 15)):224        cat_p, emb_p, rel_p = build_candidates(n_probe, dim, random.Random(7))225        L_p, _, _ = kernel_matrix(emb_p, rel_p, theta)226        t0 = time.perf_counter()227        a_p = greedy_map_naive(L_p, k_probe)228        t_a = (time.perf_counter() - t0) * 1000229        t0 = time.perf_counter()230        b_p = fast_greedy_map_inference(emb_p, rel_p, theta, k_probe)231        t_b = (time.perf_counter() - t0) * 1000232        print(f"{n_probe:>6} {k_probe:>4} {t_a:>10.2f} {t_b:>10.2f} "233              f"{t_a / t_b:>7.2f}x {str(a_p == b_p):>8}")234235    # 增量 d_i^2 与真实边际增益的一致性(论证递推正确)236    print("\n增量量 log(d_i^2) 与真实边际增益 log det(L_{Y∪{i}}) - log det(L_Y) 对照"237          "(已选前 3 个,取前 5 个未选候选):")238    Y = fast[:3]239    base = logdet_cholesky(L[np.ix_(Y, Y)])240    # 重放递推拿到第 4 轮的 d_i^2241    c = [[] for _ in range(n_items)]242    d2 = L.diagonal().copy()243    for step, j in enumerate(Y):244        d_j = math.sqrt(d2[j])245        cj = np.array(c[j]) if c[j] else None246        for i in range(n_items):247            if i in Y[:step + 1]:248                continue249            dot = float(cj @ np.array(c[i])) if cj is not None and c[i] else 0.0250            e_i = (L[j, i] - dot) / d_j251            c[i].append(e_i)252            d2[i] -= e_i ** 2253    probe = [i for i in range(n_items) if i not in Y][:5]254    print(f"{'item':>5} {'log(d_i^2)':>12} {'真实边际增益':>14} {'|差|':>10}")255    for i in probe:256        idx = Y + [i]257        true_gain = logdet_cholesky(L[np.ix_(idx, idx)]) - base258        print(f"{i:>5} {math.log(d2[i]):>12.6f} {true_gain:>14.6f} "259              f"{abs(math.log(d2[i]) - true_gain):>10.2e}")260261    # ---------------- (3) 几何解释:体积 ----------------262    print("\n== (3) 几何解释 det(L_Y) = Vol^2({B_i}) ==")263    print(f"{'集合':>14} {'log det(L_Y)':>13} {'Vol':>12} {'log det(S_Y)':>13}")264    top_rel = list(np.argsort(-rels)[:top_k])265    for name, picked in (("纯相关性 Top-K", top_rel), ("DPP 贪心", fast)):266        ld = logdet_cholesky(L[np.ix_(picked, picked)])267        vol = math.exp(0.5 * ld)268        print(f"{name:>14} {ld:>13.4f} {vol:>12.4e} "269              f"{logdet_cholesky(S[np.ix_(picked, picked)]):>13.4f}")270    print("→ DPP 选出的集合张成的超平行体体积更大,即物品之间更「排斥」。")271272    # ---------------- (4) theta 扫描 ----------------273    print("\n== (4) theta 扫描:相关性-多样性权衡 ==")274    base_rel = float(np.mean(rels[top_rel]))275    print(f"{'theta':>6} {'alpha':>8} {'avg_rel':>8} {'rel损失':>8} "276          f"{'log det S_Y':>12} {'ILD':>7} {'类目数':>6} {'类目熵':>7}")277    for theta in (0.0, 0.2, 0.4, 0.6, 0.8, 0.9, 0.99, 1.0):278        picked = fast_greedy_map_inference(embs, rels, theta, top_k)279        _, _, a = kernel_matrix(embs, rels, theta)280        m = set_metrics(cats, embs, rels, S, picked)281        loss = (base_rel - m["rel"]) / base_rel * 100282        print(f"{theta:>6.2f} {a:>8.3f} {m['rel']:>8.4f} {loss:>7.2f}% "283              f"{m['logdet_S']:>12.4f} {m['ild']:>7.4f} {m['n_cat']:>6d} "284              f"{m['entropy']:>7.4f}")285    print("注:theta=1 时正文代码令 alpha=0,此时 L'=S,相关性被完全抹掉,"286          "结果反而是「纯多样性」;真正的纯相关性需令 alpha→∞(见 theta=0.99)。")287288    print("\n== (5) 选中序列的类目分布 ==")289    for name, picked in (("纯相关性 Top-K", top_rel),290                         ("DPP theta=0.7", fast)):291        print(f"  {name:>14}: " + "->".join(cats[i] for i in picked))292293294if __name__ == "__main__":295    main()