Press "Enter" to skip to content

NLP重剑篇之逻辑回归与条件随机场

本站内容均来自兴趣收集,如不慎侵害的您的相关权益,请留言告知,我们将尽快删除.谢谢.

对数线性模型,是我心中认为最为经典的模型,其形式简单,计算方便,却又鲁棒性强,适用范围广。本章将会介绍最广为人知的两种对数线性模型:逻辑回归与线性链条件随机场。我们会从多个角度来理解他们,并用纯python实现两个算法。本章中,我们会从多个角度来得到模型的形式,某些方式在数学上可能不严谨,但能让我们将各种模型更好地联系在一起,也能启发我们寻找新的模型。 相关代码实现开源在: https://github.com/wellinxu/nlp_store ,更多内容关注知乎专栏(或微信公众号): NLP杂货铺。

 

 

逻辑回归的定义

 

从朴素贝叶斯到逻辑回归

 

从最大熵到逻辑回归

 

从广义线性模型到逻辑回归

 

逻辑回归学习算法

 

梯度下降法

 

 

条件随机场的维特比算法

 

前向计算

 

后向计算

 

概率计算

 

期望计算

 

从隐马尔可夫模型到条件随机场

 

从最大熵模型到条件随机场

 

条件随机场的概率计算

 

条件随机场的预测算法

 

条件随机场的学习算法

 

改进的迭代尺度法

 

拟牛顿法BFGS

 

参考

 

逻辑回归

 

在很多问题上,逻辑回归的准确率都很有竞争性,甚至不比很多复杂的神经网络方法低,而且形式简单,可以并行,容易人工干预且模型结果易于解释。

 

本小节,先给出逻辑回归的定义,然后从三个方面来解释,最后给出逻辑回归的训练算法。

 

逻辑回归的定义

 

二项逻辑回归模型是如下条件概率分布:

 

其中x是输入(),Y={0,1}是输出,w,b是参数。

 

从朴素贝叶斯到逻辑回归

 

逻辑回归与朴素贝叶斯看似两种不相干的模型,却有着很强的内在联系性,根据之前的篇章 ( NLP重剑篇之朴素贝叶斯与隐马尔可夫 ) ,我们可以知道,朴素贝叶斯主要求的是:

 

如果我们先对数化在指数化

 

贝叶斯中我们可以直接通过样本-特征的频率,直接估计出各个概率,而在逻辑回归中,我们需要迭代计算出每个特征的权重,所以可以得到:

 

其中则是的参数形式。

 

因为其结果表示一个概率,所以我们添加归一化:

 

当类别只有两个的时候,我们可以将分子分母同除:

 

这样我们就得到了类似逻辑回归定义的形式,只是特征通过得到,属于{0,1}特征。如此看来,逻辑回归也在计算各个类的概率,朴素贝叶斯通过频率估计概率,而逻辑回归则试图寻找一个更合适的概率,而且朴素贝叶斯的特征主要是离散型,而逻辑回归可以是离散的也可以是连续的,这样看来,逻辑回归要比朴素贝叶斯更通用一点,且似乎有着不弱于朴素贝叶斯的准确度。

 

从最大熵到逻辑回归

 

熵本来是物理化学中的一个概念,原用来描述一个系统的混乱程度,混乱程度越高熵值越高,而且系统总是往混乱变高的方向发展,一般称为熵增过程。我们把这一概念用到模型中来,希望在已有的条件下,让结果的熵最大的模型,就是我们所需要的。

 

要找到最大熵,首先先定义下熵的计算方式,我们对x的概率分布p的熵计算方式为:

 

最大熵模型里的特征比较独特,我们用n个不同的特征函数来获取特征,第i个特征表示为,特征函数是自定义的任意关于x和y的函数,一般情况下,需要有意义,在实际情况中可能会将特征函数的值域定义为{0,1},这是为了更方便计算与简化问题。

 

根据数据集,我们可以计算出与两个经验分布,基于此,我们可以计算每个特征关于的期望值:

 

