Python实现线性回归2，梯度下降算法

接上篇

4.梯度下降算法

《斯坦福大学公开课 ：机器学习课程》吴恩达讲解第二课时，是直接从梯度下降开始讲解，最后采用向量和矩阵的方式推导了解析解，国内很多培训视频是先讲解析解后讲梯度下降，个人认为梯度下降算法更为重要，它是很多算法（逻辑回归、神经网络）都可以共用的，线性回归仅仅是凑巧有全局最小值的解法，有解析解，其他大部分算法都需要遍历数据寻找全局（其实本质上是局部最优）最优解的过程。
解析解直接向量运算存在的问题：（1）矩阵必须是满秩的，如果不是，python也能模糊处理。（2）运算性能，矩阵扩大化后，量级增大对性能影响明显。
梯度下降不仅限于线性回归，非线性和神经网络同样适用。

在线性回归中，θ通过α学习速率在每个J( θ) 减小的维度上进行一个 不定式 θ参数的递减，不断更新 θ{i} 的值直到 J(θ) 收敛，达到最小值。类比于下山，每次尝试找比这个位置低的位置，一步一步到达山底（即我们求的最小值，详细请看参考资料8）

4.1批量梯度下降算法

老师常常教导我们梯度下降算法这么经典的算法，一定要自己动手推导公式，自己用程序来实现以下，我一直嗤之以鼻，线性回归这么简单的算法还需要自己敲代码实现一下吗？最近工作闲暇时，我自己用Python实现了下该算法，不敲不知道，一敲吓一大跳，太多坑在里面了，整整坑在里面2天时间，检查多遍代码最后发现是数据问题，学习率α需要取非常非常小的值才能收敛，否则程序一直处于震荡状态无法找到近似最小值。说了那么多废话，下面直接贴代码，代码里备注说明很完善，很明了，很清晰了，大家应该能看懂吧o(￣︶￣)o

def linearRegBgd2(x,y,alpha=0.005,lamba=0.005,loop_max=1000):#alpha = 0.0000001      #学习率步长，learning rate#lamba=0.005     #惩罚系数 '''参考链接：https://blog.csdn.net/mango_badnot/article/details/52328740 已将原作者计算过程进行了大量修改，就相当于是我原创了吧，将原计算的三重for循环重写为一重for循环，将运算修改为了矩阵运算原作者构造的X特征是按行存储的，与常见数据不太一致，X中第一行为特征1，第二行为特征2此函数将使用正常数据，X特征按列存储，每行是一条记录批量梯度下降算法公式：theta=theta + alpha * sum( (y-y_hat) * x)'''np.set_printoptions(precision=4) # 设置numpy输出4位小数n=len(x[0]) # 取第一行数据，查看数据列数theta=np.zeros(n)      # 初始化theta数组，默认全为0for times in range(loop_max):y_hat=np.dot(x,theta.T).reshape((-1,1))loss= (y_hat-y)  #此处是y_hat-y，对应的theta求解处增加了一个负号loss_n=np.array([]) # 为方便直接用矩阵实现sum( (y-y_hat) * x)，产生多列loss乘以loss后求和if n==1: # 判断x矩阵有几列，列为1的时候单独处理loss_n=losselse:for i in range(1,n):# 判断x矩阵有几列，列为大于2的时候，产生n列loss#print(i)if i==1:loss_n = np.column_stack((loss, loss))elif i>1:loss_n = np.column_stack((loss_n, loss))theta_old=theta # 记录迭代前的theta#tmp=alpha*(loss_n*x).sum(axis=0)theta=theta - (alpha*(x*loss_n)).sum(axis=0) #+lamba*theta   #使用矩阵求解theta，注意此处的负号，loss*x后按特征列求和作为theta，注意上方本来想实现惩罚项的，但没有想明白怎么实现if (theta-theta_old).all()<0.001: # 判断前后两次theta变换情况，当变化率很小时跳出循环break#print('x:',x,'y:',y,'loss:',loss)#print('y_hat:',y_hat.T)#print('times:',times,'theta:',theta)#print('')return theta

4.2随机梯度下降算法

随机梯度下降就是每一个样本更新一次权值，超过样本范围就循环从第一个样本再循环，同等数据量情况下（数据量不是特别巨大时）训练速度与批量梯度下降要慢一些，但是适合在线学习，来一条数据迭代训练一次。

