最大熵模型

---

前言

最大熵模型是由最大熵原理推导实现。最大熵原理是概率模型学习的一个准则，该原理认为在已知约束条件下，且没有其他的先验知识来限制模型时，选择具有最大熵的模型，可以保证不会引入偏见。

---

一、最大熵模型是什么？

最大熵原理应用到分类得到最大熵模型。
给定训练集$T=(x_1,y_1),(x_2,y_2),...,(x_N,y_N)$，联合分布P(X,Y)以及边缘分布P(X)的经验分布都可以由训练数据得到：
$$\begin{gathered} \tilde{P}(X=x,Y=y)=\frac{count(X=x,Y=y)}{N} \\ \tilde{P}(X=x)=\frac{count(X=x)}{N} \end{gathered}$$
用特征函数f(x,y)描述输入x和输出y之间的某一个事实，特征函数是一个二值函数，当x与y满足某一事实时取1，否则取0。例如，可以令特征x与标签y在训练集出现过时取1，否则取0。
特征函数f(x,y)关于经验分布$\tilde{P}(X=x,Y=y)$的期望值为：
$$E_{\tilde{P}}(f)=\sum _{x,y}\tilde{P}(x,y)f(x,y)$$
特征函数f(x,y)关于模型P(Y|X)与经验分布$\tilde{P}(x)$的期望值为：
$$E_P(f)=\sum _{x,y}\tilde{P}(x)P(y|x)f(x,y)$$
如果模型能够获取训练数据中的信息，那么就可以假设这两个期望值相等，即：
$$\sum _{x,y}\tilde{P}(x,y)f(x,y)=\sum _{x,y}\tilde{P}(x)P(y|x)f(x,y)$$
将上式作为模型学习的约束条件，条件数量对应特征函数个数，设所有满足约束条件的模型集合为：
C = { P ∣ ∑ x , y P ~ ( x , y ) f i ( x , y ) = ∑ x , y P ~ ( x ) P ( y ∣ x ) f i ( x , y ) , i = 1 , 2 , . . . , n } C=\{P|\sum_{x,y}\widetilde{P}(x,y)f_i(x,y)=\sum_{x,y}\widetilde{P}(x)P(y|x)f_i(x,y),\quad i=1,2,...,n\} C={P∣x,y∑​P (x,y)fi​(x,y)=x,y∑​P (x)P(y∣x)fi​(x,y),i=1,2,...,n}
其中n为特征函数个数。
定义在条件概率分布P(Y|X)上的条件概率熵为：
$$H(P)=-\sum _{x,y}\tilde{P}(x)P(y|x)\ln P(y|x)$$
模型集合C中条件熵H§最大的模型称为最大熵模型。

二、极大似然估计

最大熵模型的学习过程就是求解最大熵模型的过程，等价于求解以下最优化问题：
$$\begin{gathered} \max H(P)=-\sum _{x,y}\tilde{P}(x)P(y|x)\ln P(y|x) \\ s.t.\sum _{x,y}\tilde{P}(x,y)f_i(x,y)=\sum _{x,y}\tilde{P}(x)P(y|x)f_i(x,y),i=1,2,...,n \\ \sum _{x,y}P(y|x)\tilde{P}(x)=1 \end{gathered}$$
按照最优化问题的习惯，求解与上述问题等价的$\min -H(p)$。

