【机器学习基础】线性回归

在这里插入图片描述

【作者主页】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 xRd,那么该线性映射关系可以写为 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_0=1 x0=1 即可。这样上式的形式没有变化,只是向量的维度由 d d d变为 d + 1 d+1 d+1。为了表达的简洁,下面不再将常数项单独拆开表示和分析。图1展示了输入特征是一维和二维情况下的数据点和线性回归模型拟合的结果。在 d d d维输入特征和一维输出特征的情况下,线性回归模型共有 d + 1 d+1 d+1 个参数,从而给出了 d + 1 d+1 d+1 维空间中的一个 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_{i=1}^N\mathcal{L}(y_i,f_{\boldsymbol\theta}(\boldsymbol{x}_i)) J(θ)=N1i=1NL(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(yifθ(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_{i=1}^N(y_i-f_{\boldsymbol\theta}(\boldsymbol{x}_i))^2 J(θ)=2N1i=1N(yifθ(xi))2 这一总损失函数称为均方误差(mean squared error,MSE),是最常用的损失函数之一。因此,线性回归问题的优化目标为 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_{i=1}^N(y_i-f_{\boldsymbol\theta}(\boldsymbol{x}_i))^2 θminJ(θ)=θmin2N1i=1N(yifθ(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= x1TxNT   y= y1yN 其中, 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(yXθ)T(yXθ) 这里的损失函数事实上是所有样本平方误差的和,与MSE相差一步取平均值,即除以样本总数 N N N。但常数系数不影响优化得到的最终结果,为了形式简洁,我们通常在矩阵形式下省略这一系数。但是在实际计算时,为了降低样本规模的影响,绝大多数情况下还是以平均后的值作为实际损失。

  为了求函数 J ( θ ) J(\boldsymbol{\theta}) J(θ)的最小值,我们寻找其对 θ \boldsymbol\theta θ导数为零的点,即 ∂ J ∂ θ = 0 \frac{\partial J}{\partial\boldsymbol{\theta}}=\boldsymbol0 θJ=0 计算偏导数,得到
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(yTyyTXθ(Xθ)Ty+(Xθ)TXθ)=21(θyTyθyTXθθ(Xθ)Ty+θ(Xθ)TXθ)=XTy+XTXθ 通过上式得到 θ \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θ            θ=XTy+XTXθ=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.56
79248.646.006.733.0940173.071505890.91
61287.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=',', dtype='str')
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 error,RMSE)。对于真实值 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^N,RMSE为 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_{i=1}^N(y_i-\hat{y}_i)^2} LRMSE(y,y^)=N1i=1N(yiy^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, alpha=0.6, edgecolors='w', linewidth=0.5)
plt.xlabel('真实值')
plt.ylabel('预测值')
plt.title('真实值 vs 预测值')
plt.gca().set_aspect('equal', adjustable='box')
plt.plot([min(y_test), max(y_test)], [min(y_test), max(y_test)], 'r--', lw=2)# 预测误差分布
plt.subplot(1, 2, 2)
errors = y_test - y_pred
plt.hist(errors, bins=30, edgecolor='k', alpha=0.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} XRN×d y ∈ R N \boldsymbol{y}\in\mathbb{R}^{N} yRN。按照式 θ = ( 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^2+d^3) O(Nd2+d3)。虽然我们可以通过矩阵运算技巧进行优化,但时间开销仍然较大。此外,当样本很多时,存储矩阵 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_{i=1}^N(y_i-f_{\boldsymbol\theta}(\boldsymbol{x}_i))^2 J(θ)=2N1i=1N(yifθ(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_{i=1}^N(y_i-f_{\boldsymbol\theta}(\boldsymbol{x}_i))^2\right) \\[3ex] &=\boldsymbol\theta-\frac{\eta}{N}\sum_{i=1}^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_{i=1}^N(f_{\boldsymbol\theta}(\boldsymbol{x}_i)-y_i)\boldsymbol{x}_i \end{aligned} θθηθ(2N1i=1N(yifθ(xi))2)=θNηi=1N(fθ(xi)yi)θfθ(xi)=θNηi=1N(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(yXθ)T(yXθ))=θNη(XTy+XTXθ)=θ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)=θη(θTxkyk)xk 由于每次只计算一个样本,梯度下降的时间复杂度大大降低了,这样的算法就称为随机梯度下降法(stochastic gradient decent,SGD)。然而,随机梯度下降的不稳定性很高。由于单个样本计算出的梯度方向可能与所有样本算出的真正梯度方向不同,如果我们要优化的函数不是凸函数,SGD 算法就可能从原定路线偏离,收敛到其他极小值点去。因此,为了在稳定性与时间复杂度之间取得平衡,我们一般使用小批量梯度下降(mini-batch gradient decent,MBGD)算法,将样本随机分成许多大小较小的批量(batch)。每次迭代时,选取一个批量来计算函数梯度,以此估计用完整样本计算的结果。假如批量大小为 B ≪ N B\ll N BN,第 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 B=1 B=1 时就退化为SGD算法,当 B = N B=N B=N 时就退化为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, shuffle=True):# 批量计数器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(size=X.shape[1])# 随机梯度下降# 为了观察迭代过程,我们记录每一次迭代后在训练集和测试集上的均方根误差train_losses = []test_losses = []for i in range(num_epoch):# 初始化批量生成器batch_g = batch_generator(X, y_train, batch_size, shuffle=True)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, color='blue', label='train loss')
plt.plot(np.arange(num_epoch), test_losses, color='red', ls='--', label='test loss')
# 由于epoch是整数,这里把图中的横坐标也设置为整数
# 该步骤也可以省略
plt.gca().xaxis.set_major_locator(MaxNLocator(integer=True))
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} 106),就终止迭代。请实现该控制方案,并思考它和固定迭代次数之间的利弊。能不能将这两种方案同时使用呢?

