""" ============================================================================= 第 5 讲 · 代码 1：线性价值函数逼近——梯度 MC 与半梯度 TD ============================================================================= 对应讲义：lectures/lecture05_function_approximation_dqn.md 第 1–4 节 ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- 前几讲用一张表给每个状态存一个价值。状态有 1000 个时表还存得下， 但如果每个状态只被访问几次，每个格子的估计都很不准；而且相邻状态的价值明明很相似， 表格却把它们当成毫无关系的东西——**没有泛化**。 函数逼近的想法：用一个带参数 w 的函数来表示价值： V̂(s; w) = w · x(s) = Σ_i w_i x_i(s) （线性函数逼近） 其中 x(s) 是状态 s 的**特征向量**。参数的个数（比如 10 个）远少于状态数（1000 个）， 更新一个状态的估计会同时影响"特征相似"的其他状态——这就是泛化。 ----------------------------------------------------------------------------- 【环境：1000 状态随机游走（Sutton & Barto 例 9.1）】 ----------------------------------------------------------------------------- - 状态 1~1000（代码中 0~999），从中间（500）出发 - 每一步：等概率向左或向右，跳跃距离在 1~100 之间均匀随机 - 跳出左边界：回合结束，奖励 -1；跳出右边界：回合结束，奖励 +1；其他奖励 0 - γ = 1 真实价值可以用第 1 讲的解析法精确算出（解一个 1000 × 1000 的线性方程组）。 ----------------------------------------------------------------------------- 【本文件的实验】 ----------------------------------------------------------------------------- 1. 状态聚合（把 1000 个状态分成 10 组，每组共享一个参数）+ 梯度 MC 2. 同样的特征 + 半梯度 TD(0)：TD 收敛到一个**不同的、误差更大**的解（TD 不动点） 3. 特征的选择：多项式基、傅里叶基 4. "表格"只是线性函数逼近的特例（独热特征） ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python 01_linear_vfa_random_walk.py （约 30 秒） """ import numpy as np N_STATES = 1000 START = 500 MAX_JUMP = 100 # ============================================================================= # 第 1 部分：环境与真实价值 # ============================================================================= def generate_episode(rng): """ 生成一个回合。返回 (states, final_reward)：途经的状态列表，以及最终奖励（±1）。 【为什么只返回最终奖励？】这个环境中间的奖励都是 0，γ = 1， 所以回合中每个状态的回报 G_t 都等于最终奖励。 """ s = START states = [] while True: states.append(s) jump = rng.integers(1, MAX_JUMP + 1) # 1~100 s = s + jump if rng.random() < 0.5 else s - jump if s < 0: return states, -1.0 if s >= N_STATES: return states, 1.0 def build_model(): """ 用解析法（第 1 讲 8.7 节）精确计算真实价值和状态分布，只用于**评价**学习结果。 【思路】对每个状态 s，向左或向右各跳 1~100 步，每种情况概率 1/200。 跳出去的部分贡献奖励（期望奖励 R(s)），跳到内部的部分贡献转移概率 P(s'|s)。 真实价值：V = (I - P)^{-1} R。 【状态分布 μ】 从起点出发，每个状态在一个回合中被访问的期望次数 η 满足 η = e_start + P^T η （"来到 s 的次数" = 起点贡献 + 从其他状态转移过来的次数） 解出 η 后归一化：μ(s) = η(s) / Σ η。这就是策略下"各状态被访问的比例"。 【返回】(V_true, mu) """ P = np.zeros((N_STATES, N_STATES)) R = np.zeros(N_STATES) for s in range(N_STATES): for jump in range(1, MAX_JUMP + 1): for s_next in (s + jump, s - jump): if s_next < 0: R[s] += -1.0 / (2 * MAX_JUMP) elif s_next >= N_STATES: R[s] += 1.0 / (2 * MAX_JUMP) else: P[s, s_next] += 1.0 / (2 * MAX_JUMP) I = np.eye(N_STATES) V_true = np.linalg.solve(I - P, R) eta = np.linalg.solve((I - P).T, I[START]) return V_true, eta / eta.sum(), P, R def best_linear_fit(X, V_true, mu): """ 在给定特征下，**理论上最好**的线性逼近：最小化 Σ_s μ(s) [V_true(s) - w·x(s)]² 的 w。 【思路】这是一个加权最小二乘问题。把每一行乘以 sqrt(μ(s)) 就变成普通最小二乘， 用 np.linalg.lstsq 求解（比直接解正规方程数值上更稳定）。 【注意】它需要知道真实价值 V_true，所以不是一个学习算法—— 它回答的是"这组特征**最多**能做到多好"（表示能力的上限）。 """ sq = np.sqrt(mu) return np.linalg.lstsq(X * sq[:, None], V_true * sq, rcond=None)[0] def td_fixed_point(X, P, R, mu): """ 线性半梯度 TD(0) 的收敛点（"TD 不动点"），可以解析求出（讲义 4.3 节）： A w = b， A = X^T D (I - γP) X， b = X^T D R， D = diag(μ)，γ = 1 【直觉】TD 的期望更新量 E[(r + γ w·x(s') - w·x(s)) x(s)] = b - A w，令它为 0 就是不动点。 """ D = np.diag(mu) A = X.T @ D @ (np.eye(N_STATES) - P) @ X b = X.T @ D @ R return np.linalg.solve(A, b) # ============================================================================= # 第 2 部分：特征 # ============================================================================= def aggregation_features(n_groups=10): """ 状态聚合：把 1000 个状态按顺序分成 n_groups 组，特征是"属于哪一组"的独热向量。 【返回】形状 (N_STATES, n_groups) 的特征矩阵 X，X[s] 就是 x(s)。 【直觉】同一组内的所有状态共享同一个价值估计 w_g，价值函数变成"台阶"形状。 【边界】n_groups = N_STATES 时，每个状态自成一组——就是表格方法。 """ X = np.zeros((N_STATES, n_groups)) group = np.arange(N_STATES) * n_groups // N_STATES # 状态 s 属于第几组 X[np.arange(N_STATES), group] = 1.0 return X def polynomial_features(order): """ 多项式基：x(s) = [1, u, u², ..., u^order]，其中 u = s / 999 ∈ [0, 1]。 【为什么要把 s 缩放到 [0, 1]？】如果直接用 s（0~999），s^5 会大到 10^15， 不同特征的数值尺度相差巨大，梯度下降几乎无法工作。**把输入缩放到合适的范围** 是使用函数逼近时最重要的工程细节之一（神经网络也一样）。 """ u = np.arange(N_STATES) / (N_STATES - 1) return np.stack([u ** i for i in range(order + 1)], axis=1) # 【Python 小知识】np.stack(列表, axis=1)：把若干个一维数组按列并排，组成二维矩阵。 def fourier_features(order): """ 傅里叶基：x_i(s) = cos(i π u)，i = 0, 1, ..., order，u = s / 999。 【直觉】用不同频率的余弦波叠加出价值函数的形状。低频项刻画整体趋势，高频项刻画细节。 傅里叶基在实践中通常比多项式基表现更好（数值上更稳定）。 """ u = np.arange(N_STATES) / (N_STATES - 1) return np.stack([np.cos(i * np.pi * u) for i in range(order + 1)], axis=1) # ============================================================================= # 第 3 部分：梯度 MC 与半梯度 TD # ============================================================================= def gradient_mc(X, n_episodes, alpha, rng): """ 梯度蒙特卡洛：对每个回合中的每个状态 s_t，用回报 G_t 当作目标，做一步随机梯度下降。 【推导】（讲义 3 节） 想让 V̂(s; w) 接近 V^π(s)，目标函数是均方误差 J(w) = E[(V^π(s) - V̂(s; w))²]。 V^π(s) 不知道，用它的无偏样本 G_t 代替。对 w 求梯度并沿负梯度走一小步： w ← w + α [G_t - V̂(s_t; w)] ∇_w V̂(s_t; w) 线性情况下 ∇_w V̂(s; w) = x(s)，所以： w ← w + α [G_t - w·x(s_t)] x(s_t) 【和表格方法的关系】如果 x(s) 是独热向量，上式只改变 s 自己那个参数， 退化为第 3 讲的 V(s) ← V(s) + α(G - V(s))。 【参数】 X : 特征矩阵 (N_STATES, d)。 n_episodes : 回合数。 alpha : 学习率。**特征越"大"（比如特征值之和越大），α 就要越小**，否则会发散。 【返回】 (w, visit_counts)：学到的参数，以及每个状态被访问的次数（用来估计状态分布 μ）。 """ w = np.zeros(X.shape[1]) visits = np.zeros(N_STATES) for _ in range(n_episodes): states, G = generate_episode(rng) for s in states: x = X[s] w += alpha * (G - w @ x) * x visits[s] += 1 return w, visits def semi_gradient_td0(X, n_episodes, alpha, rng): """ 半梯度 TD(0)：用 TD 目标 r + γ V̂(s'; w) 代替 G_t。 w ← w + α [r + γ w·x(s') - w·x(s)] x(s) 【为什么叫"半梯度"？】 TD 目标里也含有 w（因为 V̂(s'; w) 依赖 w），但我们求梯度时**假装它是常数**， 只对 V̂(s; w) 求导。所以这不是任何目标函数的真正梯度，只是"一半"的梯度。 后果：它不再保证收敛到均方误差最小的解，而是收敛到"TD 不动点"（实验 2）。 更糟的是，在某些情况下它会发散（03 文件的"致命三角"）。 【终止】终止状态的价值为 0，TD 目标就是最终奖励。 """ w = np.zeros(X.shape[1]) for _ in range(n_episodes): states, final_reward = generate_episode(rng) for t, s in enumerate(states): x = X[s] if t + 1 < len(states): target = 0.0 + w @ X[states[t + 1]] # 中间奖励为 0，γ = 1 else: target = final_reward # 最后一步：跳出边界 w += alpha * (target - w @ x) * x return w def msve(V_hat, V_true, mu): """ 均方价值误差 VE(w) = Σ_s μ(s) [V_true(s) - V_hat(s)]²，开根号后返回（RMS 形式）。 【为什么要按 μ(s) 加权？】参数少于状态数时，不可能让所有状态都完全准确， 必须"取舍"。自然的取舍是：**经常访问的状态更重要**。μ(s) 就是策略下 各状态被访问的比例（第 1 讲 7.6 节的"状态分布"）。 """ return float(np.sqrt(np.sum(mu * (V_true - V_hat) ** 2))) def print_curve(label, V, step=100): """每隔 step 个状态打印一个价值，粗略地展示价值函数的形状。""" print(f" {label:<18s}" + "".join(f"{V[s]:6.2f}" for s in range(0, N_STATES, step)) + f"{V[-1]:6.2f}") # ============================================================================= # 第 4 部分：演示 # ============================================================================= def main(): rng = np.random.default_rng(0) print("正在精确计算 1000 个状态的真实价值和状态分布……") V_true, mu, P, R = build_model() print(" 状态 " + "".join(f"{s + 1:6d}" for s in range(0, N_STATES, 100)) + f"{N_STATES:6d}") print_curve("真实价值", V_true) # --------------------------------------------------------------------- print("\n" + "=" * 78) print("实验 1：状态聚合（10 组）+ 梯度 MC（20000 个回合，α = 5e-4）") print("=" * 78) X_agg = aggregation_features(10) w_mc, _ = gradient_mc(X_agg, 20_000, alpha=5e-4, rng=rng) V_mc = X_agg @ w_mc w_best = best_linear_fit(X_agg, V_true, mu) print_curve("梯度 MC（学到的）", V_mc) print_curve("最优聚合解", X_agg @ w_best) print(f"\n RMS 价值误差：梯度 MC {msve(V_mc, V_true, mu):.4f}；该特征下的理论最优 {msve(X_agg @ w_best, V_true, mu):.4f}") print(" → 价值函数变成了\"台阶\"：同一组的 100 个状态共享一个值。") print(" 梯度 MC 的收敛点就是均方误差最小的解（每组取组内真实价值按 μ 加权的平均）；") print(" 学到的结果还差一点，是因为常数步长带来的噪声（第 3 讲 2.3 节）：α 越小、回合越多，越接近最优。") print(f" 各组的访问频率 μ：{np.round([mu[g * 100:(g + 1) * 100].sum() for g in range(10)], 3)}") print(" （中间的组访问得多，两端的组访问得少——误差的\"权重\"也就不同。）") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("实验 2：同样的特征 + 半梯度 TD(0)（20000 个回合，α = 1e-3）") print("=" * 78) w_td = semi_gradient_td0(X_agg, 20_000, alpha=1e-3, rng=rng) V_td = X_agg @ w_td w_fix = td_fixed_point(X_agg, P, R, mu) print_curve("半梯度 TD（学到的）", V_td) print_curve("TD 不动点（解析）", X_agg @ w_fix) print_curve("最优聚合解", X_agg @ w_best) print(f"\n RMS 价值误差：半梯度 TD {msve(V_td, V_true, mu):.4f}；TD 不动点 {msve(X_agg @ w_fix, V_true, mu):.4f}；" f"最优 {msve(X_agg @ w_best, V_true, mu):.4f}") print(" → 半梯度 TD 收敛到的不是均方误差最小的解，而是\"TD 不动点\"，它在两端\"更平\"，") print(" 误差约为最优解的 2 倍。这是\"自举 + 函数逼近\"的固有偏差，不是没训练够：") print(" 学到的结果已经非常接近解析算出的不动点。") print(" （线性、同策略时可以证明：TD 不动点的误差 ≤ 最优误差 / (1-γ)。γ = 1 时这个界没有意义。）") print(" TD 的优势在别处：更新早、方差小、学得快（第 3 讲）。这里我们看的是充分训练后的收敛点。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("实验 3：特征的选择——\"能表示多好\" vs \"能学得多好\"（梯度 MC，10000 个回合，3 个种子平均）") print("=" * 78) print(" 特征 参数数 理论最优误差 学到的误差 条件数 cond(XᵀDX)") for name, make, alpha in [("多项式", polynomial_features, 5e-5), ("傅里叶", fourier_features, 2e-5)]: for order in [5, 10, 20]: X = make(order) best = msve(X @ best_linear_fit(X, V_true, mu), V_true, mu) learned = np.mean([msve(X @ gradient_mc(X, 10_000, alpha, np.random.default_rng(seed))[0], V_true, mu) for seed in range(3)]) cond = np.linalg.cond(X.T @ (mu[:, None] * X)) print(f" {name}基 {order:2d} 阶 {order + 1:3d} {best:.4f} {learned:.4f} {cond:10.1e}") print(" 观察：") print(" - \"理论最优误差\"一栏：两种特征都能把价值函数表示得非常准（阶数越高越准），") print(" 多项式基甚至比傅里叶基更准。**表示能力不是问题。**") print(" - \"学到的误差\"一栏：多项式基学得很差，比理论最优差几百倍，而且阶数越高也没有改善。") print(" - 原因在最后一栏：多项式特征（u, u², ..., u^20）彼此高度相似，XᵀDX 的条件数高达 10^7～10^17，") print(" 梯度下降在这种\"又窄又长的山谷\"里几乎走不动。傅里叶特征彼此接近正交，条件数只有 10 左右。") print(" - 傅里叶基学到的误差（约 0.03）与理论最优的差距，主要来自常数步长的噪声，而不是优化困难。") print(" - 教训：选特征时，除了\"能不能表示\"，还要看\"好不好学\"。把输入缩放到合适范围、") print(" 让特征尽量不相关，是函数逼近中最重要的工程细节之一。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("实验 4：表格方法 = 独热特征的线性函数逼近") print("=" * 78) X_tab = aggregation_features(N_STATES) # 1000 组，每组一个状态 print(f" 特征矩阵是单位矩阵吗？{np.array_equal(X_tab, np.eye(N_STATES))}") w_tab, _ = gradient_mc(X_tab, 5_000, alpha=0.01, rng=rng) w_agg5k, _ = gradient_mc(X_agg, 5_000, alpha=5e-4, rng=rng) print(f" 表格型 MC（1000 个参数，5000 回合）RMS 价值误差 = {msve(X_tab @ w_tab, V_true, mu):.4f}") print(f" 状态聚合（10 个参数， 5000 回合）RMS 价值误差 = {msve(X_agg @ w_agg5k, V_true, mu):.4f}") print(" → 表格方法有 1000 个参数、没有近似误差，但每个参数只能从\"自己那个状态\"的数据中学习，") print(" 数据有限时估计噪声很大；10 个参数的聚合方法有近似误差（台阶），但每个参数能利用") print(" 100 个状态的数据。**泛化用\"偏差\"换来了\"方差\"的大幅降低**——又一次偏差–方差权衡。") if __name__ == "__main__": main()