引入拉格朗日乘子$w_0,w_1,...,w_n$将上述带约束条件的最优化问题转化为无约束的最优化问题，定义拉格朗日函数：
$$\begin{gathered} L(P,W)=-H(P)+w_0(1-\sum _{x,y}P(y|x)\tilde{P}(x))+\sum _{i=1}^n{w}_i(\sum _{x,y}\tilde{P}(x,y)f_i(x,y)-\sum _{x,y}\tilde{P}(x)P(y|x)f_i(x,y)) \\ =\sum _{x,y}\tilde{P}(x)P(y|x)\ln P(y|x)+w_0(1-\sum _{x,y}P(y|x)\tilde{P}(x))+\sum _{i=1}^n{w}_i(\sum _{x,y}\tilde{P}(x,y)f_i(x,y)-\sum _{x,y}\tilde{P}(x)P(y|x)f_i(x,y)) \end{gathered}$$
这里最优化的原始问题是：
$$\underset{P\in C}{\min }\underset{w}{\max }L(P,w)$$
对偶问题是：
$$\underset{w}{\max }\underset{P\in C}{\min }L(P,w)$$
在上式内部的最小化部分固定w，此时L(P,w)是关于P的函数，用$\theta (P)$表示，对P(y|x)求偏导数得：
$$\begin{gathered} \frac{\partial \theta (P)}{\partial P(y|x)}=\tilde{P}(x)(1+\ln P(y|x))-\tilde{P}(x)w_0-\tilde{P}(x)\sum _{i=1}^n{w}_i{f}_i(x,y) \\ =\tilde{P}(x)(1+\ln P(y|x)-w_0-\sum _{i=1}^n{w}_i{f}_i(x,y)) \end{gathered}$$
令偏导数等于0，因为$\tilde{P}(x)>0$，所以$1+\ln P(y|x)-w_0-{\sum }_{i=1}^n{w}_i{f}_i(x,y)=0$，解得：
$$P(y|x)=\frac{e^{{\sum }_{i=1}^n{w}_i{f}_i(x,y)}}{e^{1-w_0}}$$
由于${\sum }_{y}P(y|x)=1$，得到P(y|x)关于w的表达式为：
$$\begin{gathered} P_w(y|x)=\frac{1}{Z_w(x)}{e}^{{\sum }_{i=1}^n{w}_i{f}_i(x,y)} \\ {Z}_w(x)=\sum _y{e}^{{\sum }_{i=1}^n{w}_i{f}_i(x,y)} \end{gathered}$$
由上式表示的模型$P_w=P_w(y|x)$就是最大熵模型，w是最大熵模型的参数向量，每一维度为对应特征函数的权重。
令$\mathbf{w}=[w_1,w_2,...,w_n{]}^T,\mathbf{f}(\mathbf{x},\mathbf{y})=[f_1(x,y),f_2(x,y),...,f_n(x,y)]$，则：
$$\sum _{i=1}^n{w}_i{f}_i(x,y)={\mathbf{w}}^T\mathbf{f}(\mathbf{x},\mathbf{y})$$
将$P_w(y|x)$用w表示代入$L(P,W)$的表达式并结合${\sum }_{x,y}P(y|x)\tilde{P}(x)=1$得：
$$\begin{gathered} L(P,W)=\sum _{x,y}\tilde{P}(x)P(y|x)({\mathbf{w}}^T\mathbf{f}(\mathbf{x},\mathbf{y})-\ln Z_w(x))+ \\ \sum _{x,y}\tilde{P}(x,y){\mathbf{w}}^T\mathbf{f}(\mathbf{x},\mathbf{y})-\sum _{x,y}\tilde{P}(x)P(y|x){\mathbf{w}}^T\mathbf{f}(\mathbf{x},\mathbf{y}) \\ =\sum _{x,y}\tilde{P}(x,y){\mathbf{w}}^T\mathbf{f}(\mathbf{x},\mathbf{y})-\sum _{x,y}\tilde{P}(x)P(y|x)\ln Z_w(x) \end{gathered}$$
因为：
$$\sum _y\tilde{P}(x)P(y|x)\ln Z_w(x)=\tilde{P}(x)\ln Z_w(x)$$
所以：
$$L(P,W)=\sum _{x,y}\tilde{P}(x,y){\mathbf{w}}^T\mathbf{f}(\mathbf{x},\mathbf{y})-\sum _x\tilde{P}(x)\ln Z_w(x)$$
所以最大熵模型可以转化为极大化上述函数。接下来考虑$P_w(x,y)$的对数似然函数，即
$$\begin{gathered} L_{\tilde{P}}(P_w)=\ln \prod _{x,y}P(y|x{)}^{\tilde{P}(x,y)}=\sum _{x,y}\tilde{P}(x,y)\ln P(y|x) \\ =\sum _{x,y}\tilde{P}(x,y){\mathbf{w}}^T\mathbf{f}(\mathbf{x},\mathbf{y})-\sum _{x,y}\tilde{P}(x,y)\ln Z_w(x) \end{gathered}$$
因为：
$$\sum _y\tilde{P}(x,y)\ln Z_w(x)=\ln Z_w(x)\sum _y\tilde{P}(x,y)=\ln Z_w(x)\tilde{P}(x)$$
所以：
$$L_{\tilde{P}}(P_w)=\sum _{x,y}\tilde{P}(x,y){\mathbf{w}}^T\mathbf{f}(\mathbf{x},\mathbf{y})-\sum _x\tilde{P}(x)\ln Z_w(x)$$
综上，极大化$L(P,w)$等价于最大熵模型的极大似然估计。

三、模型学习的最优化算法

