逆强化学习示例：从专家轨迹恢复奖励函数的最大熵实现

简介：面向逆向强化学习研究者和使用Java语言的开发人员，这是一套用于逆向强化学习教程的示例代码工程。项目在BURLAP强化学习代码库基础上进行了定制化改造，补充了逆向强化学习框架所需的核心模块，可直接导入集成开发环境运行，适合用于算法复现、教学演示和二次扩展。压缩包内共包含四百九十七个文件，其中绝大多数为Java源文件，数量达四百八十四个，另有少量XML配置、Maven构建描述、打包依赖的Jar包以及版权说明和项目说明文档，整体体积约二点二九兆字节，结构清晰，便于按模块查阅。该资源已有一千零四十六人浏览学习。示例中涵盖对抗式学徒学习、栅格世界等典型实验场景，同时提供多智能体性能绘图等辅助工具，可帮助读者理解逆向强化学习的整体流程、掌握BURLAP框架的改造思路，并在此基础上结合自己的课题进行算法设计与实验验证。

1. 逆强化学习示例代码在解决什么：从专家轨迹里恢复奖励函数

你在做模仿学习的时候，可能撞上过一个怪现象：行为克隆的策略在训练环境里跑得挺顺，换一张地图、换一段路况，整个策略就崩了。原因是它只学会了“照抄动作”，没拿到驱动这些动作的那个奖励函数。逆强化学习（Inverse Reinforcement Learning，IRL）做的正是这件事——从专家轨迹里反推奖励函数。IRLTutorial 这套示例代码，就是为这个目标准备的最小可运行实现。

我最早读最大熵IRL那篇论文时，卡住的不是公式，而是“去哪找一份能跑的代码”。论文里的算法和实际训练循环之间，隔着特征设计、前向规划器、稳定数值计算这一堆工程细节。IRLTutorial 的思路是：用一个小网格世界、一个线性奖励、一套最大熵IRL求解器，把从生成演示数据到训练、到可视化的一条命令链路全打通。适合刚读完论文没处下手的新研究者，也适合要在决策规划里做奖励设计、想先跑通基线的工程师。

这套示例值得花一晚上跑通的原因很简单：IRL 的算法选型、参数设置、收敛判断，全部依赖具体实现，光看公式看不出梯度方向、熵正则、特征期望这些环节是怎么互相咬合的。代码在手，你才能快速试错，知道哪个模块影响了最终恢复出来的奖励长什么样。

2. 为什么示例默认选最大熵IRL：三条路线对比与一个前置检查

2.1 先理解“歧义”：同一批轨迹对应无数个奖励函数

把问题压到最简：你有一组专家轨迹 $\mathcal{D} = {\tau_1, \dots, \tau_N}$，每条轨迹是一串状态-动作对。假设专家在优化某个未知奖励函数 $R(s,a)$，IRL 要做的就是把 $R$ 猜出来。听起来像普通的回归问题，但这里有个绕不开的约束：如果一条轨迹在奖励 $R_1$ 下是最优的，那么在 $R_1$ 基础上给所有轨迹途经的状态统一加一个常数，或者给没被轨迹访问的状态统一减一个惩罚，新奖励 $R_2$ 仍然能让这批轨迹最优。

这就是 IRL 的“歧义性”——满足“专家轨迹更优”这个不等式约束的奖励有无数个。算法比拼的关键，就是用哪种方式从这个解空间里挑出一个具体答案。世代沿用的三条路分别是：最大边际规划（MMP）挑“让专家轨迹和其他策略之间间隔最大”的那个；Abbeel 和 Ng 的学徒学习挑“特征期望匹配且间隔最大”的那个；最大熵IRL则换了个思路——在所有与专家轨迹特征期望一致的轨迹分布里，选熵最大的那个。

为什么最大熵在直觉上最稳？因为“匹配特征期望”这个约束只保证了轨迹分布的一阶统计量一致，没有理由假设专家轨迹分布比其他分布更集中。最大熵分布是对未知信息介入最少的分布。它导出的轨迹概率恰好是指数族形式 $p(\tau) \propto \exp(\sum_{t} \theta^\top f(s_t, a_t))$，对应策略自动变成 softmax 形式，后续梯度计算就有闭式结构，这也是示例代码里温度系数（熵正则）的来源。

