逻辑回归实战，一文把理论和实践结合

原创

Yerik

修改于 2020-06-16 10:04:45

7070

修改于 2020-06-16 10:04:45

文章被收录于专栏：烹饪一朵云烹饪一朵云

背景

解决的二分类问题,如手写识别0-9

目标：逻辑回归返回一个概率值[0-1]

逻辑回归的特点：快、效果好、容易实时在线预测、利于分析

方法：定义一个条件概率，如p(Y|X)相当于用模型来捕获输入X和输出Y之间的关系，如

p(Y|X)=w^Tx+b

推导

对于二分类问题，由于p(Y|X)的值域在[-∞,+∞]，为了令其至于范围压缩到[-1,1]之间，故推荐使用sigmoid函数，故得

p(y=1|x,w)=\frac{1}{1+e^{-(w^Tx+b)}} \\ p(y=0|x,w)=\frac{e^{-(w^Tx+b)}}{1+e^{-(w^Tx+b)}}

两式子合并，可得

p(y|x,w)=p(y=1|x,w)^y[1-p(y=1|x,w)^{1-y}]

关于使用sigmoid

前面说道使用该函数是由于其定义域是[-∞,+∞]，且值域是[-1,1]故选择该函数，这样的说法其实并不能服众，毕竟probit也具有相同的性质。事实上选用sigmoid的意义在于，指数族分布exponential family所具有的最佳(即maximum entropy 最大熵)性质

熵原本是信息论中的概念，用在概率分布上可以表示这个分布中所包含的不确定度，熵越大不确定度越大。所以大家可以想象到，均匀分布熵最大，因为基本新数据是任何值的概率都均等。

而我们现在关心的是，给定某些假设之后，熵最大的分布。也就是说这个分布应该在满足我假设的前提下越均匀越好。比如大家熟知的正态分布，正是假设已知mean和variance后熵最大的分布。此时回过来看logistic regression，这里假设了什么呢？

在建模预测 Y|X，并认为 Y|X 服从贝努力分布(bernoulli distribution)，所以我们只需要知道 P(Y|X)
需要一个线性模型，所以 P(Y|X) = f(wx+b)。接下来我们就只需要知道 f 是什么就行了。而我们可以通过最大熵原则推出的这个 f，就是sigmoid。

bernoulli的exponential family形式，如下

\begin{aligned} p(y_i;\phi)&=\phi^y(1-\phi)^{1-y}\\ &=e^{y\log \phi+(1-y)\log (1-\phi)}\\ &=e^{y\log \frac{\phi}{1-\phi}+\log(1-\phi)} \end{aligned}

决策边界 decision boundary

所谓决策边界就是能够把样本正确分类的一条边界，在这个分类点上，大于这个边界点，分类为1，小于这个边界点，主要有线性决策边界(linear decision boundaries)和非线性决策边界(non-linear decision boundaries)。注意：决策边界是假设函数的属性，由参数决定，而不是由数据集的特征决定。

通过该指标是用来判断是不是线性分类器，LR是线性分类器，证明如下

由： p(y|x,w)=p(y=1|x,w)^y[1-p(y=1|x,w)^{1-y}]\\ 当p(y=1|x,w)=p(y=0|x,w)=0.5(即该点恰好在边界上)时，\\得\frac{p(y=1|x,w)}{p(y=0|x,w)}=1\\ 即\log \frac{p(y=1|x,w)}{p(y=0|x,w)}=\log 1\\ \Rightarrow -(w^Tx+b)=0 \Rightarrow w^Tx+b=0

注意到w^T x+b=0也代表一条直线，不妨假设X有两个特征x1,x2，那么有w1x1+w2x2+b=0

通过决策边界，我们将构造目标函数，其意义在于，最大化似然。即最大化我们见到的是样本数据概率，反过来思考，相当于我们见到的样本和数据，就是我们未来得到的模型产生的。LR的目标函数是凸函数

比如说，我们拥有数据集：

D=\{(x_i,y_i)\}_{i=1}^n\ x_i\in R^d,y_i\in \{0,1\}D={(xi,yi)}i=1n xi∈Rd,yi∈{0,1}

前面我们得到

p(y|x,w)=p(y=1|x,w)^y[1-p(y=1|x,w)^{1-y}]p(y∣x,w)=p(y=1∣x,w)y[1−p(y=1∣x,w)1−y]

可以得到最大化目标函数：

\widehat{w_{MLE}},\widehat{b_{MLE}}=argmax_{w,b}\prod_{i=1}^np(y_i|x_i,,w,b)wMLE,bMLE=argmaxw,bi=1∏np(yi∣xi,,w,b)

其计算过程如下：

\begin{aligned} \widehat{w_{MLE}},\widehat{b_{MLE}}&=argmax_{w,b}\prod_{i=1}^np(y_i|x_i,,w,b)\\ &=argmax_{w,b}\sum_{i=1}^n\log(\prod_{i=1}^np(y_i|x_i,,w,b))，其中引入log是为了解决类型精度可能越界的问题\\ &=argmin_{w,b}-\sum_{i=1}^n\log(\prod_{i=1}^np(y_i|x_i,,w,b))，找到改了最小值得w和b\\ &=argmin_{w,b}-\sum_{i=1}^n\log(p(u_i=1|x_i,w,b)^{y_i}[1-p(y_i=1|x_i,w,b)]^{1-y_i})\\ &=argmin_{w,b}-\sum_{i=1}^n\log(p(u_i=1|x_i,w,b)+(1-y_i)\log[1-p(y_i=1|x_i,w,b)]) \end{aligned}wMLE,bMLE=argmaxw,bi=1∏np(yi∣xi,,w,b)=argmaxw,bi=1∑nlog(i=1∏np(yi∣xi,,w,b))，其中引入log是为了解决类型精度可能越界的问题=argminw,b−i=1∑nlog(i=1∏np(yi∣xi,,w,b))，找到改了最小值得w和b=argminw,b−i=1∑nlog(p(ui=1∣xi,w,b)yi[1−p(yi=1∣xi,w,b)]1−yi)=argminw,b−i=1∑nlog(p(ui=1∣xi,w,b)+(1−yi)log[1−p(yi=1∣xi,w,b)])

