机器学习算法解析—逻辑回归分类

逻辑回归模型属于广义线性回归模型，也是一种常用的二分类模型，模型假设输入向量 x 与分类概率 P(y | x; θ) 符合逻辑线性关系，并可用极大似然或极小化损失函数来估计模型的参数θ。

[欢迎访问我的博客：http://blog.csdn.net/jacobzende ]

---

分类公式:

1. 给定线性函数：

2. 给定逻辑函数：


将 z∈（-∞, +∞）值映射到［0，1］
3. 定义复合函数：

4. 定义分类公式：

合并写为：

5. 概率分割点为0.5=h(x)=g(z)，对应z＝0的线性超平面即为模型的分割面。分类公式可以很好的将点到分割面的距离转换为概率，距离越远属于相应类别的概率越大，按指数衰减。

---

参数估计：

1、参数估计方法

介绍逻辑回归模型参数估计的两种方法：极大似然估计、极小损失函数估计。
（1）极大似然估计：
给定m个训练样本，θ的似然函数等于各个样本于标注类别的概率的积：

极大似然估计即为求解逻辑似然函数的最优化问题：

（2）极小损失函数估计
逻辑回归采用log损失函数：

函数图如下，显然，预测与实际结果越不一致，损失函数的值越大：

线性回归模型采用的差方损失函数（OLS/最小二乘法）对逻辑回归模型不适用，由于逻辑回归公式是非线性的，故差方损失函数为非凸函数（非线性的二次型），最优化求解较为困难。log损失函数常用于概率输出类模型，并可以证明上述逻辑回归模型的log损失函数一定是凸函数。

由于 y 取值｛0，1｝，可将上述函数合并写为：

扩展到m个训练样本，该损失函数正好是逻辑似然函数负数的平均值，也为log损失函数提供了概率一致性的解释：

极小损失函数最优化问题即为求解：


如果定义 y 取值｛－1，＋1｝，可以得到更为紧凑的损失函数形式：

上述推导用到 g(z) = 1 - g(-z)，这个很容易证明，可以自己推导一下。

2、正则化（regularization）

为损失函数增加参数的L1或L2范数惩罚项来消去不明显的或冗余的特征以避免过拟合，L1及L2正则化惩罚项分别定义为：

L1正则化用参数的1-范数乘以系数 λ 作为惩罚项，L2正则化用参数的2-范数的平方乘以系数 λ / 2 作为惩罚项（此处除于2只有数学公式上的意义）。

另外由于截距项 theta0 没有对应的特征项(x0=1)，其目的是为了适配数据的整体偏移，因做为保留项特殊处理，不需要参与正则化。

对于 y 取值 {-1, +1}的形式，L1及L2正则化的损失函数分别为：

对上述损失函数求极小化：

第二步等式中乘以一个标量常数因子不改变最优化问题的解，最终将优化公式转化为常见的形式。

截距项不参与正则化，截距项不同于普通的特征项（可以理解为截距项的特征值＝1），截距项的作用是拟合训练样本的整体偏差（有利于其他特征项参数的正则化），而正则化主要是解决过拟合问题，因此对截距项惩罚有害无益。

---

算法：随机梯度下降

相比拟牛顿算法，随机梯度下降(SGD)算法是一种简单有效的凸函数非限制优化算法，适用于逻辑回归损失函数的最小化参数估计。

import numpy as np# 假设给定训练样本X = np.array([[-1, -1], [-2, -1], [1, 1], [2, 1]])Y = np.array([1, 1, -1, -1])# 测试样本NewX = np.array([-0.8,-1])

算法的最优参数求解过程如下：
- 初始化模型的参数θ为任意指定的值（通常为0）
- 向使损失函数J(θ)下降最深的方向改变θ的值，使得J(θ)变小
- 对每个样本执行一次上述下降过程（即上步J(θ)是单个样本的损失函数）
- 反复迭代整个样本集(设迭代次数为 iter)，直到J(θ)得到最小值

θ的更新规则为：

因函数的负的梯度（即导数）方向是函数值下降最快的方向，其中θ的n个维度是同步更新的，alpha为学习率。

# 迭代循环处理每个样本iter = 5n_samples = X.shape[0]for nIt in range(iter):    for i in range(n_samples):        #更新权重参数θ

1、计算J(θ)的导数
也即计算对于θ各个分量的偏导数，根据前面的公式，J(θ)的导数等于训练样本的成本函数的导数的均值，对于单个样本，即是求成本函数的偏导数：

（工程实现中，通常对上式分母 1 + exp(yθx) 做一些做近似处理，即可避免数值溢出也能加快计算速度，如：当 yθx > 20，则exp(yθx)达到极大值，则分母中的 1 可忽略不计，上式简化为 －yx exp(-yθx)；当 yθx < -20，则exp(yθx)趋近于0，则分母可近似＝1）

L2惩罚项的导数：

L1惩罚项无法直接求导（绝对值函数为非连续函数），可写为分段函数分别求导：

综上，对于单个样本，即J(θ)函数中 m＝1，代入上述求导公式，得到具体的参数更新公式：

（1）无正则化

# file_reg_none.pyimport numpy as np# 假设给定训练样本X = np.array([[-1, -1], [-2, -1], [1, 1], [2, 1]])Y = np.array([1, 1, -1, -1])# 测试样本NewX = np.array([-0.8,-1])# 初始化权重参数θ为0theta_intercept ＝ np.zeros(1)theta_features = np.zeros(2)# 初始化学习率alpha = 10# 迭代循环处理每个样本iter = 5n_samples = X.shape[0]for nIt in range(iter):    for i in range(n_samples):        #更新权重参数θ        p = np.dot(theta_features, X[i]) + theta_intercept        theta_features  -= alpha * ((-Y[i]*X[i]) / (1.0 + np.exp(Y[i]*p)))        theta_intercept -= alpha * ((-Y[i]*1) / (1.0 + np.exp(Y[i]*p)))# 输出权重参数print(theta_intercept, theta_features)# (array([ 4.75452057]), array([-5.24672358, -5.24550084]))# 判别测试score = np.dot(theta_features, NewX) + theta_interceptprint score, "1" if score>0 else "0"# [ 14.19740028] 1# 判别测试(概率方式：prob为正类的概率，1-prob为负类的概率）prob = np.reciprocal(1 + np.exp(-1*score))  #reciprocal: 1/xprint np.array([1 - prob, prob]).T, "1" if prob>=0.5 else "0"# [[  6.82569859e-07   9.99999317e-01]] 1

（2）L2正则化

L2正则化的更新公式直观上解释，就是每次更新前对theta j做一个极小的缩小惩罚（例如：1-αλ = 0.99)，如果 λ＝0 则与无正则化惩罚的更新公式相同。