增加一个控制代码:

if abs(train_loss - prev_train_loss) < tol:patience_count += 1if patience_count >= patience:print(f'达到终止条件,迭代停止,迭代次数:{i+1}')breakelse:patience_count = 0prev_train_loss = train_loss

tol参数用于设置损失函数值的最小变化量,patience参数用于设置连续迭代次数的容忍度。当连续patience次迭代中损失函数值变化量小于tol时,即终止迭代。

混合策略:先确保一定执行N个epoch,后面再开始用控制代码

六、学习率对迭代的影响

  在梯度下降算法中,学习率是一个非常关键的参数。我们调整上面设置的学习率,观察训练结果的变化。

_, loss1, _ = SGD(num_epoch=num_epoch, learning_rate=0.1, batch_size=batch_size)
_, loss2, _ = SGD(num_epoch=num_epoch, learning_rate=0.001, batch_size=batch_size)
plt.plot(np.arange(num_epoch), loss1, color='blue', label='lr=0.1')
plt.plot(np.arange(num_epoch), train_losses, color='red', ls='--', label='lr=0.01')
plt.plot(np.arange(num_epoch), loss2, color='green',ls='-.', label='lr=0.001')
plt.xlabel('Epoch')
plt.ylabel('RMSE')
plt.gca().xaxis.set_major_locator(MaxNLocator(integer=True))
plt.legend()
plt.show()

在这里插入图片描述

  可以看出,随着学习率增大,算法的收敛速度明显加快。那么,学习率是不是越大越好呢?我们将学习率继续上调到1.5观察结果。

_, loss3, _ = SGD(num_epoch=num_epoch, learning_rate=1.5, batch_size=batch_size)
print('最终损失:', loss3[-1])
plt.plot(np.arange(num_epoch), np.log(loss3), color='blue', label='lr=1.5')
plt.xlabel('Epoch')
plt.ylabel('log RMSE')
plt.gca().xaxis.set_major_locator(MaxNLocator(integer=True))
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 xxηJ(x)=x2ηx=(12η)x 显然,该函数在 x = 0 x=0 x=0 时可以取到唯一的全局最小值 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=(12η)kx0。考虑学习率的大小对迭代过程的影响。

  • 0 < η ≤ 0.5 0<\eta\le0.5 0<η0.5 时, 0 ≤ 1 − 2 η < 1 0\le 1-2\eta<1 012η<1,有 lim ⁡ k → + ∞ x k = 0 \lim\limits_{k\rightarrow+\infty}x_k=0 k+limxk=0 x x x在迭代过程中始终与 x 0 x_0 x0符号相同,且迭代过程可以收敛。

  • 0.5 < η < 1 0.5<\eta<1 0.5<η<1 时, − 1 < 1 − 2 η < 0 -1<1-2\eta<0 1<12η<0,同样有 lim ⁡ k → + ∞ x k = 0 \lim\limits_{k\rightarrow+\infty}x_k=0 k+limxk=0。但每次迭代 x x x都会改变符号,迭代过程可以收敛。

  • η = 1 \eta=1 η=1 时, 1 − 2 η = − 1 1-2\eta=-1 12η=1,迭代过程变为 x k + 1 = − x k x_{k+1}=-x_k xk+1=xk。因此 x k x_k xk始终在 ± x 0 \pm x_0 ±x0间变化,迭代不收敛。

  • η > 1 \eta>1 η>1 时, 1 − 2 η < − 1 1-2\eta<-1 12η<1,有 lim ⁡ k → + ∞ x k = ∞ \lim\limits_{k\rightarrow+\infty}x_k=\infty k+limxk=,每次迭代 x x x的符号会改变,且向无穷大发散。

  图6展示了在初始点 x 0 = 1 x_0=1 x0=1,学习率 η \eta η分别为0.15、0.3、0.75、1.0和1.05时,前3次迭代中 x x x的变化。可以看出,学习率在一定范围内增大时可以加速算法收敛,但学习率过大时,算法也会出现不稳定、甚至发散的情况。因此,梯度下降算法的学习率往往需要多次调整,才能找到合适的值。