为了高效求极值这个时候就需要使用梯度下降

梯度下降

\\θ←θ−η⋅∇J(θ)

其中 ∇J(θ) 是参数的梯度，根据计算目标函数采用数据量的不同，梯度下降算法又可以分为批量梯度下降算法Batch Gradient Descent，随机梯度下降算法Stochastic Gradient Descent和小批量梯度下降算法Mini-batch Gradient Descent。对于批量梯度下降算法，其 J(θ) 是在整个训练集上计算的，如果数据集比较大，可能会面临内存不足问题，而且其收敛速度一般比较慢。随机梯度下降算法是另外一个极端， J(θ) 是针对训练集中的一个训练样本计算的，又称为在线学习，即得到了一个样本，就可以执行一次参数更新。所以其收敛速度会快一些，但是有可能出现目标函数值震荡现象，因为高频率的参数更新导致了高方差。小批量梯度下降算法是折中方案，选取训练集中一个小批量样本（一般是2的倍数，如32,64,128等）计算，这样可以保证训练过程更稳定，而且采用批量训练方法也可以利用矩阵计算的优势。这是目前最常用的梯度下降算法。

导数与梯度

下面是计算函数梯度的一个例子：

\nabla(x^3+2xy-y^2)=(3x^2+2y,2x-2y)∇(x3+2xy−y2)=(3x2+2y,2x−2y)

其中∇称为梯度算子，它作用于一个多元函数，得到一个向量。

一维度的纯量x的梯度，通常用f'(x)表示。多维度的向量x的梯度，通常用∇f(x)表示。也就是说一维的纯量x的梯度就是算f(x)对x的微分，多维的向量x的梯度就是算f(x)对x所有元素的偏微分

例如：

x = \begin{bmatrix} x_1 \\ x_2 \\ · \\ · \\ · \\ x_d \end{bmatrix},\nabla f(x)=\begin{bmatrix} \frac{\partial f(x)}{\partial x_1} \\ \frac{\partial f(x)}{\partial x_2} \\ · \\ · \\ · \\ \frac{\partial f(x)}{\partial x_d} \end{bmatrix}

假设我们的x有两个维度的参数，梯度就分別需要对不同维度的参数做偏微分

x = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}, a= \begin{bmatrix} 20 \\ 1 \end{bmatrix}, b=\begin{bmatrix} 5 & 4 \\ 3 & 2 \end{bmatrix} \\ f_1(x)=\alpha^Tx+1=\begin{bmatrix} 20 \\ 1 \end{bmatrix}^T\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}+1=20x_1+x_2+1 \\ \nabla f_1(x)=\begin{bmatrix} \frac{\partial f_1(x)}{\partial x_1} \\ \frac{\partial f_1(x)}{\partial x_2} \end{bmatrix} = \begin{bmatrix} \frac{\partial (20x_1+x_2+1)}{\partial x_1} \\ \frac{\partial (20x_1+x_2+1)}{\partial x_2} \end{bmatrix} = \begin{bmatrix} 20 \\ 1 \end{bmatrix} \\ f_2(x)=x^Tbx+\alpha^Tx+1=\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}^T\begin{bmatrix} 5 & 4 \\ 3 & 2 \end{bmatrix}\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}+\begin{bmatrix} 20 \\ 1 \end{bmatrix}^T\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}+1 =5x_1^2+2x_2^2+(4+3)x_1x_2+20x_1+x_2+1=5x_1^2+2x_2^2+7x_1x_2+20x_1+x_2+1 \\ \nabla f_2(x)=\begin{bmatrix} \frac{\partial f_2(x)}{\partial x_1} \\ \frac{\partial f_2(x)}{\partial x_2} \end{bmatrix} = \begin{bmatrix} \frac{\partial (5x_1^2+2x_2^2+(4+3)x_1x_2+20x_1+x_2+1)}{\partial x_1} \\ \frac{\partial (5x_1^2+2x_2^2+(4+3)x_1x_2+20x_1+x_2+1)}{\partial x_2} \end{bmatrix}=\begin{bmatrix} 10x_1+7x_2+20 \\ 4x_2+7x_1+1 \end{bmatrix}

可导函数在某一点处取得极值的必要条件是梯度为0，梯度为0的点称为函数的驻点，这是疑似极值点。需要注意的是，梯度为0只是函数取极值的必要条件而不是充分条件，即梯度为0的点可能不是极值点。至于是极大值还是极小值，要看二阶导数

注意求导

在这里我们可能会问：直接求函数的导数(梯度)，然后令导数(梯度)为0，解方程，问题不就解决了吗？事实上没这么简单，因为这个方程可能很难解。比如下面的函数：

f(x,y)=x^3-2x^2+e^{xy}-y^3+10y^2+100\sin(xy)f(x,y)=x3−2x2+exy−y3+10y2+100sin(xy)

我们分别对x和y求偏导数，并令它们为0，得到下面的方程组：

\begin{cases} 3x^2-4x+ye^{xy}+100y\cos*(xy)=0 \\ xe^{xy}-3y^2+20y+x\cos(xy)=0 \end{cases}