我们还可以计算特征关于模型P(Y|X)与的期望值:

 

如果模型可以获取足够的信息,则。

 

根据熵的计算方式,以及上诉条件,我们可以将模型约束为:

 

得到的形式上,跟支持向量机有点像,一个最大化,几个约束条件,而且现在的约束条件都是等式条件,类似的解法,使用拉格朗日乘子法,将约束最优化问题转为无约束最优化问题,再通过对偶问题求得最终解。

 

首先我们将原始问题转化为等价形式:

 

然后引进拉格朗日乘子,得到拉格朗日函数:

 

此时原始问题转化为:

 

其对偶问题为:

 

求L对P的偏导得【1】:

 

令偏导数等于0得:

 

根据概率和为1得:

 

如果我们将问题定义为二分类,将特征函数定义为:

 

那幺最大熵模型变成:

 

这样我们又得到了跟逻辑回归非常类似的形式,当然它们之间也有区别。从公式中我们可以看出,逻辑回归比最大熵模型多了个参数b,最大熵模型在所有类别下截距b的值一样大,所以可以简化去除参数b了,而根据上一小节,参数b是与当前类别的先验概率相关,而最大熵模型的原则则是:没有信息的情况下,认为各个值的概率相等,所以最大熵模型的参数b都是一样的。但是,如果我们再添加一个特征函数:

 

此时最大熵可以表示为:

 

此时跟逻辑回归就更类似了。在二分类的时候,且特征函数定义为上述形式,最大熵模型跟逻辑回归最为相似;而在多分类的时候,最大熵模型则表现出跟逻辑回归的明显差别:

 

逻辑回归的特征只有一套,但对不同的类别,其会跟不同的参数做内积运算,也就是说,不变,而随着类别不同会改变;

 

而最大熵模型则是参数只有一套,在不同类别时,其会跟不同的特征做内积运算,也就是说,不变,随着类别不同会有改变。

 

因为概率和为1的限制,所以尽管有n个类别,也只要做n-1次计算,所以当n=2为二分类的时候,逻辑回归与最大熵模型的参数与特征都只有一套,所以此时最为相似。

 

从广义线性模型到逻辑回归

 

在介绍广义线性模型之前,我们先介绍下指数分布,如果一个分布可以写成:

 

那幺这个分布就称之为指数分布。其中是标准参数,T(y)是充分统计量,多数情况T(y)=y。

 

广义线性模型有三个假设:

 

1、y的分布符合参数为的指数分布。

 

2、模型结果是给定x的T(y)的期望:E(T(y)|x)。

 

3、

 

如果我们处理二分类问题,假设y符合伯努利分布:

 

可以看出,伯努利分布也是一种指数族分布,且:

 

可以得到:

 

根据第二个假设,模型的结果为:

 

根据第三个假设得到:

 

最终得到:

 

这样我们就得到了逻辑回归模型,此时我们也更好理解,为什幺逻辑归回的值可以认为是该类别的概率。这里我们看出,逻辑回归模型只是广义线性模型的一个特例,实际上线性回归模型也是广义线性模型的特例,也许这就是逻辑回归模型为什幺叫回归模型的原因。这里我们证明了伯努利分布属于指数族分布,实际上正态分布、泊松分布、多项式分布等也都属于指数族分布,这幺看来,广义线性模型是一个很美妙的模型。

 

逻辑回归学习算法

 

在具体使用某种算法学习之前,我们先来确定下逻辑回归学习的目标函数或者损失函数。

 

我们可以根据似然函数来获取目标函数:

 

取对数得:

 

这样我们的目标就是,最大化。

 

另一方面,根据交叉熵,我们可以得到:

 

此时我们的目标是最小化loss函数,根据上面可以看出,似然函数与交叉熵的目标是等价的,所以下面我们就用梯度下降法(求目标函数的最小值)来最小化交叉熵损失函数。

 

梯度下降法

 

