机器学习实战——最大熵模型，实战最大熵模型

机器学习实战——最大熵模型，实战最大熵模型

信息熵

香农与1948年提出的信息论以其中的信息熵这一基本概念为基础，用来表征人们对客观事件的不确定性的度量，与物理学中的熵表征物质的混乱程度有相似之处。
当处理不确定事件时，最常用的方式就是用概率方式描述，一般假定不确定的事件A每种可能的状态都有一个概率与之对应：P(Ai)s.t.∑i=1nP(Ai)=1P(Ai)≥0
而总共可能的状态如果是有限个（n个），那么就是离散的事件，用离散的随机变量表示；如果可能的状态是无穷多个，那么就是连续的事件，用连续的随机变量表示。本文所有说明仅以离散事件说明（连续事件只需将求和符合换成积分即可）。
人们对不确定的事件的信息知道得越多，信息熵就越小；知道得越少，信息熵就越大。信息熵表征了人们对不确定事件知道的信息的多少。从上述分析可以看出，不确定事件的每个状态都包含一定的信息，该信息量与概率值成反比，香农使用二进制比特来定义信息量的多少：H(Ai)=log(1P(Ai))=−logP(Ai)
log(⋅)是以2为底的对数。针对一个状态所含的信息量使用上述定义式得到，如果该状态的概率越大，则信息量越小。概率为1的状态是确定的就不携带任何信息，为0；反之概率越小的，信息量越大，但是，当概率为0时理论上携带了无穷多的信息（因为我们对其不知道任何信息）。
信息熵是表征不确定事件所有状态携带信息量的期望值：
H(A)=∑i=1npilog(1pi)=E(log(1pi))
式中以pi代替P(Ai)，因此信息熵表征了一个不确定事件A的所有可能状态所提供的平均信息量，信息熵越大，表明携带的平均信息量越大，不确定性也就越大；反之携带的平均信息量越少，不确定性越小。

最大熵原理

从上述定义可以看出当考虑所有状态时，信息熵是一个n维函数，是通过在每个维度(也就是每个状态)的信息量求加权和得到，因此求解上述定义式的最大值就是求解如下的约束最优化问题：
maxH(A)=−∑i=1npilog(pi)s.t.∑i=1npi=1pi≥0
使用Language乘数法求解得到：
L(p1,p2...pn,λ)=−∑i=1npilog(pi)+λ(∑i=1npi−1)∂L∂pi=−log(pi)−ln2+λ=0∂L∂λ=∑i=1npi−1=0
求解上述方程得到信息熵取得最大值时的条件：
p1=p2=p3....pn=1n,λ=ln2−log(n)
因此，对于未知的每个状态赋予同样的概率，使用均匀分布得到的信息熵是最大的。当人们对不确定事件进行模拟时，在一直部分知识的前提下，对未知的部分进行最合理的推断应该是最符合自然状态（最不确定状态）的推断，也就是熵最大的推断。这也与物理学中的熵增原理符合：物质总是趋于熵增大的方向变化，也就是最自然的状态变化。这种采取的保留全部的不确定性，使用最符合自然状态的推断，从满足约束条件的模型集合中选取熵最大的模型的方法就是最大熵原理。
从上述优化的角度可以看出，对于不确定的状态使用等可能的概率来寻找模型不易操作，而使用信息熵这一个数值指标去度量，通过寻找其最大值来表示等可能性在数学论证和编程实现上获得了可能。

最大熵模型