这个方程非常难以求解，对于有指数函数，对数函数，三角函数的方程，我们称为超越方程，求解的难度并不比求极值本身小。

求解

精确的求解不太可能，因此只能求近似解，这称为数值计算。工程上实现时通常采用的是迭代法，它从一个随机初始点 X0 开始，反复使用某种规则从 Xk 移动到下一个点 X k+1 ，构造这样一个数列，直到收敛到梯度为0的点处。即有下面的极限成立：

\lim_{k\rightarrow +\infty}\nabla f(x_k)=0k→+∞lim∇f(xk)=0

这些规则一般会利用一阶导数信息即梯度；或者二阶导数信息Hessian矩阵。这样迭代法的核心是得到这样的由上一个点确定下一个点的迭代公式：

x_{k+1}=h(x_k)xk+1=h(xk)

这个过程就像我们处于山上的某一位置，要到山底找水喝，因此我们必须到达最低点处：

Hessian矩阵，一个多元函数的二阶偏导数构成的方阵，描述了函数的局部曲率。黑塞矩阵最早于19世纪由德国数学家Ludwig Otto Hesse提出，并以其名字命名。黑塞矩阵常用于牛顿法解决优化问题，利用黑塞矩阵可判定多元函数的极值问题。在工程实际问题的优化设计中，所列的目标函数往往很复杂，为了使问题简化，常常将目标函数在某点邻域展开成泰勒多项式来逼近原函数，此时函数在某点泰勒展开式的矩阵形式中会涉及到黑塞矩阵。如果Hessian矩阵正定，函数有极小值如果Hessian矩阵负定，函数有极大值如果Hessian矩阵不定，则不是极值点（鞍点）这和一元函数的结果类似，Hessian矩阵可以看做是一元函数的二阶导数对多元函数的推广。一元函数的极值判别法为，假设在某点处导数等于0，则：如果二阶导数大于0，函数有极小值如果二阶导数小于0，函数有极大值如果二阶导数等于0，情况不定

存在唯一解

利用以下函数测试梯度下降法，当然了多维度也是如此

\begin{aligned} f(x)&=x^2-10x+1 \\ f'(x)&=2x-10\\ x^{(t+1)}&=x^{(t)}- \gamma f'(x^{(t)}) \\ \Rightarrow x^{(t+1)}&=x^{(t)}-\gamma(2x^{(t)}-10)=(1-2\gamma)x^{(t)}+10\gamma \\ \Rightarrow x^{(t+1)}&=(1-2\gamma)x^{(t)}+10\gamma \end{aligned}

这个例子基本上学习率可以不用太小，就可以很快就找到解，后面有跑不同学习率看几次可以得出近似解。

切线公式：y=f'(\alpha)(x-\alpha)+f(\alpha) \\ 法线公式: y=\frac{-1}{f'(\alpha)}(x-\alpha)+f(\alpha)

刚才提到我们需要先设定一个初始化的“解”，此例不妨设x(0)=20(故意跟最佳值有差距)

注意：红色线是法线，蓝色线是切线，法线和切线这两条线是垂直的，但因为x轴和y轴刻度不一样，所以看不出来它是垂直的。

推导

我们来看一元函数的泰勒展开，以便于更好的理解多元函数的泰勒展开。假设一个一元函数n阶可导，它的泰勒展开式为：

f(x+\Delta x)=f(x)+f'(x)\Delta x+\frac{1}{2}f'(x)(\Delta x)^2+...+\frac{1}{n!}f^{n}(x)(\Delta x)^{n}

如果在某一点处导数值大于0(+)，则函数在此处是增函数，加大x的值函数值会增加，减小x的值(-)函数会减小。相反的，如果在某一点处导数值小于0(-)，则函数是减函数，增加x的值函数值会减小(+)，减小x的值函数会增加。因此我们可以得出一个结论：如果x的变化很小，并且变化值与导数值反号，则函数值下降。对于一元函数，x的变化只有两个方向，要么朝左，要么朝右。

下面我们把这一结论推广到多元函数的情况。多元函数 f(x) 在x点处的泰勒展开为：

f(x+\Delta x)-f(x)=(\nabla f(x))^T\Delta x+\omicron(\Delta x)f(x+Δx)−f(x)=(∇f(x))TΔx+ο(Δx)

这里我们忽略了二次及更高的项。其中，一次项是梯度向量∇f(x)与自变量增量Δx的内积(∇f(x))^T Δx，这等价于一元函数的 f'(x0)(x-x0) 。这样，函数的增量与自变量的增量Δx 、函数梯度的关系可以表示为：

f(x+\Delta x)-f(x)=(\nabla f(x))^T \Delta x+\omicron(\Delta x)f(x+Δx)−f(x)=(∇f(x))TΔx+ο(Δx)

如果Δx足够小，在x的某一邻域内，则我们可以忽略二次及以上的项，则有：

f(x+\Delta x)-f(x) \approx(\nabla f(x))^T \Delta xf(x+Δx)−f(x)≈(∇f(x))TΔx

这里的情况比一元函数复杂多了，是一个向量，Δx有无穷多种方向，如果能保证：

(\nabla f(x))^T \Delta x < 0

则有：

f(x + \Delta x)< f(x)

令函数值递减，这就是下山的正确方向。因为有：

(\nabla f(x))^T \Delta x=\|\nabla f(x)\| \|\Delta x\|\cos \theta(∇f(x))TΔx=∥∇f(x)∥∥Δx∥cosθ

此时，||·||表示向量的模， θ 是向量 ∇f(x) 和 Δx 的夹角。因为向量的模一定大于等于0，如果