输入:样本,目标函数f(w|x,y),梯度函数,学习率,终止条件

 

输出:f(w|x,y)的极小点

 

1、取初始值,设k=0

 

2、当终止条件达到时,停止迭代,令;否则计算梯度,令,循环此步骤

 

根据梯度下降算法,我们先计算loss函数的梯度:

 

这样整个算法的代码为:

 

import random
import numpy as np

class LogisticRegress(object):
    def __init__(self, xs, ys, lr= 0.01, cycle=500):
        self.xs = xs
        self.ys = ys
        self.lr = lr    # 学习率
        self.cyclt = cycle    # 最大迭代次数
        self.params = [random.random() for i in xs[0]]
        # self.build2()
        self.build()
    def build2(self):
        # 梯度下降
        for i in range(self.cyclt):
            deltas = []
            for x, y in zip(self.xs, self.ys):
                y_hat = self.prob(x)
                delta = (y - y_hat)*x
                deltas.append(delta)
            delta = np.mean(deltas, axis=0)
            self.params = self.params - self.lr * delta
    def build(self):
        # 随机梯度下降
        for i in range(self.cyclt):
            for x, y in zip(self.xs, self.ys):
                y_hat = self.prob(x)
                delta = (y - y_hat)*x
                self.params = self.params - self.lr * delta
    def prob(self, x):
        x_w = -np.dot(x, self.params)
        p = 1-1/(1+np.exp(x_w))
        return p
    def predict(self, x):
        x_v = np.dot(x, self.params)
        if x_v < 0:
            return 0
        return 1

 

条件随机场

 

本文主要涉及的是线性链条件随机场,同样的,我们直接给出其定义与参数化形式:

 

定义:设均为线性链表示的随机变量序列,若在给定随机变量序列X的条件下,随机变量序列Y的条件概率分布P(Y|X)满足马尔科夫性:

 

则称P(Y|X)为线性链条件随机场。其参数化形式可以表示为:

 

其中是特征函数(含义参考最大熵模型),是权重。

 

从隐马尔可夫模型到条件随机场

 

根据之前章节,我们可以得到隐马尔可夫模型的计算目标:

 

上式可以转化为:

 

不妨假设存在,使得,则上式可转化为:

 

假设的所有取值为,的所有取值为,那幺上式可以转化为:

 

在隐马尔可夫中,上式中的概率值是通过样本频数直接计算出来的,这里我们将其视为位置参数,将其参数化,则得到:

 

其中z(x)是为了让概率和为1的归一化因子。

 

由此我们得到了跟条件随机场基本一致的形式,类似于逻辑回归与朴素贝叶斯的关系,条件随机场就像是参数化转移概率、观测概率、初始概率后的隐马尔可夫模型;当然除此之外,他们之间最大的差别就是可用的特征范围不同,由上式可以看出,隐马尔科夫模型只能用到局限于当下索引的特征(),而条件随机场则没有这个限制,可以用到不同位置上的特征,如。

 

从最大熵模型到条件随机场

 

上述小节中,在推导最大熵模型的时候,并没有假定y的类型,就是说,如果y是序列结果,我们依然可以得到最大熵模型的形式:

 

现在我们讨论x与y都是序列的情况,所以特征函数可以定义为:

 

此时模型可以写为:

 

这样我们得到了跟条件随机场一样的模型,由此可见,条件随机场可以看做是最大熵模型的一种特殊情况(预测值是序列,且特征函数符合某种类型)。相比较而言,最大熵模型可以解决的问题比较广,特征函数也比较泛,而条件随机场则更具针对性,且特征函数具有一定的指导意义。某种程度上说,最大熵模型给条件随机场提供了理论支撑,而条件随机场则为最大熵模型带来了更具体的应用。从这个角度看,我们只要构造一些列的特征函数,那幺不管y是怎样的结构(点,线,树,图),我们都能得到一个较好的模型。

 