def linearRegSgd(x,y,alpha=0.005,lamba=0.005,loop_max=10000):#alpha = 0.0000001      #学习率步长，learning rate#lamba=0.005     #惩罚系数 '''随机梯度下降算法公式：theta=theta + alpha * (y-y_hat) * xalpha=0.01lamba=0.005'''np.set_printoptions(precision=4) # 设置numpy输出4位小数n=len(x[0]) # 取第一行数据，查看数据列数theta=np.zeros(n)      # 初始化theta数组，默认全为0for times in range(loop_max):for i  in range(0,len(x)): # 取其中一条数据进行梯度下降# i=0y_hat=np.dot(x[i],theta.T).reshape((-1,1))loss= (y_hat-y[i])  #此处是y_hat-y，对应的theta求解处增加了一个负号theta_old=theta # 记录迭代前的thetatheta=theta - alpha*x[i]*loss[0] #+lamba*theta#求解theta，注意此处的负号，注意上方本来想实现惩罚项的，但没有想明白怎么实现if (theta-theta_old).all()<0.001: # 判断前后两次theta变换情况，当变化率很小时跳出循环break
#            print('x:',x,'y:',y,'loss:',loss)
#            print('y_hat:',y_hat.T)
#            print('times:',times,'theta:',theta)
#            print('')return theta

4.3 mini-batch梯度下降算法

如果不是每拿到一个样本即更改梯度，而是用若干个样本的平均梯度作为更新方向，这就是Mini-batch梯度下降算法。

def linearRegMgd(x,y,alpha=0.005,lamba=0.005,loop_max=1000,batch_size=9):#alpha = 0.0000001      #学习率步长，learning rate#lamba=0.005     #惩罚系数 '''mini-batch梯度下降算法公式：每次使用batch_size个数据进行计算for i=1 to m: theta=theta + alpha * (y-y_hat) * x'''np.set_printoptions(precision=4) # 设置numpy输出4位小数n=len(x[0]) # 取第一行数据，查看数据列数theta=np.zeros(n)      # 初始化theta数组，默认全为0for times in range(loop_max):for i in range(0,int(len(x)/batch_size)+1):#print('i:',i,x[i*batch_size:(i+1)*batch_size])x_mini=x[i*batch_size:(i+1)*batch_size]y_mini=y[i*batch_size:(i+1)*batch_size]y_hat=np.dot(x_mini,theta.T).reshape((-1,1))loss= (y_hat-y_mini)  #此处是y_hat-y，对应的theta求解处增加了一个负号loss_n=np.array([]) # 为方便直接用矩阵实现sum( (y-y_hat) * x)，产生多列loss乘以loss后求和if n==1: # 判断x矩阵有几列，列为1的时候单独处理loss_n=losselse:for i in range(1,n):# 判断x矩阵有几列，列为大于2的时候，产生n列loss#print(i)if i==1:loss_n = np.column_stack((loss, loss))elif i>1:loss_n = np.column_stack((loss_n, loss))theta_old=theta # 记录迭代前的theta#tmp=alpha*(loss_n*x).sum(axis=0)theta=theta - (alpha*(x_mini*loss_n)).sum(axis=0) #+lamba*theta   #使用矩阵求解theta，注意此处的负号，loss*x后按特征列求和作为theta，注意上方本来想实现惩罚项的，但没有想明白怎么实现if (theta-theta_old).all()<0.001: # 判断前后两次theta变换情况，当变化率很小时跳出循环break#print('x:',x,'y:',y,'loss:',loss)#print('y_hat:',y_hat.T)#print('times:',times,'theta:',theta)#print('')return theta

以上梯度下降算法完整代码如下

