【逆强化学习-2】最大熵学习（Maximum Entropy Learning）

0.引言

\qquad 本文是逆强化学习系列的第2篇，其余博客传送门如下：

逆强化学习0-Introduction
逆强化学习1-学徒学习

\qquad 最大熵学习是2008年出现的方法，原论文（链接见【逆强化学习0】的博客）使用的Reward的函数仍然是线性模型，但是优化的思想和之前谈到的学徒学习有本质差别，由于需要一些概率论和随机过程分析课程的知识，原paper的理论也十分晦涩难懂。本人凭借粗浅的理解给大家一个浅显易懂的解释。

原会议的presentation（PPT）永久免费
原paper见部分0- Introduction部分
GitHub源代码
本文点赞破百，解锁额外代码（DQN+maxEnt）

\qquad 学徒学习（APP）是最大化间隙策略（MMP）的一种扩展，通过求解满足最大化间隙的Reward来计算Reward，从而使得Lean的行为越来越趋向于Expert（但又不好于Expert），这种方法往往叫做Feature Matching。其缺点在于对于存在多种合理的Reward的函数或者Expert存在多种次优轨迹时，该方法就无能为力了。APP本质是有约束优化问题，而优化变量是feature的discount-expectation基向量的坐标$\theta$。然而对于每一个策略$\pi$而言，都可能存在多个Reward函数使其最优。当演示了次优行为时，需要多个策略混合来匹配特征计数，这就让Feature Matching这件事在Expert轨迹存在多个Feature期望值时变得非常模糊。在APP中，这是通过求平均的方式解决的，然而这明显不是一个合理的解决方案。
\qquad 最大熵学习（MaxEnt）同样是Feature Matching的方法，与APP不同的是，其采用了一种有原则的方式消除了这种匹配歧义。而这种原则就是最大熵原则，该原则基于一种假设——即专家系统轨迹生成自己的专家特征期望的策略是最优轨迹（即下文的约束条件1）.
\qquad 可以简单的理解为，在APP中作为损失函数的特征匹配，在MaxEnt中被放入了约束条件中，而MaxEnt正是在满足这个约束条件的情况下，要求以$\theta$为Reward函数参数时，轨迹概率分布$P(\zeta |\theta )$的信息熵最大。
\qquad 至于为什么要求信息熵最大，原paper中并无详细说明，只是指出这已经在reference里面有了相关研究，本人查阅相关资料，给出以下几个理由供大家参考：

1. 物理系统的稳定状态通常趋向于熵最大
2. 只有P恒为0的概率分布熵才为0，正态分布是所有概率分布中熵最大的（会议presentation里面说均匀分布的信息熵最大，确认过是个错误结论，试想一下均匀分布的分布区间有限而正态无限）
3. 熵越大，先验信息越少，最大熵估计也是统计决策理论中常用的一种估计原则

1.算法原理

