一文看懂线性回归（保姆级别含源码）

来源: AINLPer 微信公众号（每日论文干货分享！！）
编辑: ShuYini
校稿: ShuYini
时间: 2022- 9 -18

引言

最近作者网上看了很多关于线性回归的帖子，个人感觉比较乱，所以打算自己整理一版，希望能把线性回归说的更明白一些，为此我还整理了关于线性回归的脑图如下图所示，如果哪里有不对的地方，欢迎批评指正。

本文涉及源码下载：关注: AINLPer 回复：线性回归源码

线性回归

线性回归（英语：linear regression）是利用称为线性回归方程的最小二乘函数对一个或多个自变量和因变量之间关系进行建模的一种回归分析。这种函数是一个或多个称为回归系数的模型参数的线性组合。只有一个自变量的情况称为简单回归，大于一个自变量情况的叫做多元回归（multivariable linear regression）。【来自维基】

线性回归的应用

线性回归有很多实际用途。分为以下两大类：

1、如果目标是预测或者映射，线性回归可以用来对观测数据集的和X的值拟合出一个预测模型。当完成这样一个模型以后，对于一个新增的X值，在没有给定与它相配对的y的情况下，可以用这个拟合过的模型预测出一个y值。

2、给定一个变量y和一些变量 $X_1,...,X_p$ ，这些变量有可能与y相关，线性回归分析可以用来量化y与 $X_j$ 之间相关性的强度，评估出与y不相关的 $X_j$ ，并识别出哪些 $X_j$ 的子集包含了关于y的冗余信息。

我们大多数情况下是使用第一种方式做预测，即通过已知的x，y来预测模型，进而预测未来给x，输出y的值。

机器学习中的线性回归分类

1、单变量线性回归

2、多变量线性回归

3、多项式线性回归

线性回归求解方法

1、最小二乘法（x矩阵必须可逆）

$θ=(xT⋅x)−1⋅xT⋅y\theta = (x^ {T}\cdot x)^ {-1}\cdot x^{T} \cdot y$

详情可以参考：https://blog.csdn.net/u014445499/article/details/107209837/

2、梯度下降法

对于一个函数y=f(x)，我们要求解该函数的极值点，高中的时候就学过，通过解方程f′(x)=0，就可以得到函数的极值点(x0,y0)。不过对于计算机来说，它可不会解方程。但是它可以凭借强大的计算能力，一步一步的去把函数的极值点『试』出来。

首先要知道需要优化的目标函数为：
$J(θ)=12m∑i=1m(hθx(i)−y(i))2J(\theta)= \frac{1}{2m} \sum_{i=1}^{m}(h_{\theta}x^{(i)}-y^{(i)})^{2}$

其中：
$hθh_ {\theta }$ (x)= $θT\theta ^ {T}$ x= $θ0\theta _ {0}$ + $θ1\theta _ {1}$ $x_ {1}$ ，对 $J(θ)J(\theta)$ 求导数得到(该公式对于多项式函数同样适用)：

$∂∂θjJ(θ)=2m∑i=1m(θT⋅x(i)−y(i))xj(i)=2mXT⋅(X⋅θ−y)\frac{\partial}{\partial\theta_{j}}J(\theta)=\frac{2}{m}\sum_{i=1}^{m}(\theta^{T}\cdot x^{(i)}-y^{(i)})x_{j}^{(i)}=\frac{2}{m}X^{T}\cdot(X\cdot\theta-y)$

最后，梯度更新如下：

$θj=θj−α∂∂θjJ(θ)\theta_{j}=\theta_{j}-\alpha\frac{\partial}{\partial\theta_{j}}J(\theta)$

下图就是一个最经典的梯度下降示意图。

梯度下降法具体介绍可以参考：https://blog.csdn.net/u013898698/article/details/120720623

2.1、梯度下降分类

很多视频、文章都会提到，根据数据迭代的计算方式对梯度下降分类，主要可分为以下3类：

2.1.1、全批量梯度下降（BGD）

1、计算量大，参数更新慢，对内存的要求很高，不能以在线的形式训练模型，也就是运行时不能加入新样本2、可以得到全局最优解，参数更新比较稳定，收敛方向稳定。

2.1.2、随机梯度下降（SGD）：