IIS的想法是：假设最大熵模型当前的参数向量是$\mathbf{w}$，我们希望找到一个新的参数向量：$\mathbf{w}+\mathbf{\delta }$，使得模型的对数似然函数增大。不断重复这一过程直到找到对数似然函数的最大值。
更新前后的对数似然函数改变量为：
$$\begin{gathered} L_{\tilde{P}}(P_{w+\delta })-L_{\tilde{P}}(P_w)=(\sum _{x,y}\tilde{P}(x,y)(\mathbf{w}+\mathbf{\delta }{)}^T\mathbf{f}(\mathbf{x},\mathbf{y})-\sum _x\tilde{P}(x)\ln Z_{w+\delta }(x))- \\ (\sum _{x,y}\tilde{P}(x,y){\mathbf{w}}^T\mathbf{f}(\mathbf{x},\mathbf{y})-\sum _x\tilde{P}(x)\ln Z_w(x)) \\ =\sum _{x,y}\tilde{P}(x,y){\mathbf{\delta }}^T\mathbf{f}(\mathbf{x},\mathbf{y})-\sum _x\tilde{P}(x)\ln \frac{Z_{w+\delta }(x)}{Z_w(x)} \end{gathered}$$
利用不等式$-\ln \alpha \ge 1-\alpha$，得：
$$\begin{gathered} L_{\tilde{P}}(P_{w+\delta })-L_{\tilde{P}}(P_w)\ge \sum _{x,y}\tilde{P}(x,y){\mathbf{\delta }}^T\mathbf{f}(\mathbf{x},\mathbf{y})+1-\sum _x\tilde{P}(x)\frac{Z_{w+\delta }(x)}{Z_w(x)} \\ =\sum _{x,y}\tilde{P}(x,y){\mathbf{\delta }}^T\mathbf{f}(\mathbf{x},\mathbf{y})+1-\sum _x\tilde{P}(x)\sum _y{P}_w(y|x)e^{{\mathbf{\delta }}^T\mathbf{f}(\mathbf{x},\mathbf{y})} \end{gathered}$$
令$A(\mathbf{\delta }|\mathbf{w})={\sum }_{x,y}\tilde{P}(x,y){\mathbf{\delta }}^T\mathbf{f}(\mathbf{x},\mathbf{y})+1-{\sum }_x\tilde{P}(x){\sum }_y{P}_w(y|x)e^{{\mathbf{\delta }}^T\mathbf{f}(\mathbf{x},\mathbf{y})}$为对数似然函数的改变量的一个下界，因为$\mathbf{\delta }$是一个n维向量，而$A(\mathbf{\delta }|\mathbf{w})$在e的指数部分有$\mathbf{\delta }$，所以若此时求偏导会导致导数的每一个分量都有$\mathbf{\delta }$。因此，尝试寻找一个新的下界，使得指数部分只存在$\mathbf{\delta }$的单个分量。
令$f^{\#}(x,y)=1^T\mathbf{f}(\mathbf{x},\mathbf{y})$即(x,y)处特征出现的次数，则：
$$\begin{gathered} A(\mathbf{\delta }|\mathbf{w})=\sum _{x,y}\tilde{P}(x,y){\mathbf{\delta }}^T\mathbf{f}(\mathbf{x},\mathbf{y})+1-\sum _x\tilde{P}(x)\sum _y{P}_w(y|x)e^{f^{\#}(x,y)\frac{{\mathbf{\delta }}^T\mathbf{f}(\mathbf{x},\mathbf{y})}{f^{\#}(x,y)}} \\ =\sum _{x,y}\tilde{P}(x,y){\mathbf{\delta }}^T\mathbf{f}(\mathbf{x},\mathbf{y})+1-\sum _x\tilde{P}(x)\sum _y{P}_w(y|x)e^{f^{\#}(x,y){\sum }_{i=1}^n{\delta }_i\frac{f_i(x,y)}{f^{\#}(x,y)}} \end{gathered}$$
由于$\frac{f_i(x,y)}{f^{\#}(x,y)}\ge 0,{\sum }_{i=1}^n\frac{f_i(x,y)}{f^{\#}(x,y)}=1$且e的指数函数为下凸函数，则由琴生不等式得:
$$e^{f^{\#}(x,y){\sum }_{i=1}^n{\delta }_i\frac{f_i(x,y)}{f^{\#}(x,y)}}\le \sum _{i=1}^n\frac{f_i(x,y)}{f^{\#}(x,y)}{e}^{f^{\#}(x,y){\delta }_i}$$
所以：
$$A(\mathbf{\delta }|\mathbf{w})\ge su{m}_{x,y}\tilde{P}(x,y){\mathbf{\delta }}^T\mathbf{f}(\mathbf{x},\mathbf{y})+1-\sum _x\tilde{P}(x)\sum _y{P}_w(y|x)\sum _{i=1}^n\frac{f_i(x,y)}{f^{\#}(x,y)}{e}^{f^{\#}(x,y){\delta }^i}$$
令$B(\mathbf{\delta }|\mathbf{w})={\sum }_{x,y}\tilde{P}(x,y){\mathbf{\delta }}^T\mathbf{f}(\mathbf{x},\mathbf{y})+1-{\sum }_x\tilde{P}(x){\sum }_y{P}_w(y|x){\sum }_{i=1}^n\frac{f_i(x,y)}{f^{\#}(x,y)}{e}^{f^{\#}(x,y){\delta }^i}$为对数似然函数改变量的一个新下界，将其对$\mathbf{\delta }$的任一分量求导得：
$$\frac{\partial B(\mathbf{\delta }|\mathbf{w})}{\partial {\delta }_i}=\sum _{x,y}\tilde{P}(x,y)f_i(x,y)-\sum _{x,y}\tilde{P}(x)P_w(y|x)f_i(x,y)e^{f^{\#}(x,y){\delta }_i}$$
令偏导数等于0得到:
$$\sum _{x,y}\tilde{P}(x)P_w(y|x)f_i(x,y)e^{f^{\#}(x,y){\delta }_i}=E_{\tilde{P}}(f_i)$$
依次对$\mathbf{\delta }$的分量求上式的解即可解出$\mathbf{\delta }$。

为简化求解过程，一般假设$f^{\#}(x,y){\delta }_i=M$为一常数，即假设对任意的(x,y)，特征出现的次数相同。那么${\delta }_i$可显式得出：
$${\delta }_i=\frac{1}{M}\ln \frac{E_{\tilde{P}}(f_i)}{E_P(f_i)}$$
若$f^{\#}(x,y){\delta }_i$不为常数，则需要通过数值计算求解，如使用牛顿法迭代求解。

四、代码实现

代码如下：

import numpy as np
from collections import defaultdict
from sklearn.datasets import load_digits
from time import sleep
from tqdm import tqdm


def load_data():
    '''
    加载sklearn自带的手写数字识别数据集
    返回 输入、输出
    '''

    digits = load_digits()
    xs = digits.data
    ys = digits.target
    return xs, ys


class MaxEntropy(object):
    def __init__(self, train_xs, train_ys, test_xs, test_ys, iterations=200):
        self.train_xs = train_xs
        self.train_ys = train_ys
        self.test_xs = test_xs
        self.test_ys = test_ys
        self.iterations = iterations
        self.example_cnt, self.feature_cnt = train_xs.shape
        self.feature_label_pairs_dict, self.pair_cnt = self.get_feature_label()  # pair_cnt相当于特征函数个数
        self.weights = np.zeros((self.pair_cnt,))
        self.xy2id_dicts, self.id2xy_dict = self.get_search_dict()
        self.e_hat_f = self.get_e_hat_f()
        self.e_f = self.get_e_f()
        self.f_weight = self.get_f_weight()

    def get_feature_label(self):
        """
        获取
        每个特征的特征标签对字典
        总的特征标签对的个数
        :return:
        """
        feature_label_pairs = [defaultdict(int) for feature_id in range(self.feature_cnt)]
        pair_cnt = 0

        for example_id in range(self.example_cnt):
            for feature_id in range(self.feature_cnt):

                feature_label_pairs[feature_id][(self.train_xs[example_id][feature_id], self.train_ys[example_id])] += 1
        for pair_dict in feature_label_pairs:
            pair_cnt += len(pair_dict)

        return feature_label_pairs, pair_cnt

    def get_search_dict(self):
        """
        保存特征标签对xy对应id的字典数组
        保存特征标签对id对应xy的的字典
        :return:
        """
        idx = 0
        xy2id_dicts = [dict() for feature_id in range(self.feature_cnt)]
        id2xy_dict = dict()
        for feature_id in range(self.feature_cnt):
            for (x, y) in self.feature_label_pairs_dict[feature_id].keys():
                xy2id_dicts[feature_id][(x, y)] = idx
                id2xy_dict[idx] = (x, y)
                idx += 1

        return xy2id_dicts, id2xy_dict

    def get_e_hat_f(self):
        """
        特征函数f(x,y)关于经验分布P_hat(X,Y)的期望值
        :return:
        """
        e_hat_f = np.zeros((self.pair_cnt,))

        for feature_id in range(self.feature_cnt):
            for (x, y) in self.feature_label_pairs_dict[feature_id]:
                idx = self.xy2id_dicts[feature_id][(x, y)]
                e_hat_f[idx] = self.feature_label_pairs_dict[feature_id][(x, y)] / self.example_cnt

        return e_hat_f

    def get_pw_xy(self, x):
        """
        求解条件概率分布P_w(y|x)
        :param x:
        :return:
        """
        uniques = np.unique(self.train_ys)
        numerators = np.zeros((max(uniques)+1, ))

        for feature_id in range(self.feature_cnt):
            for label in uniques:
                if (x[feature_id], label) in self.feature_label_pairs_dict[feature_id]:
                    idx = self.xy2id_dicts[feature_id][(x[feature_id], label)]
                    numerators[label] += self.weights[idx]

        numerators = np.exp(numerators)
        z_w = np.sum(numerators)

        return numerators / z_w

    def predict(self, x):
        '''
        预测标签
        :param X:要预测的样本
        :return: 预测值
        '''
        return np.argmax(self.get_pw_xy(x))

    def test(self):
        '''
        对测试集进行测试
        :return:
        '''
        # 错误值计数
        error_cnt = 0
        # 对测试集中所有样本进行遍历
        for i in range(len(self.test_xs)):
            # 预测该样本对应的标签
            result = self.predict(self.test_xs[i])
            # 如果错误，计数值加1
            if result != self.test_ys[i]:
                error_cnt += 1
        # 返回准确率
        return 1 - error_cnt / len(self.test_xs)

    def get_e_f(self):
        """
        计算特征函数f(x,y)关于模型P(X｜Y）与经验分布P_hat(X)的期望
        :return:
        """
        e_f = np.zeros((self.pair_cnt,))

        for example_id in range(self.example_cnt):
            pw_xy = self.get_pw_xy(self.train_xs[example_id])
            for feature_id in range(self.feature_cnt):
                if (self.train_xs[example_id][feature_id], self.train_ys[example_id]) \
                        in self.feature_label_pairs_dict[feature_id]:
                    idx = self.xy2id_dicts[feature_id][self.train_xs[example_id][feature_id], self.train_ys[example_id]]
                    e_f[idx] += (1 / self.example_cnt) * pw_xy[self.train_ys[example_id]]

        return e_f

    def get_f_weight(self):
        """
        计算iis引进的一个量1 / f_#(x,y)
        :return:
        """
        f_weight = np.zeros((self.pair_cnt,))

        for feature_id_1 in range(self.feature_cnt):
            for x_y_pair_1 in self.feature_label_pairs_dict[feature_id_1]:
                idx = self.xy2id_dicts[feature_id_1][x_y_pair_1]
                similar_x_y_cnt = 0
                x_y_count = self.feature_label_pairs_dict[feature_id_1][x_y_pair_1]
                for feature_id_2 in range(self.feature_cnt):
                    for x_y_pair_2 in self.feature_label_pairs_dict[feature_id_2]:
                        if x_y_pair_2 == x_y_pair_1:
                            similar_x_y_cnt += self.feature_label_pairs_dict[feature_id_2][x_y_pair_2]
                f_weight[idx] = x_y_count / similar_x_y_cnt

        return f_weight

    def iis_train(self):
        """
        采用改进尺度迭代法进行训练
        :return:
        """
        for i in tqdm(range(self.iterations)):
            # sigma_array = self.f_weight * np.log(self.e_hat_f / self.e_f)
            sigma_array = (1 / self.feature_cnt) * np.log(self.e_hat_f / self.e_f)
            self.weights = self.weights + sigma_array
            if (i+1) % 5 == 0:
                accuracy = self.test()
                print('the accuracy is:%.4f' % accuracy)
            sleep(0.03)


if __name__ == '__main__':
    features, targets = load_data()

    train_count = int(len(features) * 0.8)

    train_xs, train_ys = features[:train_count], targets[:train_count]
    test_xs, test_ys = features[train_count:], targets[train_count:]

    # 初始化最大熵类
    maxEnt = MaxEntropy(train_xs, train_ys, test_xs, test_ys)
    # 开始训练
    print('start to train')
    maxEnt.iis_train()
    # 开始测试
    print('start to test')
    accuracy = maxEnt.test()
    print('the accuracy is:%.4f' % accuracy)

---