接上篇
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数组，默认全为0
for 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=loss
else:
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) * x
alpha=0.01
lamba=0.005
'''
np.set_printoptions(precision=4) # 设置numpy输出4位小数
n=len(x[0]) # 取第一行数据，查看数据列数
theta=np.zeros(n) # 初始化theta数组，默认全为0
for times in range(loop_max):
for i in range(0,len(x)): # 取其中一条数据进行梯度下降
# i=0
y_hat=np.dot(x[i],theta.T).reshape((-1,1))
loss= (y_hat-y[i]) #此处是y_hat-y，对应的theta求解处增加了一个负号
theta_old=theta # 记录迭代前的theta
theta=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数组，默认全为0
for 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=loss
else:
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 LinearRegression
def 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):
# 最小二乘法直接求解theta
xtx = np.dot(x.T, x)
if np.linalg.det(xtx) == 0.0: # 判断xtx行列式是否等于0，奇异矩阵不能求逆
print('Can not resolve the problem')
return
theta_lsq = np.dot(np.dot(np.linalg.inv(np.dot(x.T, x)), x.T), y)
return theta_lsq
def 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数组，默认全为0
for times in range(loop_max):
y_hat=np.dot(x.T,theta)
loss= y-y_hat
for 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=theta
theta=theta+(alpha*x.T*loss_n.T).sum(axis=0) #+lamba*theta
if (theta-theta_old).all()<0.001:
break
#print('x:',x,'y:',y,'loss:',loss)
#print('times:',times,'theta:',theta)
#print('')
return theta
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数组，默认全为0
for 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=loss
else:
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
def 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) * x
alpha=0.01
lamba=0.005
'''
np.set_printoptions(precision=4) # 设置numpy输出4位小数
n=len(x[0]) # 取第一行数据，查看数据列数
theta=np.zeros(n) # 初始化theta数组，默认全为0
for times in range(loop_max):
for i in range(0,len(x)): # 取其中一条数据进行梯度下降
# i=0
y_hat=np.dot(x[i],theta.T).reshape((-1,1))
loss= (y_hat-y[i]) #此处是y_hat-y，对应的theta求解处增加了一个负号
theta_old=theta # 记录迭代前的theta
theta=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
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数组，默认全为0
for 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=loss
else:
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
if __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对应列的数据，其他列作为x
y = 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、【机器学习】线性回归的梯度下降法