\cos \theta \geqslant -1cosθ⩾−1

只有当：

\theta = \pi

时 cos θ 有极小值-1，此时梯度和 Δx 反向，即夹角为180度。因此当向量 Δx 的模大小一定时，当：

\Delta x=- \alpha \nabla f(x)Δx=−α∇f(x)

即在梯度相反的方向函数值下降的最快。此时有：

\cos \theta = -1

函数的下降值为：

只要梯度不为0，往梯度的反方向走函数值一定是下降的。直接用 Δx=-∇f(x) 可能会有问题，因为 x+ Δx 可能会超出x的邻域范围之外，此时是不能忽略泰勒展开中的二次及以上的项的，因此步伐不能太大。一般设：

f(x)Δx=−α∇f(x)

其中 α 为一个接近于0的正数，称为步长，由人工设定，用于保证 x+ Δx 在x的邻域内，从而可以忽略泰勒展开中二次及更高的项，则有：

(\nabla f(x))^T \Delta x=-\alpha(\nabla f(x))^T(\nabla f(x)) \leqslant 0(∇f(x))TΔx=−α(∇f(x))T(∇f(x))⩽0

从初始点 X0 开始，使用如下迭代公式：

X_{k+1}=x_k-\alpha \nabla f(x_k)

只要没有到达梯度为0的点，则函数值会沿着序列 Xk 递减，最终会收敛到梯度为0的点，这就是梯度下降法。迭代终止的条件是函数的梯度值为0（实际实现时是接近于0），此时认为已经达到极值点。注意我们找到的是梯度为0的点，这不一定就是极值点，后面会说明。梯度下降法只需要计算函数在某些点处的梯度，实现简单，计算量小。

实现细节问题

在训练的过程中不免会遇到一些细节的问题，最基本的就是调参了

初始值的设定

一般的，对于不带约束条件的优化问题，我们可以将初始值设置为0，或者设置为随机数，对于神经网络的训练，一般设置为随机数，这对算法的收敛至关重要

学习率的设定

学习率设置为多少，也是实现时需要考虑的问题。最简单的，我们可以将学习率设置为一个很小的正数，如0.001。另外，可以采用更复杂的策略，在迭代的过程中动态的调整学习率的值。比如前1万次迭代为0.001，接下来1万次迭代时设置为0.0001。

面临的问题

这个是其实是要基于业务来考虑的，因为在现实层面上，局部极值不可避免

局部极小值

有些函数可能有多个局部极小值点，下面是一个例子：

存在局部最小值

设计一个有局部极小值和全域极小值的函数，到四次方，如令x=10，函数的值就很大

\begin{aligned} f(x)&=x^4-50x^3-x+1 \\ f'(x)&=4x^3-150x^2-1 \\ x^{(t+1)}&=x^{(t)}- \gamma f'(x^{(t)}) \\ \Rightarrow x^{(t+1)}&=x^{(t)}-\gamma(4x^{(t)^3}-150x^{(t)^2}-1) \\ \Rightarrow x^{(t+1)}&=-4\gamma x^{(t)^3}+150\gamma x^{(t)^2}+x^{(t)}+\gamma \end{aligned}

我们需要先设定一个初始化的“解”，此例我设置x(0)=-20(故意跟最佳值有差距)

可以看到学习率太小，初始值不好，解就会掉到局部极小值。

当学习率为0.0004对此例子来说，虽然步伐够大跳出了局部极值，但到全域极值时，因为步伐太大，所以走不到最好的值。

这个学习率(0.0003)对此例子来说就够了，可以走到全域极值。

鞍点

鞍点是指梯度为0，Hessian矩阵既不是正定也不是负定，即不定的点。下面是鞍点的一个例子，假设有函数：

x^2-y^2x2−y2

显然在(0, 0)这点处不是极值点，但梯度为0，下面是梯度下降法的运行结果：

在这里，梯度下降法遇到了鞍点，认为已经找到了极值点，从而终止迭代过程，而这根本不是极值点。

对于怎么逃离局部极小值点和鞍点，有一些解决方案，在这里我们暂时不细讲，以后有机会再专门写文章介绍。对于凸优化问题，不会遇到上面的局部极小值与鞍点问题，即梯度下降法一定能找到全局最优解。凸优化的概念将在SIGAI后续的文章中介绍。

没有极值

虽然说微分可以找极值，但很多函数既无最大值，也无最小值，因为函数的长相弯弯曲曲很多次，有局部极值或鞍部，所以一次微分等于0求得的可能是极值，也可以是相对极值。比如这个例子

f(x)=5x_1^2+2x_2^2+7x_1x_2+10x_1+x_2+1\\ \nabla f(x)=\begin{bmatrix} 10x_1+7x_2+10 \\ 4x_2+7x_1+1 \end{bmatrix}\\ \nabla f(x)==\begin{bmatrix} 10x_1+7x_2+10 \\ 4x_2+7x_1+1 \end{bmatrix}=\begin{bmatrix} 0 \\ 0 \end{bmatrix}解联立方程得\\ \Rightarrow \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}=\begin{bmatrix} \frac{11}{3} \\ \frac{-20}{3} \end{bmatrix}

这个方程式可以找到极值“解”让f(x)最小即f(x)=16，但这个值真的是最小吗? 不妨找个点随便代入，如

\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}=\begin{bmatrix} -20 \\ 20 \end{bmatrix}，f(\begin{bmatrix} x_1 \\ x_2 \end{bmatrix})=-179[x1x2]=[−2020]，f([x1x2])=−179

这个值比微分的最佳解还要小，所以可以得知微分等于0找到的不一定是最佳解，所以用梯度下降法，可以找到更好的解。

