一、前言
本这会是激动人心的一章，因为我们将首次接触到最优化算法。仔细想想就会发现，其实我们日常生活中遇到过很多最优化问题，比如如何在最短时间内从A点到达B点？如何投入最少工作量却获得最大的效益？如何设计发动机使得油耗最少而功率最大？可见，最优化的作用十分强大。接下来，我们介绍几个最优化算法，并利用它们训练出一个非线性函数用于分类。读者不熟悉回归也没关系，
假设现在有一些数据点，我们用一条直线对这些点进行拟合（该线称为最佳拟合直线)，这个拟合过程就称作回归。利用Logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。这里的“回归”一词源于最佳拟合，表示要找到最佳拟合参数集。训练分类器时的做法就是寻找最佳拟合参数，使用的是最优化算法。接下来介绍这个二值型输出分类器的数学原理。
二、Logistic回归与梯度上升算法
Logistic回归是众多分类算法中的一员。通常，Logistic回归用于二分类问题，当然它也可以用于多分类问题，让我们来了解一下，什么是Logistic回归。
1、Logistic回归
假设现在有一些数据点，我们利用一条直线对这些点进行拟合(该线称为最佳拟合直线)，这个拟合过程就称作为回归，如下图所示：
Logistic回归是分类方法，它利用的是Sigmoid函数阈值在[0,1]这个特性。Logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。其实，Logistic本质上是一个基于条件概率的判别模型(Discriminative Model)。
所以要想了解Logistic回归，我们必须先看一看Sigmoid函数 ，我们也可以称它为Logistic函数。它的公式如下：
下面这张图片，为我们展示了Sigmoid函数的样子。
z是一个矩阵，θ是参数列向量(要求解的)，x是样本列向量(给定的数据集)。θ^T表示θ的转置。g(z)函数实现了任意实数到[0,1]的映射，这样我们的数据集([x0,x1,...,xn])，不管是大于1或者小于0，都可以映射到[0,1]区间进行分类。hθ(x)给出了输出为1的概率。比如当hθ(x)=0.7，那么说明有70%的概率输出为1。输出为0的概率是输出为1的补集，也就是30%。
如果我们有合适的参数列向量θ([θ0,θ1,...θn]^T)，以及样本列向量x([x0,x1,...,xn])，那么我们对样本x分类就可以通过上述公式计算出一个概率，如果这个概率大于0.5，我们就可以说样本是正样本，否则样本是负样本。
接下来就要介绍逻辑回归代价函数了，什么是代价函数，代价函数可以理解为是回归直线上的值与所有真实点之间距离平方的平均差值，代价函数越小，拟合直线拟合的越成功，所以目标变为求代价函数最小值。下面介绍代价函数。
根据sigmoid函数的特性，我们可以做出如下的假设：
式即为在已知样本x和参数θ的情况下，样本x属性正样本(y=1)和负样本(y=0)的条件概率。理想状态下，根据上述公式，求出各个点的概率均为1，也就是完全分类都正确。但是考虑到实际情况，样本点的概率越接近于1，其分类效果越好。比如一个样本属于正样本的概率为0.51，那么我们就可以说明这个样本属于正样本。另一个样本属于正样本的概率为0.99，那么我们也可以说明这个样本属于正样本。但是显然，第二个样本概率更高，更具说服力。我们可以把上述两个概率公式合二为一：
合并出来的Loss，我们称之为损失函数(Loss Function)。当y等于1时，(1-y)项(第二项)为0；当y等于0时，y项(第一项)为0。为s了简化问题，我们对整个表达式求对数，(将指数问题对数化是处理数学问题常见的方法)：
这个损失函数，是对于一个样本而言的。给定一个样本，我们就可以通过这个损失函数求出，样本所属类别的概率，而这个概率越大越好，所以也就是求解这个损失函数的最大值。既然概率出来了，那么最大似然估计也该出场了。假定样本与样本之间相互独立，那么整个样本集生成的概率即为所有样本生成概率的乘积，便可得到如下公式：
其中，m为样本的总数，y(i)表示第i个样本的类别，x(i)表示第i个样本，需要注意的是θ是多维向量，x(i)也是多维向量。
综上所述，满足代价函数J(θ)的最大的θ值即是我们需要求解的模型。但是我们一般在公式加上 -1/m,所以变为梯度下降求最小值。
怎么求解使J(θ)最小的θ值呢？因为是求最小值，所以我们需要使用梯度下降算法。当然如果使J(θ) := -J(θ)，那么问题就从求极小值转换成求极大值了，梯度上升算法与梯度下降算法本质一样的，它们的思想都是相同的，学会其一，就也会了另一个。
2、梯度下降算法
说了半天，梯度上升算法又是啥？J(θ)太复杂，我们先看个简单的求极大值的例子。一个看了就会想到高中生活的函数：
举个栗子，存在以下函数
求极值，先求函数的导数：
令导数为0，可求出x=2即取得函数f(x)的极小值。极小值等于f(2)=-4
此时我们就可以用迭代的方法来做。就像爬坡一样，一点一点逼近极值。这种寻找最佳拟合参数的方法，就是最优化算法。爬坡这个动作用数学公式表达即为：
其中，α为步长，也就是学习速率，控制更新的幅度。也就是下降时步幅。比如从(0,0)开始，迭代路径就是1->2->3->4->...->n，直到求出的x为函数极大值的近似值，停止迭代。我们可以编写代码，来实现这一过程：
def Gradient_Ascent_test():
def f_prime(x_old): #f(x)的导数
return 2 * x_old - 4
x_old = -1 #初始值，给一个小于x_new的值
x_new = 0 #梯度上升算法初始值，即从(0,0)开始
alpha = 0.01 #步长，也就是学习速率，控制更新的幅度
learningrate = 0.00000001 #精度，也就是更新阈值
while abs(x_new - x_old) > learningrate:
x_old = x_new
x_new = x_old - alpha * f_prime(x_old) #上面提到的公式
print(x_new) #打印最终求解的极值近似值
if __name__ == '__main__':
Gradient_Ascent_test()
结果很显然，已经非常接近我们的真实极值2了。这一过程，就是梯度下降算法。那么同理，J(θ)这个函数的极值，也可以这么求解。公式可以这么写：
由上小节可知J(θ)为：
那么，现在我只要求出J(θ)的偏导，就可以利用梯度下降算法，求解J(θ)的极小值了。
那么现在开始求解J(θ)对θ的偏导，求解如下：
其中：
再由：
可得：
接下来，就剩下第三部分：
综上所述：
因此，梯度上升迭代公式为：
但是我们习惯加负号。于是公式为：
你发现其实这两个公式一样的，就是提个负号就行，知道了，梯度下降迭代公式，我们就可以自己编写代码，计算最佳拟合参数了。
三、Python3实战
1、数据准备
数据集已经为大家准备好，下载地址：数据集下载
先看下数据集有哪些数据：
这个数据有两维特征，因此可以将数据在一个二维平面上展示出来。我们可以将第一列数据(X1)看作x轴上的值，第二列数据(X2)看作y轴上的值。而最后一列数据即为分类标签。根据标签的不同，对这些点进行分类。
那么，先让我们编写代码，看下数据集的分布情况：
从上图可以看出数据的分布情况。假设Sigmoid函数的输入记为z，那么z=w0x0 + w1x1 + w2x2，即可将数据分割开。或者可以理解为z=wx+b其中，wx分别代表w1w2，x1x2，x0为全是1的向量,w0=b。我们要求的就是未知的参数为w0，w1，w2，为什么时，代价函数的值最小，也就是我们需要求的回归系数(最优参数)。
2、训练算法
在编写代码之前，让我们回顾下梯度上升迭代公式：
将上述公式矢量化（矩阵相乘一定注意行列相等）：
根据矢量化的公式，编写代码如下： #创建标签列表
import matplotlib.pyplot as plt
import numpy as np
import random
#梯度上升算法
def loadDataSet():
dataMat = []
labelMat = []
fr = open('F:\machinelearning\logisitic\\testSet.txt')
for line in fr.readlines():
lineArr = line.strip().split()
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelMat.append(int(lineArr[2]))
fr.close()
return dataMat, labelMat
"""
函数说明:绘制数据集
"""
def plotDataSet():
dataMat, labelMat = loadDataSet() #加载数据集
dataArr = np.array(dataMat) #转换成numpy的array数组
n = np.shape(dataMat)[0] #数据个数
xcord1 = []; ycord1 = [] #正样本
xcord2 = []; ycord2 = [] #负样本
for i in range(n): #根据数据集标签进行分类
if int(labelMat[i]) == 1:
xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2]) #1为正样本
else:
xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2]) #0为负样本
fig = plt.figure()
ax = fig.add_subplot(111) #添加subplot
ax.scatter(xcord1, ycord1, s = 20, c = 'red', marker = 's',alpha=.5)#绘制正样本
ax.scatter(xcord2, ycord2, s = 20, c = 'green',alpha=.5) #绘制负样本
plt.title('DataSet') #绘制title
plt.xlabel('x'); plt.ylabel('y') #绘制label
plt.show()
#函数说明:sigmoid函数
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
"""
函数说明:梯度下降算法
Parameters:
dataMatIn - 数据集
classLabels - 数据标签
Returns:
weights.getA() - 求得的权重数组(最优参数)
"""
def gradAscent(dataMatIn, classLabels):
dataMatrix = np.mat(dataMatIn) #转换成numpy的mat
labelMat = np.mat(classLabels).transpose() #转换成numpy的mat,并进行转置
m, n = np.shape(dataMatrix) #返回dataMatrix的大小。m为行数,n为列数。
alpha = 0.01 #移动步长,也就是学习速率,控制更新的幅度。
maxCycles = 500 #最大迭代次数
weights = np.ones((n,1))
for k in range(maxCycles):
h = sigmoid(dataMatrix * weights) #梯度下降矢量化公式
error = h - labelMat
weights = weights - alpha * dataMatrix.transpose() * error
return weights.getA() #将矩阵转换为数组，返回权重数组
if __name__ == '__main__':
dataMat, labelMat = loadDataSet()
print(gradAscent(dataMat, labelMat))
我们已经已经求出权重参数，接下来绘制决策边界。
3、绘制决策边界
已知z=w0x0 + w1x1 + w2x2，令z=0，则0=w0 + w1x1 + w2x2。横坐标为x1，纵坐标为x2。这个方程未知的参数为w0，w1，w2，也就是我们需要求的回归系数(最优参数)。x2 =( -w0 - w1x1)/w2,编写代码。
import matplotlib.pyplot as plt
import numpy as np
def loadDataSet():
dataMat = []
labelMat = []
fr = open('F:\machinelearning\logisitic\\testSet.txt')
for line in fr.readlines():
lineArr = line.strip().split()
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelMat.append(int(lineArr[2]))
fr.close()
return dataMat, labelMat
#函数说明:sigmoid函数
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
"""
函数说明:梯度下降算法
Parameters:
dataMatIn - 数据集
classLabels - 数据标签
Returns:
weights.getA() - 求得的权重数组(最优参数)
"""
def gradAscent(dataMatIn, classLabels):
dataMatrix = np.mat(dataMatIn) #转换成numpy的mat
labelMat = np.mat(classLabels).transpose() #转换成numpy的mat,并进行转置
m, n = np.shape(dataMatrix) #返回dataMatrix的大小。m为行数,n为列数。
alpha = 0.01 #移动步长,也就是学习速率,控制更新的幅度。
maxCycles = 500 #最大迭代次数
weights = np.ones((n,1))
for k in range(maxCycles):
h = sigmoid(dataMatrix * weights) #梯度下降矢量化公式
error = h - labelMat
weights = weights - alpha * dataMatrix.transpose() * error
return weights.getA() #将矩阵转换为数组，返回权重数组
#绘制决策边界
def plotBestFit(weights):
dataMat, labelMat = loadDataSet() #加载数据集
dataArr = np.array(dataMat) #转换成numpy的array数组
n = np.shape(dataMat)[0] #数据个数
xcord1 = []; ycord1 = [] #正样本
xcord2 = []; ycord2 = [] #负样本
for i in range(n): #根据数据集标签进行分类
if int(labelMat[i]) == 1:
xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2]) #1为正样本
else:
xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2]) #0为负样本
fig = plt.figure()
ax = fig.add_subplot(111) #添加subplot
ax.scatter(xcord1, ycord1, s = 20, c = 'red', marker = 's',alpha=.5)#绘制正样本
ax.scatter(xcord2, ycord2, s = 20, c = 'green',alpha=.5) #绘制负样本
x = np.arange(-3.0, 3.0, 0.1)
#令z=0，求x2 关于x1的函数
y = (-weights[0] - weights[1] *x) / weights[2]
plt.plot(x, y)
plt.title('Bestfit') #绘制title
plt.xlabel('X1'); plt.ylabel('X2') #绘制label
plt.show()
if __name__ == '__main__':
dataMat, labelMat = loadDataSet()
#梯度下降
weights = gradAscent(dataMat, labelMat)
plotBestFit(weights)
可以看到结果还不错，但是梯度上升算法在每次更新回归系数时都需要遍历整个数据集，该方法在处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升算法。
4、改进的随机梯度上升算法
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
m,n = np.shape(dataMatrix) #返回dataMatrix的大小。m为行数,n为列数。
weights = np.ones(n) #参数初始化
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
alpha = 4/(1.0+j+i)+0.01 #降低alpha的大小，每次减小1/(j+i)。
randIndex = int(random.uniform(0,len(dataIndex))) #随机选取样本
h = sigmoid(sum(dataMatrix[dataIndex[randIndex]]*weights)) #选择随机选取的一个样本，计算h
error = h - classLabels[dataIndex[randIndex]] #计算误差
weights = weights - alpha * error * dataMatrix[dataIndex[randIndex]] #更新回归系数
del(dataIndex[randIndex]) #删除已经使用的样本
return weights
该算法第一个改进之处在于，alpha在每次迭代的时候都会调整，alpha会随着迭代次数不断减小，但永远不会减小到0，因为这里还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果需要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。另一点值得注意的是，在降低alpha的函数中，alpha每次减少1/(j+i)，其中j是迭代次数，i是样本点的下标。第二个改进的地方在于更新回归系数(最优参数)时，只使用一个样本点，并且选择的样本点是随机的，每次迭代不使用已经用过的样本点。这样的方法，就有效地减少了计算量，并保证了回归效果。
编写代码如下:
import matplotlib.pyplot as plt
import numpy as np
def loadDataSet():
dataMat = []
labelMat = []
fr = open('F:\machinelearning\logisitic\\testSet.txt')
for line in fr.readlines():
lineArr = line.strip().split()
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelMat.append(int(lineArr[2]))
fr.close()
return dataMat, labelMat
#函数说明:sigmoid函数
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
#随机梯度算法
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
m,n = np.shape(dataMatrix) #返回dataMatrix的大小。m为行数,n为列数。
weights = np.ones(n) #参数初始化
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
alpha = 4/(1.0+j+i)+0.01 #降低alpha的大小，每次减小1/(j+i)。
randIndex = int(random.uniform(0,len(dataIndex))) #随机选取样本
h = sigmoid(sum(dataMatrix[dataIndex[randIndex]]*weights)) #选择随机选取的一个样本，计算h
error = h - classLabels[dataIndex[randIndex]] #计算误差
weights = weights - alpha * error * dataMatrix[dataIndex[randIndex]] #更新回归系数
del(dataIndex[randIndex]) #删除已经使用的样本
return weights #返回
#绘制决策边界
def plotBestFit(weights):
dataMat, labelMat = loadDataSet() #加载数据集
dataArr = np.array(dataMat) #转换成numpy的array数组
n = np.shape(dataMat)[0] #数据个数
xcord1 = []; ycord1 = [] #正样本
xcord2 = []; ycord2 = [] #负样本
for i in range(n): #根据数据集标签进行分类
if int(labelMat[i]) == 1:
xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2]) #1为正样本
else:
xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2]) #0为负样本
fig = plt.figure()
ax = fig.add_subplot(111) #添加subplot
ax.scatter(xcord1, ycord1, s = 20, c = 'red', marker = 's',alpha=.5)#绘制正样本
ax.scatter(xcord2, ycord2, s = 20, c = 'green',alpha=.5) #绘制负样本
x = np.arange(-3.0, 3.0, 0.1)
#令z=0，求x2 关于x1的函数
y = (-weights[0] - weights[1] *x) / weights[2]
plt.plot(x, y)
plt.title('Bestfit') #绘制title
plt.xlabel('X1'); plt.ylabel('X2') #绘制label
plt.show()
if __name__ == '__main__':
dataMat, labelMat = loadDataSet()
#随机梯度下降
weights = stocGradAscent1(np.array(dataMat), labelMat)
plotBestFit(weights)
5、回归系数与迭代次数的关系
可以看到分类效果也是不错的。不过，从这个分类结果中，我们不好看出迭代次数和回归系数的关系，也就不能直观的看到每个回归方法的收敛情况。因此，我们编写程序，绘制出回归系数和迭代次数的关系曲线：
现在只能英文了，如果看过前面的KNN文章应该知道可以设置中文，python升级3.9了，matplotlib有点问题，以后在解决了。
由于改进的随机梯度上升算法，随机选取样本点，所以每次的运行结果是不同的。但是大体趋势是一样的。我们改进的随机梯度上升算法收敛效果更好。为什么这么说呢？让我们分析一下。我们一共有100个样本点，改进的随机梯度上升算法迭代次数为150。而上图显示15000次迭代次数的原因是，使用一次样本就更新一下回归系数。因此，迭代150次，相当于更新回归系数150*100=15000次。简而言之，迭代150次，更新1.5万次回归参数。从上图左侧的改进随机梯度上升算法回归效果中可以看出，其实在更新2000次回归系数的时候，已经收敛了。相当于遍历整个数据集20次的时候，回归系数已收敛。训练已完成。
再让我们看看上图右侧的梯度上升算法回归效果，梯度上升算法每次更新回归系数都要遍历整个数据集。从图中可以看出，当迭代次数为300多次的时候，回归系数才收敛。凑个整，就当它在遍历整个数据集300次的时候已经收敛好了。
改进的随机梯度上升算法，在遍历数据集的第20次开始收敛。而梯度上升算法，在遍历数据集的第300次才开始收敛。
评论