百度网站开发,wordpress 一栏主题,域名卖给别人有风险吗,做流量的网站【作者主页】Francek Chen 【专栏介绍】 ⌈ ⌈ ⌈Python机器学习 ⌋ ⌋ ⌋ 机器学习是一门人工智能的分支学科#xff0c;通过算法和模型让计算机从数据中学习#xff0c;进行模型训练和优化#xff0c;做出预测、分类和决策支持。Python成为机器学习的首选语言#xff0c;… 【作者主页】Francek Chen 【专栏介绍】 ⌈ ⌈ ⌈Python机器学习 ⌋ ⌋ ⌋ 机器学习是一门人工智能的分支学科通过算法和模型让计算机从数据中学习进行模型训练和优化做出预测、分类和决策支持。Python成为机器学习的首选语言依赖于强大的开源库如Scikit-learn、TensorFlow和PyTorch。本专栏介绍机器学习的相关算法以及基于Python的算法实现。 【GitCode】专栏资源保存在我的GitCode仓库https://gitcode.com/Morse_Chen/Python_machine_learning。 文章目录 一、线性回归的映射形式和学习目标二、线性回归的解析方法三、动手实现线性回归的解析方法四、使用sklearn中的线性模型五、梯度下降算法六、学习率对迭代的影响 本文将逐步引入一些数学工具讲解另一个较为简单的机器学习算法——线性回归linear regression。与上一篇文章介绍的k近邻算法不同线性回归是一种基于数学模型的算法其首先假设数据集中的样本与标签之间存在线性关系再建立线性模型求解该关系中的各个参数。在实际生活中线性回归算法因为其简单易算在统计学、经济学、天文学、物理学等领域中都有着广泛应用。下面我们从线性回归的数学描述开始讲解线性回归的原理和实践。
一、线性回归的映射形式和学习目标 顾名思义在“线性”回归问题中我们假设输入与输出成线性关系。设输入 x ∈ R d \boldsymbol{x}\in\mathbb{R}^d x∈Rd那么该线性映射关系可以写为 f θ ( x ) θ T x θ 1 x 1 ⋯ θ d x d f_{\boldsymbol{\theta}}(\boldsymbol{x})\boldsymbol{\theta}^{\mathrm{T}}\boldsymbol{x}\theta_1x_1\cdots\theta_dx_d fθ(x)θTxθ1x1⋯θdxd 其中 θ ∈ R d \boldsymbol{\theta}\in\mathbb{R}^d θ∈Rd 是模型的参数我们通常以脚标 f θ f_{\boldsymbol{\theta}} fθ的形式来表示 θ \boldsymbol{\theta} θ是模型 f f f的参数。如果线性回归需要包含常数项我们只需再添加一维参数 θ 0 \theta_0 θ0以及对应的常数特征 x 0 1 x_01 x01 即可。这样上式的形式没有变化只是向量的维度由 d d d变为 d 1 d1 d1。为了表达的简洁下面不再将常数项单独拆开表示和分析。图1展示了输入特征是一维和二维情况下的数据点和线性回归模型拟合的结果。在 d d d维输入特征和一维输出特征的情况下线性回归模型共有 d 1 d1 d1 个参数从而给出了 d 1 d1 d1 维空间中的一个 d d d维超平面。 图1 一维和二维输入特征的线性模型 在机器学习中我们一般先设计损失函数由模型预测的标签与真实标签之间的误差计算损失的值并通过最小化损失函数来训练模型调整模型参数。设共有 N N N个输入数据 x 1 , ⋯ , x N \boldsymbol{x}_1,\cdots,\boldsymbol{x}_N x1,⋯,xN其对应的标签分别是 y 1 , ⋯ , y N y_1,\cdots,y_N y1,⋯,yN那么模型的总损失为 J ( θ ) 1 N ∑ i 1 N L ( y i , f θ ( x i ) ) J(\boldsymbol{\theta})\frac{1}{N}\sum_{i1}^N\mathcal{L}(y_i,f_{\boldsymbol\theta}(\boldsymbol{x}_i)) J(θ)N1i1∑NL(yi,fθ(xi)) 其中 L ( y i , f θ ( x i ) ) \mathcal{L}(y_i,f_{\boldsymbol\theta}(\boldsymbol{x}_i)) L(yi,fθ(xi)) 为单个样本的损失函数用来衡量真实标签与预测标签之间的距离。我们定义损失函数为 L ( y i , f θ ( x i ) ) 1 2 ( y i − f θ ( x i ) ) 2 \mathcal{L}(y_i,f_{\boldsymbol\theta}(\boldsymbol{x}_i))\frac{1}{2}(y_i-f_{\boldsymbol\theta}(\boldsymbol{x}_i))^2 L(yi,fθ(xi))21(yi−fθ(xi))2 式中的系数 1 2 \frac{1}{2} 21是为了求导后与系数2抵消便于计算。均方误差的形状如图2所示。我们可以看到在 y i y_i yi与 f θ ( x i ) f_{\boldsymbol\theta}(\boldsymbol{x}_i) fθ(xi)距离较小的情况下损失也较小并且变化不大而随着两者的距离增大其损失以二次速度迅速增长。这样的损失函数的设计可以让模型倾向忽略预测已经很精准的数据而重点关注预测和标签差距较大的数据。要知道特征和标签数据的采集经常会带上些许的偏差或者噪声例如我们做物理实验时测量一个物体的尺寸往往需要测量多次后取平均来降低测量噪声带来的不确定度。因此当预测和标签已经足够接近时没有必要将精力放在进一步消除最后一点损失上。 图2 均方误差示意 将该损失函数代入总误差可得 J ( θ ) 1 2 N ∑ i 1 N ( y i − f θ ( x i ) ) 2 J(\boldsymbol{\theta})\frac{1}{2N}\sum_{i1}^N(y_i-f_{\boldsymbol\theta}(\boldsymbol{x}_i))^2 J(θ)2N1i1∑N(yi−fθ(xi))2 这一总损失函数称为均方误差mean squared errorMSE是最常用的损失函数之一。因此线性回归问题的优化目标为 min θ J ( θ ) min θ 1 2 N ∑ i 1 N ( y i − f θ ( x i ) ) 2 \min_{\theta} J(\boldsymbol{\theta})\min_{\theta}\frac{1}{2N}\sum_{i1}^N(y_i-f_{\boldsymbol\theta}(\boldsymbol{x}_i))^2 θminJ(θ)θmin2N1i1∑N(yi−fθ(xi))2
二、线性回归的解析方法 我们使用正规的解方程法即最小二乘法。为了使表达式更简洁我们进一步将数据聚合把输入向量和标签组合成矩阵 X ( x 1 T ⋮ x N T ) y ( y 1 ⋮ y N ) \boldsymbol{X} \begin{pmatrix} \boldsymbol{x}_1^\mathrm{T} \\ \vdots \\ \boldsymbol{x}_N^\mathrm{T} \end{pmatrix}\ \ y\begin{pmatrix} y_1 \\ \vdots \\ y_N \end{pmatrix} X x1T⋮xNT y y1⋮yN 其中 X \boldsymbol{X} X的每一行对应一个数据实例的特征向量每一列对应了一个具体特征在各个数据实例上的取值。这样以向量的平方作为损失函数可将总损失写为 J ( θ ) 1 2 N ( y − X θ ) T ( y − X θ ) J(\boldsymbol{\theta})\frac{1}{2N}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})^{\mathrm{T}}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}) J(θ)2N1(y−Xθ)T(y−Xθ) 这里的损失函数事实上是所有样本平方误差的和与MSE相差一步取平均值即除以样本总数 N N N。但常数系数不影响优化得到的最终结果为了形式简洁我们通常在矩阵形式下省略这一系数。但是在实际计算时为了降低样本规模的影响绝大多数情况下还是以平均后的值作为实际损失。 为了求函数 J ( θ ) J(\boldsymbol{\theta}) J(θ)的最小值我们寻找其对 θ \boldsymbol\theta θ导数为零的点即 ∂ J ∂ θ 0 \frac{\partial J}{\partial\boldsymbol{\theta}}\boldsymbol0 ∂θ∂J0 计算偏导数得到 0 ∂ J ∂ θ ∂ ∂ θ 1 2 ( y T y − y T X θ − ( X θ ) T y ( X θ ) T X θ ) 1 2 ( ∂ y T y ∂ θ − ∂ y T X θ ∂ θ − ∂ ( X θ ) T y ∂ θ ∂ ( X θ ) T X θ ∂ θ ) − X T y X T X θ \begin{aligned} \boldsymbol0\frac{\partial J}{\partial\boldsymbol{\theta}}\frac{\partial}{\partial\boldsymbol{\theta}}\frac{1}{2}(\boldsymbol{y}^{\mathrm{T}}\boldsymbol{y}-\boldsymbol{y}^{\mathrm{T}}\boldsymbol{X}\boldsymbol{\theta}-(\boldsymbol{X}\boldsymbol{\theta})^{\mathrm{T}}\boldsymbol{y}(\boldsymbol{X}\boldsymbol{\theta})^{\mathrm{T}}\boldsymbol{X}\boldsymbol{\theta}) \\[2ex] \frac{1}{2}\left(\frac{\partial\boldsymbol{y}^{\mathrm{T}}\boldsymbol{y}}{\partial\boldsymbol{\theta}}-\frac{\partial\boldsymbol{y}^{\mathrm{T}}\boldsymbol{X}\boldsymbol{\theta}}{\partial\boldsymbol{\theta}}-\frac{\partial(\boldsymbol{X}\boldsymbol{\theta})^{\mathrm{T}}\boldsymbol{y}}{\partial\boldsymbol{\theta}}\frac{\partial(\boldsymbol{X}\boldsymbol{\theta})^{\mathrm{T}}\boldsymbol{X}\boldsymbol{\theta}}{\partial\boldsymbol{\theta}}\right) \\[2ex] -\boldsymbol{X}^{\mathrm{T}}\boldsymbol{y}\boldsymbol{X}^{\mathrm{T}}\boldsymbol{X}\boldsymbol\theta \end{aligned} 0∂θ∂J∂θ∂21(yTy−yTXθ−(Xθ)Ty(Xθ)TXθ)21(∂θ∂yTy−∂θ∂yTXθ−∂θ∂(Xθ)Ty∂θ∂(Xθ)TXθ)−XTyXTXθ 通过上式得到 θ \boldsymbol\theta θ的解析解为 0 − X T y X T X θ ⇒ X T X θ X T y ⇒ θ ( X T X ) − 1 X T y \begin{aligned} \boldsymbol0 -\boldsymbol{X}^{\mathrm{T}}\boldsymbol{y}\boldsymbol{X}^{\mathrm{T}}\boldsymbol{X}\boldsymbol\theta \\ \Rightarrow \ \ \boldsymbol{X}^{\mathrm{T}}\boldsymbol{X}\boldsymbol\theta \boldsymbol{X}^{\mathrm{T}}\boldsymbol{y} \\ \Rightarrow \ \ \ \ \ \ \ \ \ \ \ \ \boldsymbol\theta (\boldsymbol{X}^{\mathrm{T}}\boldsymbol{X})^{-1}\boldsymbol{X}^{\mathrm{T}}\boldsymbol{y} \end{aligned} 0⇒ XTXθ⇒ θ−XTyXTXθXTy(XTX)−1XTy 于是算法学到的模型对训练数据的预测为 f θ ( X ) X θ X ( X T X ) − 1 X T y f_{\boldsymbol{\theta}}(\boldsymbol{X})\boldsymbol{X}\boldsymbol\theta\boldsymbol{X}(\boldsymbol{X}^{\mathrm{T}}\boldsymbol{X})^{-1}\boldsymbol{X}^{\mathrm{T}}\boldsymbol{y} fθ(X)XθX(XTX)−1XTy 最小二乘法的缺点当 X T X \boldsymbol{X}^{\mathrm{T}}\boldsymbol{X} XTX不可逆时无法求解即使可逆逆矩阵求解可能计算很复杂求得的权系数向量 θ \boldsymbol\theta θ可能不稳定即样本数据的微小变化可能导致 θ \boldsymbol\theta θ的巨大变化从而使得回归模型不稳定缺乏泛化能力。 下面我们用NumPy库中的线性代数相关工具直接用解析解来计算线性回归模型。
三、动手实现线性回归的解析方法 采用的数据集由房屋信息与房屋售价组成。其中房屋信息包含所在区域平均收入、区域平均房屋年龄、区域平均房间数、区域平均卧室数、人口等。我们希望根据某一区域中房屋的整体信息用线性模型预测该区域中房屋的平均售价。表1展示了其中的3条数据。 表1 房屋信息数据示例 区域平均收入区域平均房屋年龄区域平均房间数区域平均卧室数区域人口房屋售价79545.465.687.014.0923086.801059033.5679248.646.006.733.0940173.071505890.9161287.065.868.515.1336882.151058987.98 我们首先读入并处理数据并且划分训练集与测试集为后续算法实现做准备。这里我们会用到sklearn中的数据处理工具包preprocessing中的StandardScalar类。该类的fit函数可以根据输入的数据计算平均值和方差并用计算结果将数据标准化使其均值为 0、方差为 1。例如数组[0, 1, 2, 3, 4]经过标准化就变为[-1.41, -0.71, 0.00, 0.71, 1.41]。对输入数据进行标准化可以避免不同特征的数据之间数量级差距过大导致的问题。以上面列出的数据条目为例区域平均收入在 1 0 4 10^4 104数量级而区域平均卧室数在 1 0 0 10^0 100数量级如果直接用原始数据进行训练假如区域平均收入在运算中产生了 0.1 0.1 0.1的误差约为 1 0 1 10^1 101就几乎足够掩盖区域平均卧室数带来的影响。因此通常来说我们在训练前将不同特征的数据放缩到同一量级上。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from sklearn.preprocessing import StandardScaler# 从源文件加载数据并输出查看数据的各项特征
lines np.loadtxt(USA_Housing.csv, delimiter,, dtypestr)
header lines[0]
lines lines[1:].astype(float)
print(数据特征, , .join(header[:-1]))
print(数据标签, header[-1])
print(数据总条数, len(lines))# 划分训练集与测试集
ratio 0.8
split int(len(lines) * ratio)
np.random.seed(0)
lines np.random.permutation(lines)
train, test lines[:split], lines[split:]# 数据归一化
scaler StandardScaler()
scaler.fit(train) # 只使用训练集的数据计算均值和方差
train scaler.transform(train)
test scaler.transform(test)# 划分输入和标签
x_train, y_train train[:, :-1], train[:, -1].flatten()
x_test, y_test test[:, :-1], test[:, -1].flatten()我们按照线性回归的解析方法的推导利用NumPy库中的工具直接进行矩阵运算并输出预测值与真实值的误差。衡量误差的标准也有很多这里我们采用均方根误差rooted mean squared errorRMSE。对于真实值 y 1 , y 2 , ⋯ , y N y_1,y_2,\cdots,y_N y1,y2,⋯,yN和预测值 y ^ 1 , y ^ 2 , ⋯ , y ^ N \hat{y}_1,\hat{y}_2,\cdots,\hat{y}_N y^1,y^2,⋯,y^NRMSE为 L R M S E ( y , y ^ ) 1 N ∑ i 1 N ( y i − y ^ i ) 2 \mathcal{L}_\mathrm{RMSE}(\boldsymbol{y},\hat{\boldsymbol{y}})\sqrt{\frac{1}{N}\sum_{i1}^N(y_i-\hat{y}_i)^2} LRMSE(y,y^)N1i1∑N(yi−y^i)2 RMSE与MSE非常接近但是平方再开方的操作使得RMSE应当与 y \boldsymbol{y} y具有相同的量纲从直观上易于比较。我们可以简单认为对于任意样本 x \boldsymbol{x} x模型预测的标签 y ^ \hat{y} y^与真实值 y y y之间的偏差大致就等于RMSE的值。而MSE由于含有平方其量纲和数量级相对来说不够直观但其更容易求导。因此我们常将MSE作为训练时的损失函数而用RMSE作为模型的评价指标。
# 在X矩阵最后添加一列1代表常数项
X np.concatenate([x_train, np.ones((len(x_train), 1))], axis-1)
# 表示矩阵相乘X.T表示矩阵X的转置np.linalg.inv函数可以计算矩阵的逆
theta np.linalg.inv(X.T X) X.T y_train
print(回归系数, theta)# 在测试集上使用回归系数进行预测
X_test np.concatenate([x_test, np.ones((len(x_test), 1))], axis-1)
y_pred X_test theta# 计算预测值和真实值之间的RMSE
rmse_loss np.sqrt(np.square(y_test - y_pred).mean())
print(RMSE, rmse_loss)四、使用sklearn中的线性模型 接下来我们使用sklearn中已有的工具LinearRegression来实现线性回归模型。可以看出该工具计算得到的回归系数与RMSE都和我们用解析方式计算的结果相同。
from sklearn.linear_model import LinearRegression# 初始化线性模型
linreg LinearRegression()
# LinearRegression的方法中已经考虑了线性回归的常数项所以无须再拼接1
linreg.fit(x_train, y_train)# coef_是训练得到的回归系数intercept_是常数项
print(回归系数, linreg.coef_, linreg.intercept_)
y_pred linreg.predict(x_test)# 计算预测值和真实值之间的RMSE
rmse_loss np.sqrt(np.square(y_test - y_pred).mean())
print(RMSE, rmse_loss)对数据进行可视化并观察预测回归直线和预测误差分布。
import matplotlib.pyplot as plt# 可视化
plt.figure(figsize(12, 6))# 真实值与预测值对比
plt.subplot(1, 2, 1)
plt.scatter(y_test, y_pred, alpha0.6, edgecolorsw, linewidth0.5)
plt.xlabel(真实值)
plt.ylabel(预测值)
plt.title(真实值 vs 预测值)
plt.gca().set_aspect(equal, adjustablebox)
plt.plot([min(y_test), max(y_test)], [min(y_test), max(y_test)], r--, lw2)# 预测误差分布
plt.subplot(1, 2, 2)
errors y_test - y_pred
plt.hist(errors, bins30, edgecolork, alpha0.7)
plt.xlabel(预测误差)
plt.ylabel(频数)
plt.title(预测误差分布)plt.rcParams[font.sans-serif] [SimHei]
plt.rcParams[axes.unicode_minus] Falseplt.tight_layout()
plt.show()五、梯度下降算法 虽然对于线性回归问题我们在选取平方损失函数后可以通过数学推导得到问题的解析解。但是这样的做法有一些严重的缺陷。第一解析解中涉及大量的矩阵运算非常耗费时间和空间。假设样本数目为 N N N特征维度为 d d d那么 X ∈ R N × d \boldsymbol{X}\in\mathbb{R}^{N\times d} X∈RN×d y ∈ R N \boldsymbol{y}\in\mathbb{R}^{N} y∈RN。按照式 θ ( X T X ) − 1 X T y \boldsymbol\theta (\boldsymbol{X}^{\mathrm{T}}\boldsymbol{X})^{-1}\boldsymbol{X}^{\mathrm{T}}\boldsymbol{y} θ(XTX)−1XTy 进行计算的时间复杂度大约是 O ( N d 2 d 3 ) O(Nd^2d^3) O(Nd2d3)。虽然我们可以通过矩阵运算技巧进行优化但时间开销仍然较大。此外当样本很多时存储矩阵 X \boldsymbol{X} X也会占用大量空间。第二在更广泛的机器学习模型中大多数情况下我们都无法得到解析解或求解析解非常困难。因此我们通常会采用数值模拟的方法避开复杂的计算经过一定次数的迭代得到与解析解误差很小的数值解。本节继续以平方损失函数的线性回归为例介绍机器学习中非常常用的数值计算方法梯度下降gradient decent, GD算法。 回顾梯度的意义我们可以发现梯度的方向就是函数值上升最快的方向。那么反过来说梯度的反方向就是函数值下降最快的方向。如果我们将参数不断沿梯度的反方向调整就可以使函数值以最快的速度减小。当函数值几乎不再改变时我们就找到了函数的一个局部极小值。而对于部分较为特殊的函数其局部极小值就是其全局最小值。在此情况下梯度下降算法最后可以得到全局最优解。我们暂时不考虑具体哪些函数满足这样的条件而是按照直观的思路进行简单的推导。设模型参数为 θ \boldsymbol\theta θ损失函数为 J ( θ ) J(\boldsymbol\theta) J(θ)。那么梯度下降的公式为 θ ← θ − η ∇ θ ( θ ) \boldsymbol\theta\leftarrow\boldsymbol\theta-\eta\nabla_{\boldsymbol\theta}(\boldsymbol\theta) θ←θ−η∇θ(θ) 其中 ← \leftarrow ←表示令左边的变量等于右边表达式的值 η \eta η是参数更新的步长称为学习率learning rate。我们将上节推导的带平均的线性回归损失函数 J ( θ ) 1 2 N ∑ i 1 N ( y i − f θ ( x i ) ) 2 J(\boldsymbol{\theta})\frac{1}{2N}\sum\limits_{i1}^N(y_i-f_{\boldsymbol\theta}(\boldsymbol{x}_i))^2 J(θ)2N1i1∑N(yi−fθ(xi))2 代入进去就得到 θ ← θ − η ∇ θ ( 1 2 N ∑ i 1 N ( y i − f θ ( x i ) ) 2 ) θ − η N ∑ i 1 N ( f θ ( x i ) − y i ) ∇ θ f θ ( x i ) θ − η N ∑ i 1 N ( f θ ( x i ) − y i ) x i \begin{aligned} \boldsymbol\theta\leftarrow\boldsymbol\theta-\eta\nabla_{\boldsymbol\theta}\left(\frac{1}{2N}\sum_{i1}^N(y_i-f_{\boldsymbol\theta}(\boldsymbol{x}_i))^2\right) \\[3ex] \boldsymbol\theta-\frac{\eta}{N}\sum_{i1}^N(f_{\boldsymbol\theta}(\boldsymbol{x}_i)-y_i)\nabla_{\boldsymbol\theta}f_{\boldsymbol\theta}(\boldsymbol{x}_i) \\[3ex] \boldsymbol\theta-\frac{\eta}{N}\sum_{i1}^N(f_{\boldsymbol\theta}(\boldsymbol{x}_i)-y_i)\boldsymbol{x}_i \end{aligned} θ←θ−η∇θ(2N1i1∑N(yi−fθ(xi))2)θ−Nηi1∑N(fθ(xi)−yi)∇θfθ(xi)θ−Nηi1∑N(fθ(xi)−yi)xi 如果写成矩阵形式上式就等价于 θ ← θ − η ∇ θ ( 1 2 N ( y − X θ ) T ( y − X θ ) ) θ − η N ( − X T y X T X θ ) θ − η N X T ( f θ ( X ) − y ) \begin{aligned} \boldsymbol\theta\leftarrow\boldsymbol\theta-\eta\nabla_{\boldsymbol\theta}\left(\frac{1}{2N}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})^{\mathrm{T}}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})\right) \\[2ex] \boldsymbol\theta-\frac{\eta}{N}(-\boldsymbol{X}^{\mathrm{T}}\boldsymbol{y}\boldsymbol{X}^{\mathrm{T}}\boldsymbol{X}\boldsymbol\theta) \\[2ex] \boldsymbol\theta-\frac{\eta}{N}\boldsymbol{X}^{\mathrm{T}}(f_{\boldsymbol{\theta}}(\boldsymbol{X})-\boldsymbol{y}) \end{aligned} θ←θ−η∇θ(2N1(y−Xθ)T(y−Xθ))θ−Nη(−XTyXTXθ)θ−NηXT(fθ(X)−y) 图3展示了在二维情况下MSE损失函数梯度下降的迭代过程。图中的椭圆线条是损失函数的等值线颜色表示损失函数值的大小越偏蓝色的地方损失越小越偏红色的地方损失越大。3条曲线代表了从3个不同的初值进行梯度下降的过程中参数值的变化情况。可以看出在MSE损失的情况下无论初始值如何参数都不断沿箭头所指的梯度反方向变化最终到达损失函数的最小值点 θ ∗ \boldsymbol{\theta}^* θ∗。这是因为MSE损失函数关于参数是凸函数所以无论起点如何沿梯度方向都可以到达是损失函数最小的地方。 图3 二维MSE损失函数的梯度下降示意 如果需要优化的损失函数是非凸的梯度下降就可能陷入局部极小值无法达到全局最优。如图4所示损失函数在空间上有两个局部极小值点 θ 1 ∗ \boldsymbol{\theta}_1^* θ1∗和 θ 2 ∗ \boldsymbol{\theta}_2^* θ2∗其中 θ 1 ∗ \boldsymbol{\theta}_1^* θ1∗上方的是全局最小值两条曲线分别代表从不同的起始参数出发进行梯度下降的结果。可以看出如果参数的初始值比较靠下梯度下降算法就只能收敛到较差的解 θ 2 ∗ \boldsymbol{\theta}_2^* θ2∗上。图中的两个起始参数位置其实非常接近在这样的情况下模型得到的结果受初始的随机参数以及计算误差的影响非常大从而非常不稳定。当模型较为复杂时我们通常无法直观判断损失函数是否为凸函数也无法先验地得知函数有几个局部极小值点、其中哪个才是全局最优还很难控制随机生成的初始参数的位置。因此在现代的深度学习中非凸函数的优化仍然是一个重要的研究课题。 图4 有两个局部极小值的损失函数 梯度下降的公式已经不含有矩阵求逆和矩阵相乘等时间复杂度很高的运算但当样本量很大时计算矩阵与向量相乘仍然很耗时矩阵的存储问题也没有解决。因此我们可以每次只随机选一个样本计算其梯度并进行梯度下降。设选取的样本为 x k \boldsymbol{x}_k xk其参数更新可以写为 θ ← θ − η ∇ θ ( 1 2 ( y k − θ T x k ) 2 ) θ − η ( θ T x k − y k ) x k \begin{aligned} \boldsymbol\theta\leftarrow\boldsymbol\theta-\eta\nabla_{\boldsymbol\theta}\left(\frac{1}{2}(y_k-\boldsymbol{\theta}^{\mathrm{T}}\boldsymbol{x}_k)^2\right) \\[2ex] \boldsymbol\theta-\eta(\boldsymbol{\theta}^{\mathrm{T}}\boldsymbol{x}_k-y_k)\boldsymbol{x}_k \end{aligned} θ←θ−η∇θ(21(yk−θTxk)2)θ−η(θTxk−yk)xk 由于每次只计算一个样本梯度下降的时间复杂度大大降低了这样的算法就称为随机梯度下降法stochastic gradient decentSGD。然而随机梯度下降的不稳定性很高。由于单个样本计算出的梯度方向可能与所有样本算出的真正梯度方向不同如果我们要优化的函数不是凸函数SGD 算法就可能从原定路线偏离收敛到其他极小值点去。因此为了在稳定性与时间复杂度之间取得平衡我们一般使用小批量梯度下降mini-batch gradient decentMBGD算法将样本随机分成许多大小较小的批量batch。每次迭代时选取一个批量来计算函数梯度以此估计用完整样本计算的结果。假如批量大小为 B ≪ N B\ll N B≪N第 i i i个小批量中的数据为 X ( i ) \boldsymbol{X}_{(i)} X(i)和 y ( i ) \boldsymbol{y}_{(i)} y(i)那么相应的梯度下降公式变为 θ ← θ − η B X ( i ) T ( f θ ( X ( i ) T ) − y ( i ) ) \boldsymbol\theta\leftarrow\boldsymbol\theta-\frac{\eta}{B}\boldsymbol{X}_{(i)}^{\mathrm{T}}(f_{\boldsymbol{\theta}}(\boldsymbol{X}_{(i)}^{\mathrm{T}})-\boldsymbol{y}_{(i)}) θ←θ−BηX(i)T(fθ(X(i)T)−y(i)) 可以看出MBGD算法当 B 1 B1 B1 时就退化为SGD算法当 B N BN BN 时就退化为GD算法。对每个小批量来说用来计算梯度的矩阵 X \boldsymbol{X} X的大小就从 N × d N\times d N×d 下降到 B × d B\times d B×d时间和空间复杂度同样大大降低了。对于MBGD来说反复随机抽样进行迭代大概率可以在平均意义上消除梯度估计的偏差并且还可以通过调整 B B B来控制随机性的大小。图5中展示了在MSE损失下进行梯度下降时参数的轨迹绿色实线表示从起点进行梯度下降蓝色虚线表示从MBGD红色虚线表示SGD。可以看出全样本的GD每次计算出的都是精确的梯度值下降轨迹完全沿梯度方向MBGD的轨迹存在一定的振荡但是始终与实线偏离不远并且最后也可以收敛到最优解的位置SGD的轨迹则震荡很大并且在最优解附近时由于其随机性较大其轨迹会在周围反复抖动很难真正收敛到最优解。虽然MBGD也有在最优解附近抖动的情况但是其抖动幅度较小在合适的情况下可以认为它得到了最优解很好的近似。 图5 几种梯度下降算法的参数变化轨迹 虽然SGD与MBGD理论上是不同的算法但是作为SGD的扩展在现代深度学习中SGD已经成为了MBGD的代名词在不少代码库中SGD就是MBGD。因此不再区分这两个算法统一称为SGD。下面我们来动手实现SGD算法完成相同的线性回归任务并观察SGD算法的表现。首先我们实现随机划分数据集和产生批量的函数。
# 该函数每次返回大小为batch_size的批量
# x和y分别为输入和标签
# 若shuffle True则每次遍历时会将数据重新随机划分
def batch_generator(x, y, batch_size, shuffleTrue):# 批量计数器batch_count 0if shuffle:# 随机生成0到len(x)-1的下标idx np.random.permutation(len(x))x x[idx]y y[idx]while True:start batch_count * batch_sizeend min(start batch_size, len(x))if start end:# 已经遍历一遍结束生成breakbatch_count 1yield x[start: end], y[start: end]接下来是算法的主体部分。我们提前设置好迭代次数、学习率和批量大小并用上面的梯度下降公式 θ ← θ − η B X ( i ) T ( f θ ( X ( i ) T ) − y ( i ) ) \boldsymbol\theta\leftarrow\boldsymbol\theta-\frac{\eta}{B}\boldsymbol{X}_{(i)}^{\mathrm{T}}(f_{\boldsymbol{\theta}}(\boldsymbol{X}_{(i)}^{\mathrm{T}})-\boldsymbol{y}_{(i)}) θ←θ−BηX(i)T(fθ(X(i)T)−y(i)) 不断迭代最后将迭代过程中RMSE的变化曲线绘制出来。可以看出最终得到的结果和上面精确计算的结果虽然有差别但已经十分接近RMSE也在可以接受的范围内。另外在模型优化中我们一般将一次参数更新称为一步step例如进行一次梯度下降而将遍历一次所有训练数据称为一轮epoch。
def SGD(num_epoch, learning_rate, batch_size):# 拼接原始矩阵X np.concatenate([x_train, np.ones((len(x_train), 1))], axis-1)X_test np.concatenate([x_test, np.ones((len(x_test), 1))], axis-1)# 随机初始化参数theta np.random.normal(sizeX.shape[1])# 随机梯度下降# 为了观察迭代过程我们记录每一次迭代后在训练集和测试集上的均方根误差train_losses []test_losses []for i in range(num_epoch):# 初始化批量生成器batch_g batch_generator(X, y_train, batch_size, shuffleTrue)train_loss 0for x_batch, y_batch in batch_g:# 计算梯度grad x_batch.T (x_batch theta - y_batch)# 更新参数theta theta - learning_rate * grad / len(x_batch)# 累加平方误差train_loss np.square(x_batch theta - y_batch).sum()# 计算训练和测试误差train_loss np.sqrt(train_loss / len(X))train_losses.append(train_loss)test_loss np.sqrt(np.square(X_test theta - y_test).mean())test_losses.append(test_loss)# 输出结果绘制训练曲线print(回归系数, theta)return theta, train_losses, test_losses# 设置迭代次数学习率与批量大小
num_epoch 20
learning_rate 0.01
batch_size 32
# 设置随机种子
np.random.seed(0)_, train_losses, test_losses SGD(num_epoch, learning_rate, batch_size)# 将损失函数关于运行次数的关系制图可以看到损失函数先一直保持下降之后趋于平稳
plt.plot(np.arange(num_epoch), train_losses, colorblue, labeltrain loss)
plt.plot(np.arange(num_epoch), test_losses, colorred, ls--, labeltest loss)
# 由于epoch是整数这里把图中的横坐标也设置为整数
# 该步骤也可以省略
plt.gca().xaxis.set_major_locator(MaxNLocator(integerTrue))
plt.xlabel(Epoch)
plt.ylabel(RMSE)
plt.legend()
plt.show()解析法与梯度下降法比较 解析方法与梯度下降方法各有优劣。前者能直接得到精确的解但计算的时空开销较大其表达式一般情况下也难以计算后者通过数值近似方法用较小的时空复杂度得到了与精确解接近的结果但是往往需要手动调整学习率和迭代次数。整体而言梯度下降方法更为常用它也衍生出了许多更有效、更快速的优化方法是现代深度学习的基础之一。
思考
1调整SGD算法中batch_size的大小观察结果的变化。对较大规模的数据集batch_size过小和过大分别有什么缺点 过小的批量大小训练过程中需要更多的参数更新步骤因此训练速度可能会变慢。较小的批量可能无法充分利用计算资源导致训练过程中的并行性较低无法充分利用GPU等加速计算资源。梯度的估计可能会更加不稳定因为每个批量中的样本可能并不代表整个数据集的分布情况这可能会导致训练过程中的震荡。 过大的批量大小训练过程中每个参数更新所使用的样本数目较多因此训练速度可能会变慢特别是在大规模数据集上。较大的批量可能需要较大的内存来存储计算过程中的梯度信息和中间结果这可能会导致内存不足或者性能下降。过大的批量可能会导致训练过程中陷入局部极小值因为每次参数更新所使用的样本数目较多可能会限制参数空间的搜索范围。
2在SGD算法的代码中我们采用了固定迭代次数的方式但是这样无法保证程序执行完毕时迭代已经收敛也有可能迭代早已收敛而程序还在运行。另一种方案是如果损失函数值连续 M M M次迭代都没有减小或者减小的量小于某个预设精度 ϵ \epsilon ϵ例如 1 0 − 6 10^{-6} 10−6就终止迭代。请实现该控制方案并思考它和固定迭代次数之间的利弊。能不能将这两种方案同时使用呢
增加一个控制代码
if abs(train_loss - prev_train_loss) tol:patience_count 1if patience_count patience:print(f达到终止条件迭代停止迭代次数{i1})breakelse:patience_count 0prev_train_loss train_losstol参数用于设置损失函数值的最小变化量patience参数用于设置连续迭代次数的容忍度。当连续patience次迭代中损失函数值变化量小于tol时即终止迭代。
混合策略先确保一定执行N个epoch后面再开始用控制代码
六、学习率对迭代的影响 在梯度下降算法中学习率是一个非常关键的参数。我们调整上面设置的学习率观察训练结果的变化。
_, loss1, _ SGD(num_epochnum_epoch, learning_rate0.1, batch_sizebatch_size)
_, loss2, _ SGD(num_epochnum_epoch, learning_rate0.001, batch_sizebatch_size)
plt.plot(np.arange(num_epoch), loss1, colorblue, labellr0.1)
plt.plot(np.arange(num_epoch), train_losses, colorred, ls--, labellr0.01)
plt.plot(np.arange(num_epoch), loss2, colorgreen,ls-., labellr0.001)
plt.xlabel(Epoch)
plt.ylabel(RMSE)
plt.gca().xaxis.set_major_locator(MaxNLocator(integerTrue))
plt.legend()
plt.show()可以看出随着学习率增大算法的收敛速度明显加快。那么学习率是不是越大越好呢我们将学习率继续上调到1.5观察结果。
_, loss3, _ SGD(num_epochnum_epoch, learning_rate1.5, batch_sizebatch_size)
print(最终损失, loss3[-1])
plt.plot(np.arange(num_epoch), np.log(loss3), colorblue, labellr1.5)
plt.xlabel(Epoch)
plt.ylabel(log RMSE)
plt.gca().xaxis.set_major_locator(MaxNLocator(integerTrue))
plt.legend()
plt.show()然而算法的RMSE随着迭代不但没有减小反而发散了。注意由于原始的损失数值过大我们将纵轴改为了损失的对数因此图中的直线对应于实际的指数增长。上图最后的损失约为 1 0 77 10^{77} 1077。为了进一步说明这一现象产生的原因我们用一个简单的例子来详细说明。假设我们要优化的目标函数是 J ( x ) x 2 J(x)x^2 J(x)x2那么梯度下降的迭代公式为 x ← x − η ∇ J ( x ) x − 2 η x ( 1 − 2 η ) x x\leftarrow x-\eta\nabla J(x)x-2\eta x(1-2\eta)x x←x−η∇J(x)x−2ηx(1−2η)x 显然该函数在 x 0 x0 x0 时可以取到唯一的全局最小值 J ( 0 ) 0 J(0)0 J(0)0。设 x x x的初始值为 x 0 x_0 x0经过 k k k次迭代后 x x x变为 x k ( 1 − 2 η ) k x 0 x_k(1-2\eta)^kx_0 xk(1−2η)kx0。考虑学习率的大小对迭代过程的影响。 当 0 η ≤ 0.5 0\eta\le0.5 0η≤0.5 时 0 ≤ 1 − 2 η 1 0\le 1-2\eta1 0≤1−2η1有 lim k → ∞ x k 0 \lim\limits_{k\rightarrow\infty}x_k0 k→∞limxk0。 x x x在迭代过程中始终与 x 0 x_0 x0符号相同且迭代过程可以收敛。 当 0.5 η 1 0.5\eta1 0.5η1 时 − 1 1 − 2 η 0 -11-2\eta0 −11−2η0同样有 lim k → ∞ x k 0 \lim\limits_{k\rightarrow\infty}x_k0 k→∞limxk0。但每次迭代 x x x都会改变符号迭代过程可以收敛。 当 η 1 \eta1 η1 时 1 − 2 η − 1 1-2\eta-1 1−2η−1迭代过程变为 x k 1 − x k x_{k1}-x_k xk1−xk。因此 x k x_k xk始终在 ± x 0 \pm x_0 ±x0间变化迭代不收敛。 当 η 1 \eta1 η1 时 1 − 2 η − 1 1-2\eta-1 1−2η−1有 lim k → ∞ x k ∞ \lim\limits_{k\rightarrow\infty}x_k\infty k→∞limxk∞每次迭代 x x x的符号会改变且向无穷大发散。 图6展示了在初始点 x 0 1 x_01 x01学习率 η \eta η分别为0.15、0.3、0.75、1.0和1.05时前3次迭代中 x x x的变化。可以看出学习率在一定范围内增大时可以加速算法收敛但学习率过大时算法也会出现不稳定、甚至发散的情况。因此梯度下降算法的学习率往往需要多次调整才能找到合适的值。 图6 不同学习率对梯度下降的影响 附以上文中的数据集及相关资源下载地址 链接https://pan.quark.cn/s/49692450efbe 提取码QE8x