在这里插入图片描述

图6 不同学习率对梯度下降的影响

:以上文中的数据集及相关资源下载地址:
链接:https://pan.quark.cn/s/49692450efbe
提取码:QE8x

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/395996.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

集成视触觉传感器的机器人操作学习

强化学习是一种仿人学习的方法&#xff0c;其在不断与环境交互试错的过程中进行学习&#xff0c;提高自身的认知。其具有如下的优点&#xff0c;首先是数据依赖性低&#xff0c;强化学习通过与环境的交互来学习&#xff0c;减少了对标记数据的依赖性&#xff0c;可以大量的减少…

Linux 系统框架分析(一)

一、linux内核结构框图 对内核结构框图有个总体的把握&#xff0c;有助于理解为什么驱动要这样写&#xff0c;为什么写的应用程序所用的C库接口能够产生这么多的事情。 框图可以看出来&#xff0c;linux系统&#xff0c;包括五个系统 一、Linux内核结构介绍 Linux 内核是操作…

Spring及相关框架的重要的问题

Java框架 问题一&#xff1a;Spring框架中的单例bean是线程安全的吗&#xff1f; 看下图&#xff0c;不能被修改的成员变量就是无状态的类&#xff0c;无状态的类没有线程安全问题&#xff0c;所以在开发中尽量避免可修改的成员变量。 回答&#xff1a;不是线程安全的&#xf…

Oracle一对多(一主多备)的DG环境如何进行switchover切换?

本文主要分享Oracle一对多(一主多备)的DG环境的switchover切换&#xff0c;如何进行主从切换&#xff0c;切换后怎么恢复正常同步&#xff1f; 1、环境说明 本文的环境为一主两备&#xff0c;数据库版本为11.2.0.4&#xff0c;主要信息如下&#xff1a; 数据库IPdb_unique_n…

Github 2024-08-09 开源项目日报 Top10

根据Github Trendings的统计,今日(2024-08-09统计)共有10个项目上榜。根据开发语言中项目的数量,汇总情况如下: 开发语言项目数量Python项目6TypeScript项目4Jupyter Notebook项目1Cuda项目1Sentry:开发者优先的错误跟踪和性能监控平台 创建周期:5093 天开发语言:Python,…

android系统中data下的xml乱码无法查看问题剖析及解决方法

背景&#xff1a; Android12高版本以后系统生成的很多data路径下的xml都变成了二进制类型&#xff0c;根本没办法看xml的内容具体如下&#xff1a; 比如想要看当前系统的widget的相关数据 ./system/users/0/appwidgets.xml 以前老版本都是可以直接看的&#xff0c;这些syste…

旅游出行必备商城小程序的设计

管理员账户功能包括&#xff1a;系统首页&#xff0c;个人中心&#xff0c;用户管理&#xff0c;新闻类型管理&#xff0c;新闻资讯管理&#xff0c;商品类型管理&#xff0c;旅游商品管理&#xff0c;旅游景点&#xff0c;景点分类&#xff0c;系统管理 微信端账号功能包括&am…

GitHub的常用操作

目录 GitHub GitHub加速 克隆GitHub上的项目到本地 克隆GitHub上指定分支的项目 把本地项目上传到GitHub上管理 删除分支里的内容 单个仓库管理多个项目 上传项目到新建的分支 目前正在逐步熟悉GitHub&#xff0c;打算把整理好的代码上传到GitHub上&#xff0c;建立属…

C++ 类与对象

面向对象程序设计基本特点 特点&#xff1a; 抽象&#xff08;数据抽象&#xff0c;行为抽象&#xff09; 数据抽象&#xff1a;int hour,int minute.....,车&#xff1a;长&#xff0c;宽&#xff0c;高.... 功能抽象&#xff1a;showTime(),setTime() .....车&#xff1a;刹车…