# file_reg_L2.pyimport numpy as np# 假设给定训练样本X = np.array([[-1, -1], [-2, -1], [1, 1], [2, 1]])Y = np.array([1, 1, -1, -1])# 测试样本NewX = np.array([-0.8,-1])# 初始化权重参数θ为0theta_intercept ＝ np.zeros(1)theta_features = np.zeros(2)# 初始化学习率alpha = 10#初始化正则化系数 λlumda = 0.0001# 迭代循环处理每个样本iter = 5n_samples = X.shape[0]for nIt in range(iter):    for i in range(n_samples):        #更新权重参数θ        p = np.dot(theta_features, X[i]) + theta_intercept        r = max(0.0,  (1.0 - alpha*lumda))        theta_features  = theta_features  * r - alpha * ((-Y[i]*X[i]) / (1.0 + np.exp(Y[i]*p)))        theta_intercept -= alpha * ((-Y[i]*1) / (1.0 + np.exp(Y[i]*p)))# 输出权重参数print(theta_intercept, theta_features)# (array([ 4.73725143]), array([-5.16735701, -5.16601054]))# 判别测试score = np.dot(theta_features, NewX) + theta_interceptprint score, "1" if score>0 else "0"# [ 14.03714758] 1# 判别测试(概率方式：prob为正类的概率，1-prob为负类的概率）prob = np.reciprocal(1 + np.exp(-1*score))  #reciprocal: 1/xprint np.array([1 - prob, prob]).T, "1" if prob>=0.5 else "0"# [[  8.01205493e-07   9.99999199e-01]] 1

（3）L1正则化
定义符号函数sign(x)：当x>0时sign(x)=1，当x<0时sign(x)=—1，当x=0时sign(x)=0。代入L1惩罚项的梯度公式，L1参数更新公式可写为：