在开始训练之前，我建议先做一个前置检查：计算专家轨迹的特征期望，再计算均匀随机策略的特征期望，看它们之间的距离。这个检查能提前暴露“演示数据没有信息量”的问题，省得把锅甩给算法。

import numpy as np demo_exp = np.mean(demo_features, axis=0) # 从轨迹统计专家特征计数 uniform_policy = np.ones((S, A)) / A # 均匀策略 rand_exp = compute_expected_feature(uniform_policy, sa_features, transitions) dist = np.linalg.norm(demo_exp - rand_exp) print("feature gap:", dist, "demo norm:", np.linalg.norm(demo_exp))

如果dist比demo_exp的范数小一个数量级，说明专家轨迹和随机探索的访问分布差别不大，轨迹里根本没携带偏好信息。这时候换更复杂的 IRL 算法也恢复不出有意义的奖励，该做的是重新设计特征、增加演示轨迹数量、或者调大专家轨迹里的行为差异性。

2.2 三条算法路线的选型对比：为什么最大熵是默认值

对于要在教程示例里落地的主算法，我更看重实现复杂度和对噪声演示的容忍度之间的平衡。下面这张表是我自己选型时的对比基准：

从工程角度解释一下这张表。MMP 用 0-1 间隔做优化目标，对轨迹里的异常点非常敏感：一条混进来的“错误轨迹”可能让间隔约束被严重违背，恢复出的奖励权重直接跑偏。学徒学习用特征期望做匹配，相当于把轨迹压缩成特征均值，抗噪能力好一些，但它每轮仍然要解一个硬策略，这个“内层强化学习”一旦收敛慢，整体训练时间会成倍上涨。

最大熵IRL 多出来的工程量只是一次软价值迭代和一次期望特征计算，得到的策略是概率分布形式。这样即使演示轨迹不是严格最优，算法也能用概率去拟合，而不是强行要求轨迹满足硬约束。梯度形式也有记忆点：最大熵 IRL 的参数更新方向是“当前策略的特征期望 − 专家轨迹的特征期望”，这一点在第 4 章的代码里会反复出现。

选择最大熵IRL做示例默认算法，还有一个实际原因：它的超参数和调试信号足够直观。温度系数控制策略平滑度，梯度范数能直接反映特征匹配进度，特征期望差值能可视化地展示在状态空间上。这些特性让教程代码不会变成黑匣子，新手能一眼看出“现在学到哪一步了”。

3. 把IRLTutorial跑通：环境、演示数据、训练与可视化的完整链路

3.1 拿到源码先看目录：四个模块各管一件事

我拿到一份示例代码，习惯先看目录结构，而不是直接跑训练脚本。IRLTutorial 的常见组织方式是把四个关注点拆开：环境构建、算法求解、脚本入口、可视化工具。具体分工大致是：

environment/ # 网格世界定义、转移表构建、状态特征生成 algorithms/ # 软价值迭代、期望特征计算、最大熵IRL梯度更新 scripts/ # generate_demos.py / train_irl.py / visualize_reward.py visualization/ # 奖励热力图、轨迹覆盖图绘制

这样拆的好处是环境与前向规划器解耦。后面你要迁移到自己的任务，只需要改environment/下的接口，算法代码基本不用动。示例代码的依赖也很轻——核心只需要 numpy，画图才用 matplotlib，不需要 GPU，CPU 上几十秒就能跑完整个训练流程。

创建环境时我会用 conda 隔离，避免污染系统 Python：

conda create -n irl python=3.9 conda activate irl pip install numpy matplotlib

如果本地已经有 Python 环境，跳过 conda 也没关系，关键是 numpy 版本不要过旧，否则np.einsum和np.linalg.solve在部分场景下行为不一致。没有特殊需求就不要上 PyTorch，示例代码用的是查表式的网格世界，用深度学习框架反而是给调试增加负担。

3.2 生成专家演示数据：先搞清楚“专家”是谁

IRL 的数据是轨迹，不是带回报标签的样本。示例代码里的演示数据生成脚本，通常的做法是：先写死一份真实奖励（比如在网格里布置正负信号），用普通价值迭代求出最优策略，再从策略里按一个小概率随机扰动采样轨迹。扰动的作用是让轨迹不显得“过分完美”，因为真实专家演示总是有噪声的。

