浏览本文须要的背景知识点:线性回归、最大似然预计、一丢丢编程常识
一、引言
后面几节咱们学习了规范线性回归,而后介绍了三种正则化的办法 - 岭回归、Lasso回归、弹性网络回归,这些线性模型解决的都是回归的问题。最开始还介绍了两种简略的算法-PLA与口袋算法,他们解决的是分类问题。
那么咱们能应用回归的形式来解决分类问题么,答案是必定的,这就是上面要介绍的模型 - 对数几率回归算法1(Logistic Regression Algorithm),也有被直译为逻辑回归。
二、模型介绍
对数几率回归的模型函数
既然要通过回归的形式来解决分类的问题,能够通过先进行回归剖析,而后通过一个函数将间断的后果映射成离散的分类后果,例如上面的函数表达式:
$$y=\left\{\begin{array}{ll}0 & \left(f\left(w^{T} x\right) \leq C\right) \\1 & \left(f\left(w^{T} x\right)>C\right)\end{array}\right.$$
在介绍对数几率回归模型之前,先来看看一类被称为 S 函数2(Sigmoid function)的函数,这类函数的图像成 S 型,其中最常见的一种函数为逻辑函数3(Logistic function),函数图像如下图所示:
通过图像能够看到,这个函数当自变量越大时,函数值越趋近于1,当自变量最小时,函数值越趋近于0,当自变量为0时,函数值为0.5。当下面表达式中的 C 等于 0.5 时,能够看作将间断的后果映射成离散的后果。
逻辑函数的函数表达式为:
$$f(z)=\frac{1}{1+e^{-z}}$$
将线性方程带入逻辑函数中,就失去了对数几率回归的函数方程:
$$f(x)=\frac{1}{1+e^{-w^Tx}}$$
对数几率回归的代价函数
咱们须要一个代价函数来示意数据拟合的状况,这时很容易想到的是应用与线性回归一样的均方差4(mean-square error/MSE)的办法来做为代价函数
$$\operatorname{Cost}(w)=\sum_{i=1}^{N}(y-\hat{y})^{2}$$
但对数几率回归是应用似然函数5(likelihood function)的对数模式来作为其代价函数,前面会阐明为什么应用这种形式比MSE更适合。
最大似然预计:
思考一个抛硬币的例子。假如这个硬币侧面跟背面轻重不同。咱们把这个硬币抛80次(侧面记为H,背面记为T)。并把抛出一个侧面的概率记为 p,抛出一个背面的概率记为 1 - p。假如咱们抛出了49个侧面,31个背面,即49次H,31次T。假如这个硬币是咱们从一个装了三个硬币的盒子外头取出的。这三个硬币抛出侧面的概率别离为 1/3、1/2、2/3 。这些硬币没有标记,所以咱们无奈晓得哪个是哪个。应用最大似然预计,基于二项分布中的概率品质函数公式,通过这些试验数据(即采样数据),咱们能够计算出哪个硬币的可能性最大。咱们能够看到当 p = 2/3 时,似然函数获得值最大。
$$\begin{array}{l}L\left(p=\frac{1}{3} \mid H=49, T=31\right)=C_{80}^{49}\left(\frac{1}{3}\right)^{49}\left(1-\frac{1}{3}\right)^{31} \approx 0.000 \\L\left(p=\frac{1}{2} \mid H=49, T=31\right)=C_{80}^{49}\left(\frac{1}{2}\right)^{49}\left(1-\frac{1}{2}\right)^{31} \approx 0.012 \\L\left(p=\frac{2}{3} \mid H=49, T=31\right)=C_{80}^{49}\left(\frac{2}{3}\right)^{49}\left(1-\frac{2}{3}\right)^{31} \approx 0.054\end{array}$$
第一种示意形式:
在二元分类的状况下,例如 y 取 0 或者 1,将对数几率回归的函数方程的后果看作概率,能够写作下式:
$$\left\{\begin{array}{ll}P(y=1 \mid x, w)= & h(x) \\P(y=0 \mid x, w)= & 1-h(x)\end{array}\right.$$
因为 y 只能取 0 或者 1,能够将下面两个式写成一个式子:
$$P(y \mid x, w)=h(x)^{y}+[1-h(x)]^{1-y}$$
写出似然函数,其中 N 为样本总数量,将每一个概率累乘就失去了似然函数:
$$L(w)=\prod_{i=1}^{N} h\left(X_{i}\right)^{y_{i}}\left[1-h\left(X_{i}\right)\right]^{1-y_{i}}$$
在 w 的所有可能取值中找一个使得似然函数取到最大值,这时求得的解就是 w 的最大似然预计6(maximum likelihood estimation/MLE)
$$w=\underset{w}{\operatorname{argmax}}(L(w))=\underset{w}{\operatorname{argmax}}\left(\prod_{i=1}^{N} h\left(X_{i}\right)^{y_{i}}\left[1-h\left(X_{i}\right)\right]^{1-y_{i}}\right)$$
因为带连乘运算的代价函数不好优化,这里咱们对似然函数取自然对数并且取反,S 函数的取值为(0,1),似然函数的取值也是(0,1),对其取对数不影响其枯燥性。这样就失去了对数几率回归的代价函数:
$$\begin{aligned}\operatorname{Cost}(w) &=-\ln L(w) \\&=-\sum_{i=1}^{N}\left(y_{i} \ln h\left(X_{i}\right)+\left(1-y_{i}\right) \ln \left(1-h\left(X_{i}\right)\right)\right)\end{aligned}$$
因为加了一步取反的操作,这是就不是求最大值,而是将其改为求最小值:
$$w=\underset{w}{\operatorname{argmin}}\left(-\sum_{i=1}^{N}\left(y_{i} \ln h\left(X_{i}\right)+\left(1-y_{i}\right) \ln \left(1-h\left(X_{i}\right)\right)\right)\right)$$
第二种示意形式:
咱们先来看下 S 函数的一个性质:
$$\begin{aligned}f(z) &=\frac{1}{1+e^{-z}}=\frac{e^{z}}{1+e^{z}} \\f(-z) &=\frac{1}{1+e^{z}}=1-f(z)\end{aligned}$$
在二元分类的状况下, 当 y 的值取 -1 或者 1时,将对数几率回归的函数方程的后果看作概率,能够写作下式
$$\left\{\begin{array}{l}P(y=1 \mid x, w)=h(x)=\frac{1}{1+e^{-w^{T} x}} \\P(y=-1 \mid x, w)=1-h(x)=h(-x)=\frac{1}{1+e^{w^{T} x}}\end{array}\right.$$
因为 y 只能取 -1 或者 1,能够将下面两个式写成一个式子:
$$P(y \mid x, w)=\frac{1}{1+e^{-y w^{T} x}}$$
写出似然函数,其中 N 为样本总数量,将每一个概率累乘就失去了似然函数:
$$L(w)=\prod_{i=1}^{N} \frac{1}{1+e^{-y_{i} w^{T} X_{i}}}$$
还是求最大似然预计:
$$w=\underset{w}{\operatorname{argmax}}\left(\prod_{i=1}^{N} \frac{1}{1+e^{-y_{i} w^{T} X_{i}}}\right)$$
这里与第一种形式一样,咱们对似然函数取自然对数并且取反,就失去了对数几率回归的代价函数:
$$\begin{aligned}\operatorname{Cost}(w) &=-\ln L(w) \\&=-\sum_{i=1}^{N} \ln \frac{1}{1+e^{-y_{i} w^{T} X_{i}}} \\&=-\sum_{i=1}^{N} \ln 1-\ln \left(1+e^{-y_{i} w^{T} X_{i}}\right) \\&=\sum_{i=1}^{N} \ln \left(1+e^{-y_{i} w^{T} X_{i}}\right)\end{aligned}$$
也是求代价函数的最小值:
$$w=\underset{w}{\operatorname{argmin}}\left(\sum_{i=1}^{N} \ln \left(1+e^{-y_{i} w^{T} X_{i}}\right)\right)$$
在sklearn中应用的就是下面这个代价函数。
对数几率回归的正则化:
与线性回归一样,对数几率回归的代价函数也能够加上正则化的惩办项,有两种形式来增加正则项,一个是给正则项增加系数来管制大小,另一个是给代价函数增加系数来管制大小,实质作用是一样的,在 sklearn 中应用的是(2)式中的模式。
L1 正则化:
$$\begin{aligned}\operatorname{Cost}(w)_{L 1} &=\operatorname{Cost}(w)+\lambda\|w\|_{1} & (1) \\&=C \cdot \operatorname{Cost}(w)+\|w\|_{1} & (2)\end{aligned}$$
L2 正则化:
$$\begin{aligned}\operatorname{Cost}(w)_{L2} &=\operatorname{Cost}(w)+\lambda\|w\|_{2}^2 & (1) \\&=C \cdot \operatorname{Cost}(w)+\|w\|_{2}^2 & (2)\end{aligned}$$
弹性网络正则化:
$$\begin{aligned}\operatorname{Cost}(w)_{E N} &=\operatorname{Cost}(w)+\lambda \rho\|w\|_{1}+\frac{\lambda(1-\rho)}{2}\|w\|_{2}^{2} \\&=C \cdot \operatorname{Cost}(w)+\rho\|w\|_{1}+\frac{(1-\rho)}{2}\|w\|_{2}^{2}\end{aligned}$$
对数几率回归的代价函数最小化的优化办法有多种,例如后面几节介绍过的坐标降落法。无正则化或者L2正则化的状况下能够用梯度降落法7 (Gradient descent)、牛顿法8 (Newton's method)来等等。
在 scikit-learn 中能够看到更多的优化办法,其中大多是后面提到算法的减速版本或是变体,例如随机均匀梯度降落法(Stochastic Average Gradient/SAG)、随机均匀梯度降落减速法(SAGA)、L-BFGS算法(Limited-memory Broyden–Fletcher–Goldfarb–Shanno/L-BFGS),前面会一一介绍这些算法。
三、算法步骤
梯度降落法
梯度降落法与坐标降落法的思路一样,都是通过一次一次更新权重系数 w,一步一步的迫近代价函数的最小值。坐标降落法是通过固定某一个坐标轴,求出部分最优解,而后更新 w 的值。梯度降落法令是求出函数的梯度后,沿着梯度的反方向再找到一个适合的步长系数迭代更新 w 的值。
能够了解为碗中的一个小球从某一点自在滚落,必然会沿着变动量最大的方向挪动,最初达到最低点。
<center>坐标降落法 VS 梯度降落法</center>
具体步骤:
(1)初始化权重系数 w,例如初始化为零向量
(2)计算降落方向,梯度降落法将梯度的反方向作为降落方向
$$\Delta w = -\frac{\partial \operatorname{Cost}(w) }{\partial w} $$
(3)应用线搜寻9 办法来计算降落的步长 ,以使每一步迭代都满足Wolfe条件10
(4)更新权重系数 w
$$w=w+\alpha \Delta w$$
(5)反复步骤(2)~(4)直到梯度的大小足够小
牛顿法
牛顿法与梯度降落法一样,也是一种降落办法,两个算法的步骤大致相同,区别就在于抉择的降落方向不一样。
<center>绿色为梯度降落,红色为牛顿法,牛顿法的门路更加间接</center>
具体步骤:
(1)初始化权重系数 w,例如初始化为零向量
(2)计算降落方向,牛顿法应用二阶泰勒开展,将黑塞矩阵的逆矩阵与梯度的乘积的反方向作为降落方向
$$\Delta w=-H^{-1} \frac{\partial \operatorname{Cost}(w)}{\partial w}$$
(3)应用线搜寻9 办法来计算降落的步长 ,以使每一步迭代都满足Wolfe条件10
(4)更新权重系数 w
$$w=w+\alpha \Delta w$$
(5)反复步骤(2)~(4)直到梯度的大小足够小
四、原理证实
对数几率回归的代价函数是凸函数
凸函数除了通过后面几节中介绍的办法判断外,如果函数是二阶可导的,当函数的黑塞矩阵是半正定的,该函数就为凸函数。
多元二次可微的连续函数在凸集上是凸的,当且仅当它的黑塞矩阵在凸集的外部是半正定的。
- 第一种代价函数:
(1)将代价函数连加运算去掉失去一个新函数
(2)带入 h(x)
(3)将 e 的指数符号变换一下
(4)对数除法开展为对数的减法
(5)开展括号
(6)第二项与第四项一样,化简失去
$$\begin{aligned}f(w) &=-(y \ln h(x)+(1-y) \ln (1-h(x))) & (1)\\&=-\left(y \ln \frac{1}{1+e^{-w^{T} x}}+(1-y) \ln \left(1-\frac{1}{1+e^{-w^{T} x}}\right)\right) & (2)\\&=-\left(y \ln \frac{e^{w^{T} x}}{1+e^{w^{T} x}}+(1-y) \ln \left(\frac{1}{1+e^{w^{T} x}}\right)\right) & (3)\\&=-y\left(w^{T} x-\ln \left(1+e^{w^{T} x}\right)\right)-(1-y)\left(0-\ln \left(1+e^{w^{T} x}\right)\right) & (4) \\&=-y w^{T} x+y\left(1+e^{w^{T} x}\right)+\ln \left(1+e^{w^{T} x}\right)-y\left(1+e^{w^{T} x}\right) & (5)\\&=\ln \left(1+e^{w^{T} x}\right)-y w^{T} x & (6)\end{aligned}$$
(1)下面失去的 f(w)
(2)依据求导公式求一阶导数
(3)利用S函数的性质,将分子化简一下并合并同类项
(4)依据求导公式求二阶导数
(5)整顿一下后果,能够看到二阶导数是一个负数乘以一个由 x 向量与 x 向量组成的矩阵
$$\begin{aligned}f(w) &=\ln \left(1+e^{w^{T} x}\right)-y w^{T} x & (1)\\\frac{\partial f(w)}{\partial w} &=\frac{1}{1+e^{w^{T} x}} e^{w^{T} x} x-y x & (2)\\&=\left(\frac{1}{1+e^{-w^{T} x}}-y\right) x & (3)\\\frac{\partial^{2} f(w)}{\partial w \partial w^{T}} &=-\frac{1}{\left(1+e^{-w^{T} x}\right)^{2}} e^{-w^{T} x}x(-x^{T}) & (4)\\&=\frac{e^{-w^{T} x}}{\left(1+e^{-w^{T} x}\right)^{2}} x x^{T} & (5)\end{aligned}$$
- 第二种代价函数:
(1)将代价函数连加运算去掉失去一个新函数
(2)依据求导公式求一阶导数
(3)利用S函数的性质,将分子化简一下
(4)依据求导公式求二阶导数
(5)y的平方为1,整顿一下后果,能够看到二阶导数也是一个负数乘以一个由 x 向量与 x 向量组成的矩阵
$$\begin{aligned}f(w) &=\ln \left(1+e^{-y w^{T} x}\right) & (1)\\\frac{\partial f(w)}{\partial w} &=\frac{1}{1+e^{-y w^{T} x}} e^{-y w^{T} x}\left(-y x\right) & (2)\\&=-\frac{y}{1+e^{y w^{T} x}} x & (3)\\\frac{\partial^{2} f(w)}{\partial w \partial w^{T}} &=\frac{y}{\left(1+e^{y w^{T} x}\right)^{2}} e^{y w^{T} x} x \left(yx^{T}\right) & (4)\\&=\frac{e^{y w^{T} x}}{\left(1+e^{y w^{T} x}\right)^{2}} x x^{T} & (5)\end{aligned}$$
(1)令二阶导数为 H 矩阵
(2)矩阵正定的判断办法,带入二阶导数,两种二阶导数的常数局部都用 C 示意,并且 C>0
(3)利用转置的性质改写一下
(4)能够看到最初的后果大于等于 0,阐明 H 为半正定矩阵
$$\begin{aligned}H &=\frac{\partial^{2} f(w)}{\partial w \partial w^{T}} & (1) \\v^{T} H v &=C \cdot v^{T} x x^{T} v \quad(C>0) & (2) \\&=C \cdot\left(x^{T} v\right)^{T}\left(x^{T} v\right) & (3) \\&=C \cdot\left|x^{T} v\right|^{2} \geq 0 & (4)\end{aligned}$$
两种新函数的黑塞矩阵都是半正定的,那么新函数是凸函数,多个凸函数相加的函数仍然是凸函数,则阐明代价函数是一个凸函数,证毕。
为什么不必MSE作为代价函数11
- MSE 作为代价函数对谬误数据的惩办力度有余
当预测数据齐全谬误时(例如理论值为 0 预测值为 1,或理论值为 1 预测值为 0)
用MSE作为代价函数:
$$\text { Cost }=(1-0)^{2}=1$$
应用对数似然函数作为代价函数(第一种):
$$\text { Cost }=-(1 \cdot \ln 0+(1-1) \cdot \ln (1-0))=\infty$$
应用对数似然函数作为代价函数(第二种):
(1)理论值为 1 预测值为 -1
(2)预测值为 -1,其对应的h(x) = 0
(3)这时线性组合的值为负无穷大
(4)这时的代价为无穷大
$$\begin{array}{l}y=1, \hat{y}=-1 & (1)\\h(x)=0 & (2)\\w^{T} x=-\infty & (3)\\\text { Cost }=\ln \left(1+e^{-1 \cdot(-\infty)}\right)=\infty & (4)\end{array}$$
能够看到 MSE 的代价值远远小于对数似然函数的代价值,阐明 MSE 的代价函数不会重重地惩办谬误分类,哪怕是齐全谬误的分类。
- MSE 作为代价函数不是一个凸函数
先来看下下面 S函数导函数的一个性质,上面会用到这个性质:
$$\begin{aligned}h(z) &=\frac{1}{1+e^{-z}} \\\frac{\partial h(z)}{\partial z} &=-\frac{1}{\left(1+e^{-z}\right)^{2}} e^{-z}(-1) \\&=\frac{e^{-z}}{\left(1+e^{-z}\right)^{2}} \\&=\frac{1}{1+e^{-z}} \frac{e^{-z}}{1+e^{-z}} \\&=h(z)(1-h(z))\end{aligned}$$
(1)将 MSE 作为代价函数
(2)求一阶导数,因为间接对 w 求导较为简单,所以拆解为两步求导
(3)先对 y^ 求导,前面再对 w 求导,用到了下面说的S函数的性质,留神最初还有一个x
(4)开展括号
(5)求二阶导数,也是拆解为两步求导
(6)先对 y^ 求导,前面再对 w 求导
$$\begin{aligned}\operatorname{Cost}(w) &=(y-\hat{y})^{2} & (1)\\\frac{\partial \operatorname{Cost}(w)}{\partial w} &=\frac{\partial \operatorname{Cost}(w)}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial w} & (2)\\&=2(y-\hat{y})(-1) \hat{y}(1-\hat{y}) x & (3)\\&=-2\left(y \hat{y}-\hat{y}^{2}-y \hat{y}^{2}+\hat{y}^{3}\right) x & (4)\\\frac{\partial^{2} \operatorname{Cost}(w)}{\partial w\partial w^T} &=\frac{\partial \frac{\partial \operatorname{Cost}(w)}{\partial w}}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial w^T} & (5)\\&=-2\left(y-2 \hat{y}-2 y \hat{y}+3 \hat{y}^{2}\right) \hat{y}(1-\hat{y}) xx^{T} & (6)\end{aligned}$$
同下面证实一样,当初只须要看常数项是否是大于等于 0
(1)当 y = 1 时,带入二阶导函数中
(2)因式分解
(3)能够看到,当 y^ 取不同范畴的时候,常数项不是始终都是大于等于0的
$$\begin{aligned}g(\hat{y}) &=-2\left(1-4 \hat{y}+3 \hat{y}^{2}\right) \quad(y=1) & (1) \\&=-2(3 \hat{y}-1)(\hat{y}-1) & (2) \\& \Rightarrow\left\{\begin{array}{ll}\hat{y} \in\left[0, \frac{1}{3}\right) & g(\hat{y})<0 \\\hat{y} \in\left[\frac{1}{3}, 1\right] & g(\hat{y}) \geq 0\end{array}\right. & (3)\end{aligned}$$
(1)当 y = 0 时,带入二阶导函数中
(2)因式分解
(3)能够看到,当 y^ 取不同范畴的时候,常数项也不是始终都是大于等于0的
$$\begin{array}{rlr}g(\hat{y}) & =-2\left(3 \hat{y}^{2}-2 \hat{y}\right) \quad (y=0) & (1)\\& =-2 \hat{y}(3 \hat{y}-2) & (2)\\& \Rightarrow\left\{\begin{aligned}\hat{y} \in\left[0, \frac{2}{3}\right] & g(\hat{y}) \geq 0 \\\hat{y} \in\left(\frac{2}{3}, 1\right] & g(\hat{y})<0\end{aligned}\right. & (3)\end{array}$$
因为应用MSE作为代价函数的黑塞矩阵不能保障是半正定的,阐明其代价函数不是一个凸函数。
梯度降落法的降落方向
先来看下一元的泰勒公式12 ,其中o(...) 示意括号内式子的高阶无穷小,多元泰勒公式就是后面的一、二阶导数改写成雅可比矩阵和黑塞矩阵的模式:
$$f(x)=f(a)+\frac{f^{\prime}(a)}{1 !}(x-a)+\frac{f^{(2)}(a)}{2 !}(x-a)^{2}+\cdots+\frac{f^{(n)}(a)}{n !}(x-a)^{n}+o\left((x-a)^{n}\right)$$
我的指标是每次迭代后的代价函数值比上次都小:
$$\operatorname{Cost}(w_{n+1}) \lt \operatorname{Cost}(w_{n})$$
将对数几率回归的代价函数用多元的一阶泰勒公式近似,雅可比矩阵为代价函数的梯度:
$$\operatorname{Cost}(w) \approx \operatorname{Cost}(A)+\nabla \operatorname{Cost}(A)^{T}(w-A)$$
(1)将 w_n 带入 A,w_n+1 带入 w
(2)w_n+1 - w_n = w,移项能够看到须要保障梯度乘 w 小于 0
(3)当 w 为负的梯度时,能够保障梯度乘 w 小于 0
(4)这样就失去了梯度降落的降落方向
$$\begin{aligned}\operatorname{Cost}\left(w_{n+1}\right) & \approx \operatorname{Cost}\left(w_{n}\right)+\nabla \operatorname{Cost}\left(w_{n}\right)^{T}\left(w_{n+1}-w_{n}\right) & (1) \\\operatorname{Cost}\left(w_{n+1}\right)-\operatorname{Cost}\left(w_{n}\right) & \approx \nabla \operatorname{Cost}\left(w_{n}\right)^{T} \Delta w<0 & (2) \\\nabla \operatorname{Cost}\left(w_{n}\right)^{T} \Delta w &=-\nabla \operatorname{Cost}\left(w_{n}\right)^{T} \nabla \operatorname{Cost}\left(w_{n}\right)<0 & (3) \\\Delta w &=-\nabla \operatorname{Cost}\left(w_{n}\right) & (4) \end{aligned}$$
牛顿法的降落方向
将对数几率回归的代价函数用多元的二阶泰勒公式近似,雅可比矩阵为代价函数的梯度,H为代价函数的黑塞矩阵:
$$\operatorname{Cost}(w) \approx \operatorname{Cost}(A)+\nabla \operatorname{Cost}(A)^{T}(w-A)+\frac{1}{2}(w-A)^{T} H(w-A)$$
(1)将 w_n 带入 A,并对近似函数求梯度,另梯度等于 0 向量,这时获得以后的极小值
(2)这样就失去了牛顿法的降落方向
$$\begin{aligned}\nabla \operatorname{Cost}(w) &=\nabla \operatorname{Cost}\left(w_{n}\right)+H\left(w-w_{n}\right)=0 & (1) \\w &=w_{n}-H^{-1} \nabla \operatorname{Cost}\left(w_{n}\right) & (2)\end{aligned}$$
五、代码实现
应用 Python 实现对数几率回归算法(梯度降落法):
import numpy as npc_1 = 1e-4c_2 = 0.9def cost(X, y, w): """ 对数几率回归的代价函数 args: X - 训练数据集 y - 指标标签值 w - 权重系数 return: 代价函数值 """ power = -np.multiply(y, X.dot(w)) p1 = power[power <= 0] p2 = -power[-power < 0] # 解决 python 计算 e 的指数幂溢出的问题 return np.sum(np.log(1 + np.exp(p1))) + np.sum(np.log(1 + np.exp(p2)) - p2)def dcost(X, y, W): """ 对数几率回归的代价函数的梯度 args: X - 训练数据集 y - 指标标签值 w - 权重系数 return: 代价函数的梯度 """ return X.T.dot(np.multiply(-y, 1 / (1 + np.exp(np.multiply(y, X.dot(w))))))def direction(d): """ 更新的方向 args: d - 梯度 return: 更新的方向 """ return -ddef sufficientDecrease(X, y, w, step): """ 判断是否满足充沛降落条件(sufficient decrease condition) args: X - 训练数据集 y - 指标标签值 w - 权重系数 step - 步长 return: 是否满足充沛降落条件 """ d = dcost(X, y, w) p = direction(d) return cost(X, y, w + step * p) <= cost(X, y, w) + c_1 * step * p.T.dot(d)def curvature(X, y, w, step): """ 判断是否满足曲率条件(curvature condition) args: X - 训练数据集 y - 指标标签值 w - 权重系数 step - 步长 return: 是否满足曲率条件 """ d = dcost(X, y, w) p = direction(d) return -p.T.dot(dcost(X, y, w + step * p)) <= -c_2 * p.T.dot(d)def select(step_low, step_high): """ 在范畴内抉择一个步长,间接取中值 args: step_low - 步长范畴开始值 step_high - 步长范畴完结值 return: 步长 """ return (step_low + step_high) / 2def lineSearch(X, y, w, step_init, step_max): """ 线搜寻步长,使其满足 Wolfe 条件 args: X - 训练数据集 y - 指标标签值 w - 权重系数 step_init - 步长初始值 step_max - 步长最大值 return: 步长 """ step_i = step_init step_low = step_init step_high = step_max i = 1 d = dcost(X, y, w) p = direction(d) while (True): # 不满足充沛降落条件或者前面的代价函数值大于前一个代价函数值 if (not sufficientDecrease(X, y, w, step_i) or (cost(X, y, w + step_i * p) >= cost(X, y, w + step_low * p) and i > 1)): # 将以后步长作为搜寻的右边界 # return search(X, y, w, step_prev, step_i) step_high = step_i else: # 满足充沛降落条件并且满足曲率条件,即曾经满足了 Wolfe 条件 if (curvature(X, y, w, step_i)): # 间接返回以后步长 return step_i step_low = step_i # 抉择下一个步长 step_i = select(step_low, step_high) i = i + 1def logisticRegressionGd(X, y, max_iter=1000, tol=1e-4, step_init=0, step_max=10): """ 对数几率回归,应用梯度降落法(gradient descent) args: X - 训练数据集 y - 指标标签值 max_iter - 最大迭代次数 tol - 变动量容忍值 step_init - 步长初始值 step_max - 步长最大值 return: w - 权重系数 """ # 初始化 w 为零向量 w = np.zeros(X.shape[1]) # 开始迭代 for it in range(max_iter): # 计算梯度 d = dcost(X, y, w) # 当梯度足够小时,完结迭代 if np.linalg.norm(x=d, ord=1) <= tol: break # 应用线搜寻计算步长 step = lineSearch(X, y, w, step_init, step_max) # 更新权重系数 w w = w + step * direction(d) return w
应用 Python 实现对数几率回归算法(牛顿法):
import numpy as npc_1 = 1e-4c_2 = 0.9def cost(X, y, w): """ 对数几率回归的代价函数 args: X - 训练数据集 y - 指标标签值 w - 权重系数 return: 代价函数值 """ power = -np.multiply(y, X.dot(w)) p1 = power[power <= 0] p2 = -power[-power < 0] # 解决 python 计算 e 的指数幂溢出的问题 return np.sum(np.log(1 + np.exp(p1))) + np.sum(np.log(1 + np.exp(p2)) - p2)def dcost(X, y, w): """ 对数几率回归的代价函数的梯度 args: X - 训练数据集 y - 指标标签值 w - 权重系数 return: 代价函数的梯度 """ return X.T.dot(np.multiply(-y, 1 / (1 + np.exp(np.multiply(y, X.dot(w))))))def ddcost(X, y, w): """ 对数几率回归的代价函数的黑塞矩阵 args: X - 训练数据集 y - 指标标签值 w - 权重系数 return: 代价函数的黑塞矩阵 """ exp = np.exp(np.multiply(y, X.dot(w))) result = np.multiply(exp, 1 / np.square(1 + exp)) X_r = np.zeros(X.shape) for i in range(X.shape[1]): X_r[:, i] = np.multiply(result, X[:, i]) return X_r.T.dot(X)def direction(d, H): """ 更新的方向 args: d - 梯度 H - 黑塞矩阵 return: 更新的方向 """ return - np.linalg.inv(H).dot(d)def sufficientDecrease(X, y, w, step): """ 判断是否满足充沛降落条件(sufficient decrease condition) args: X - 训练数据集 y - 指标标签值 w - 权重系数 step - 步长 return: 是否满足充沛降落条件 """ d = dcost(X, y, w) H = ddcost(X, y, w) p = direction(d, H) return cost(X, y, w + step * p) <= cost(X, y, w) + c_1 * step * p.T.dot(d)def curvature(X, y, w, step): """ 判断是否满足曲率条件(curvature condition) args: X - 训练数据集 y - 指标标签值 w - 权重系数 step - 步长 return: 是否满足曲率条件 """ d = dcost(X, y, w) H = ddcost(X, y, w) p = direction(d, H) return -p.T.dot(dcost(X, y, w + step * p)) <= -c_2 * p.T.dot(d)def select(step_low, step_high): """ 在范畴内抉择一个步长 args: step_low - 步长范畴开始值 step_high - 步长范畴完结值 return: 步长 """ return (step_low + step_high) / 2def lineSearch(X, y, w, step_init, step_max): """ 线搜寻步长,使其满足 Wolfe 条件 args: X - 训练数据集 y - 指标标签值 w - 权重系数 step_init - 步长初始值 step_max - 步长最大值 return: 步长 """ step_i = step_init step_low = step_init step_high = step_max i = 1 d = dcost(X, y, w) H = ddcost(X, y, w) p = direction(d, H) while (True): # 不满足充沛降落条件或者前面的代价函数值大于前一个代价函数值 if (not sufficientDecrease(X, y, w, step_i) or (cost(X, y, w + step_i * p) >= cost(X, y, w + step_low * p) and i > 1)): # 将以后步长作为搜寻的右边界 # return search(X, y, w, step_prev, step_i) step_high = step_i else: # 满足充沛降落条件并且满足曲率条件,即曾经满足了 Wolfe 条件 if (curvature(X, y, w, step_i)): # 间接返回以后步长 return step_i step_low = step_i # 抉择下一个步长 step_i = select(step_low, step_high) i = i + 1def logisticRegressionNewton(X, y, max_iter=1000, tol=1e-4, step_init=0, step_max=10): """ 对数几率回归,应用牛顿法(newton's method) args: X - 训练数据集 y - 指标标签值 max_iter - 最大迭代次数 tol - 变动量容忍值 step_init - 步长初始值 step_max - 步长最大值 return: w - 权重系数 """ # 初始化 w 为零向量 w = np.zeros(X.shape[1]) # 开始迭代 for it in range(max_iter): # 计算梯度 d = dcost(X, y, w) # 计算黑塞矩阵 H = ddcost(X, y, w) # 当梯度足够小时,完结迭代 if np.linalg.norm(d) <= tol: break # 应用线搜寻计算步长 step = lineSearch(X, y, w, step_init, step_max) # 更新权重系数 w w = w + step * direction(d, H) return w
六、第三方库实现
scikit-learn13 实现无正则化的对数几率回归:
from sklearn.linear_model import LogisticRegression# 初始化对数几率回归器,无正则化reg = LogisticRegression(penalty="none")# 拟合线性模型reg.fit(X, y)# 权重系数w = reg.coef_# 截距b = reg.intercept_
scikit-learn13 实现L1正则化的对数几率回归:
from sklearn.linear_model import LogisticRegression# 初始化对数几率回归器,L1正则化,应用坐标降落法reg = LogisticRegression(penalty="l1", C=10, solver="liblinear")# 拟合线性模型reg.fit(X, y)# 权重系数w = reg.coef_# 截距b = reg.intercept_
scikit-learn13 实现L2正则化的对数几率回归:
from sklearn.linear_model import LogisticRegression# 初始化对数几率回归器,L2正则化reg = LogisticRegression(penalty="l2", C=10)# 拟合线性模型reg.fit(X, y)# 权重系数w = reg.coef_# 截距b = reg.intercept_
scikit-learn13 实现弹性网络正则化的对数几率回归:
from sklearn.linear_model import LogisticRegression# 初始化对数几率回归器,弹性网络正则化reg = LogisticRegression(penalty="elasticnet", C=10, l1_ratio=0.5, solver="saga")# 拟合线性模型reg.fit(X, y)# 权重系数w = reg.coef_# 截距b = reg.intercept_
scikit-learn 每种优化办法所反对的正则项:
七、动画演示
下图展现了演示数据,其中红色示意标签值为1、蓝色示意标签值为-1:
下图为应用梯度降落法拟合数据的后果,其中浅红色示意拟合后依据权重系数计算出预测值为1的局部,浅蓝色示意拟合后依据权重系数计算出预测值为-1的局部:
下图比照了不同的正则化形式对后果的影响,正则化惩办项系数 C = 0.01,弹性网络中 = 0.5。
八、思维导图
九、参考文献
- https://en.wikipedia.org/wiki...
- https://en.wikipedia.org/wiki...
- https://en.wikipedia.org/wiki...
- https://en.wikipedia.org/wiki...
- https://en.wikipedia.org/wiki...
- https://en.wikipedia.org/wiki...
- https://en.wikipedia.org/wiki...
- https://en.wikipedia.org/wiki...
- https://en.wikipedia.org/wiki...
- https://en.wikipedia.org/wiki...
- https://towardsdatascience.co...
- https://en.wikipedia.org/wiki...
- https://scikit-learn.org/stab...
残缺演示请点击这里
注:本文力求精确并通俗易懂,但因为笔者也是初学者,程度无限,如文中存在谬误或脱漏之处,恳请读者通过留言的形式批评指正
本文首发于——AI导图,欢送关注