""" ============================================================================= 第 13 讲 · 代码 4：偏好贝叶斯优化 —— 只能问"哪个更好"时，怎样找到最好的参数 ============================================================================= 对应讲义：lectures/lecture13_dueling_pbo.md 第 5 节 ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- 外骨骼的步态参数 x ∈ [0, 1]（例如步频），每个使用者觉得最舒服的 x 不同。人说不出"这个步态值 0.73 分"， 只能走两种步态后说"这个更舒服"（第 10 讲第 6.3 节，Tucker 等 2020）。每次比较都要人真的走一遍， 所以只能问几十次。怎样尽快找到最舒服的 x？ 偏好贝叶斯优化 = 第 6 讲的**高斯过程偏好学习**（效用 f(x) 的先验是 GP，比较的似然是 BT） + 第 12 讲的**采集函数**（但现在每次要挑**一对** x 去比较）。 几种"挑哪一对"的方法： · 随机一对 · 最好 vs 最不确定：一个是当前后验均值最高的 x，另一个是后验标准差最大的 x · 最好 vs 乐观：一个是当前最好的 x，另一个是"和最好的比，乐观地看能好多少"最大的 x（第 6 讲代码 04 演示 3） · 自我陪练（双 Thompson 采样）：从后验抽两条效用函数，各取最大值点（Sui 等 2017；Tucker 等 2020 的 CoSpar） · EUBO：挑让 E[max(f(x₁), f(x₂))] 最大的一对 —— "这一对里较好的那个，期望有多好"（Astudillo 等 2023） 真实的效用：x ≈ 0.7 处一个高峰，x ≈ 0.2 处一个矮峰。人的回答有噪声：P(a 更舒服) = σ(3 (f(a) − f(b)))。 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 04_preferential_bo.py # 约 30 秒 需要 numpy。依赖第 6 讲的 bayes_tools.py（GP 偏好学习；它会再导入第 5 讲的 bt_numpy.py）。 """ import math import os import sys import numpy as np _HERE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, os.path.join(_HERE, "..", "lecture06")) from bayes_tools import gp_pref_laplace # noqa: E402 from duel_tools import pad, print_header # noqa: E402 GRID = np.linspace(0.0, 1.0, 61) KAPPA = 3.0 METHODS = {"random": "随机一对", "best_var": "最好 vs 最不确定", "best_ucb": "最好 vs 乐观", "sparring": "自我陪练（双 Thompson）", "eubo": "EUBO"} def utility(x): """真实的舒适度：x ≈ 0.7 处高峰（1.0），x ≈ 0.2 处矮峰（0.5）。只在模拟里用。""" return 1.0 * np.exp(-((x - 0.7) / 0.1) ** 2) + 0.5 * np.exp(-((x - 0.2) / 0.15) ** 2) F = utility(GRID) def norm_cdf(z): """标准正态的累积分布函数（数组）。""" return 0.5 * (1.0 + np.vectorize(math.erf)(np.asarray(z) / math.sqrt(2.0))) def norm_pdf(z): """标准正态的密度。""" return np.exp(-0.5 * np.asarray(z) ** 2) / math.sqrt(2 * math.pi) def eubo_pair(mean, cov): """ 在所有网格点对 (i, j) 上计算 EUBO = E[max(f_i, f_j)]，返回最大的一对。 【公式】(f_i, f_j) 服从二元正态。记 θ = sqrt(σ_i² + σ_j² − 2 Cov_ij)，α = (μ_i − μ_j)/θ： E[max(f_i, f_j)] = μ_i Φ(α) + μ_j Φ(−α) + θ φ(α) （两个正态变量最大值的期望，Clark 1961） 【直觉】这一对里较好的那个，期望有多好。两个都有希望、而且彼此不太相关（θ 大）的一对，EUBO 最高 —— "比较之后，至少有一个是好的"。 """ sd2 = np.diag(cov) i, j = np.triu_indices(len(mean), k=1) theta = np.sqrt(np.clip(sd2[i] + sd2[j] - 2 * cov[i, j], 1e-12, None)) alpha = (mean[i] - mean[j]) / theta val = mean[i] * norm_cdf(alpha) + mean[j] * norm_cdf(-alpha) + theta * norm_pdf(alpha) k = int(np.argmax(val)) return int(i[k]), int(j[k]) def choose_pair(method, fit, rng): """按方法挑下一对要比较的网格下标。""" G = len(GRID) mean, cov = fit.mean, fit.cov best = int(np.argmax(mean)) if method == "random": a, b = rng.choice(G, size=2, replace=False) return int(a), int(b) if method == "best_var": sd = np.sqrt(np.clip(np.diag(cov), 0, None)) sd[best] = -1 return best, int(np.argmax(sd)) if method == "best_ucb": var_diff = np.clip(np.diag(cov) + cov[best, best] - 2 * cov[:, best], 0, None) ucb = mean - mean[best] + 2 * np.sqrt(var_diff) ucb[best] = -np.inf return best, int(np.argmax(ucb)) if method == "sparring": L = np.linalg.cholesky(cov + 1e-8 * np.eye(G)) a = int(np.argmax(mean + L @ rng.normal(size=G))) b = int(np.argmax(mean + L @ rng.normal(size=G))) if a == b: # 两个样本推荐了同一个点：再抽一次 b = int(np.argmax(mean + L @ rng.normal(size=G))) if a == b: b = int(rng.integers(G)) return a, b if method == "eubo": return eubo_pair(mean, cov) raise ValueError(method) def run(method, rep, budget=40, init=3): """ 一次偏好贝叶斯优化：先随机问 init 对，再按方法问 budget 对。 【返回】dict：问题数 → 推荐点（后验均值最大的点）的真实舒适度和最大值之差 """ rng = np.random.default_rng(rep) I, J, Y = [], [], [] for _ in range(init): a, b = rng.choice(len(GRID), size=2, replace=False) I.append(a), J.append(b), Y.append(int(rng.random() < 1 / (1 + math.exp(-KAPPA * (F[a] - F[b]))))) out = {} fit = gp_pref_laplace(GRID, np.array(I), np.array(J), np.array(Y), ell=0.1, sf=1.0) for t in range(1, budget + 1): a, b = choose_pair(method, fit, rng) I.append(a), J.append(b) Y.append(int(rng.random() < 1 / (1 + math.exp(-KAPPA * (F[a] - F[b]))))) fit = gp_pref_laplace(GRID, np.array(I), np.array(J), np.array(Y), ell=0.1, sf=1.0) if t in (5, 10, 20, 40): out[t] = F.max() - F[int(np.argmax(fit.mean))] return out def main(reps=40): print_header("偏好贝叶斯优化：五种挑选比较对的方法") print(f"真实舒适度：x ≈ 0.7 处高峰（1.0），x ≈ 0.2 处矮峰（0.5）；人的回答 P(a 更舒服) = σ({KAPPA:.0f}(f(a) − f(b)))。") print(f"先随机问 3 对，再问 40 对。{reps} 次实验的平均：“推荐的步态”（后验均值最高的 x）的舒适度离最大值还差多少，") print("以及差距 < 0.05（找到了高峰）的实验比例：\n") checks = [5, 10, 20, 40] print(f" {pad('方法', 24)}" + "".join(pad(f"{c} 个问题后", 12, "right") for c in checks) + pad("找到高峰", 10, "right")) for method, name in METHODS.items(): res = [run(method, rep) for rep in range(reps)] gaps = {c: np.mean([r[c] for r in res]) for c in checks} found = np.mean([r[40] < 0.05 for r in res]) print(f" {pad(name, 24)}" + "".join(f"{gaps[c]:>12.3f}" for c in checks) + f"{found:>10.2f}") print(""" 观察： · EUBO 和“最好 vs 乐观”表现最好（40 个问题后约 85% 的实验找到了高峰）。两者都把比较集中在 “最好的可能在哪里”上：EUBO 挑“两个都可能好、又不太相关”的一对；“最好 vs 乐观”拿当前最好的去挑战 最有希望超过它的点 —— 即使结果常常一边倒，它也很快排除了没希望的区域。 · “最好 vs 最不确定”最差：另一个点总是被选在最没把握的地方（常常是边缘），而那里多半不好， 答案几乎总是“最好的那个赢”，信息很少（第 7 讲的 p(1 − p)）。 · 自我陪练（双 Thompson 采样）在这里和随机一对差不多：问题很少时，后验很宽，两个样本的最大值点 几乎可以落在任何地方 —— 和第 12 讲冷启动里一样，Thompson 类方法的随机探索在预算很小时代价大。 · 随机一对在一维、61 个候选点的问题上也不算太差；维数一高，随机就不行了。 · 实际的外骨骼研究（Tucker 等 2020，CoSpar）用的是“自我陪练”的思路，并加上“让人也可以说两个差不多”、 “按人的反馈调整探索范围”等改进；结论依赖于问题的规模和噪声，这个一维玩具只说明方法之间的差异确实存在。 · 第 6 讲代码 04 演示 3 用的就是“最好 vs 乐观”：当时我们看到它的不确定性带收得慢（比较一边倒）， 这里看到它找最优点却很快 —— 学准整个函数和找到最好的点，是两个不同的目标。""") if __name__ == "__main__": main()