""" ============================================================================= 第 12 讲 · 代码 3：上下文老虎机——LinUCB 与线性汤普森采样 ============================================================================= 对应讲义：lectures/lecture12_thompson_sampling.md 第 6 节 ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- 到目前为止，每个臂只有一个固定的成功率 μ_a。可在推荐系统里，**不同的用户喜欢不同的东西**： 体育迷爱点体育新闻，科技迷爱点科技新闻。只找"平均最好的一篇文章"远远不够。 **上下文老虎机（contextual bandit）**：每一步先看到一个**上下文** x_t（比如用户的兴趣特征）， 再选一个臂，奖励的期望依赖于 (x_t, 臂)。目标是学会"对什么样的用户推荐什么"。 这里用最常见的**线性模型**：臂 a 对用户 x 的点击率 = θ_a · x（θ_a 未知）。 - LinUCB（Li et al., 2010）：用岭回归估计 θ_a，再加上一个"不确定性"奖励——线性版的 UCB。 - 线性汤普森采样（Agrawal & Goyal, 2013）：从 θ_a 的（近似）后验中抽样——线性版的汤普森采样。 本文件的环境（每次实验随机生成一个）： - 5 个兴趣话题；每个用户的兴趣 x 是这 5 个话题上的一个比例（加起来等于 1），大多数用户集中在一两个话题上； - 8 篇文章，每篇主要属于一个话题：对该话题的"纯粹"读者点击率 30%～60%，对其他读者 0%～10%； - 用户 x 点击文章 a 的概率 = θ_a · x（恰好是线性的）。 本文件： 1. 上下文有多重要：不看用户、只推荐"平均最好"的文章，最多能达到多少点击率？ 2. 各种算法的遗憾：不看上下文的 UCB1、按用户类型分组的 UCB1、ε-贪心、LinUCB、线性汤普森采样 3. 探索强度的影响：不探索（贪心）同样会卡住；理论上的探索系数偏保守 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- cd code/lecture12 python 03_contextual_linucb.py （约 10 秒） """ import numpy as np from bayes_bandits import pad np.set_printoptions(precision=3, suppress=True) D = 5 # 上下文的维数（兴趣话题数） K = 8 # 臂数（文章数） # ============================================================================= # 第 1 部分：环境 # ============================================================================= def make_environments(n_runs, T, seed): """ 为 n_runs 次独立实验各生成一个环境，并**预先**生成所有的用户和随机数。 【返回】字典： "theta" : 形状 (n_runs, K, D)，每篇文章对每个"纯粹话题读者"的点击率 "X" : 形状 (n_runs, T, D)，每一步到来的用户的兴趣向量 "U" : 形状 (n_runs, T)，每一步用来决定"是否点击"的均匀随机数 【为什么预先生成？】所有算法用**同一批用户、同一批随机数**（"公共随机数"）： 如果某个算法在第 t 步选了文章 a，它是否被点击由 U[t] < θ_a·x_t 决定。 这样不同算法之间的差异只来自它们的选择，而不是运气，比较更公平、所需重复次数更少。 【Python 小知识】rng.dirichlet(alpha, size) 生成"加起来等于 1"的随机比例向量； alpha 越小，生成的向量越集中在少数几个分量上。 """ rng = np.random.default_rng(seed) theta = rng.uniform(0.0, 0.1, size=(n_runs, K, D)) main_topic = rng.integers(D, size=(n_runs, K)) # 每篇文章的主话题 r_idx, k_idx = np.meshgrid(np.arange(n_runs), np.arange(K), indexing="ij") theta[r_idx, k_idx, main_topic] = rng.uniform(0.3, 0.6, size=(n_runs, K)) X = rng.dirichlet(0.3 * np.ones(D), size=(n_runs, T)) U = rng.random((n_runs, T)) return {"theta": theta, "X": X, "U": U} def click_probs(theta, x): """每篇文章对当前用户的点击率 p[r, a] = θ_{r,a} · x_r。形状 (n_runs, K)。""" return np.einsum("rkd,rd->rk", theta, x) # 【Python 小知识】np.einsum 用下标字符串描述"哪些维度相乘、哪些维度求和"： # "rkd,rd->rk" 表示：对每个 r、k，把 theta[r,k,:] 和 x[r,:] 逐元素相乘后在 d 上求和。 # ============================================================================= # 第 2 部分：算法（全部向量化：同时处理 n_runs 次实验） # ============================================================================= class LinearModels: """ 为每次实验、每个臂各维护一个岭回归模型： A_a = λI + Σ x xᵀ（只累加选了臂 a 的那些步）， b_a = Σ r·x， θ̂_a = A_a⁻¹ b_a 【为什么直接维护 A⁻¹？】每一步都要用 A⁻¹，如果每次都重新求逆太慢。 每加一个数据点，A 变成 A + x xᵀ，它的逆可以用 **Sherman–Morrison 公式** 在 O(d²) 内更新： (A + x xᵀ)⁻¹ = A⁻¹ - (A⁻¹ x)(A⁻¹ x)ᵀ / (1 + xᵀ A⁻¹ x) 【λ 的作用】正则化：没有数据时 A = λI 仍然可逆；λ 越大，θ̂ 越被"拉向 0"。 【两个关键量】对当前用户 x： mean_a = θ̂_a · x ——预测的点击率 width_a = sqrt(xᵀ A_a⁻¹ x) ——预测的不确定性：x 所在的方向上数据越少，它越大 """ def __init__(self, n_runs, lam=1.0): self.A_inv = np.tile(np.eye(D) / lam, (n_runs, K, 1, 1)) # 形状 (n_runs, K, D, D) self.b = np.zeros((n_runs, K, D)) self.rows = np.arange(n_runs) def predict(self, x): """返回 (mean, width)，形状都是 (n_runs, K)。""" theta_hat = np.einsum("rkij,rkj->rki", self.A_inv, self.b) mean = np.einsum("rki,ri->rk", theta_hat, x) width = np.sqrt(np.maximum(np.einsum("ri,rkij,rj->rk", x, self.A_inv, x), 0.0)) return mean, width def update(self, a, x, r): """把 (x, r) 加入每次实验中被选中的臂 a 的模型。""" Ainv = self.A_inv[self.rows, a] # (n_runs, D, D) u = np.einsum("rij,rj->ri", Ainv, x) # A⁻¹ x denom = 1.0 + np.einsum("ri,ri->r", x, u) self.A_inv[self.rows, a] = Ainv - np.einsum("ri,rj->rij", u, u) / denom[:, None, None] self.b[self.rows, a] += r[:, None] * x def run(algo, env, T, param=None, seed=0, checkpoints=(1000, 10000)): """ 在 env 描述的 n_runs 个环境上运行一种算法 T 步，返回在 checkpoints 时刻的平均累积遗憾， 以及最后 1000 步的平均每步遗憾。 【参数】algo 取值： "ucb1" : 不看上下文的 UCB1（把每篇文章当成一个普通的臂），param = c "cluster_ucb" : 先把用户按"最感兴趣的话题"分成 D 类，每类单独跑一个 UCB1，param = c "egreedy" : 岭回归 + ε-贪心，param = ε（ε = 0 即贪心） "linucb" : LinUCB，param = α（不确定性奖励的系数） "lints" : 线性汤普森采样，param = v（后验的尺度） 【遗憾】相对于"知道 θ、对每个用户都选点击率最高的文章"的理想策略。 """ rng = np.random.default_rng(seed) theta, X, U = env["theta"], env["X"], env["U"] n_runs = theta.shape[0] rows = np.arange(n_runs) models = LinearModels(n_runs) counts = np.zeros((n_runs, D, K)) # 给 ucb1 / cluster_ucb 用的计数（ucb1 只用第 0 类） sums = np.zeros((n_runs, D, K)) cum = np.zeros(n_runs) record, last = [], np.zeros(n_runs) for t in range(T): x = X[:, t] p = click_probs(theta, x) if algo in ("ucb1", "cluster_ucb"): g = np.argmax(x, axis=1) if algo == "cluster_ucb" else np.zeros(n_runs, dtype=int) n, s = counts[rows, g], sums[rows, g] # (n_runs, K) tot = np.maximum(n.sum(axis=1, keepdims=True), 1) score = np.where(n > 0, s / np.maximum(n, 1) + np.sqrt(param * np.log(tot) / np.maximum(n, 1)), np.inf) score = score + 1e-9 * rng.random(score.shape) # 随机打破平局 else: mean, width = models.predict(x) if algo == "egreedy": score = mean + 1e-9 * rng.random(mean.shape) elif algo == "linucb": score = mean + param * width elif algo == "lints": # θ̃_a ~ N(θ̂_a, v² A_a⁻¹) ⇒ θ̃_a·x ~ N(mean_a, v² width_a²)：只需要对"得分"抽样 score = mean + param * width * rng.standard_normal(mean.shape) else: raise ValueError(algo) a = np.argmax(score, axis=1) if algo == "egreedy": explore = rng.random(n_runs) < param a = np.where(explore, rng.integers(K, size=n_runs), a) r = (U[:, t] < p[rows, a]).astype(float) step_regret = p.max(axis=1) - p[rows, a] cum += step_regret if t >= T - 1000: last += step_regret if algo in ("ucb1", "cluster_ucb"): counts[rows, g, a] += 1 sums[rows, g, a] += r else: models.update(a, x, r) if t + 1 in checkpoints: record.append(cum.mean()) return record, last.mean() / 1000 def main(): T, n_runs = 10000, 100 env = make_environments(n_runs, T, seed=0) theta, X = env["theta"], env["X"] # --------------------------------------------------------------------- print("=" * 78) print("演示 1：上下文有多重要？（100 个随机环境的平均）") print("=" * 78) P = np.einsum("rkd,rtd->rtk", theta, X) # 每一步、每篇文章的点击率 oracle = P.max(axis=2).mean() best_fixed = P.mean(axis=1).max(axis=1).mean() # 每个环境中"平均点击率最高的那一篇" random_ = P.mean() print(f" 随机推荐 平均点击率 {random_:.4f}") print(f" 永远推荐\"平均最好\"的那一篇文章 平均点击率 {best_fixed:.4f} ← 任何不看用户的算法的上限") print(f" 按用户推荐最合适的文章（理想） 平均点击率 {oracle:.4f}") print(f" → 看上下文能把点击率从 {best_fixed:.3f} 提高到 {oracle:.3f}（{oracle / best_fixed:.1f} 倍）。") print(" 不看上下文的算法，即使完美地找到了\"平均最好\"的文章，每一步也要比理想策略少") print(f" {oracle - best_fixed:.3f} 的点击率——相对于理想策略，它的遗憾是**线性**的。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print(f"演示 2：各种算法的遗憾（T = {T}，100 个随机环境的平均；遗憾相对于理想的按用户推荐）") print("=" * 78) algos = [("UCB1（不看用户，c = 0.5）", "ucb1", 0.5), ("按用户类型分组的 UCB1", "cluster_ucb", 0.5), ("岭回归 + ε-贪心（ε = 0.05）", "egreedy", 0.05), ("岭回归 + 贪心（ε = 0）", "egreedy", 0.0), ("LinUCB（α = 0.5）", "linucb", 0.5), ("线性汤普森采样（v = 0.2）", "lints", 0.2)] print(" " + pad("算法", 30) + " R(1000) R(10000) 最后 1000 步每步遗憾") for name, algo, param in algos: rec, last = run(algo, env, T, param, seed=1) print(" " + pad(name, 30) + f"{rec[0]:8.1f} {rec[1]:8.1f} {last:.4f}") print(" 观察：") print(" - 不看用户的 UCB1：遗憾线性增长，最后每步遗憾停在演示 1 算出的\"上下文的价值\"附近——它找到了") print(" 平均最好的文章，但这对个性化推荐远远不够。") print(" - 按用户类型分组：把每类用户当成一个独立的老虎机，好了很多；但它把\"兴趣混合\"的用户粗暴地归为一类，") print(" 而且各类之间不共享任何信息。") print(" - 线性模型：所有用户共享同一组参数 θ_a，一个用户的反馈能帮助预测**所有**用户的点击率，") print(" 连从没见过的兴趣组合也能预测——这就是\"泛化\"。LinUCB 和线性汤普森采样的遗憾最小。") print(" - 但线性模型本身不解决探索问题：岭回归 + 贪心的遗憾是 LinUCB 的 6 倍多，比粗糙的\"分组 UCB1\"还差得多，") print(" 最后每步遗憾停在 0.1 左右——它没有持续进步，而是卡住了。") print(" 原因和普通老虎机中的贪心一样：一开始所有文章的预测都是 0，某篇文章先被点击，它的预测变成正数，") print(" 之后这类用户就一直被推荐它；其他文章从来没被推荐过，预测永远是 0，也就永远没有机会。") print(" - ε-贪心能避免卡住，但探索是盲目的：已经很确定的用户和文章，也在以 ε 的概率随机推荐。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 3：探索强度的影响（T = 10000，100 个随机环境的平均）") print("=" * 78) print(" " + pad("算法", 30) + " R(1000) R(10000) 最后 1000 步每步遗憾") for name, algo, param in ([(f"LinUCB（α = {a}）", "linucb", a) for a in [0.0, 0.1, 0.5, 1.0, 2.0]] + [(f"线性汤普森采样（v = {v}）", "lints", v) for v in [0.05, 0.2, 0.5, 1.0]]): rec, last = run(algo, env, T, param, seed=1) print(" " + pad(name, 30) + f"{rec[0]:8.1f} {rec[1]:8.1f} {last:.4f}") print(" 观察：") print(" - α = 0 就是贪心（卡住）；α 太小，探索不够；α 太大，探索过多。这与第 11 讲 UCB 的 c 完全一样。") print(" - LinUCB 论文中理论推荐的 α = 1 + sqrt(ln(2/δ)/2)，δ = 0.05 时约为 2.4——实际表现偏保守，") print(" 较小的 α（0.5～1）更好。线性汤普森采样的理论尺度 v 更大（与维数、置信度有关），实践中也要小得多。") print(" - 不确定性奖励 α·sqrt(xᵀA⁻¹x) 的含义：在当前用户的兴趣方向上，这篇文章的数据有多少。") print(" 一个从没给体育迷推荐过的文章，对体育迷的不确定性很大，即使它已经被科技迷点击过很多次。") if __name__ == "__main__": main()