1、运算速度很快，同时能够在线学习。2、随机梯度下降参数更新的过程震荡很大，目标函数波动剧烈，参数更新方向有很大的波动。3、其较大的波动可能收敛到比批量梯度下降更小的局部极小值，因为会从一个极小值跳出来。

2.1.3、小批量梯度下降（MBGD）：

1、降低了更新参数的方差，使得收敛过程更加的稳定2、能够利用高度优化的矩阵运算，很高效的求得每小批数据的梯度

梯度下降法和最小二乘法求解对比

梯度下降法：

– 缺点

---- 1、需要选择学习率 alpha

---- 2、需要多次的迭代

---- 3、特征值取值范围比较大需要进行特征缩放（归一化）

– 优点

---- 1、当特征数n很大的时候，效果比较好

最小二乘法

– 缺点

---- 1、不需要选择学习率 alpha

---- 2、没有多次的迭代

---- 3、没有特征缩放（归一化）

– 优点

---- 1、当特征数n很大的时候，运算的比较慢，求逆矩阵的复杂度为: $n^3$

两种算法选择：

在深度学习中，因为涉及参数特别多，动不动涉及的参数高达几十几百万，所以基本上全部都是通过梯度下降的方法求解。后面会有这种方法的代码实现，其实这块主要帮助大家加深对梯度下降的理解，在实际项目中，主要还是通过mini-batch的思想来更新梯度。以此衍生出来的方法也有很多，例如：SGD、SGDM（SGD+Momentum）、Adagrad、RMSProp、Adma(2015年)、swats（2017）、AMSGrad(2018)、AdaBound(2019)、Cycle-LR、OneCycle-LR，目前常用的主要有：SGDM、Adma、SWATS。这里就不展开来讨论了，后面会专门出一期方法对比的文章。

线性回归代码演示

单变量线性回归

# 【数据生成】根据函数y=4.7x+2 来构造的x、y数据，并画出来
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
x=2*np.random.rand(100,1)
y=2+4.7*x+np.random.rand(100,1)
plt.plot(x,y,'b.')
plt.xlabel('x')
plt.ylabel('y')
plt.axis([0,2,0,15])
plt.show()

在这里插入图片描述

# 【最小二乘法求解】通过回归公式来计算权重和偏置项theta，并画出求解得到的值
x_b=np.c_[np.ones((100,1)),x]
theta=np.linalg.inv(x_b.T.dot(x_b)).dot(x_b.T).dot(y)# 两个点画一条线
xx=np.array([[0],[2]])
xx_tmp=np.c_[np.ones((2,1)),xx]
yy=xx_tmp.dot(theta)
print('方程直接求解w、b的值分别为{0},{1}'.format(theta[1],theta[0]))
plt.plot(xx,yy,'r--')
plt.plot(x,y,'b.')
plt.xlabel('x')
plt.ylabel('y')
plt.axis([0,2,0,15])
plt.show()

在这里插入图片描述

# 【批量梯度下降求解】用梯度下降的方式求解出y、x之间的关系,用所有的样本进行训练
# 从下面结果可以看到，计算出来的theta值和上面直接计算出来的值是差不多的。
m=100
lr=0.05  # 学习率
iterate=1000
theta=np.random.rand(2,1)
x_b=np.c_[np.ones((100,1)),x]
for i in range(iterate):gradent=2/m*x_b.T.dot(x_b.dot(theta)-y)theta=theta-lr*gradent
print("theta:",theta)

【批量梯度下降（BGD）】主要是在整个数据集，通过对所有的样本的计算来求解梯度的方向。
【学习率对梯度下降的影响】主要观察不同学习率对实际训练结果的影响，这里会设置的学习率有四个：0.01、0.05、0.1、0.5

theta_for_bgd=[]
# 批量梯度求解函数
def diff_lr_for_gradent(iterations,lr):j_save=[]x_b=np.c_[np.ones((100,1)),x]m=len(x_b)theta=np.random.rand(2,1)plt.plot(x,y,'b.')for iteration in range(iterations):y_pred=x_b.dot(theta)# 计算损失函数并存放到[]里面J=1/(2*m)*(y-y_pred).T.dot(y-y_pred)j_save.append(float(J[0]))   # 画出拟合直线plt.plot(x,y_pred,'r-')# 梯度计算及更新gradent=2/m*x_b.T.dot(x_b.dot(theta)-y)theta=theta-lr*gradentif lr==0.05:theta_for_bgd.append(theta)# 画出损失函数和训练次数的图像 plt.xlabel('X')plt.ylabel("Y")plt.axis([0,2,0,15])plt.title('lr={}'.format(lr)) return  j_save