模型学习的是条件概率P(Y|X)，属于判别式模型。
对于给定的训练数据集，首先可以得出P(X,Y)和P(X)的经验分布P^(X,Y)和P^(X)：
P^(X,Y)=Count(X=x,Y=y)NP^(X)=Count(X=x)N
式中，N为训练样本的个数，Count(⋅)表示满足条件的样本个数。
其次从样本中选取或者设计n个特征函数：
fi(x,y)={1,0,if x,y 满足预定事实otherwise，i=1,2,...n
这些特征在建模过程中对输出有影响，所建立的模型应该满足所有这些特征的约束，所有符合这些约束的模型产生一个集合C。
特征函数关于经验分布P^(X,Y)的期望：
EP^(f)=∑x,yP^(X,Y)f(x,y)
特征函数关于理论模型P(X,Y)的期望：
EP(f)=∑x,yP(X,Y)f(x,y)=∑x,yP(X)P(Y|X)f(x,y)
理论模型应该与训练样本中的分布一致，因此P(X)=P^(X)，则
EP(f)=∑x,yP^(X)P(Y|X)f(x,y)
如果模型P(Y|X)能够从训练数据中获取携带的信息，那么就需要满足上述两个期望相等，即有：
EP^(f)=EP(f)
也就是：
∑x,yP^(X,Y)f(x,y)=∑x,yP^(X)P(Y|X)f(x,y)
上述就是模型P(Y|X)学习时必须满足的约束条件，一共选取的是n个特征函数，则对应有n个约束条件。
满足所有约束条件的模型集合：
C={P|EP^(fi)=EP(fi),i=1,2...n}
从中选取条件熵最大的模型p^(y|x)作为判别x属于y的概率大小的判定函数，即有
p^(y|x)=argmaxP∈CH(P)
定义在模型P(Y|X)上的条件熵为：
H(P)=H(Y|X)=∑xP(x)H(Y|X=x)=−∑xP(x)∑yP(y|x)log(P(y|x))=−∑x,yP(x)P(y|x)log(P(y|x))
其中建立的模型的分布P(x)应该与训练数据的经验分布P^(x)相等，故有
p^(y|x)=argmaxP∈CH(P)=−∑x,yP^(x)P(y|x)log(P(y|x))
从上述满足特征函数定义的约束条件下的模型集合中选取条件熵最大的模型p^(y|x)就得到了最大熵模型，就可以用它来进行分类或判别。

模型学习

通过上述分析可以得出，最大熵模型的学习就是求解如下的约束最优化问题：
argminP∈C−H(P)=∑x,yP^(x)P(y|x)log(P(y|x))s.t.EP^(fi)=EP(fi),i=1,2...n∑yP(y|x)=1
使用Language乘数法求解，n个约束共有n个乘子w1,w2...wn，同时外加归一化约束条件的乘子w0，以p表示模型变量P(y|x)则
L(p,w)=∑x,yP^(x)plog(p)+∑i=1nwi(∑x,yP^(x,y)fi(x,y)−∑x,yP^(x)pf(x,y))+w0(∑yp−1)
由于上述函数为变量p的凸函数，因此原始问题与其对偶问题的解是等价的，原始问题为
argminP∈CmaxwL(p,w)
对偶问题为
argmaxwminP∈CL(p,w)
首先求解内层最小值，以p为自变量，求导得：
∂L∂p=∑x,yP^(x)(log(p)+1+w0−∑i=1nwifi(x,y))=0p=exp(∑i=1nwifi(x,y)−1−w0)=exp(∑i=1nwifi(x,y))/exp(1+w0)
由归一化条件得
∑yp=∑yexp(∑i=1nwifi(x,y))/exp(1+w0)=1exp(1+w0)=∑yexp(∑i=1nwifi(x,y))
故有
p^(y|x)=exp(∑ni=1wifi(x,y))∑yexp(∑ni=1wifi(x,y))
这便是最大熵原理选择出来的最大熵模型，从这个结果可以看出，最大熵模型与Lgistic回归模型类似，都是在对数空间下的线性模型。模型的学习就是要确定各个特征的权重wi，这一步利用上述对偶问题的外部最优化问题求解，就是求解已知利用权重和特征表示的最大熵模型p^条件下关于权重w的最大值
argmaxwΨ(w)=L(p^(y|x),w)=∑x,yP^(x,y)∑i=1nwifi(x,y)−∑x,yP^(x)P(y|x)log(Zw(x))=∑x,yP^(x,y)∑i=1nwifi(x,y)−∑xP^(x)log(Zw(x))
上述最后一步利用P(y|x)对y求和的归一化条件得出。然后上述问题就变为无约束优化问题，可以利用GIS、IIS、拟牛顿法或者其他最优化方法求解即可得到最优解w∗。那么将学习到的最优解w∗带入最大熵模型p^中记得到了最终的分类决策函数。

决策函数的解释

根据学习模型为判别式模型，使用学习的策略为使用对数损失函数：
L(Y,f(X))=−log(P(Y|X))
以优化经验风险最小化为策略进行：
argminfRemp(f)=∑i=1NP^(x,y)L(Y,f(X))=argmaxf∑i=1NP^(xi,yi)(∑lnwlfl(xi,yi)−log(Zw(xi)))=argmaxf∑x,yP^(x,y)∑lnwlfl(x,y)−∑xP^(x)log(Zw(x))
最后一步利用的是联合分布求和得到边缘分布的方式，即∑yP^(x,y)=P^(x)。因此可以看出上式就是前面优化问题的对偶形式Ψ(w)的形式相同，也就得出了目标函数的最优化问题与经验风险最小化策略等价。

具体实现

特征函数的设定

当使用最大熵模型进行训练时，核心问题就是如何选定或者设计特征函数，由于是判别式模型，选取的特征函数在结果中并没有什么可解释性，因此这个过程是需要针对具体问题具体分析，同时可能需要结合具体领域的背景和先验知识来确定，因此此处本人不做太多讨论，这里为了实验的顺利进行，我使用的垃圾邮件分类的数据，数据一共有57个维度，分别从特征词的频率、大写字母的平均长度、总数等方面统计，因此我使用了如下的方式定义了八个特征函数：分别是特征词频率之和、大写字母平均长度、总数和最长片段数目进行设计。实现如下

FEATURE_FUNCTION = 8
def featureFunc(x, y):
    feats = zeros(FEATURE_FUNCTION)
    if (sum(x[0,:54]) > 5.4 and y == 0) or (sum(x[0,:54]) < 5.4 and y == 1):
        feats[0] = 1
    else:
        feats[1] = 1
    if (x[0,54] > 3 and y == 0) or (x[0,54] < 3 and y == 1):
        feats[2] = 1
    else:
        feats[3] = 1
    if (x[0,55] > 100 and y == 0) or (x[0,55] < 50 and y == 1):
        feats[4] = 1
    else:
        feats[5] = 1
    if (x[0,56] > 500 and y == 0) or (x[0,56] < 300 and y == 1):
        feats[6] = 1
    else:
        feats[7] = 1
    return feats

经验分布的求解

从前面的模型分析可以得出，针对训练数据需要计算经验分布P^(x,y)和P^(x)，这两个分布分布定义了如下的两个函数

def EmpProbX(ds, x):
    '''Build the empirical distribution of X'''
    N   = ds.shape[0]
    Nx  = mat(ones(N)).transpose() * x
    res = sum(ds - Nx, 1)
    return len(res[res == 0]) * 1.0 / N

def EmpProbXY(ds, labels, x, y):
    '''Build the empirical distribution of X,Y'''
    data = mat(ones((ds.shape[0], ds.shape[1] + 1)))
    data[:, :-1] = ds
    data[:, -1]  = labels
    line = mat(zeros(x.size + 1))
    line[0, :-1] = x
    line[0, -1]  = y
    return EmpProbX(data, line)

P^(x,y)的计算过程就是，给定训练数据和类标记，以及某一个样本(x,y)，找出训练集合中与之相同的样本个数，然后除以训练集合总数即可，P^(x)的计算与之类似，但是只需要给定样本数据和计算样本x。详细过程是使用Numpy提供的矩阵的运算实现的，构建两个矩阵之后进行元素与元素求差，最后统计按行求和为0的个数。

拟牛顿法中的梯度求解

训练过程使用拟牛顿法求解，需要将上述求解目标函数Ψ(w)的最大值变为最优化方法的标准形式——求解最小值，只需要将相减的两项互换即可。因此互换之后的梯度
gwi=−∂Ψ(w)∂wi=∑x,yP^(x)∑ni=1wifi(x,y)Zw(x)fi(x,y)−∑x,yP^(x,y)fi(x,y)=∑x,yP^(x)P(y|x)fi(x,y)−∑x,yP^(x,y)fi(x,y)=Ep(fi)−EP^(fi)
按照上述方式构造的求解梯度的函数如下：

def ComputeGrad(ds, y, w):
    '''Compute the gradient'''
    g = mat(zeros(FEATURE_FUNCTION)).transpose()
    featsList = [] ##Compute feature function
    for m in range(ds.shape[0]):
        f = featureFunc(ds[m,:], y[m,0])
        featsList.append(f)

    for k in range(FEATURE_FUNCTION):
        ##Compute expectation of empirical prob of X,Y
        eep = 0.0
        for m in range(ds.shape[0]):
            p = EmpProbXY(ds, y, ds[m,:], y[m,0])
            eep += p * featsList[m][k]
        ##Compute expectation of theoritical prob of X,Y
        etp = 0.0; 
        for m in range(ds.shape[0]):
            p = EmpProbX(ds, ds[m,:])
            f = featsList[m][k]
            t = tackleExpModel(ds[m,:], w, sum(featsList[m] * w))
            etp += p * f * t
        g[k, 0] = etp - eep
    return g

其中计算exp(∑ni=1wifi(x,y))/∑yexp(∑ni=1wifi(x,y))时涉及到的指数次方计算过程，可能会出现数值溢出，使用了如下的方式进行特殊处理：

def tackleExpModel(x, w, xi):
    '''Use the exp function may encounter numeric overflow'''
    Z = [sum(featureFunc(x, y) * w) for y in Y]
    m = max(Z)
    return exp(xi - m - log(sum([exp(z-m) for z in Z])))

首先是按每个分量y求∑ni=1wifi(x,y)保存在列表Z中，然后用Z中的最大值提出来，对整个先去对数，然后再求指数的形式进行优化，最后求和。形式化如下（以zi表示Z中第i个元素）
ezi/∑i=1nezi=exp{log(ezi)−log(em∑i=1nezi−m)}=exp{zi−m−log∑i=1nezi−m}
其中m是Z中元素的最大值，这样所有元素减去最大值之后再进行指数计算就能避免溢出。

拟牛顿法求解

求解上述每个梯度分量之后，利用拟牛顿法就可以求解了，有关拟牛顿法的内容不再赘述，可以参考有关最优化方法的书籍。详细实现就是按照拟牛顿法的步骤，在限定的最大迭代次数内进行迭代，每次迭代使用上面的方法计算梯度g、矩阵H和权重w，只到梯度的范数达到给定的要求或者梯度开始增大或者达到最大迭代次数就跳出迭代。（具体代码在此不予给出）

预测

使用学习到的权重，在给定的特征函数下依次计算每个类别下的概率值，得出最大的概率值就判定为相应类别。

def MaxEntClassify(w, predict):
    res = [-1.0, 0]
    s   = 0.0
    for y in Y:
        tmp = exp(sum(featureFunc(predict, y) * w))
        s += tmp
        if tmp > res[0]:
            res[0] = tmp
            res[1] = y
    res[0] = res[0] / s
    return res

程序返回的是一个列表，第一项是判定的所在类别的概率，第二项是判定的类别的标号。
最后，使用训练数据进行训练后，进行了简单的测试，结果如下图所示：