本节并没有从经典的马尔科夫随机场、概率无向图、势函数等来解释条件随机场;粗浅地认为,那样解释条件随机场需要引用大量的定义、定理等,如果讲清楚不仅需要大量篇幅,而且需要非常深的理论知识,如果一略而过,又会有类似于【这个所谓的势函数又是哪里蹦出来】的种种疑问;所以本节先直接给出定义,然后从最大熵模型的角度来得到条件随机场,这样只要引入熵的定义即可。

 

根据定义我们可以根据crf相关输入获取其特征情况:

 

import numpy as np
import math

def exp_dot(w, f):
    """
    :param w: array-like, 1*n
    :param f: array-like, 1*n
    :return: float
    """
    return math.exp(np.dot(w, f))

class CrfInput(object):
    def __init__(self, x, y, label_num, features_fun, predict=False):
        self.x = x
        self.y = y
        self.start_y = -1    # 初始状态
        self.labels = range(label_num)
        self.label_num = label_num
        self.features_fun = features_fun  # 特征函数集合
        self.predict = predict    # 是否是预测
        self.all_features_by_x = []    # [{j:{i:feats}}]
        self.features_by_xy = []    # [[]]
        self.F = []    # 各特征值在序列上的求和, [], 1*len(features)
        self.f_sharp = 1    # F#的值, 所有特征的和
        self.get_x_features()  # 获取x所有情况下的特征
        if not self.predict:
            self.get_xy_features()  # 获取xy的特征
            self.F = np.sum(self.features_by_xy, 0)
            self.f_sharp = np.sum(self.F)
    def get_x_features(self):
        """基于x获取所有可能的特征"""
        for i, xi in enumerate(self.x):
            tem_fs = {}    # {j:{i:feats}}
            yis = self.labels
            if i == 0:
                yis = [self.start_y]
            for yi in yis:
                tem_fs_yi = {}    # {i:feats}
                for yj in self.labels:
                    tem_fs_yi[yj] = [f(x, yi, yj, i) for f in self.features_fun]
                tem_fs[yi] = tem_fs_yi
            self.all_features_by_x.append(tem_fs)
    def get_xy_features(self):
        """基于xy,获取当前样本特征"""
        for i, xi in enumerate(self.x):
            if i == 0:
                tem_fs = [f(x, self.start_y, self.y[i], i) for f in self.features_fun]
            else:
                tem_fs = [f(x, self.y[i-1], self.y[i], i) for f in self.features_fun]
            self.features_by_xy.append(tem_fs)

 

条件随机场的概率计算

 

与隐马尔科夫模型类似,我们需要计算条件随机场的一些概率问题,我们先将公式做些简单变形:

 

其中

 

基于此公式,我们也定义出条件随机场的前向,后向概率计算公式。

 

前向计算

 

由此可以得到代码:

 

class CrfCal(object):
    def __init__(self, crf_input: CrfInput, w):
        self.crf_input = crf_input
        self.start_y = -1    # 初始状态
        self.w = w    # 特征权重
        self.z_w = 1    # 分母
        self.alphas = []
        self.betas = []
        self.cal()
    def cal(self):
        self.alpha()
        self.beta()
    def alpha(self):
        """
        前向计算
        :param fs: 特征,[{j:{i:feats}}]
        :return:z_w, alphas:[t[i]]
        """
        for t, tem_f in enumerate(self.crf_input.all_features_by_x):
            if t == 0:
                alpha = [exp_dot(self.w, tem_f[-1][i]) for i in self.crf_input.labels]
            else:
                alpha = []
                for i in self.crf_input.labels:
                    new_alpha_i = sum([self.alphas[t-1][j] * exp_dot(self.w, tem_f[j][i]) for j in self.crf_input.labels])
                    alpha.append(new_alpha_i)
            self.alphas.append(alpha)
        self.z_w = sum(self.alphas[-1])

 

后向计算

 

其中后向计算中的计算与李航老师的统计学习方法中不同,应该是书中错误。由此可得到代码:

 