根据下方画图（红色线是拟合线、绿色线表示了迭代次数和损失函数的关系）：

可以看出来lr越小，其拟合的直线向中间逼近的比较慢（损失函数收敛慢），随着lr的增加，逼近的速度也越快（损失函数收敛快）。

但是当lr很大的时候，反而不能趋向实际要拟合的直线（损失函数会变得发散），这说明了在训练过程中要选择合适的lr。

iterations=100
plt.figure(figsize=(15,4))plt.subplot(241)
j_save=diff_lr_for_gradent(iterations,0.01)
tmp=np.arange(iterations)
plt.subplot(245)
plt.plot(tmp,j_save,'g-')plt.subplot(242)
j_save=diff_lr_for_gradent(iterations,0.05)
tmp=np.arange(iterations)
plt.subplot(246)
plt.plot(tmp,j_save,'g-')plt.subplot(243)
j_save=diff_lr_for_gradent(iterations,0.1)
tmp=np.arange(iterations)
plt.subplot(247)
plt.plot(tmp,j_save,'g-')plt.subplot(244)
j_save=diff_lr_for_gradent(iterations,0.5)
tmp=np.arange(iterations)
plt.subplot(248)
plt.plot(tmp,j_save,'g-')
plt.show()

在这里插入图片描述
【随机梯度下降（SGD）】随机是指每个batch的数据是对整个数据空间的随机采样，只要采样够随机，随机梯度的期望和真实梯度是相等的。放在模型训练里，哪怕一个batch就一个数据点，只要多迭代几个batch，最终效果一定是一样的。

【自适应学习率】学习率我们希望是一开始的时候步长比较长，后来慢慢越走越小，进而达到损失函数收敛的最佳参数。

Tip：这里会有两个超参数：epoch、batch

Epoch：定义了学习算法在整个训练数据集中的工作次数，Batch：定义在更新内部模型参数之前要处理的样本数。

详情可以参考：https://blog.csdn.net/weixin_42137700/article/details/84302045

theta_for_sgd=[]
x_b=np.c_[np.ones((100,1)),x]
m=len(x_b)
epoches=50 #训练次数
theta=np.random.rand(2,1)
#batch=m # 假设batch的大小和训练数据的大小一样# lr学习率更新
def lr_update(iterations):return 5/(50+iterations)for epoch in range(epoches):for i in range(m):# 在训练数据中随机的选择x，ydata_index=np.random.randint(m)xi=x_b[data_index:data_index+1]yi=y[data_index:data_index+1]# 依据该条数据更新梯度gradent=2*xi.T.dot(xi.dot(theta)-yi)lr=lr_update(epoch*m+i)theta=theta-lr*gradenttheta_for_sgd.append(theta)# 画出部分图用来展示if epoch<5 and i< 10 :y_pred=x_b.dot(theta)plt.plot(x,y_pred,'r-')plt.plot(x,y,'b.')
plt.xlabel('X')
plt.ylabel('Y')
plt.axis([0,2,0,15])
plt.title('SGD')
plt.show()

在这里插入图片描述
【小批量梯度下降（MBGD）】把数据分为若干个批，按批来更新参数，利用每一批的数据来计算梯度的方向，下降起来就不容易跑偏，减少了随机性。

在实际操作的时候，需要重新对训练数据做索引打乱操作（shuffle）。

下面是代码展示：

theta_for_mbgd=[]
epoches=50  # 50次迭代训练
batch=5 # 每次训练批次为5
m=len(x_b)
# lr学习率更新函数
def lr_update(t):return 5/(50+t)
t=0
# 初始化theta
theta=np.random.rand(2,1)for epoch in range(epoches):# 对训练数据打乱--洗牌shuffle_indexes=np.random.permutation(m) x_b_shuffle=x_b[shuffle_indexes]y_shuffle=y[shuffle_indexes]for i in range(0,m,batch):t=t+1     x_b_tmp=x_b_shuffle[i:i+batch]y_tmp=y_shuffle[i:i+batch]# 计算梯度gradent=2/batch*x_b_tmp.T.dot(x_b_tmp.dot(theta)-y_tmp)theta=theta-lr_update(t)*gradenttheta_for_mbgd.append(theta)# 画出部分图像if epoch<5 and t< 10 :y_pred=x_b.dot(theta)plt.plot(x,y_pred,'r-')plt.plot(x,y,'b.')plt.axis([0,2,0,15])
plt.xlabel('X')
plt.ylabel('Y')
plt.title('MBGD')
plt.show()

