梯度下降的线性回归用python_python实现线性回归的梯度下降法

在机器学习中，有时候需要对原始的模型构造损失函数，然后通过优划算法对损失函数进行优划，从而找到最优的参数使损失函数达到最小值。但损失函数一般都比较复杂，难以从函数本身找到最优的参数，因此实际应用过程中使用得较多的就是梯度下降法。通过逐渐改变参数，使损失函数逐渐收敛，最终确定参数值使损失函数的值最小。

梯度下降的方式分为三种：批量梯度下降(Batch Gradient Descent)，随机梯度下降(stochastic gradient descent)，小批量随机下降(mini-batch Gradient Descent)。批量梯度下降(BGD)针对的是整个数据集，通过对所有的样本的计算来求解梯度的方向。但是在实际应用中数据量往往会很大，从而导致使用批量梯度下降的效率变得很低，因此引入了随机批量下降法(SGD)，随机梯度下降法是每次迭代只针对单个样本计算损失函数对参数的偏导数，从而寻找一个梯度来更新参数，因此效率高于BGD，但是SGD也有一定的局限性，它的收敛效果不如BGD好，因为每次只使用单个样本进行计算。而小批量随机下降(MBGD)又针对SGD进行了改进，一次性针对一小部分样本来计算，因此在提高效率的同时也比SGD更加稳定，降低了随机性。

接下来就这三种梯度下降方法针对线性回归问题来进行python的实现

批量梯度下降法

import numpy as np

X = 2 * np.random.random(size=20000).reshape(-1, 2)

y = X[:, 0] * 2. + X[:, 1] * 3. + 3. + np.random.normal(size=10000)

temp = np.ones((len(y), 1))

X_b = np.hstack((temp, X)) #为了矩阵运算方便在X中加上全为1的一列

theta = np.zeros(X_b.shape[1]) #theta是参数，梯度下降通过不断更新theta的值使损失函数达到最小值

eta = 0.01 #eta代表是学习速率

episilon = 1e-8 #episilon用来判断损失函数是否收敛

print(X_b.shape)

print(y.shape)

print(theta.shape)

这是运行结果

(10000, 3)

(10000,)

(3,)

在示例中我们创建了一个二维数组，原始的X是只有两个特征，为了后续矩阵计算的方便，我们给X添加了一个全为1的特征，添加后的结果存在X_b中，后续计算均使用X_b进行。eta代表学习速率，如果eta太小的话，theta下降得就会很慢，但是eta太大的话，参数可能就会朝着增大损失函数的方向移动。episilon用来判断损失函数是否已经收敛。

def J(theta, X_b, y):

'''

损失函数

'''

return np.sum((y - np.dot(X_b, theta))**2) / len(y)

def dJ(theta, X_b, y):

'''

损失函数对theta的偏导数

'''

gradient = X_b.T.dot(X_b.dot(theta) - y) * 2. / len(y)

return gradient

这里定义了两个函数，J代表损失函数，dJ代表损失函数对参数theta的偏导数。

def gradient_decent(theta, X_b, y):

'''

梯度下降过程

'''

while True:

last_theta = theta

theta = theta - eta * dJ(theta, X_b, y)

if abs(J(theta, X_b, y) - J(last_theta, X_b, y)) <= episilon: #判断损失函数是否收敛，也可以限定最大迭代次数

break

return theta

这是梯度下降的整个过程，判断前后两次损失函数的差值是否小于episilon，若小于的话则判断为收敛中止循环，若大于的话则继续梯度下降过程

接下来调用该函数

rst = gradient_decent(theta, X_b, y)

print(rst)

[2.97546654 2.00406701 3.0292728 ]

可以看到rst中的三个数于我们之前设置y的那个方程的三个系数大致是一样的。rst[0]代表的是截距，即原函数中没有x的常数项3.0，rst[1]代表的是X第0个特征的系数，原函数中是2.0， rst[2]代表X的第一个特征的系数，原函数中是3.0。

随机梯度下降法

由于随机梯度下降法的两次迭代的损失函数的差值随机性很强，所以一般不使用episilon来判断收敛，而是设置一个最大的迭代次数，迭代完之后就return结果。

首先还是使用之前的数据

X = 2 * np.random.random(size=20000).reshape(-1, 2)

y = X[:, 0] * 2. + X[:, 1] * 3. + 3. + np.random.normal(size=10000)

temp = np.ones((len(y), 1))

X_b = np.hstack((temp, X)) #为了矩阵运算方便在X中加上全为1的一列

theta = np.zeros(X_b.shape[1]) #theta是参数，梯度下降通过不断更新theta的值使损失函数达到最小值

eta = 0.01 #eta代表是学习速率

episilon = 1e-8 #episilon用来判断损失函数是否收敛

def dJ_sgd(theta, X_b_i, y_i):

return X_b_i.T.dot(X_b_i.dot(theta) - y_i) * 2

def sgd(X_b_i, y, theta, n_iters):

t0 = 5

t1 = 50

def learn_rate(t):

return t0/(t + t1)

theta = theta

for cur_iter in range(n_iters):

rand_i = np.random.randint(len(X_b))

gradient = dJ_sgd(theta, X_b[rand_i], y[rand_i])

theta = theta - learn_rate(cur_iter) * gradient

return theta

print(sgd(X_b, y, theta, n_iters=len(X_b)//3))

dJ_sgd是损失函数对theta的偏导数，这里传进去的参数是X_b和y中的某一个样本的数据。在实现过程中还定义了一个learn_rate的函数，因为在随机梯度下降过程中学习速率应该要随次数递减，因为存在随机性强的特点，防止下降到局部最小值附近又发生较大的改变。

运行结果

[2.99905964 2.06011999 3.01684189]

小批量随机下降

小批量随机下降与随即梯度下降思路大致相同，只不过传参的时候需要传入一部分数据

def dJ_mbgd(theta, X_b_n, y_n, num):

return X_b_n.T.dot(X_b_n.dot(theta) - y_n) * 2 / num

def mbgd(theta, X_b, y, num, n_iters):

t0 = 5

t1 = 50

theta = theta

num = num

def learn_rate(t):

return t0/(t + t1)

for cur_iter in range(n_iters):

x_index = np.random.randint(0, len(y), num)

gradient = dJ_mbgd(theta, X_b[x_index,], y[x_index], num)

theta = theta - learn_rate(cur_iter) * gradient

return theta

print(mbgd(theta, X_b, y, num=20, n_iters=len(X_b)//3))

num代表每次下降使用多少个样本数据。

运行结果

[2.98763611 2.00527425 3.01625422]