def beta(self):
        """
        后向计算
        :param fs: 特征,[{j:{i:feats}}]
        :return:z_w, betas:[t[i]]
        """
        self.betas = [[0] for i in self.crf_input.all_features_by_x]
        T = len(self.crf_input.x)
        for t in reversed(range(T)):
            if t == T - 1:
                beta = [1 for i in self.crf_input.labels]
            else:
                tem_f = self.crf_input.all_features_by_x[t + 1]
                beta = []
                for i in self.crf_input.labels:
                    new_beta_i = sum([exp_dot(self.w, tem_f[i][j]) * self.betas[t+1][j] for j in self.crf_input.labels])
                    beta.append(new_beta_i)
            self.betas[t] = beta
        # todo beta计算z_w,书上有误
        z_w = sum([exp_dot(self.w, self.crf_input.all_features_by_x[0][-1][j]) * self.betas[0][j] for j in self.crf_input.labels])
        return z_w

 

概率计算

 

按照前向-后向计算方法,可以计算条件概率:

 

由此可得到代码:

 

def p_yi_t(self, t):
        yi = self.crf_input.y[t]
        p = self.alphas[t][yi] * self.betas[t][yi]/self.z_w
        return p
    def p_yij_t(self, t):
        yi = self.start_y if t == 0 else self.crf_input.y[t - 1]
        yj = self.crf_input.y[t]
        if t == 0:
            p = exp_dot(self.w, self.crf_input.all_features_by_x[t][yi][yj]) * self.betas[t][yj] / self.z_w
        else:
            p = self.alphas[t-1][yi] * exp_dot(self.w, self.crf_input.all_features_by_x[t][yi][yj]) * self.betas[t][yj] / self.z_w
        return p

 

期望计算

 

同样利用前向后向计算方法,可以计算特征函数关于不同分布的数学期望。特征函数关于条件分布P(Y|X)的数学期望:

 

假设经验分布为,特征函数关于联合分布P(X,Y)的数学期望是:

 

上面引用部分,是统计学习方法一书中给出,而我认为这样计算期望不太正确,应该按照下面的方式计算:

 

由此可得到代码:

 

def w_f_yx_hat(self):
        """w*f在 \hat P(x,y)下的期望, BFGS loss的一部分"""
        wf_hat = [0 for i in self.crf_input.features_fun]
        for t, fs_xy_t in enumerate(self.crf_input.features_by_xy):
            for k, fk in enumerate(fs_xy_t):
                diff = np.dot(self.w, fk)
                wf_hat[k] += diff
        return wf_hat
    def pw(self):
        """crf概率"""
        pw = exp_dot(self.w, self.crf_input.F) / self.z_w
        return pw
    def f_yx_pw(self):
        """BFGS 梯度的一部分,上面公式没有对x概率求和部分的实现"""
        pw = self.pw()
        f_pw = [0 for i in self.crf_input.features_fun]
        for t, fs_xy_t in enumerate(self.crf_input.features_by_xy):
            for k, fk in enumerate(fs_xy_t):
                diff = pw * fk
                f_pw[k] += diff
        return f_pw

 

条件随机场的预测算法

 

跟隐马尔科夫模型类似,条件随机场也使用维特比算法来预测,根据定义我们需要计算:

 

这样条件随机场的预测问题就变成了求非规范化概率最大的最优路径问题。

 

条件随机场的维特比算法

 

输入:模型特征向量F(y,x),权重W,观测序列x

 

输出:最优路径

 

1) 初始化

 

2)递推,对t=2,3,…,T

 

3)终止

 

4)路径回溯

 

由此可得到算法:

 