在这里插入图片描述
【三种梯度下降方法对比】下面画出了theta在三种梯度下降方法中的变化。如下图所示：
BGD红色曲线，沿着曲线直接到达终点，没有出现起伏波动；
SGD蓝色曲线，最终达到终点，但是起伏比较大；
MBGD绿色曲线，达到终点，有起伏，但是起伏不大。

theta_for_bgd=np.array(theta_for_bgd)
theta_for_sgd=np.array(theta_for_sgd)
theta_for_mbgd=np.array(theta_for_mbgd)
plt.figure(figsize=(8,4))
plt.subplot(1,3,1)
plt.plot(theta_for_bgd[:,0],theta_for_bgd[:,1],'r-s',linewidth=1,label='bgd')
plt.legend(loc='upper left')
plt.subplot(1,3,2)
plt.plot(theta_for_sgd[:,0],theta_for_sgd[:,1],'b-*',linewidth=1,label='sgd')
plt.legend(loc='upper left')
plt.subplot(1,3,3)
plt.plot(theta_for_mbgd[:,0],theta_for_mbgd[:,1],'g-+',linewidth=1,label='mbgd')
plt.legend(loc='upper left')plt.show()

在这里插入图片描述
三种梯度算法总结：

全批量梯度下降（BGD）

1、计算量大，参数更新慢，对内存的要求很高，不能以在线的形式训练模型，也就是运行时不能加入新样本2、可以得到全局最优解，参数更新比较稳定，收敛方向稳定。

随机梯度下降（SGD）：

1、运算速度很快，同时能够在线学习。2、随机梯度下降参数更新的过程震荡很大，目标函数波动剧烈，参数更新方向有很大的波动。3、其较大的波动可能收敛到比批量梯度下降更小的局部极小值，因为会从一个极小值跳出来。

小批量梯度下降（MBGD）：

1、降低了更新参数的方差，使得收敛过程更加的稳定2、能够利用高度优化的矩阵运算，很高效的求得每小批数据的梯度

多项式线性回归

之前的时候，我们采用的直线方程：y=4.7x+2然后加上随机变量生成测试数据，这里，我们首先通过y=2x^2+x+5来生成测试数据。生成的数据图像如下图所示：

tip：实际应用要处理的数据可能还要离散、无规律，这里主要是用来解释多项式回归问题。

# 生成y=3x^3+2x^2+x+5数据测试数据
m=100 # x的维度
x_2=4*np.random.rand(m,1)
y_2=2*x_2**2+x_2+5+2*np.random.randn(m,1) # 后面添加了一个正态分布
plt.plot(x_2,y_2,'r.')
plt.show()

在这里插入图片描述

# 按照上面数据，我们还是通过一元线性函数进行求解如下:
epoches=20
batch_size=8
x_2_b=np.c_[np.ones((100,1)),x_2]
# print(x_2_b)
theta=np.random.rand(2,1)for epoch in range(epoches):# 数据shuffleshuffle_indexes=np.random.permutation(m) x_b_shuffle=x_2_b[shuffle_indexes]y_shuffle=y_2[shuffle_indexes]for i in range(0,m,batch_size):x_b_tmp=x_b_shuffle[i:i+batch_size]y_tmp=y_shuffle[i:i+batch_size]# 计算梯度gradent=2/batch_size*x_b_tmp.T.dot(x_b_tmp.dot(theta)-y_tmp)theta=theta-lr_update(t)*gradentplt.plot(x_2,y_2,'r.')
y_pred=x_2_b.dot(theta)
plt.plot(x_2,y_pred,'b*')
plt.show()

在这里插入图片描述