下面就简单介绍一下这个熵，对于连续变量而言，信息熵通常表示为
$$Ent={\int }_{x\backsim \pi }-p(x)logp(x)$$
对于强化学习任务而言，最大化信息熵写为：
$$\begin{array}{rl}& \max \sum _{\zeta \in D}-P(\zeta |\theta )logP(\zeta |\theta )\\ s.t.& \{\begin{array}{l}{\sum }_{\zeta \in D}P(\zeta |\theta )f_{\zeta }=\tilde{f}\\ {\sum }_{\zeta \in D}P(\zeta |\theta )=1\end{array}\end{array}$$

构造拉格朗日函数
$$L(P,\lambda ,\mu )=\sum _{\zeta \in D}[P(\zeta |\theta )logP(\zeta |\theta )+\lambda (P(\zeta |\theta )f_{\zeta }-\tilde{f})+\mu (P(\zeta |\theta )-1)]$$

应用拉格朗日函数的KKT条件
$$\begin{array}{clc}\nabla L_P=& {\sum }_{\zeta \in D}logP(\zeta |\theta )+1+\lambda f_{\zeta }+\mu =0& ①\\ \nabla L_{\lambda }=& {\sum }_{\zeta \in D}P(\zeta |\theta )f_{\zeta }-\tilde{f}=0& ②\\ \nabla L_{\mu }=& {\sum }_{\zeta \in D}P(\zeta |\theta )-1=0& ③\end{array}$$
由①③式得
$$P(\zeta |\theta )=\frac{exp(-1-\mu -\lambda f_{\zeta })}{{\sum }_{\zeta \in D}exp(-1-\mu -\lambda f_{\zeta })}$$
\qquad 光靠这个式子肯定是解不出最优的$\theta$的，这就要提到原paper的另一个假设——使用$\theta$参数的Reward函数$R_{\theta }(\tau )$时，$\zeta$轨迹的概率$P(\zeta |\theta )$正比于$R_{\theta }(\zeta )$的自然指数，再加上概率归一性约束，可得专家系统策略的轨迹概率为：
$$P(\zeta |\theta )=\frac{exp(R_{\theta }(\zeta ))}{{\int }_{\tau \in D}\left[exp(R_{\theta }(\tau ))\mathrm{d}\tau \right]}$$
需要注意的是，这里的$R(\theta )$指的是累积奖赏而非单步奖赏。

上式中的$Z$在paper中又被称为partial function，原文是已知原系统的dynamic model的，即已知系统的状态转移概率。在不知道状态转移概率时$Z$无法直接求得，通常也有三种方法：

1. 拉普拉斯近似（Laplace Approximation）
2. 值函数近似（Value Function Approximation）
3. 采样近似（Sample-Based Approximation）

有读者肯定会疑问，原paper中给出的损失函数不是最大信息熵而是最大似然，这又是为什么。原paper中给出了一个让人难以理解的解释：

Maximizing the entropy of the distribution over paths subject to the feature constraints from observed data implies that we maximize the likelihood of the observed data under the maximum entropy (exponential family) distribution derived above (Jaynes 1957).

——即从观测数据上满足feature matching的约束（约束1）的条件下最大化轨迹分布的信息熵等价于在最大信息熵分布的条件下从观测数据最大化似然。本文不对此深究，感兴趣的朋友可以研究一下下面这篇论文

Jaynes, E. T. 1957. Information theory and statistical mechanics. Physical Review 106:620–630.

\qquad 而假设是专家系统是最大熵分布的，因此对专家轨迹概率使用最大似然，得到
$$L(\theta )=\sum _{\zeta \in E}logp(\zeta |\theta )$$
即轨迹概率的最大似然。代入最大熵分布下的轨迹概率公式（其中E代表Expert的轨迹空间，而D代表Agent的轨迹空间（可以认为是全部轨迹空间），E的采样空间是与损失函数直接挂钩的，而D的采样空间则用来对$Z$估计的）：
$$\begin{array}{rl}L(\theta )& =\sum _{\tau \in E}logp(\tau |\theta )\\ & =\sum _{\tau \in E}log\frac{1}{Z}exp(R_{\theta }(\tau ))\\ & =\sum _{\tau \in E}{R}_{\theta }(\tau )-MlogZ\\ & =\sum _{\tau \in E}{R}_{\theta }(\tau )-Mlog\sum _{\tau \in D}exp(R_{\theta }(\tau ))\\ {\nabla }_{\theta }L& =\sum _{\tau \in E}\frac{d{R}_{\theta }(\tau )}{d\theta }-M\frac{1}{{\sum }_{\tau \in D}exp(R_{\theta }(\tau ))}\sum _{\tau \in D}\left[exp(R_{\theta }(\tau ))\frac{d{R}_{\theta }(\tau )}{d\theta }\right]\\ & =\sum _{\tau \in E}\frac{d{R}_{\theta }(\tau )}{d\theta }-M\sum _{\tau \in D}\left[\frac{exp(R_{\theta }(\tau ))}{{\sum }_{\tau \in D}exp(R_{\theta }(\tau ))}\frac{d{R}_{\theta }(\tau )}{d\theta }\right]\\ & =\sum _{\tau \in E}\frac{d{R}_{\theta }(\tau )}{d\theta }-M\sum _{\tau \in D}\left[p(\tau |\theta )\frac{d{R}_{\theta }(\tau )}{d\theta }\right]\\ & =\sum _{\tau \in E}\frac{d{R}_{\theta }(\tau )}{d\theta }-M\sum _{s_i\in S}\left[p(s|\theta )\frac{d{r}_{\theta }(s)}{d\theta }\right]\end{array}$$
归一化的损失函数为：
$${\nabla }_{\theta }\overline{L}=\frac{1}{M}\sum _{\tau \in E}\frac{d{R}_{\theta }(\tau )}{d\theta }-\sum _{s_i\in S}\left[p(s|\theta )\frac{d{r}_{\theta }(s)}{d\theta }\right]$$
其中M是专家轨迹的条数，如果状态空间是无限的，则不能直接套用此公式（但不代表无法解决）。
对于线性Reward，轨迹的累积Reward为
$$R_{\theta }(\zeta )={\theta }^T{f}_{\zeta }=\sum _{s_j\in \zeta }{\theta }^T{f}_{s_j}$$
Expert产生的Feature Expectation为
$$\tilde{f}=\frac{1}{m}\sum _i{f}_{{\tilde{\zeta }}_i}$$
损失函数梯度可表示为：
$${\nabla }_{\theta }L=\tilde{f}-\sum _{s_i\in S}{D}_{s_i}{f}_{si}$$
其中D为状态访问频次（State Visitiation Frequency），可以通过不断与环境互动近似出。
总结一下这个公式的推导需要注意一下几点：

1. 最大熵原则是建立在Feature Matching的基础上的，而轨迹概率分布的公式则是由最大熵原则+约束推导出的
2. 最大熵原则的轨迹概率分布公式，未知配分函数项Z是在全部轨迹集上求和，因此是使用Agent的轨迹进行近似
3. 最大化专家系统的轨迹概率似然其实是一个与原问题等价的优化问题，因此损失函数导数的第一项是在Expert的Demonstrations上求和（或求积分），而不是Agent的。

2.仿真

\qquad 本文的仿真平台参照了github上的资源，并进行了略微修改（仿真环境在学徒学习那篇有了详细介绍，在此就不再赘述了）：

GitHub源代码
本文点赞破百，解锁额外代码（DQN+maxEnt）

使用方法仍然是直接运行train.py即可（注意需要在mountaincar/maxent/的目录下运行）。和学徒学习的代码一样，也是基于Q-Table的。
\qquad 源代码中对Feature没有做任何的提取，直接将每个状态（20个位置采样×20个速度采样总共400个离散状态）作为Feature。假设不同特征之间是解耦的，Feature Matrix就是对角矩阵，即因此状态访问频次×特征即特征访问频次。
\qquad 在源代码中learner_feature_expectations即特征访问频次，而归一化之后即为梯度的第二项
源代码的其中一部分如下：

	expert = expert_feature_expectations(feature_matrix, demonstrations)
    learner_feature_expectations = np.zeros(n_states)

    theta = -(np.random.uniform(size=(n_states,)))

    episodes, scores = [], []

    for episode in range(30000):
        state = env.reset()
        score = 0

        if (episode != 0 and episode == 10000) or (episode > 10000 and episode % 5000 == 0):
            learner = learner_feature_expectations / episode
            maxent_irl(expert, learner, theta, theta_learning_rate)
                
        while True:
            state_idx = idx_state(env, state)
            action = np.argmax(q_table[state_idx])
            next_state, reward, done, _ = env.step(action)
            
            irl_reward = get_reward(feature_matrix, theta, n_states, state_idx)
            next_state_idx = idx_state(env, next_state)
            update_q_table(state_idx, action, irl_reward, next_state_idx)
            
            learner_feature_expectations += feature_matrix[int(state_idx)]

看完原github的程序，本人还有一个疑点，即是maxent.py文件中的一段

def maxent_irl(expert, learner, theta, learning_rate):
    gradient = expert - learner
    theta += learning_rate * gradient
    # Clip theta
    for i in range(len(theta)):
    	if theta[i]>0:
    		theta[i]=0

\qquad 原文中的clip theta实际上是防止theta超过0，类似于深度学习中的梯度截断操作，然而这个操作在我尝试多次之后并无用处，而且也没有任何意义（因为在theta>0后，也会在下一次迭代时通过学习使得theta重新<0）。本人的建议是增加一个学习率递减的schedule，并且将clip的范围从[-inf,0]修改到[-0.5,0.5]，可以获得相对稳定的学习率曲线，下面分别是clip(-0.5,0.5)和clip(-0.5,0)的对比:

clip(-0.5,0.5)	clip(-0.5,0)

\qquad 可以发现增加一部分梯度的正向范围反而更有利于学习，这是由于expert-learn的极小值是在0处取得的，然而在学习率固定时，改函数会在离散迭代时在0附近震荡，梯度若在0处截断会导致学习率锐减为0（可能正是原作者用意？）。

下面是test.py保存的几组gif的图

Reward=-158	Reward=-138	Reward=-146

希望本文对您有帮助，谢谢阅读！