下图我将上式子画出來它的坐标跟微分解还有梯度法如何让解更新。

这边我只跑100次，因为解在无穷大的地方，但可以看到loss值不断在变小中。

总结步骤

设p(y=1|x,w)=\frac{1}{e^{-w^Tx_i+n}}=\sigma(w^Tx+b)\\ 由argmin_{w,b}-\sum_{i=1}^n\log(p(u_i=1|x_i,w,b)+(1-y_i)\log[1-p(y_i=1|x_i,w,b)])\\ =argmin_{w,b}-\sum_{i=1}^n\log(\sigma(w^Tx_i+b)+(1-y_i)\log[1-\sigma(w^Tx_i+b)])\\ 得\partial\frac{L(w,b)}{w}=-\sum_{i=1}^ny_i\frac{\sigma(w^Tx_i+b[1-\sigma(w^Tx_i+b)]}{\sigma(w^Tx_i+b)}x_i+(1-y_i)\frac{-\sigma(w^Tx_i+b)[1-\sigma(w^Tx_i+b)]}{1-\sigma(w^Tx_i+b)}x_i=\sum_{i=1}^n[\sigma(w^Tx_i+b)-y_i]x_i\\ \partial\frac{L(w,b)}{b}=-\sum_{i=1}^ny_i\frac{\sigma(w^Tx_i+b[1-\sigma(w^Tx_i+b)]}{\sigma(w^Tx_i+b)}+(1-y_i)\frac{-\sigma(w^Tx_i+b)[1-\sigma(w^Tx_i+b)]}{1-\sigma(w^Tx_i+b)}=\sum_{i=1}^n[\sigma(w^Tx_i+b)-y_i]\\ 其中\sigma(w^Tx_i+b)是预测值，y_i是真实值设p(y=1∣x,w)=e−wTxi+n1=σ(wTx+b)由argminw,b−i=1∑nlog(p(ui=1∣xi,w,b)+(1−yi)log[1−p(yi=1∣xi,w,b)])=argminw,b−i=1∑nlog(σ(wTxi+b)+(1−yi)log[1−σ(wTxi+b)])得∂wL(w,b)=−i=1∑nyiσ(wTxi+b)σ(wTxi+b[1−σ(wTxi+b)]xi+(1−yi)1−σ(wTxi+b)−σ(wTxi+b)[1−σ(wTxi+b)]xi=i=1∑n[σ(wTxi+b)−yi]xi∂bL(w,b)=−i=1∑nyiσ(wTxi+b)σ(wTxi+b[1−σ(wTxi+b)]+(1−yi)1−σ(wTxi+b)−σ(wTxi+b)[1−σ(wTxi+b)]=i=1∑n[σ(wTxi+b)−yi]其中σ(wTxi+b)是预测值，yi是真实值

步骤归纳

最小化F(w)

设置初始w，计算F(w)
计算梯度∇F(w) 下降方向：dir = -∇F(w)
尝试梯度更新 w new = w + 步长 ∗ dir 得到下降后的 w new和F(w new)
如果 F(w new)-F(W)较小，停止；否则 w=w; F(w)=F(w new)跳到第2步

随机梯度下降法

对于有些机器学习问题，我们的目标函数是对样本的损失函数。假设训练样本集有N个样本，训练时优化的目标是这个数据集上的平均损失函数：

L(w) = \frac{1}{N}\sum_{i=1}^N L(w,x_i,y_i)

其中L(w,xi,yi) 是对单个训练样本 (xi,yi) 的损失函数，w是需要学习的参数。如果训练时每次都用所有样本计算梯度并更新，成本太高，作为改进可以在每次迭代时选取一批样本，将损失函数定义在这些样本上。

批量随机梯度下降法在每次迭代中使用上面目标函数的随机逼近值，即只使用 M<-N 个样本来近似计算损失函数。在每次迭代时要优化的目标函数变为：

L(w) \approx \frac{1}{M}\sum_{i=1}^M L(w,x_i,y_i)

已经证明，随机梯度下降法在数学期望的意义下收敛，即随机采样产生的梯度的期望值是真实的梯度。因为每次迭代时的目标函数实际上是不一样的，因此随机梯度下降法并不能保证每次迭代时函数值一定下降。

每次从训练样本中抽取一个样本进行更新，每次都不用遍历所有数据集，迭代速度快，需要迭代更多次数
每次选取方向不一定是最优

初始化w',b'\\ for\ t=1,2...\\ w^{t+1}=w^t-\eta\sum_{i=1}^n[\sigma(w^Tx_i+n)-y_i]x_i\\ b^{t+1}=b^t-\eta\sum_{i=1}^n[\sigma(w^Tx_i+n)-y_i]

批量梯度下降和随机梯度下降，折中的方法：Mini-batch Gradent Descen

正则化项

从上文可知，

p(y=1|x,w)=\frac{1}{1+e^{-(w^Tx+b)}} \\ p(y=0|x,w)=\frac{e^{-(w^Tx+b)}}{1+e^{-(w^Tx+b)}}

当w变得非常大的时候，我们会发现p(y=1|x,w)≈1，p(y=0|x,w)≈0，这是非常理想的情况，而此时就出现了所谓过拟合现象

事实上在拟合中会出现过拟合和欠拟合两种现象

欠拟合
- 模型没有很好的捕捉到数据特征
- 学习不充分，导致没有很好拟合
- 方案
  - 模型复杂化
  - 增加训练轮数
过拟合
- 过度训练，加强了对噪声学习
- 过度拟合样本，泛化性差
- 方案
  - 交叉验证
  - 增加数据
  - 特征筛选
  - 正则化

带着问题出发：

正则项存在的意义是什么，为什么要使用正则项？正则项是如何防止过拟合的？
有哪几种正则项，如何表示，它们的相同点和不同点是什么？
不同正则项的使用场景是什么，如何选取正则项呢？

为什么要使用正则化项

其实正则化项是对参数的控制。其目的主要有两个

实现参数的稀疏，这样可以简化模型，避免过拟合。在一个模型中重要的特征并不是很多，如果考虑所有的特征都是有作用的，那么就会对训练集进行充分的拟合，导致在测试集的表现并不是很好，所以我们需要稀疏参数，简化模型。
尽可能保证参数小一些，这又是为啥呢？因为越是复杂的模型，它会对所有的样本点进行拟合，如果在这里包含异常的样本，就会在小区间内产生很大的波动，不同于平均水平的高点或者低点，这样的话，会导致其导数很大，我们知道在多项式导数中，只有参数非常大的时候，才会产生较大的导数，所以模型越复杂，参数值也就越大。为了避免这种过度的拟合，需要控制参数值的大小。

总而言之，避免w过大，L1正则化和L2正则化上可以看做是损失函数的惩罚

正则项的分类

正则项有三种：L0、L1、L2

L0正则化的值是模型参数中非零参数的个数。
L1正则化表示各个参数绝对值之和。
L2正则化标识各个参数的平方的和的开方值。

L0正则化

保证参数稀疏化来防止过拟合，可以用非零参数，来进行特征选择。但是L0正则化不好求，因此采用L1正则化。L1正则化是L0正则化的最优凸近似，比L0容易求解，并且可以实现稀疏的效果。

L1正则化

L1正则化可以产生稀疏权值矩阵，即产生一个稀疏模型，可以用于特征选择

L1正则化也叫lasso，它往往是替代L0正则化来防止过拟合的。为啥用L1范数，因为L1范数就是各个参数的绝对值相加，我们已知，参数的值的大小和模型的复杂度是成正比的，因此复杂模型，L1范数就会大，导致损失函数大。下面定量的分析：

在原始的代价函数后面加上一个L1正则化项，即所有权重w的绝对值的和，乘以λ/n。如下：

C = C_0 + \frac{\lambda}{n}\sum_w|w|

计算导数得：

\frac{\partial C}{\partial w}=\frac{\partial C_0}{\partial w}+ \frac{\lambda}{n}sgn(w)

上式中sgn(w)表示w的符号。那么权重w的更新规则为：

w \rightarrow w'=w-\frac{\eta \lambda}{n}sgn(w)-\eta\frac{\partial C_0}{\partial w}

现在来观察正则求导项，可知当w为正时，更新后的w变小；当w为负时，更新后的w变大。因此它的效果就是让w往0靠，使网络中的权重尽可能为0，也就相当于减小了网络复杂度，防止过拟合。另外，上面没有提到一个问题，当w为0时怎么办？当w等于0时，|w|是不可导的，所以我们只能按照原始的未经正则化的方法去更新w，这就相当于去掉ηλsgn(w)/n这一项，所以我们可以规定sgn(0)=0，这样就把w=0的情况也统一进来了。

注意到L1正则化是权值的绝对值之和，J是带有绝对值符号的函数，因此J是不完全可微的。机器学习的任务就是要通过一些方法（比如梯度下降）求出损失函数的最小值。考虑二维的情况，即只有两个权值w1和w2，此时L=|w1|+|w2|，对于梯度下降法，求解J的过程可以画出等值线，同时L1正则化的函数L也可以在w1w2的二维平面上画出来。如下图：

在图中，当J等值线与L首次相交的地方就是最优解。上图中J与L在L的一个顶点处相交，这个顶点就是最优解。注意到这个顶点的值是(w1,w2)=(0,w)。可以直观想象，因为L函数有很多突出的角（二维情况下四个，多维情况下更多），J与这些角接触的机率会远大于与L其它部位接触的机率，而在这些角上，会有很多权值等于0，这就是为什么L1正则化可以产生稀疏模型，进而可以用于特征选择。

L2正则化

L2正则化可以防止模型过拟合（overfitting）；一定程度上，L1也可以防止过拟合

L2正则化也是防止过拟合的，原因和L1一样一样的，就是形式不同。L2范数是各参数的平方和再求平方根。对于L2的每个元素都很小，但是不会为0，只是接近0，参数越小说明模型越简单，也就越不容易产生过拟合。L2正则化也叫做“岭回归”。

来让我们看看具体的例子，对于房屋价格预测我们可能有上百种特征，与刚刚所讲的多项式例子不同，我们并不知道哪些是高阶多项式的项。所以，如果我们有一百个特征，我们并不知道如何选择关联度更好的参数，如何缩小参数的数目等等。因此在正则化里，我们要做的事情，就是把减小我们的代价函数（例子中是线性回归的代价函数）所有的参数值，因为我们并不知道是哪一个或哪几个要去缩小。因此，我们需要修改代价函数，在这后面添加一项，就像我们在方括号里的这项。当我们添加一个额外的正则化项的时候，我们收缩了每个参数。

现在进行定量的分析：

L2正则化就是在代价函数后面再加上一个正则化项：

C = C_0 + \frac{\lambda}{2n}\sum_w w^2

C0代表原始的代价函数，后面那一项就是L2正则化项，它是这样来的：所有参数w的平方的和，除以训练集的样本大小n。λ就是正则项系数，权衡正则项与C0项的比重。另外还有一个系数1/2，1/2经常会看到，主要是为了后面求导的结果方便，后面那一项求导会产生一个2，与1/2相乘刚好凑整。L2正则化项是怎么避免overfitting的呢？我们推导一下看看，先求导：

\begin{aligned} \frac{\partial C}{\partial w}&=\frac{\partial C_0}{\partial w}+\frac{\lambda}{n}w\\ \frac{\partial C}{\partial b}&=\frac{\partial C_0}{\partial b} \end{aligned}

可以发现L2正则化项对b的更新没有影响，但是对于w的更新有影响:

w \rightarrow w-\eta\frac{\partial C_0}{\partial w}-\frac{\eta \lambda}{n}w=(1-\frac{\eta \lambda}{n})w-\eta\frac{\partial C_0}{\partial w}

在不使用L2正则化时，求导结果中w前系数为1，现在w前面系数为 1-ηλ/n ，因为η、λ、n都是正的，在样本量充足的时候，1-ηλ/n小于1，它的效果是减小w，这也就是权重衰减的由来。当然考虑到后面的导数项，w最终的值可能增大也可能减小。

同理，假设有如下带L2正则化的损失函数，同样可以画出他们在二维平面上的图形，如下：

L2正则化二维平面下L2正则化的函数图形是个圆，与方形相比，被磨去了棱角。因此J与L相交时使得w1或w2等于零的机率小了许多，这就是为什么L2正则化不具有稀疏性的原因。

总结

L1会趋向于产生少量的特征，而其他的特征都是0，而L2会选择更多的特征，这些特征都会接近于0。Lasso在特征选择时候非常有用，而Ridge就只是一种规则化而已。在所有特征中只有少数特征起重要作用的情况下，选择Lasso比较合适，因为它能自动选择特征。而如果所有特征中，大部分特征都能起作用，而且起的作用很平均，那么使用Ridge也许更合适。

理论联系实际

接下来我们将利用一些私有的脱敏数据进行试验，逻辑回归以及正则化项

逻辑回归

加载数据

def loaddata(file, delimeter):
    data = np.loadtxt(file, delimiter=delimeter)
    print('Dimensions: ',data.shape)
    print(data[:5,:])
    return(data)


data = loaddata('./data/sample_1.txt', ',')

整理训练样本

np.c_是按行连接两个矩阵，就是把两矩阵左右相加，要求行数相等，这么做的目的在于可以将b扩展进矩阵中

X = np.c_[np.ones((data.shape[0],1)), data[:,0:2]]
y = np.c_[data[:,2]]

而允许这么扩展的原理在于

w^Tx+b=\begin{bmatrix} w_1 \\ w_2 \\ w_3 \end{bmatrix}*\begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}+b=w_1*x_1+w_2*x_2+w_3*x_3+b\\ \Rightarrow \begin{bmatrix} w_1 \\ w_2 \\ w_3 \\ b \end{bmatrix}*\begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ 1 \end{bmatrix} =W^Tx

可视化数据

def plotData(data, label_x, label_y, label_pos, label_neg, axes=None):
    # 获得正负样本的下标(即哪些是正样本，哪些是负样本)
    neg = data[:,2] == 0
    pos = data[:,2] == 1
    
    if axes == None:
        axes = plt.gca() # 获取图的坐标信息
        
    axes.scatter(data[pos][:,0], data[pos][:,1], marker='*', c='r', s=30, linewidth=2, label=label_pos)
    axes.scatter(data[neg][:,0], data[neg][:,1], c='c', s=30, label=label_neg)
    axes.set_xlabel(label_x)
    axes.set_ylabel(label_y)
    axes.legend(frameon= True, fancybox = True);
    
plotData(data, 'English', 'Math', 'Pass', 'Fail')

接下来则是开始训练数据，首先是定义sigmoid函数

def sigmoid(z):
    return(1 / (1 + np.exp(-z)))

接下来定义损失函数，也就是误差函数,即前文我们提到的

argmin_{w,b}-\sum_{i=1}^n\log(\sigma(w^Tx_i+b)+(1-y_i)\log[1-\sigma(w^Tx_i+b)])

def lossFunction(theta, X, y):
    m = y.size
    h = sigmoid(X.dot(theta))
    J = -1.0*(1.0/m)*(np.log(h).T.dot(y)+np.log(1-h).T.dot(1-y))
               
    if np.isnan(J[0]):
        return(np.inf)
    return J[0]

接着定义梯度下降函数

def gradient(theta, X, y):
    m = y.size
    h = sigmoid(X.dot(theta.reshape(-1,1)))
    
    grad =(1.0/m)*X.T.dot(h-y)

    return(grad.flatten())

即使用X的特征和权重相乘而得到预测值h，之后预测值减去真实值再乘x，也就是前文所说推导的公式

\partial\frac{L(w,b)}{w}=\sum_{i=1}^n[\sigma(w^Tx_i+b)-y_i]x_i=grad

此时，我们已然准备好了相关的函数，现在需要初始化函数因为不知道最优参数是啥，故直接用np.zeros初始化为0

initial_theta = np.zeros(X.shape[1])
loss = lossFunction(initial_theta, X, y)
grad = gradient(initial_theta, X, y)
print('Loss: \n', loss)
print('Grad: \n', grad)

得到结果，这个相当于做了一轮的迭代

为了执行多轮迭代，我们需要使用minimize函数

res = minimize(lossFunction, initial_theta, args=(X,y), jac=gradient, options={'maxiter':400})
print(res)

得到结果

其中x: array([-25.16131634, 0.2062316 , 0.20147143])第一个参数是他的偏置，后面两个参数则是对应X的第一列和第二列的特征

最后对结果进行预测

def predict(theta, X, threshold=0.5):
    p = sigmoid(X.dot(theta.T)) >= threshold
    return(p.astype('int'))

sigmoid(np.array([1, 81, 57]).dot(res.x.T))

得到一个值，0.9537929840904646可以认为接近于1，预测是准确的，之后可视化处理

plt.scatter(81, 57, s=60, c='y', marker='v', label='(81, 57)')
plotData(data, 'English', 'Math', 'Pass', 'Fail')
x1_min, x1_max = X[:,1].min(), X[:,1].max(),
x2_min, x2_max = X[:,2].min(), X[:,2].max(),
xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max), np.linspace(x2_min, x2_max))

print(xx1.shape)

h = sigmoid(np.c_[np.ones((xx1.ravel().shape[0],1)), xx1.ravel(), xx2.ravel()].dot(res.x.T))
h = h.reshape(xx1.shape)
plt.contour(xx1, xx2, h, 1, linewidths=1, colors='b');

minimize函数

scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)