# 通过上图可以发现，通过直线的方法并不能够很好的拟合出曲线，这个时候就要选择多元函数进行拟合，这就是前面介绍的多项式回归
# 即利用y=ax^2+bx+c进行拟合：
epoches=500
batch_size=16
x_2_c=np.c_[np.ones((100,1)),x_2,x_2**2] # 这里可以通过添加x_2**3\x_2**4\....来计算不同degree的结果
# print(x_2_b)
theta=np.random.rand(3,1)for epoch in range(epoches):# 数据shuffleshuffle_indexes=np.random.permutation(m) x_b_shuffle=x_2_c[shuffle_indexes]y_shuffle=y_2[shuffle_indexes]for i in range(0,m,batch_size):x_b_tmp=x_b_shuffle[i:i+batch_size]y_tmp=y_shuffle[i:i+batch_size]# 计算梯度gradent=2/batch_size*x_b_tmp.T.dot(x_b_tmp.dot(theta)-y_tmp)theta=theta-lr_update(t)*gradentplt.plot(x_2,y_2,'r.')
x2y_pred=x_2_c.dot(theta)
# print(xy_pred)
plt.plot(x_2,x2y_pred,'b.')
plt.show()

在这里插入图片描述

多变量线性回归

当一个函数不仅仅和一参数有关系的时候（例如：吴恩达教程里面房价可能和房子面积、房子建设时间、房子位置等；李宏毅教程里面的神奇宝贝cp值可能和神奇宝贝的属性值、成长系数、防御值等），就涉及到了多项式回归。
这里以函数z=3x+6y+2创造数据集：

# 依据z=3x+6y+10加上随机数字来生成训练数据（x,y,z），后面以该组数据作为训练数据进行训练，让训练结果无线逼近z=3x+6y+10
# 画出z=3x+6y+10的视图 scatter  contour3D  scatter3D  plot_surface  plot_wireframe
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
m=50
x=np.linspace(1,20,m)
y=np.linspace(1,20,m)
x,y=np.meshgrid(x,y)
x1=x.reshape(-1,1)
y1=y.reshape(-1,1)
z1=3*x1+6*y1+2+5*np.random.randn(m*m,1)
z=z1.reshape(m,m) # 此时构造的训练数据集为：x1(m*m,1),x2(m*m,1),z1((m*m,1))# 设置图像大小
plt.figure(dpi=100)ax=plt.axes(projection='3d')
ax.plot_surface(x,y,z)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()

在这里插入图片描述
依据上面生成的数据，通过梯度下降的方法训练得到的结果为：y=3.03x+6.01y+1.48,会发现这和y=3x+6y+2是非常接近的，另外下面两张图分别也展示了实际生成的平面、以及随着迭代次数增加loss的变化曲线。在实际应用中，多变量线性回归可能不知两个特征x，y，可能还有更多的值，但是都可以通过类似的方法训练。

# 根据提取下降来计算theta,
epoches=100
batch_size=256
theta_recode=[]
theta=np.random.rand(3,1)
x_y_b=np.c_[np.ones((m*m,1)),x1,y1]
t=1
y_loss=[]# lr学习率更新函数，学习率一开始选择的比较大，后面逐渐变小。
def lr_update(t):return 5/(1000+t)
# 损失函数
def costFunction(X,y,theta):inner =np.power( X @ theta - y, 2)return np.sum(inner) / (2 * len(X))
for epoch in range(epoches):shuffle_indexes=np.random.permutation(m*m) x_y_b_shuffle=x_y_b[shuffle_indexes]z_shuffle=z1[shuffle_indexes]for i in range(0,m*m,batch_size):x_y_b_temp=x_y_b_shuffle[i:i+batch_size]z_temp=z_shuffle[i:i+batch_size]gradent=1/batch_size*x_y_b_temp.T.dot(x_y_b_temp.dot(theta)-z_temp)theta=theta-lr_update(t)*gradent t+=1loss=costFunction(x_y_b,z1,theta)y_loss.append(loss)
#     print('this is the {} epoch, the loss is {}'.format(epoch,loss))print('theta={}'.format(theta))
# plot_function(theta)
plt.figure(dpi=150)
ax=plt.subplot(1,2,1,projection='3d')
# ax=plt.axes(projection='3d')
z_end=theta[1]*x+theta[2]*y+theta[0]
# 画图对比，之前的图和现在图对比
ax.plot_surface(x,y,z)
ax.plot_surface(x,y,z_end,color='r')plt.subplot(1,2,2)
x_loss=np.linspace(1,epoches,epoches)
plt.plot(x_loss,y_loss,'g-')
plt.subplots_adjust(wspace =0.5, hspace =0)
plt.show()