""" ============================================================================= 第 4 讲 · 代码 4：识别与 Rashomon 效应 —— 哪些结论是数据真正支持的？ ============================================================================= 对应讲义：lectures/lecture04_factor_models.md 第 5 节 ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- 因子模型学出了每个用户的口味向量、每个回答的特征向量。很诱人的做法是去"解读"它们： "第一维代表详细程度""这个用户在第二维上很高"……可是这些数字有多少是数据真正决定的？ 演示 1：旋转不变性 —— 把所有用户向量和回答向量一起旋转同一个角度，预测**一点都不变**。 所以单独的坐标没有意义，只有向量之间的点积（夹角、长度）有意义。 演示 2：换一个随机种子重新训练，得到的嵌入"看起来完全不同"，预测却几乎一样； 用"普鲁克分析"（Procrustes，找最合适的旋转 / 翻转）对齐以后，两者几乎重合， 也能和真实的风格维度对上。 演示 3：Rashomon 效应 —— 换不同的 K、不同的正则化，得到一批验证集上"一样好"的模型。 它们在具体某个用户、某个回答上的预测和推荐，却常常不一样。 （名字来自黑泽明 1950 年的电影《罗生门》：几个目击者对同一件事给出了各不相同、却都说得通的叙述。） ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 04_identification_rashomon.py # 约 8 秒 只用到 Python 标准库，以及第 1 讲的 toy_world.py、本目录的 factor_models.py。 """ import math import random from factor_models import FactorModel, TasteWorld, dot, evaluate, fit_mf, log_loss, predict_prob, split from toy_world import mean, pad, print_header def rotate(vec, theta): """把二维向量逆时针旋转 θ 弧度。""" c, s = math.cos(theta), math.sin(theta) return [c * vec[0] - s * vec[1], s * vec[0] + c * vec[1]] def correlation(xs, ys): """皮尔逊相关系数。""" mx, my = mean(xs), mean(ys) sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) return sxy / (sum((x - mx) ** 2 for x in xs) * sum((y - my) ** 2 for y in ys)) ** 0.5 def procrustes_2d(source, target): """ 二维的正交普鲁克分析：找一个旋转（可能再加一次翻转）R，让 R·source 尽量接近 target。 【思路】对二维向量，最优的旋转角有封闭解：θ = atan2(Σ(s_x t_y − s_y t_x), Σ(s_x t_x + s_y t_y))。 翻转的情况：先把 source 的第二个坐标取负，再求一次最优旋转；两种里取误差小的。 （高维时用奇异值分解求解，第 5 讲起有 numpy 可以用。） 【参数】source、target：等长的二维向量列表（同一组回答，分别来自两个模型 / 真实值） 【返回】(对齐后的 source 列表, 是否翻转, 旋转角) 【边界】只允许旋转和翻转，不允许缩放 —— 两边的整体尺度不同时，对齐后仍有尺度差（可以先各自标准化）。 """ best = None for flip in (False, True): src = [[a, -b] if flip else [a, b] for a, b in source] num = sum(s[0] * t[1] - s[1] * t[0] for s, t in zip(src, target)) den = sum(s[0] * t[0] + s[1] * t[1] for s, t in zip(src, target)) theta = math.atan2(num, den) aligned = [rotate(s, theta) for s in src] err = sum((a[0] - t[0]) ** 2 + (a[1] - t[1]) ** 2 for a, t in zip(aligned, target)) if best is None or err < best[0]: best = (err, aligned, flip, theta) return best[1], best[2], best[3] # ============================================================================= # 演示 1：旋转不变性 # ============================================================================= def demo_rotation(model, world, test): print_header("演示 1：把所有向量一起旋转 —— 预测一点都不变") theta = math.radians(73) rotated = FactorModel(model.k, dict(model.a), {i: rotate(v, theta) for i, v in model.u.items()}, dict(model.z), {j: rotate(v, theta) for j, v in model.v.items()}) scaled = FactorModel(model.k, dict(model.a), {i: [2 * x for x in v] for i, v in model.u.items()}, dict(model.z), {j: [x / 2 for x in v] for j, v in model.v.items()}) diff_rot = max(abs(predict_prob(model, i, j) - predict_prob(rotated, i, j)) for i, j, _ in test) diff_scl = max(abs(predict_prob(model, i, j) - predict_prob(scaled, i, j)) for i, j, _ in test) j0 = world.items[0] print(f"一个拟合好的 K=2 模型。回答 {j0} 的向量：原来 {[round(x, 3) for x in model.v[j0]]}，" f"旋转 73° 以后 {[round(x, 3) for x in rotated.v[j0]]}") print(f" 测试集上预测概率的最大差别：旋转 {diff_rot:.1e}；用户向量 ×2、回答向量 ÷2 {diff_scl:.1e}") print(""" 观察： · 旋转矩阵 R 满足 RᵀR = I，所以 (R·U)·(R·V) = Uᵀ RᵀR V = U·V：点积不变，预测就不变。更一般地， 用户向量乘以任何可逆矩阵 A 的转置、回答向量乘以 A⁻¹，点积也不变。数据只能确定所有的点积 （“格拉姆矩阵” U Vᵀ），确定不了坐标本身 —— 这就是因子模型的**不可识别性**。 · 缩放（用户 ×2、回答 ÷2）也不改变预测。L2 正则化会消除缩放的自由度（让两边的长度“平衡”时惩罚最小）， 但消除不了旋转和翻转（旋转不改变长度，正则化项也不变）。 · 所以“第一维代表详细程度”这种解读，要先固定一个旋转（例如和某个已知的特征对齐，演示 2）才有意义。""") # ============================================================================= # 演示 2：不同的随机种子 # ============================================================================= def demo_seeds(world, train, test): print_header("演示 2：换个随机种子 —— 嵌入看起来不同，预测几乎一样") models = [fit_mf(train, world.users, world.items, 2, 3.0, random.Random(seed)) for seed in (1, 2, 3)] items = world.items print(f" {pad('种子', 10)} {pad('测试对数损失', 12, 'right')} {pad('第 1 维和种子 1 的相关', 22, 'right')} " f"{pad('预测和种子 1 的相关', 18, 'right')}") base = models[0] base_pred = [predict_prob(base, i, j) for i, j, _ in test] for seed, m in zip((1, 2, 3), models): ll = evaluate(m, test)[0] c1 = correlation([m.v[j][0] for j in items], [base.v[j][0] for j in items]) cp = correlation([predict_prob(m, i, j) for i, j, _ in test], base_pred) print(f" {pad(f'种子 {seed}', 10)} {ll:>12.4f} {c1:>22.3f} {cp:>18.4f}") print("\n 用普鲁克分析把种子 2、3 的回答向量对齐到种子 1（只允许旋转和翻转）：") for seed, m in zip((2, 3), models[1:]): aligned, flip, theta = procrustes_2d([m.v[j] for j in items], [base.v[j] for j in items]) c = correlation([a[0] for a in aligned], [base.v[j][0] for j in items]) print(f" 种子 {seed}：翻转 {'是' if flip else '否'}，旋转 {math.degrees(theta):+.0f}°，对齐后第 1 维的相关 = {c:.3f}") # 和真实的风格维度对齐：先把两边都标准化（真实风格和学到的向量尺度不同） def standardize(vecs): sd = math.sqrt(mean([v[0] ** 2 + v[1] ** 2 for v in vecs]) / 2) return [[v[0] / sd, v[1] / sd] for v in vecs] true = standardize([world.style[j] for j in items]) aligned, flip, theta = procrustes_2d(standardize([base.v[j] for j in items]), true) c_d = correlation([a[0] for a in aligned], [t[0] for t in true]) c_f = correlation([a[1] for a in aligned], [t[1] for t in true]) print(f"\n 把种子 1 的回答向量对齐到真实的风格（详细程度、正式程度）：") print(f" 对齐后第 1 维和“详细程度”的相关 = {c_d:.3f}，第 2 维和“正式程度”的相关 = {c_f:.3f}") print(""" 观察： · 三个种子的测试对数损失几乎一样，预测的相关接近 1；可是“第 1 维”的坐标在不同种子之间相关很低，甚至是负的。 如果你只看一个种子的结果就宣布“第 1 维 = 详细程度”，换个种子这个结论就变了。 · 普鲁克分析找到合适的旋转 / 翻转以后，不同种子的嵌入几乎完全重合 —— 它们是**同一个**解的不同“摆法”。 · 对齐到真实的风格维度以后，学到的两维和“详细程度”“正式程度”高度相关：模型确实学到了真实的结构， 只是它不知道这两个方向该叫什么名字、该朝哪边。真实项目里没有“真实风格”可以对齐，通常用已知的特征 （例如回答的长度）去对齐或解释隐向量。""") # ============================================================================= # 演示 3：Rashomon 效应 # ============================================================================= def demo_rashomon(world, train, valid, test, eps=0.01): print_header("演示 3：Rashomon 效应 —— 一批“一样好”的模型，具体预测却不一样") configs = [(k, l2, seed) for k in (1, 2, 3) for l2 in (2.0, 3.0, 5.0) for seed in (11, 12)] fitted = [] for k, l2, seed in configs: m = fit_mf(train, world.users, world.items, k, l2, random.Random(seed)) vprobs = [predict_prob(m, i, j) for i, j, _ in valid] vloss = log_loss(vprobs, [y for _, _, y in valid]) fitted.append((vloss, k, l2, seed, m)) best = min(f[0] for f in fitted) rashomon = [f for f in fitted if f[0] <= best * (1 + eps)] print(f"训练了 {len(fitted)} 个模型（K ∈ {{1, 2, 3}} × 正则化 ∈ {{2, 3, 5}} × 2 个种子）。") print(f"验证集对数损失最好的是 {best:.4f}；在它的 {eps:.0%} 以内的模型（“Rashomon 集合”）有 {len(rashomon)} 个：") for vloss, k, l2, seed, _ in sorted(rashomon, key=lambda f: f[0]): print(f" K = {k}，正则化 {l2:g}，种子 {seed}：验证损失 {vloss:.4f}") models = [f[4] for f in rashomon] # (1) 测试格子上的"预测分歧"：有多少格子，Rashomon 集合里的模型对"接受还是拒绝"意见不一 disagree = 0 spread = [] for i, j, _ in test: ps = [predict_prob(m, i, j) for m in models] if min(ps) < 0.5 < max(ps): disagree += 1 spread.append(max(ps) - min(ps)) # (2) 每个用户的"第一推荐"：有多少用户，不同模型推荐的第一名不一样 top_diff = 0 for i in world.users: tops = {max(world.items, key=lambda j: predict_prob(m, i, j)) for m in models} if len(tops) > 1: top_diff += 1 print(f"\n 测试集的 {len(test)} 个格子里，模型对“接受还是拒绝”意见不一的：{disagree} 个（{disagree / len(test):.1%}）") print(f" 同一个格子，不同模型给出的接受概率最多相差：平均 {mean(spread):.3f}，最大 {max(spread):.3f}") print(f" {len(world.users)} 个用户里，不同模型推荐的“第一名”不一样的：{top_diff} 个（{top_diff / len(world.users):.0%}）") print(""" 观察： · 这些模型在验证集上的表现几乎没有区别 —— 单看总体指标，你没有理由偏爱其中任何一个。 · 可是它们对**具体**的用户和回答，常常给出不同的判断：一部分格子的“接受 / 拒绝”预测会翻转， 相当一部分用户会因为“选了哪个模型”而收到不同的第一推荐。Breiman（2001）称之为 Rashomon 效应， Marx 等（2020）把它叫“预测多重性”（predictive multiplicity）。 · 和演示 1、2 的区别：旋转不变性是**同一个**模型的不同写法，预测完全相同；Rashomon 效应是**真正不同**的模型， 总体一样好，细节不同。前者是“参数不可识别”，后者是“数据不足以在一批好模型里做选择”。 · 对 RLHF 的含义：很多不同的奖励模型能同样好地解释人类反馈，但它们对某些回答的打分不同， 优化不同的奖励模型会得到不同的策略（Skalse 等 2023）。报告时最好给出“一批好模型”的范围， 而不只是一个“最好的模型”；第 6 讲的奖励模型集成就是一种做法。""") def main(): rng = random.Random(1950) # 1950：电影《罗生门》上映的年份 world = TasteWorld(rng) obs = world.sample(rng, 0.6) train, rest = split(obs, rng, test_frac=0.3) valid, test = rest[:len(rest) // 2], rest[len(rest) // 2:] model = fit_mf(train, world.users, world.items, 2, 3.0, random.Random(1)) demo_rotation(model, world, test) demo_seeds(world, train, test) demo_rashomon(world, train, valid, test) if __name__ == "__main__": main()