逻辑回归算法原理_逻辑回归与二阶优化算法

本文主要讨论二阶优化算法，逻辑回归只是一个用来帮助实现的手段而已，不会过多讨论。在之前的文章中有介绍过神经网络的基本单元——神经元

永远在你身后：Numpy实现神经网络框架(1)​zhuanlan.zhihu.com

不考虑激活函数的话，那么它就是一个线性回归的模型

是输入向量，包含样本的特征值，

权重向量。线性回归估算的是一个连续变量的

条件期望

说简单点就是，用已有的数据去拟合一条“直线”，然后通过该直线的解析式来计算新的输入的结果（条件期望）

可以看到，这里的输出对应的是具体的各个数值，是“连续的”。举个栗子，给定房屋面积，预测房价。面积的取值范围可以是从1到10000（平米）等等，房价可以是1万，1亿或者更高

但是有时候碰到一些分类任务，例如给定温湿度，风力等因素，预测第二天是否下雨，那么我们就希望获得一个概率。此时的输出的条件期望就不应该再是1,2,3...这样数值了，而应该是0~1之间的一个数来表示第二天下雨的概率

逻辑回归

这里使用

来表示给定条件下（

），第二天下雨（

）的概率，所以，此时要对线性回归的输出——

——做一个变换，符合要求的一个就是Logistic变换，公式如下

将该式两边取

的指数，化简（为方便，下面将

用

表示）：

令

，就得到了熟悉的sigmoid函数：

有的资料会将上式右边分子分母同除

：

其函数图像如下所示

可以看到其输出的范围是(0,1)，符合我们的要求，而这种对线性回归的输出进行Logistic变换的回归也就是逻辑回归。现在还需要解决的一个问题就是，如何更新参数？

极大似然估计

我们使用极大似然估计的方法来求解在观测到的结果的条件下最为可能的参数（关于极大似然估计更多的解释知乎上有很多）。然后是前面所讨论的正例的概率如下：

相应的，反例概率就是：

综合成一般形式：

然后定义似然函数

上式是连乘形式，不利于微分，所以对两边取自然对数，得到对数似然函数

上面（1.2）到（1.3）是套用了（1.1）式，而（1.3）到（1.4）是因为：

所以：

这里我们要求的是似然值最大的参数，但是使用损失函数所要求的是损失最小，所以将上面的对数似然函数加一个负号，就得到了损失函数。又因为有

个样本，所以损失再取均值：

最后，推导损失关于

每一个分量的偏导

以上就是逻辑回归的基本介绍，下面上代码，基本就是照着公式敲的，不多注释了

import numpy as np
import matplotlib.pyplot as plt

def sigmoid(z):     # 分情况处理防止溢出
    return (1.0 - 1/(1 + np.exp(z))) if z < 0 else 1/(1 + np.exp(-z))

def loss(y, z):     # 分情况处理防止溢出
    return z*(1-y) if z > 100 else (np.log(1+np.exp(z)) - y*z)

class Logistic(object):
    def __init__(self, d, lr=1e-2):
        # d: 样本的维数
        self.w = np.random.randn(d, 1)
        self.b = 0
        self.lr = lr
        self.max_niter = 1000   # 最大迭代次数
        self.eps = 1e-3         # 损失减小的阈值
        self.sigmoid = np.vectorize(sigmoid)
        self.loss = np.vectorize(loss)

    def fit(self, X, Y):
        # X: 矩阵，每一行是一个样本
        pre_loss = float('inf')
        for _ in range(self.max_niter):
            for x, y in zip(X, Y):
                gradient = y - self.sigmoid(np.dot(x, self.w) + self.b)
                self.b += self.lr * gradient
                self.w += (self.lr * gradient * x).reshape(-1, 1)

            A = np.dot(X, self.w) + self.b
            loss = self.loss(Y, A).mean()       # 计算平均损失
            if (pre_loss - loss) < self.eps:    # 上一次损失减去当前损失小于一定阈值则停止迭代
                print(_)
                break
            pre_loss = loss

    def predict(self, X):
        '''
        X: 待预测的样本矩阵，每行一个样本
        '''
        return self.sigmoid(np.dot(X, self.w) + self.b)

以上是逻辑回归的实现代码，并且对sigmoid函数和损失函数进行了处理，因为当

为一个很大的正数时，

会溢出，所以可以将

分情况进行处理，如下式：

然后是针对损失函数防止溢出，把（1.5）拿过来，改成计算单个样本损失的形式：

这样就避免了当

时，

函数会出现下溢出。并且，当

时

下面给出要给一个具体例子的代码：

def plot_decision_boundary(model, axis):    # 绘制决策边界
    x0, x1 = np.meshgrid(
        np.linspace(axis[0], axis[1], int((axis[1]-axis[0])*10)),
        np.linspace(axis[2], axis[3], int((axis[3]-axis[2])*10))
    )
    X_new = np.vstack((x0.ravel(), x1.ravel())).T
    y_predict = model.predict(X_new)
    zz = y_predict.reshape(x0.shape)

    from matplotlib.colors import ListedColormap
    custom_cmap = ListedColormap(['#ef9a9a', '#fff59d'])
    plt.contourf(x0, x1, zz, linewidth=5, cmap=custom_cmap)

if __name__ == "__main__":
    data = np.loadtxt('./dataset/Trainset.txt', np.float)
    X = data[:,:2]
    Y = data[:,-1:]
    lr = Logistic(2)
    lr.fit(X, Y)

    pos, neg = [], []
    for x, y in zip(X, Y):
        if y == 1: pos.append(x)
        else: neg.append(x)

    pos = np.array(pos)
    neg = np.array(neg)

    plt.figure()
    plot_decision_boundary(lr, axis=[-2.5, 3.5, -4.5, 16])
    plt.scatter(pos[:,0], pos[:,1], c='r')
    plt.scatter(neg[:,0], neg[:,1], c='b')
    plt.show()

