梯度下降算法代码实现

梯度下降算法的矩阵算法及代码实现:

由上期梯度下降算法中公式可看出,主要难点在于梯度的计算,即公式中$\frac{1}{m}{\sum }_{i=1}^m(h_{\theta }(xi)-yi)$ 部分的计算.
这部分的计算用矩阵形式来算的话会更容易理解一些,接下来要开始用代码矩阵来算梯度咯!
首先

导包:

import numpy as np
import os    #系统
%matplotlib inline 
import matplotlib as mpl
import matplotlib.pyplot as plt

#随机种子

np.random.seed(42)

#保存图像

DIR = '.'   # 当前目录

#定义存储图片的函数

def save_fig(fig_id,tight_layout = True):
    path = os.path.join(DIR,fig_id + '.png')    # 设置图片要保存在的路径
#path = os.path.join(DIR,'images',MODEL_ID ,fig_id +'.png')
    print('Saving figure',fig_id)               # 打印存储过程
    plt.savefig(path,format='png',dpi=300)      # 保存图片

#过滤警告信息

import warnings
warnings.filterwarnings(action='ignore',message='internal gelsd')

#生成训练数据

#(特征部分)
X = 2*np.random.rand(100,1)   
#(标签部分)
y = 4 + 3 * X + np.random.randn(100,1)  

#画原始数据的图像

plt.plot(X,y,'b.')
plt.xlabel('$x_1$',fontsize=18)
plt.ylabel('$y$',rotation=0,fontsize=18)
plt.axis([0,2,0,15])    # 横纵坐标轴刻度大小
save_fig('generated_data_plot')     # 保存图片
plt.show()

#添加新特征

X_b = np.c_[np.ones((100,1)),X]

#创建测试数据

X_new = np.array([[0],[2]])
X_new_b = np.c_[np.ones((2,1)),X_new]

1. 用现成包做线性回归模型

#从sklearn包里导入线性回归模型
from sklearn.linear_model import LinearRegression
#创建线性回归对象
lin_reg = LinearRegression()
#拟合训练数据
lin_reg.fit(X,y)
#输出截距,斜率
lin_reg.intercept_,lin_reg.coef_

#对策测试集进行预测
lin_reg.predict(X_new)

预测结果为:
array([[4.21509616],
[9.75532293]])

2.手写代码

①.批量梯度下降

eta = 0.1   # 步长,学习率
n_iterations = 1000       #迭代次数
m = 100              #数据量
theta = np.random.randn(2,1)    #随机取θ

for iteration in range(n_iterations):
    gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
    theta = theta - eta * gradients

theta_path_bgd = []
def plot_gradient_descent(theta,eta,theta_path=None):
    m = len(X_b)
    plt.plot(X,y,'b.')
    n_iterations = 1000
    for iteration in range(n_iterations):
        if iteration < 10:
            y_predict = X_new_b.dot(theta)
            style = 'b-'
            plt.plot(X_new,y_predict,style)
        gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
        theta = theta - eta * gradients
        if theta_path is not None:
            theta_path.append(theta)
    plt.xlabel('$x_1$',fontsize = 18)
    plt.axis([0,2,0,15])
    plt.title(r'$\eta = {}$'.format(eta),fontsize = 16)

np.random.seed(42)
theta = np.random.randn(2,1)

plt.figure(figsize = (10,4))
plt.subplot(131);plot_gradient_descent(theta,eta = 0.02)
plt.ylabel('$y$',rotation = 0,fontsize = 18)
plt.subplot(132);plot_gradient_descent(theta,eta = 0.1,theta_path = theta_path_bgd)
plt.subplot(133);plot_gradient_descent(theta,eta = 0.5)

save_fig('gradient_descent_plot')
plt.show()

上述代码中gradients的求法用的便是矩阵算法:

其中
X_b.dot(theta) 是100行2列的矩阵*2行1列的矩阵,求出的值即为$$h_{\theta }={\theta }_0+{\theta }_1\cdot x$$
y的值是一100行1列的矩阵,则 X_b.dot(theta) - y 的最终值为100行1列的矩阵.
X_b.T 是100行两列的X_b矩阵的转置,转制为2行,100列的矩阵,最终与上述 X_b.dot(theta) - y 最终的100行1列值相乘,结果为两个值,两列一行,求出的值即为 $${\theta }_0,{\theta }_1$$

依次循环,则gradients的值在不断更新迭代,最终绘出线性回归模型为:

②.随机梯度下降

theta_path_sgd = []
m = len(X_b)
np.random.seed(42)

n_epochs = 5

theta = np.random.randn(2,1)

for epoch in range(n_epochs):
    for i in range(m):
        if epoch == 0 and i < 20:
            y_predict = X_new_b.dot(theta)
            style = 'b-'
            plt.plot(X_new,y_predict,style)
        random_index = np.random.randint(m)
        xi = X_b[random_index:random_index + 1]
        yi = y[random_index:random_index + 1]
        gradients = 2 * xi.T.dot(xi.dot(theta) - yi)
        eta = 0.1
        theta = theta - eta * gradients
        theta_path_sgd.append(theta)
        
plt.plot(X,y,'b.')
plt.xlabel('$x_1$',fontsize = 18)
plt.ylabel('$y$',rotation=0,fontsize=18)
plt.axis([0,2,0,15])
save_fig('sgd_plot')
plt.show()

③.小批量梯度下降

theta_path_mgd = []

n_iterations = 50
minibatch_size = 20

np.random.seed(42)
theta = np.random.randn(2,1)

for epoch in range(n_iterations):
    # 乱序数据
    shuffled_indices = np.random.permutation(m)
    X_b_shuffled = X_b[shuffled_indices]
    y_shuffled = y[shuffled_indices]
    for i in range(0,m,minibatch_size):
        xi = X_b_shuffled[i:i + minibatch_size]
        yi = y_shuffled[i:i + minibatch_size]
        gradients = 2/minibatch_size * xi.T.dot(xi.dot(theta) - yi)
        eta = 0.1
        theta = theta - eta * gradients
        theta_path_mgd.append(theta)

我们再来看一下三种方法的绘制过程

theta_path_bgd = np.array(theta_path_bgd)
theta_path_sgd = np.array(theta_path_sgd)
theta_path_mgd = np.array(theta_path_mgd)

plt.figure(figsize = (7,4))
plt.plot(theta_path_sgd[:,0],theta_path_sgd[:,1],'r-s',linewidth = 1,label = 'Stochastic')
plt.plot(theta_path_bgd[:,0],theta_path_bgd[:,1],'b-o',linewidth = 3,label = 'Batch')
plt.plot(theta_path_mgd[:,0],theta_path_mgd[:,1],'g-+',linewidth = 2,label = 'Mini-batch')
plt.legend(loc = 'upper left',fontsize = 16)
plt.xlabel(r'$\theta_0$',fontsize=20)
plt.ylabel(r'$\theta_1$',fontsize=20,rotation = 0)
plt.axis([2.5,4.5,2.3,3.9])
save_fig('gradient_descent_paths_plot')
plt.show()

运行的结果图为

由图中可以看出,
批量梯度算法是蓝色那条线,目标明确但需要计算量很大且速度会比较慢;
随机梯度下降为红色线条,目标比较混乱,由随机值进行,计算量比较小但因为目标不明确,其行进的方向非常曲折;
绿色线段为小批量梯度下降,小批量梯度下降的速度与目标明确程度介于批量和随机梯度之间,相比较而言是较优选择.

今天的内容就要到此结束啦,一起前行在研究的小伙伴们或大神们,上述理解有错误或疑问的话欢迎讨论哦!~