# -*- coding: utf-8 -*-
"""@Time    : 2018/10/22 17:23@Author  : hanzi5@Email   : **@163.com@File    : linear_regression_bgd.py@Software: PyCharm
"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegressiondef load_dataset(n): # 原作者产生数据的函数noise = np.random.rand(n)X = [[x, 1.] for x in range(n)]y = [(0.5 * X[i][0]  + 1. + noise[i]) for i in range(n)]return np.array(X).T, np.array(y).T # 注意X，W，y的维数def linearRegLsq(x,y):# 最小二乘法直接求解thetaxtx = np.dot(x.T, x)if np.linalg.det(xtx) == 0.0: # 判断xtx行列式是否等于0，奇异矩阵不能求逆print('Can not resolve the problem')returntheta_lsq = np.dot(np.dot(np.linalg.inv(np.dot(x.T, x)), x.T), y)return theta_lsqdef linearRegBgd1(x,y,alpha=0.005,lamba=0.005,loop_max=1000):#alpha = 0.0000001      #学习率步长，learning rate#lamba=0.005     #惩罚系数 '''已将原作者计算过程进行了大量修改，已将原作者计算过程进行了大量修改，就相当于是我原创了吧，将原计算的三重for循环重写为一重for循环参考链接：https://blog.csdn.net/mango_badnot/article/details/52328740 原作者构造的X特征是按行存储的，与常见数据不太一致，X中第一行为特征1，第二行为特征2批量梯度下降算法公式：theta=theta + alpha * sum( (y-y_hat) * x)2018-10-22 17:18 测试通过，此函数不做多的备注了，详细备注放在了linearRegBgd2里'''np.set_printoptions(precision=4) # 设置numpy输出4位小数n=len(x) # 取第一行数据，查看数据列数theta=np.zeros(n)      # 初始化theta数组，默认全为0for times in range(loop_max):y_hat=np.dot(x.T,theta)loss= y-y_hatfor i in range(1,n):if i==1:loss_n = np.row_stack((loss, loss))elif i>1:loss_n = np.row_stack((loss_n, loss))theta_old=thetatheta=theta+(alpha*x.T*loss_n.T).sum(axis=0) #+lamba*thetaif (theta-theta_old).all()<0.001:break#print('x:',x,'y:',y,'loss:',loss)#print('times:',times,'theta:',theta)#print('')return thetadef linearRegBgd2(x,y,alpha=0.005,lamba=0.005,loop_max=1000):#alpha = 0.0000001      #学习率步长，learning rate#lamba=0.005     #惩罚系数 '''参考链接：https://blog.csdn.net/mango_badnot/article/details/52328740 已将原作者计算过程进行了大量修改，就相当于是我原创了吧，将原计算的三重for循环重写为一重for循环，将运算修改为了矩阵运算原作者构造的X特征是按行存储的，与常见数据不太一致，X中第一行为特征1，第二行为特征2此函数将使用正常数据，X特征按列存储，每行是一条记录批量梯度下降算法公式：theta=theta + alpha * sum( (y-y_hat) * x)'''np.set_printoptions(precision=4) # 设置numpy输出4位小数n=len(x[0]) # 取第一行数据，查看数据列数theta=np.zeros(n)      # 初始化theta数组，默认全为0for times in range(loop_max):y_hat=np.dot(x,theta.T).reshape((-1,1))loss= (y_hat-y)  #此处是y_hat-y，对应的theta求解处增加了一个负号loss_n=np.array([]) # 为方便直接用矩阵实现sum( (y-y_hat) * x)，产生多列loss乘以loss后求和if n==1: # 判断x矩阵有几列，列为1的时候单独处理loss_n=losselse:for i in range(1,n):# 判断x矩阵有几列，列为大于2的时候，产生n列loss#print(i)if i==1:loss_n = np.column_stack((loss, loss))elif i>1:loss_n = np.column_stack((loss_n, loss))theta_old=theta # 记录迭代前的theta#tmp=alpha*(loss_n*x).sum(axis=0)theta=theta - (alpha*(x*loss_n)).sum(axis=0) #+lamba*theta   #使用矩阵求解theta，注意此处的负号，loss*x后按特征列求和作为theta，注意上方本来想实现惩罚项的，但没有想明白怎么实现if (theta-theta_old).all()<0.001: # 判断前后两次theta变换情况，当变化率很小时跳出循环break#print('x:',x,'y:',y,'loss:',loss)#print('y_hat:',y_hat.T)#print('times:',times,'theta:',theta)#print('')return thetadef linearRegSgd(x,y,alpha=0.005,lamba=0.005,loop_max=10000):#alpha = 0.0000001      #学习率步长，learning rate#lamba=0.005     #惩罚系数 '''随机梯度下降算法公式：for i=1 to m: theta=theta + alpha * (y-y_hat) * xalpha=0.01lamba=0.005'''np.set_printoptions(precision=4) # 设置numpy输出4位小数n=len(x[0]) # 取第一行数据，查看数据列数theta=np.zeros(n)      # 初始化theta数组，默认全为0for times in range(loop_max):for i  in range(0,len(x)): # 取其中一条数据进行梯度下降# i=0y_hat=np.dot(x[i],theta.T).reshape((-1,1))loss= (y_hat-y[i])  #此处是y_hat-y，对应的theta求解处增加了一个负号theta_old=theta # 记录迭代前的thetatheta=theta - alpha*x[i]*loss[0] #+lamba*theta#求解theta，注意此处的负号，注意上方本来想实现惩罚项的，但没有想明白怎么实现if (theta-theta_old).all()<0.001: # 判断前后两次theta变换情况，当变化率很小时跳出循环break
#            print('x:',x,'y:',y,'loss:',loss)
#            print('y_hat:',y_hat.T)
#            print('times:',times,'theta:',theta)
#            print('')return thetadef linearRegMgd(x,y,alpha=0.005,lamba=0.005,loop_max=1000,batch_size=9):#alpha = 0.0000001      #学习率步长，learning rate#lamba=0.005     #惩罚系数 '''mini-batch梯度下降算法公式：每次使用batch_size个数据进行计算for i=1 to m: theta=theta + alpha * (y-y_hat) * x'''np.set_printoptions(precision=4) # 设置numpy输出4位小数n=len(x[0]) # 取第一行数据，查看数据列数theta=np.zeros(n)      # 初始化theta数组，默认全为0for times in range(loop_max):for i in range(0,int(len(x)/batch_size)+1):#print('i:',i,x[i*batch_size:(i+1)*batch_size])x_mini=x[i*batch_size:(i+1)*batch_size]y_mini=y[i*batch_size:(i+1)*batch_size]y_hat=np.dot(x_mini,theta.T).reshape((-1,1))loss= (y_hat-y_mini)  #此处是y_hat-y，对应的theta求解处增加了一个负号loss_n=np.array([]) # 为方便直接用矩阵实现sum( (y-y_hat) * x)，产生多列loss乘以loss后求和if n==1: # 判断x矩阵有几列，列为1的时候单独处理loss_n=losselse:for i in range(1,n):# 判断x矩阵有几列，列为大于2的时候，产生n列loss#print(i)if i==1:loss_n = np.column_stack((loss, loss))elif i>1:loss_n = np.column_stack((loss_n, loss))theta_old=theta # 记录迭代前的theta#tmp=alpha*(loss_n*x).sum(axis=0)theta=theta - (alpha*(x_mini*loss_n)).sum(axis=0) #+lamba*theta   #使用矩阵求解theta，注意此处的负号，loss*x后按特征列求和作为theta，注意上方本来想实现惩罚项的，但没有想明白怎么实现if (theta-theta_old).all()<0.001: # 判断前后两次theta变换情况，当变化率很小时跳出循环break#print('x:',x,'y:',y,'loss:',loss)#print('y_hat:',y_hat.T)#print('times:',times,'theta:',theta)#print('')return thetaif __name__ == "__main__":path = 'D:/python_data/8.Advertising.csv'data = pd.read_csv(path)  # TV、Radio、Newspaper、Sales#data_matrix = data.as_matrix(columns=['TV', 'Radio', 'Newspaper', 'Sales'])  # 转换数据类型data_array = data.values[:,1:]  # 转换数据类型，去除第一列序号列x = data_array[:, :-1] #去除y对应列的数据，其他列作为xy = data_array[:, -1].reshape((-1,1))# 0、绘制图像，查看数据分布情况plt.plot(data['TV'], y, 'ro', label='TV')plt.plot(data['Radio'], y, 'g^', label='Radio')plt.plot(data['Newspaper'], y, 'mv', label='Newspaer')plt.legend(loc='lower right')plt.grid()plt.show()# 1、使用sklearn LinearRegression包求解θlinreg = LinearRegression()model = linreg.fit(x, y)print('')print('1、sklearn LinearRegression包求解θ:','coef:', linreg.coef_[0],',intercept:', linreg.intercept_)# 2、最小二乘法，直接求解析解θtheta_lsq = linearRegLsq(x,y)print('')print('2、最小二乘法，theta解析解：',theta_lsq.reshape(1,3)[0])# 3、批量梯度下降求解theta# 注意下面两个函数alpha都是非常小的值，取过大的值时，不收敛#x1, y1 = load_dataset(10)#theta1=linearRegBgd1(x1, y1)theta1=linearRegBgd1(x.T, y.T,alpha = 0.0000001)print('')print('3.1、批量梯度下降，linearRegBgd1函数，theta：',theta1)theta2=linearRegBgd2(x, y,alpha = 0.0000001)print('')print('3.2、批量梯度下降，linearRegBgd2函数，theta：',theta2)theta3=linearRegSgd(x, y,alpha = 0.000001,loop_max=10000)print('')print('3.3、随机梯度下降，linearRegSgd函数，theta：',theta3)theta4=linearRegMgd(x, y,alpha = 0.000001,loop_max=10000,batch_size=11)print('')print('3.4、mini-batch梯度下降，linearRegMgd函数，theta：',theta4)

程序运行计算结果：

数据见以下链接，复制到txt文档中，另存为CSV格式即可。
数据

参考资料：
1、python实现线性回归
2、机器学习：线性回归与Python代码实现
3、python实现简单线性回归
4、Machine-Learning-With-Python
5、《机器学习》西瓜书，周志华
6、机器学习视频，邹博
7、 斯坦福大学公开课 ：机器学习课程
8、马同学高等数学 如何直观形象的理解方向导数与梯度以及它们之间的关系？
9、【机器学习】线性回归的梯度下降法