解释

fun: 求最小值的目标函数
x0:变量的初始猜测值，如果有多个变量，需要给每个变量一个初始猜测值。minimize是局部最优的解法.
args:常数值，fun中没有数字，都以变量的形式表示，对于常数项，需要在这里给值
method:求极值的方法，官方文档给了很多种。一般使用默认。该参数代表采用的方式，默认是BFGS, L-BFGS-B, SLSQP中的一种，可选TNC
jac：该参数就是计算梯度的函数，和fun参数类似，第一个必须为theta且其shape必须为(n,)即一维数组,最后返回的梯度也必须为一个一维数组
constraints:约束条件，针对fun中为参数的部分进行约束限制
options可以设置最大迭代轮数，以字典的形式来进行设置，例如：options={'maxiter':400}

正则化项

准备数据

data = loaddata('./data/sample_2.txt', ',')

y = np.c_[data[:,2]]
X = data[:,0:2]

plotData(data, 'x', 'y', 'y = 1', 'y = 0')

我们的目标就是将这批数据进行分类，比如画一个圈，圈内和圈外泾渭分明？

poly = PolynomialFeatures(10) # 多项式特征，最高10维，例如[a,b]的多项式交互式输出[1,a,b,ab]
XX = poly.fit_transform(data[:,0:2])
# 看看形状(特征映射后x有多少维了)
print(XX.shape)