class Crf(object):
    def __init__(self, xs, ys, label_num, features_fun):
        self.xs = xs
        self.ys = ys
        self.labels = range(label_num)
        self.label_num = label_num
        self.features_fun = features_fun    # 特征函数集合
        self.start_label = -1
        self.w = [0 for i in self.features_fun]
        self.crf_inputs = [CrfInput(x, y, self.label_num, self.features_fun) for x, y in zip(self.xs, self.ys)]
    def predict(self, x):
        """维特比预测算法"""
        deltas = []
        phis = []
        result_y = [1 for i in x]
        crf_input = CrfInput(x, result_y, self.label_num, self.features_fun, True)
        for i, xi in enumerate(x):
            delta, phi = [0 for i in self.labels], {}
            if i == 0:
                for yi in self.labels:
                    delta[yi] = np.dot(self.w, crf_input.all_features_by_x[i][self.start_label][yi])
                    phi[yi] = 0
            else:
                for yi in self.labels:
                    tem_delta = [deltas[i-1][yj] + np.dot(self.w, crf_input.all_features_by_x[i][yj][yi]) for yj in self.labels]
                    max_index = int(np.argmax(tem_delta))
                    delta[yi] = tem_delta[max_index]
                    phi[yi] = max_index
            deltas.append(delta)
            phis.append(phi)
        y_t = int(np.argmax(deltas[-1]))
        result_y[-1] = y_t
        for t in reversed(range(len(x)-1)):
            y_t1 = result_y[t+1]
            result_y[t] = phis[t+1][y_t1]
        return result_y

 

条件随机场的学习算法

 

根据模型,首先计算对数似然概率:

 

改进的迭代尺度法

 

改进的迭代尺度法(IIS)的基本思想是:假设当前参数向量是w,希望找到一个新的参数向量,使得模型的对数似然函数值增加,改进的迭代尺度法就是找到这样一种方法,不断更新参数,且每次参数更新能够一个变量一个变量更新,即一次只更新。

 

为了达到目的,我们先计算:

 

根据不等式:,我们得到:

 

为了能一次只优化一个变量,我们希望的下界对的导数只包含,所以此时的下界还不满足,我幺需要继续缩小下界。我们设,可以得到:

 

根据jensen不等式得:

 

所以得到:

 

易证是凸函数,则根据其导数等于0可以计算B的最大值,其导数是:

 

令其等于0得:

 

由此得到具体算法:

 

输入:特征函数,数据X,Y

 

输出:最优参数

 

1、根据X,Y计算经验分布

 

2、初始化

 

3、for i in [1,2,…,n]

 

令是方程

 

的解(可用牛顿法求解),更新

 

4、如果所有都收敛,则返回,否则重复3。

 

另一方面,如果,M是常数,推导卡壳了,先放着。

 

由此可得代码:

 

def e_f_yx_hat(self):
        Fs = [crf_input.F for crf_input in self.crf_inputs]
        ef = np.mean(Fs, 0)
        return ef
    def newton_method(self, fx, gx):
        """
        牛顿法求fx=0的解
        :param fx: 方程
        :param gx: 梯度
        :return:
        """
        x = 0
        while True:
            fx_i = fx(x)
            if abs(fx_i) < 0.001:
                break
            gx_i = gx(x)
            x = x - fx_i/gx_i
        return x
    def iis(self):
        """改进的迭代尺度学习算法"""
        ef_hat = self.e_f_yx_hat()
        f_sharps = [crf_input.f_sharp for crf_input in self.crf_inputs]
        done = False
        cycle = 0
        while not done:
            cycle += 1
            done = True
            crf_cals = [CrfCal(crf_input, self.w) for crf_input in self.crf_inputs]
            pws = [crf_cal.pw() for crf_cal in crf_cals]
            Fs = [crf_input.F for crf_input in self.crf_inputs]
            for k in range(len(self.features_fun)):
                f_dk = lambda dk: np.mean([pw * F[k] * np.exp(dk * f_sharp)
                                       for pw, F, f_sharp in zip(pws, Fs, f_sharps)], 0) - ef_hat[k]
                g_dk = lambda dk: np.mean([pw * F[k] * np.exp(dk * f_sharp) * f_sharp
                                       for pw, F, f_sharp in zip(pws, Fs, f_sharps)], 0)
                dk = self.newton_method(f_dk, g_dk)
                # dk2 = ef_hat[k]/(sum([pw * np.mean([tf[k] for tf in fs_xy]) for pw, fs_xy in zip(pws, fs_xys)], 0))
                # dk2 /= 7
                # print(dk, dk2)
                self.w[k] += dk
                if abs(dk) > 0.0001:
                    done = False
            if cycle > 160:
                break

 