Trainset.txt的内容如下：

-0.386323   3.989286    1
0.556921    8.294984    1
1.224863    11.587360   0
-1.347803   -2.406051   1
1.196604    4.951851    1
0.275221    9.543647    0
0.470575    9.332488    0
-1.889567   9.542662    0
-1.527893   12.150579   0
-1.185247   11.309318   0
-0.445678   3.297303    1
1.042222    6.105155    1
-0.618787   10.320986   0
1.152083    0.548467    1
0.828534    2.676045    1
-1.237728   10.549033   0
-0.683565   -2.166125   1
0.229456    5.921938    1
-0.959885   11.555336   0
0.492911    10.993324   0
0.184992    8.721488    0
-0.355715   10.325976   0
-0.397822   8.058397    0
0.824839    13.730343   0
1.507278    5.027866    1
0.099671    6.835839    1
-0.344008   10.717485   0
1.785928    7.718645    1
-0.918801   11.560217   0
-0.364009   4.747300    1
-0.841722   4.119083    1
0.490426    1.960539    1
-0.007194   9.075792    0
0.356107    12.447863   0
0.342578    12.281162   0
-0.810823   -1.466018   1
2.530777    6.476801    1
1.296683    11.607559   0
0.475487    12.040035   0
-0.783277   11.009725   0
0.074798    11.023650   0
-1.337472   0.468339    1
-0.102781   13.763651   0
-0.147324   2.874846    1
0.518389    9.887035    0
1.015399    7.571882    0
-1.658086   -0.027255   1
1.319944    2.171228    1
2.056216    5.019981    1
-0.851633   4.375691    1
-1.510047   6.061992    0
-1.076637   -3.181888   1
1.821096    10.283990   0
3.010150    8.401766    1
-1.099458   1.688274    1
-0.834872   -1.733869   1
-0.846637   3.849075    1
1.400102    12.628781   0
1.752842    5.468166    1
0.078557    0.059736    1
0.089392    -0.715300   1
1.825662    12.693808   0
0.197445    9.744638    0
0.126117    0.922311    1
-0.679797   1.220530    1
0.677983    2.556666    1
0.761349    10.693862   0
-2.168791   0.143632    1
1.388610    9.341997    0
0.317029    14.739025   0

最终，训练结束后，决策边界绘制效果如下：

下面讨论使用二阶优化算法（牛顿法）结合讨论

牛顿法

首先，令损失函数表示如下

设损失函数连续、可导，那么在损失最小出的点，参数的梯度（导数）必然为0：

也就是所说找到梯度为0点，就找到了使损失最小的点。基于这一点，我们再来看牛顿法（也叫牛顿迭代法），它的迭代公式可以从泰勒展开式中的前几项中所推出，不过为了更加直观，所以还是先从函数图像的栗子入手

现有一函数

，要求出

的解，假设该函数图像如下

因为没有方法直接算出该函数值为0处的解析解，所以只能使用其他的方法来求近似解。如上图所示，首先在

轴上随便（不要太随便）找一个

，并计算对应的函数值

，得到坐标点

过该点做函数的切线，可以知道该切线会交

轴于

（只要切线的斜率不为0）。此时我们可以得到一个三角形，如下图所示

在该三角形中，对边比邻边就是函数在

处的斜率

通过整理上式，可以得到：

这样，我们就可以准确的得到

的具体值，从而继续进行迭代

把前面的公式整理一下，得到迭代公式

观察该式，可以发现，当：

此时

，即可停止迭代。不过因为分母不能为0，所以也就等价于

，也就是一开始题目所求。那么这和优化算法有什么关系呢？

现在假设有一函数

，该函数连续可导，现求

，套用迭代公式就是

然后，令

（导函数也是一个函数对吧），可得：

这样一直迭代到满足停止条件时，所得

就可以满足

。也就是我们优化算法最初的目的了。另外，从泰勒展开式也可推导出来，不感兴趣的可以直接略过

---

泰勒展开式：

取其二阶展开

两边对求

导

我们要求

，所以令

，得到

---

再就是具体栗子，将牛顿法应用于逻辑回归。损失关于参数的一阶偏导前面已经得到了：

然后再求二阶偏导

（2.2）到（2.3）式套用了sigmoid函数的导数公式：

其中

，再将一、二阶偏导代入（2.1）式，得到

上面的计算是针对参数的每个分量单独处理，分别计算每一个参数的一、二阶导数。由于

是一个向量，所以也可以直接使用向量来计算，首先将（2.1）式改为：

之所以把二阶导数写成倒数形式，是为了接下来的替换。由于

是向量，所以将一阶导数写成梯度形式：

再将二阶导数换成海森矩阵（详细介绍见知乎）

令

就得到了多元函数的牛顿迭代公式

此处要求海森矩阵的逆矩阵，其是否可逆不说，而且一旦样本的特征数量很大，计算量就上天了，所以需要寻求一种近似的替代方法

再回顾一下泰勒展开式的二阶展开求导后的式子

当

是向量时，上式变成

化简

令

式（2.4）变为

或

式（2.5）（2.6）称为拟牛顿条件，满足这两个条件之一的矩阵，即可用来替代

或

关于拟牛顿法的具体算法有DFP算法、BFGS算法（L-BFGS）、Broyden类算法，具体如操作在《统计学习方法》中有详细的介绍，手头没有书的可以知乎百度，这里也顺手贴个链接：统计学习方法，拟牛顿法