python scripts/generate_demos.py \ --env-name gridworld-v0 \ --num-demos 12 \ --expert-noise 0.1 \ --seed 42

--num-demos 12是经验值。少于 6 条时特征期望的方差太大，训练出来的权重反复横跳；多于 30 条对这个小网格世界收益递减，训练时间反而变长。--expert-noise 0.1表示每一步有约 10% 的概率随机替换动作，它让轨迹里保留下“探索”信号，这对后面最大熵IRL 估计策略熵很有用。输出是一个轨迹列表文件，每条轨迹由一串状态编号和动作编号组成，文件体积很小，可以直观打开检查。

注意一个关键点：脚本内部保存了真实奖励R_true，这是为了训练结束后做对照验证。不要把R_true传给训练脚本，否则就是作弊。IRL 的难度正在于只有轨迹，没有回报。如果你要用自己的演示数据，需要先把轨迹里的状态、动作映射成和transition表一致的下标，脚本里通常会留一个load_own_demos.py模板。

3.3 训练与可视化的完整命令：一条命令看到奖励长什么样

生成好演示数据后，训练脚本做三件事：加载轨迹、跑最大熵IRL 更新奖励权重、保存结果到 checkpoint 文件。训练参数先照抄默认值，跑通后再调：

python scripts/train_irl.py \ --algo maxent \ --epochs 150 \ --lr 0.02 \ --entropy-reg 0.3 \ --save checkpoints/maxent_gridworld_v0.npz python scripts/visualize_reward.py \ --ckpt checkpoints/maxent_gridworld_v0.npz \ --savefig figures/reward_heatmap.png

--epochs 150对约 7×7=49 个状态的小网格足够稳定收敛，一般 50 轮后权重方向已经固定，剩下的是让数值更稳定。--lr 0.02是针对线性奖励的经验值：学习率大于 0.1 时权重后期会在最优点附近震荡，小于 0.005 时收敛过慢，你会在训练日志里看到梯度范数一直降不下去。--entropy-reg 0.3是温度系数，比论文默认的 1.0 略低，目地是让恢复出的策略更贴专家行为，而不是过度平滑。

训练结束后，可视化脚本会输出两张图：一张是恢复奖励的热力图，另一张是专家轨迹的覆盖位置。如果算法正常工作，高奖励区域应该和轨迹覆盖密集区域基本重合。这是最快的反馈信号，也是检查“奖励是否恢复对了”的第一道关卡。

4. 核心实现拆到函数级：特征映射、软价值迭代与权重更新

4.1 特征是怎么喂给奖励函数的：线性奖励的设计与代码

最大熵IRL 的一个常见实现是线性奖励：$R(s,a) = \theta^\top f(s, a)$，其中 $f$ 是特征向量，$\theta$ 是待学习的权重。选择线性形式的核心原因是可解释：训练结束后，你能直接读出每个语义特征对奖励的贡献，这比奖励网络的黑匣子友好太多。另一个原因是从梯度角度看，线性奖励让更新公式变成简单的“特征期望之差”，实现代码短且不易出错。

对于网格世界，特征要由用户预先设计。示例代码里常用三维特征：是否到达终点、是否踩到危险区域、到目标的曼哈顿距离（取负值）。构造代码如下：

def build_features(grid): S = grid.size F = 3 features = np.zeros((S, F)) goals = grid.goal_positions for s, (x, y) in enumerate(grid.states): features[s, 0] = 1.0 if (x, y) in goals else 0.0 features[s, 1] = 1.0 if grid.hazard[x, y] else 0.0 features[s, 2] = -abs(x - goals[0][0]) - abs(y - goals[0][1]) features[:, 2] /= np.abs(features[:, 2]).max() return features