拟牛顿法BFGS

 

拟牛顿法是牛顿法的近似,为了减少海赛矩阵的逆运算而得到的。根据上面我们已经得到目标函数:

 

其梯度是:

 

基于此,我们得到具体算法:

 

输入:特征函数,数据集X,Y,目标函数 loss(w) ,梯度g(w),精度

 

输出:最优参数

 

0、根据X,Y计算经验分布

 

1、选定初始点,取为正定对称矩阵,设k=0

 

2、,若,则停止计算,;否则转3

 

3、由,求出

 

4、一维搜索:求使得:

 

5、

 

6、计算,若,则停止计算,;否则计算:

 

实际只要计算(使用两次Sherman-Morrison公式):

 

其中:

 

7、k=k+1,转3

 

其中的更新推导【2】如下图所示(图中H表示B):

 

由此可得到代码:

 

def loss(self, w):
        """BFGS loss函数"""
        n = len(self.xs)
        crf_cals = [CrfCal(crf_input, w) for crf_input in self.crf_inputs]
        loss = (np.sum([math.log(crf_cal.z_w) for crf_cal in crf_cals])
                -np.sum([crf_cal.w_f_yx_hat() for crf_cal in crf_cals]))/n
        return loss
    def gradient(self, w):
        """BFGS loss的梯度函数"""
        n = len(self.xs)
        crf_cals = [CrfCal(crf_input, w) for crf_input in self.crf_inputs]
        g = (np.sum([crf_cal.f_yx_pw() for crf_cal in crf_cals], 0)
             - np.sum([crf_input.F for crf_input in self.crf_inputs], 0)) / n
        return g
    def bfgs(self):
        """拟牛顿法BFGS学习算法"""
        loss_w = self.loss
        g_w = self.gradient
        I = np.identity(len(self.w))
        D = np.identity(len(self.w))
        g_k = g_w(self.w)
        cycle = 0
        while True:
            cycle += 1
            print(self.loss(self.w))
            if np.linalg.norm(g_k) < 0.0001:
                return
            pk = -np.dot(D, g_k)
            lambda_k = self.linear_search(loss_w, pk)
            delta = lambda_k * pk
            self.w += delta
            if cycle > 20:
                return
            g_k_1 = g_w(self.w)
            if np.linalg.norm(g_k_1) < 0.0001:
                return
            # 更新D
            y = g_k_1 - g_k
            n = np.dot(delta, y)    # 分母
            m = np.outer(delta, y)    # 分子
            D = np.dot(
                np.dot(I-m/n, D),
                np.transpose(I - m/n)
                )+np.outer(delta, delta)/n
            g_k = g_k_1
    def linear_search(self, fx, x):
        """
        tree.py有类似实现,非凸时不一定最优
        黄金分割搜索法,假设区间在0.00001与0.1之间
        """
        a = 0.00001
        b = 0.1
        e = 0.00001
        while b - a > e:
            loss_a1 = fx(self.w + a * x)
            loss_b1 = fx(self.w + b * x)
            if loss_a1 < loss_b1:
                b = a * 0.382 + b * 0.618
            else:
                a = a * 0.618 + b * 0.382
        return (a + b)/2

 

参考

 

【1】李航《统计学习方法》最大熵模型p(y|x)推导的正确过程

 

【2】BFGS中B的逆矩阵迭代公式推导:附录A

Be First to Comment

发表评论

电子邮件地址不会被公开。 必填项已用*标注