定义损失函数

def lossFunctionReg(theta, reg, *args):
    m = y.size
    h = sigmoid(XX.dot(theta))
    
    J = -1.0*(1.0/m)*(np.log(h).T.dot(y)+np.log(1-h).T.dot(1-y)) + (reg/(2.0*m))*np.sum(np.square(theta[1:]))
    #  (reg/(2.0*m))*np.sum(np.square(theta[1:]))  即L2正则，开平方加和
    
    if np.isnan(J[0]):
        return(np.inf)
    return(J[0])

定义gd

def gradientReg(theta, reg, *args):
    m = y.size
    h = sigmoid(XX.dot(theta.reshape(-1,1)))
      
    grad = (1.0/m)*XX.T.dot(h-y) + (reg/m)*np.r_[[[0]],theta[1:].reshape(-1,1)]
    # (reg/m)*np.r_[[[0]],theta[1:].reshape(-1,1)] 这个是L2正则的求导
        
    return(grad.flatten())

初始化参数

init_theta = np.zeros(XX.shape[1])
# 初始化权重
lossFunctionReg(init_theta, 1, XX, y)

进行正则化项

fig, axes = plt.subplots(1,4, sharey = True, figsize=(17,5))

for i, C in enumerate([0.0, 1.0, 10.0, 100.0]):
    # 最优化 costFunctionReg
    res2 = minimize(lossFunctionReg, init_theta, args=(C, XX, y), jac=gradientReg, options={'maxiter':10000})
    
    # 准确率
    accuracy = 100.0*sum(predict(res2.x, XX) == y.ravel())/y.size    

    # 对X,y的散列绘图
    plotData(data, 'x', 'y', 'y = 1', 'y = 0', axes.flatten()[i])
    
    # 画出决策边界
    x1_min, x1_max = X[:,0].min(), X[:,0].max(),
    x2_min, x2_max = X[:,1].min(), X[:,1].max(),
    xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max), np.linspace(x2_min, x2_max))
    h = sigmoid(poly.fit_transform(np.c_[xx1.ravel(), xx2.ravel()]).dot(res2.x))
    h = h.reshape(xx1.shape)
    axes.flatten()[i].contour(xx1, xx2, h, [0.5], linewidths=1, colors='g');       
    axes.flatten()[i].set_title('Train acc {}% with lambda = {}'.format(np.round(accuracy, decimals=2), C))

运行结果如图所示

决策边界，我们不妨分别来看看正则化系数lambda太大太小分别会出现什么情况：

Lambda = 0 : 就是没有正则化，这样的话，就过拟合咯
Lambda = 1||10 : 这才是正确的打开方式
Lambda = 100 : 正则化项太激进，导致基本就没拟合出决策边界

PolynomialFeatures

可以理解为专门生成多项式特征，并且多项式包含的是相互影响的特征集，从低维空间往高维空间进行映射，比如：一个输入样本是２维的。形式如[a,b] ,则二阶多项式的特征集如下[1,a,b,a2,ab,b2]。通过这种形式，可以轻松的将x扩展为X向量

sklearn.preprocessing.PolynomialFeatures(degree=2, *, interaction_only=False, include_bias=True, order='C')

参数：

degree: integer，多项式阶数，默认为2；
interaction_only : boolean, default = False，如果值为true(默认是false),则会产生相互影响的特征集；
include_bias : boolean，是否包含偏差列。

参考内容

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

正则表达式

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

正则表达式