使用Cisco进行模拟配置OSPF路由协议

OSPF路由协议 1.实验目的 1&#xff09;理解OSPF 2&#xff09;掌握OSPF的配置方法 3&#xff09;掌握查看OSPF的相关信息 2.实验流程 开始 → 布置拓扑 → 配置IP地址 → 配置OSPF路由并验证PC路由的连通性 → 查看路由器路由信息 → 查看路由协议配置与统计信息 → 查看O…

【从零开始一步步学习VSOA开发】VSOA命令行工具vcx

VSOA命令行工具vcx vcx 介绍 vcx 是一个使用 VSOA RPC 客户端功能执行器&#xff0c;支持 RPC SET/GET 调用。 [rootsylixos:/root]# [rootsylixos:/root]# vcx -help USAGE: vcx [options] url -h : Show help message. -v : Show vcx version. -z …

[MRCTF2020]PYWebsite-1

打开以后查看源码信息 看到flag.php试着打开 提示看到&#xff0c;需要后端审计代码&#xff0c;而且应该要改ip&#xff0c;改成自己本地&#xff0c;burp抓包看一下 改X-Forwarded-For:127.0.0.1 得到flag flag{74242eb7-844f-4638-8aae-9ec37870d585}

通过LLM大模型将「白雪公主的故事」转为图数据存储

&#x1f4a1; 本次将使用LLM大模型将「白雪公主的故事」转为图数据存储于neo4j数据库中&#xff0c;并展示图数据部分的效果 故事内容 很久很久以前&#xff0c;在一个遥远的王国里&#xff0c;有一位美丽的**王后**生下了一个皮肤像雪一样白皙、嘴唇像血一样鲜红的**女儿**…

网页设计模板范例

随着互联网的发展&#xff0c;网页设计变得越来越重要。一个吸引人的网页设计可以吸引更多的用户&#xff0c;提升用户体验&#xff0c;并且使网站内容更加易于浏览和理解。在这篇文章中&#xff0c;我将为大家介绍一个网页设计模板范例。 1. 选择合适的颜色和字体&#xff1a;…

golang for range time.Ticker 和 time.Timer时间通道使用示例 - 每隔指定时间执行一次,执行指定时长后退出执行

golang中的 ticker和timer时间通道除了可以使用for select case语句来执行外&#xff0c; 还可以使用 for range语句来执行ticker或者timer时间通道。 for range time.Ticker 和 time.Timer时间通道使用示例 下面的示例演示了time.Ticker 和 time.Timer的区别和使用演示。 Ti…

JAVA分布式CAP原则

分布式CAP原则主要是使用SpringCloud框架的时候会涉及到该部分知识点 CAP原则指的是在一个分布式系统中&#xff0c;一致性&#xff0c;可用性&#xff0c;分区容错性 实际的项目开发中&#xff0c;这三者往往是不可兼顾的。 AP&#xff1a;牺牲一致性&#xff0c;保证可用性…

【学习笔记】Matlab和python双语言的学习(动态规划)

文章目录 前言一、动态规划动态规划的基本步骤示例1示例2 三、代码实现----Matlab1.示例12.示例2 四、代码实现----python1.示例12.示例2 总结 前言 通过模型算法&#xff0c;熟练对Matlab和python的应用。 学习视频链接&#xff1a; https://www.bilibili.com/video/BV1EK411…

Linux驱动开发—平台总线模型详解

文章目录 1.平台总线介绍1.1平台总线模型的组成部分1.2平台总线模型的优势 2.使用平台总线模型开发驱动2.1注册platform设备2.2注册platform驱动2.3效果演示 1.平台总线介绍 Linux 平台总线模型&#xff08;Platform Bus Model&#xff09;是一种设备驱动框架&#xff0c;用于…

Debezium报错处理系列之第116篇:Caused by: java.lang.NumberFormatException: null

Debezium报错处理系列之第116篇:Caused by: java.lang.NumberFormatException: null 一、完整报错二、错误原因三、解决方法Debezium从入门到精通系列之:研究Debezium技术遇到的各种错误解决方法汇总: Debezium从入门到精通系列之:百篇系列文章汇总之研究Debezium技术遇到的…

【从零开始一步步学习VSOA开发】VSOA命令行工具vcl

VSOA命令行工具vcl 介绍 vcl 是一个 VSOA 命令监听器&#xff0c;支持订阅来自特定服务器发布的所有数据。 [rootsylixos:/root]# [rootsylixos:/root]# vcl -help USAGE: vcl [options] host [topics] -h : Show help message. -v : Show vcl version. …