虽然上述公式中梯度计算本身没有问题，但是SGD算法本身的近似性(每次更新用随机性较强的单个样本的梯度来近似替代整体样本的梯度），会导致任何 θj 都无法如期望收敛到0（会在0附近震荡），无法实现L1正则化的降维优化效果。

一种解决方案就是对上式中最后的L1惩罚项进行截取限制：当前两项的计算结果大于零的时候，减去αλ，如果 θj < 0 则 赋值 θj＝0；反之当前两项计算结果小于零时，加上αλ，如果 θj > 0 则 赋值 θj＝0，更新公式可写为：

该方案通过避免正负震荡，增加了 θj 收敛到 0 的可能性，但是由于单个样本梯度值的波动性非常大（整体样本的梯度应该是逐步收敛到0）， θj 可能重新远离 0 值，L1惩罚需要通过缓慢的过程才能再次将其拉回 0 值，这必定导致训练结束前重新离开 0 值的 θj 无法收敛到0，一个改进的解决办法就是把被截取的惩罚项的值累计起来加到下一次更新的惩罚项中，以保证把从 0 值离开的 θj 以更快的速度拉回 0 值（远离 0 值的 θj 没有被截取过，因此不受这一个过程的影响），更新公式可写为：

上市中u代表累计可以生效的惩罚，q表示累计已经生效的惩罚，初值均赋0。

# file_reg_L1.pyimport numpy as np# 假设给定训练样本X = np.array([[-1, -1], [-2, -1], [1, 1], [2, 1]])Y = np.array([1, 1, -1, -1])# 测试样本NewX = np.array([-0.8,-1])# 初始化权重参数θ为0theta_intercept ＝ np.zeros(1)theta_features = np.zeros(2)# 初始化学习率alpha = 10#初始化正则化系数 λlumda = 0.0001#初始化累计变量u = 0q = np.zeros(2)# 迭代循环处理每个样本iter = 5n_samples = X.shape[0]for nIt in range(iter):    for i in range(n_samples):        #更新权重参数θ        u += alpha * lumda        p = np.dot(theta_features, X[i]) + theta_intercept        theta_features  -= alpha * ((-Y[i]*X[i]) / (1.0 + np.exp(Y[i]*p)))        theta_intercept -= alpha * ((-Y[i]*1) / (1.0 + np.exp(Y[i]*p)))        for j in range(2):            z = theta_features[j]            if theta_features[j] > 0.0:                theta_features[j] = max(0.0, theta_features[j] - ((u + q[j])))            elif theta_features[j] < 0.0:                theta_features[j] = min(0.0, theta_features[j] + ((u - q[j])))            q[j] += theta_features[j] - z# 输出权重参数print(theta_intercept, theta_features)# (array([ 4.75086043]), array([-5.23041214, -5.22916138]))# 判别测试score = np.dot(theta_features, NewX) + theta_interceptprint score, "1" if score>0 else "0"# [ 14.16435152] 1# 判别测试(概率方式：prob为正类的概率，1-prob为负类的概率）prob = np.reciprocal(1 + np.exp(-1*score))  #reciprocal: 1/xprint np.array([1 - prob, prob]).T, "1" if prob>=0.5 else "0"# [[  7.05504827e-07   9.99999294e-01]] 1

2、学习率 α 取值策略
（1）固定值
这种方案优点是实现简单，缺点主要是：步长取小了收敛非常缓慢，步长取大了会出现震荡，影响收敛的精度，如果取一个不慢不快的步长，则在快速下降的区域显得缓慢、在平缓下降的底部区域又可能引起震荡。
（2）逐步衰减
包括线性衰减和非线性衰减两大类。线性衰减策略为：每执行一定次数的更新（如：1000次）就将学习率缩小一定的比例（如：乘以0.9）。非线性的衰减会采用一些启发式的经验公式来保证学习率更好的适应当前位置的曲率（如：二阶导数等），此处不再展开讨论。

[我的更多博客：http://blog.csdn.net/jacobzende ]
3、补充说明
1. ElasticNet（即：组合使用L1、L2正则化）的更新公式＝L2的占比＊L2的更新量 ＋ （1—L2的占比）＊L1的更新量。
2. 正则化系数以及迭代的次数的设置会显著的影响模型的效果，每次迭代前重新打乱训练样本的顺序有利于模型的收敛。
3. 对于稀疏类型的训练样本，由于截距项的更新频次远高于各特征项，因此截距项应该使用一个更小的学习率（如：乘以0.01），以平衡两者的更新强度。
4. 如果各个特征量级不同，对特征做归一化(Scaling)可以显著的加快算法的收敛速度。
5. SGD算法的优点是高效、易于实现和调优，相比起计算复杂度较高的lbfgs、newton-cg算法，其非常适用于大样本高维度的应用场景。
6. 凸函数的最大参数估计问题可采用随机梯度上升算法求解（参数更新的方向与SGD相反）。
7. LR可以扩展为多分类模型（如：CRF、MaxEnt、Softmax模型等），可以通过多个二分类的LR来组合求解（OVA：one versus all），不过直接求解多分类模型效果会更好，我们在相关的章节再具体介绍.