这里有一个常见选择：特征取的是“下一状态”的特征，而不是当前状态。因为在网格世界里，动作的价值体现在它把智能体带到哪个状态，$R(s,a)$ 具体实现时用 $f(s')$ 代替 $f(s,a)$，依靠transition表来查。构造状态-动作特征等价于对所有动作做一次索引：

def build_sa_features(features, transitions): S, A = transitions.shape F = features.shape[1] sa_features = np.zeros((S, A, F)) for s in range(S): for a in range(A): ns = transitions[s, a] sa_features[s, a] = features[ns] return sa_features

逻辑说明：transitions[s, a]是执行动作 $a$ 后到达的状态编号，features[ns]取出该状态的特征向量，赋值给当前状态-动作对。参数上要留意两点。第一，特征维度不能太大，网格世界用 3～5 维足够；维度多了但演示轨迹有限，必然出现某些列长期为 0，权重语义无法解释。第二，特征列要尽量做归一化，否则量纲差异悬赏时，同一学习率下部分权重飞了，部分权重还没动。如果你换到连续状态空间，做法类似，只是把“查表取特征”换成径向基函数（RBF）或者一个小型网络编码器，线性权重结构可以先保留。

4.2 前向规划器：软价值迭代的实现与两个关键参数

前向规划器是最大熵IRL 的“内层循环”，它负责回答一个问题：给定当前奖励权重 $\theta$，智能体在环境中会表现出什么样的策略分布？最大熵假设下，这个策略不是硬 max，而是 softmax 形式，所以要用软价值迭代而非普通价值迭代，即把 Bellman 最优方程里的 max 替换成 log-sum-exp。

def soft_value_iteration(sa_features, theta, transitions, gamma=0.95, tau=0.3, max_iters=200, tol=1e-5): S, A, F = sa_features.shape V = np.zeros(S) for _ in range(max_iters): rewards = sa_features.dot(theta) # (S, A) 当前奖励 Q = rewards + gamma * V[transitions] # (S, A) 动作价值 m = np.max(Q, axis=1, keepdims=True) # 数值稳定 V_new = tau * (m + np.log(np.sum(np.exp((Q - m) / tau), axis=1))) if np.max(np.abs(V_new - V)) < tol: V = V_new break V = V_new Q = rewards + gamma * V[transitions] logits = Q / tau logits -= logits.max(axis=1, keepdims=True) # 稳定 softmax：先减最大值再指数 policy = np.exp(logits) policy /= policy.sum(axis=1, keepdims=True) return V, policy

这段代码有三个地方值得细看。第一，Q = rewards + gamma * V[transitions]是向量化写法：transitions 形状为 (S, A)，V[transitions]直接得到每个 (s, a) 的下一状态价值，不必写 for 循环。第二，V_new的 log-sum-exp 必须先减最大值m再取指数，否则当某个动作优势远大于其他动作时，exp(Q / tau)在 tau 很小时会出现 float 溢出，训练日志里突然全是 NaN。第三，最后的策略计算用稳定 softmax：先减最大值、再指数、再归一化。数值稳定性是这节最容易翻车的地方，我见过太多实现因为这里少写一行，把问题归咎于算法不收敛。

参数上，tau是熵正则温度：越大策略越均匀，模拟专家“犯错”的概率越高；越小策略越贪心，越接近硬最优。示例代码的网格世界用 0.2～0.5 比较合适。gamma是折扣因子，0.9～0.95 就够；如果拉到 0.99，状态访问频率方程会变得病态，线性方程组求解的精度会下降。判断收敛不能只靠固定迭代次数，输出 Bellman 残差会更直观：如果残差停在某个非零值上不下降，优先检查 transitions 表是否有越界下标。

4.3 反向更新与收敛判断：梯度方向是特征期望之差

有了软策略之后，反向更新的核心就一句话：当前策略的特征期望减专家轨迹的特征期望，沿负梯度方向更新权重。

def compute_expected_feature(policy, sa_features, transitions, gamma=0.95): S, A = policy.shape P_pi = np.zeros((S, S)) for s in range(S): for a in range(A): ns = transitions[s, a] P_pi[s, ns] += policy[s, a] A_mat = np.eye(S) - gamma * P_pi.T b = np.full(S, (1.0 - gamma) / S) # 均匀初始分布 d = np.linalg.solve(A_mat, b) # 折扣状态访问频率 model_exp = np.zeros(sa_features.shape[2]) for s in range(S): model_exp += d[s] * np.einsum('a,af->f', policy[s], sa_features[s]) return model_exp

逻辑说明：先把策略概率写成状态转移矩阵 $P_\pi$，然后解线性方程组 $(I - \gamma P_\pi^\top) d = (1-\gamma)/S \cdot \mathbf{1}$，得到折扣状态访问频率 $d$。最后把每个状态下的期望特征按访问频率加权求和。这里没有用蒙特卡洛采样，是因为小状态空间下解析解的精度远高于采样，且不会引入随机噪声让训练曲线抖动。

主训练循环很短，却包含了最大熵IRL 的全部信息：

theta = np.zeros(F) demo_exp = compute_demo_feature(demos, sa_features) # 统计专家轨迹特征均值 for epoch in range(150): V, policy = soft_value_iteration(sa_features, theta, transitions) model_exp = compute_expected_feature(policy, sa_features, transitions) grad = model_exp - demo_exp theta -= 0.02 * grad theta -= 0.01 * theta # L2 权重衰减 if epoch % 10 == 0: print(f"epoch {epoch:3d} | grad_norm {np.linalg.norm(grad):.4f} | theta {theta}")

理解这个梯度方向很关键：它不是在提升“当前轨迹的回报”，而是在缩小“当前策略访问分布”与“专家轨迹访问分布”的特征差距。当两者特征期望一致时，梯度为零，训练收敛。我见过有人误把这段理解成“让奖励输出变大”，结果在评估时发现奖励函数完全偏离真实语义。theta -= 0.01 * theta是 L2 衰减，作用是抑制某些特征维在样本不足时无约束增长，增强权重的稳定性。

收敛判断有一个很实用的观察矩阵：

5. 避坑：权重发散、轨迹污染、收敛误判的常见原因

5.1 权重变成 NaN：先把数值稳定性补上再谈算法

我最早跑类似实现时，训练到第 40 轮左右 theta 突然全是nan。当时第一反应是学习率太大，调小后照样崩。后来逐行检查才发现问题出在两处：一是np.exp(Q / tau)在 tau 偏小时溢出，二是 softmax 里没有减最大值，中间计算结果先一步超过 float64 上限。这不是算法问题，是数值实现问题。

解决办法就是第 4.2 节的稳定 logsumexp 和稳定 softmax。如果你看到 NaN 出现的位置集中在某个特征维度，再去检查该特征列在演示轨迹里是否始终为 0。零列的反向梯度本身含有除零风险，训练前打印demo_exp能提前排除这个坑。另外把学习率从 0.1 降到 0.02，权重衰减从 0 调到 0.01，这一组组合在大多数网格示例里都能稳住训练过程。

5.2 专家轨迹太少、太干净或太脏，权重都会横跳

当--num-demos低于 6 时，训练日志里的 theta 几乎每 20 轮换一次符号，看起来像“不收敛”。原因是特征期望的统计量在少量样本下方差太大，梯度方向本身就在随机游走。另一个极端是演示轨迹完全没有噪声，所有轨迹都严格贴着最优路线走，此时最大熵的熵正则作用被削弱，恢复出的策略过于自信，在新状态上反而偏离。

处理方式分三条：把演示数量抬到 12 条以上；在生成轨迹时加 0.05～0.1 的随机扰动，保留行为多样性；温度 tau 不要设到 1.0 以上，否则策略平滑到几乎没有信息。我一般会用脚本固定随机种子，并保存演示数据的噪声参数，否则后面复现时这些细节会全部变成玄学。

5.3 软价值迭代不收敛：先查转移表，再调折扣因子

如果软价值迭代的 Bellman 残差在几十轮后停在某个非零值上，我先怀疑的不是迭代次数，而是转移表。transitions数组里只要有一格下标越界，V[transitions]就会悄悄取到错误状态，价值迭代自然是乱的。便宜的检查方法是一行断言：

assert transitions.min() >= 0 and transitions.max() < S, "transition 表越界"

这行断言应该放在每次前向计算之前，成本极低，却能拦住一大类环境定义错误。另一个常见原因是 gamma 设成 0.99。在状态访问频率方程里，gamma 越接近 1，矩阵 $I - \gamma P_\pi^\top$ 越接近奇异，np.linalg.solve的解会变得极其敏感，任何微小策略变化都会让访问频率大幅波动。网格示例用 0.9～0.95 即可，gamma 过高带来的收益远小于它引发的数值副作用。

5.4 loss 在下降但行为翻车：别只盯着训练日志

最迷惑人的情况是：梯度范数在降，特征期望差在缩小，但把训练好的奖励重新放回环境采样，生成轨迹只覆盖了演示轨迹的前半段，后半段完全跑偏。原因是最大熵IRL 的训练目标是轨迹似然，不是行为匹配；当 tau 偏大时，策略分布变得平滑，每个动作的概率都被平均掉，轨迹多样性靠“乱走”来凑，特征期望虽然匹配了，但具体路径结构已经完全变形。

应对方法是加一个与训练目标无关的外部验证指标。我用的是行为距离，即状态访问频率分布的绝对差：

def behavior_distance(demo_trajs, gen_trajs, S): def freq(trajs): f = np.zeros(S) for t in trajs: for s, a in t: f[s] += 1 return f / f.sum() return np.abs(freq(demo_trajs) - freq(gen_trajs)).sum() / 2

这个值越小越好，超过 0.3 基本说明策略结构和专家行为差异很大。训练完不要只看 checkpoint 里的权重数值，必须做一次重放验证。这步验证也是把 IRL 落地到实际任务前最便宜的安全网。

5.5 特征稀疏导致的“假零权重”：没访问过不等于不重要

如果专家轨迹里某个特征对应的状态从未出现，demo_exp的那一维就是 0，梯度对该维几乎没有信号，权重会被 L2 衰减慢慢压到 0。这个 0 并不代表“真实奖励中该特征不重要”，只能说明“当前数据里没有证据”。这会让特征分析变得误导：你看着权重表，以为某个维度无用，实际上只是轨迹没覆盖到。

训练前打印demo_exp，统计全零列的占比。如果超过总维度的三分之一，优先合并相似特征、减少特征维度，或者补充覆盖这些状态的演示轨迹。千万不要在特征稀疏时用全局归一化强行激活该列，那只会在反向传播时引入一个始终有梯度的噪声源。这个检查只要一行打印，却能在你换新任务时节省大量排查时间。

6. 把示例迁移到自己的任务：验证奖励恢复的三个技巧

6.1 验证先用行为距离，再用奖励对齐，最后再用权重符号

从网格世界迁移到自己的任务时，第一件事不是改代码，而是定验证指标。我推荐按顺序做三层验证：先用第 5.4 节的行为距离做快速筛选，行为距离大的模型直接淘汰，不浪费精力看权重；通过后再做奖励对齐——把学习到的奖励和真实奖励（生成演示时用的R_true）做归一化相关计算，相关系数大于 0.7 说明恢复出了正确的相对排序；最后才看权重符号，符号正确足以支撑“哪些区域该给正奖励、哪些该给负奖励”这类定性判断，但不要强求数值完全一致。IRL 本身受歧义性限制，权重数值不是唯一解，过度追求数值对齐反而会误导后续调参。

6.2 迁移时只动三个接口：特征、前向规划器、演示数据来源

把示例从网格世界搬到连续控制或导航任务，不需要改动最大熵IRL 的主循环，只需要替换三个接口。特征方面：查表特征换成 RBF 特征或小网络编码器，保持线性权重结构，可解释性仍然保留；如果特征换完梯度崩了，先回到语义标签（比如目标方向、障碍距离、速度）这类低维特征，不要一上来就上高维传感器原始值。前向规划器方面：软价值迭代换成你的仿真器或 MPC 的软化版本，核心要求是能输出当前策略分布和期望特征，其他都无所谓。演示数据来源方面：把脚本里伪造的R_true换成人或历史日志的轨迹，同时注意状态编号和动作空间维度的映射一致性。

前阵子我把这套示例往一个二维导航任务上迁移，最直接的血泪经验是：直接用原始传感器特征跑了整整两天，权重一直在震荡，后来换成“目标相对方向 + 最近障碍距离 + 速度大小”三个语义特征，主循环一行没改，训练就开始正常收敛。那之后我养成了一个习惯：训练前一定打印demo_exp分布，确认每条特征在演示里有足够覆盖，再做任何调参。这套示例的边界和长处在网格世界看得很清楚：线性奖励 + 解析期望特征，能让你把 IRL 的主循环理解到极致，再往里塞网络、塞连续控制，你才有底气判断问题到底出在算法还是出在特征。希望这一篇能帮你少踩几个我用时间换来的坑。

本文还有配套的精品资源，点击获取