[机器学习] - 岭回归与Lasso回归

注：在学习机器学习中回归算法时，随时都会接触最小二乘法原理和求解过程，最小二乘法可能对于理工科的学生低头不见抬头见的知识点，特点是在我学习《数值分析》课程中，老师讲了些最小二乘法的历史[…吧啦吧啦啥高斯、啥勒让德…]听的很有趣，感觉自己啥都懂，实际上上手编程又不会。借着对机器学习的回归模型的学习，提前介绍下我对它情有独钟的理解，以及它求解过程。

最开始使用最小二乘法是在数学建模中建立的时间序列预测模型，几乎对每一个多项式确定系数的求解问题，都有它的影子。同时最开始接触到了岭回归和Lasso回归，但是一直也无法对其基本原理有一个透彻、直观的理解，也不知道它归属于机器学习中的范畴。直到最近再次接触到这个概念，经过一番苦思冥想后终于有了我自己的理解。

1 普通最小二乘法

百度一下： 最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。
用一句话概括：最小化损失函数来确定系数。

最小二乘法的理论：
　　在我们研究特征/输入变量xix_ixi与输出变量y之间的相互关系时，通常可以得到一系列成对的一组数据{x1x_1x1，x2x_2x2，… xnx_nxn，y}；

如果输入特征只有一个x时，也就是直线方程：y=w0+w1xy={w_0}+{w_1}xy=w0+w1x
将这些数据描绘在x -y直角坐标系中，若发现这些点在一条直线附近，可以令这条直线方程如：
　　　　　　　　　　　　　　　　　　　　　　　 y=w0+w1xy = {w_0} + {w_1}xy=w0+w1x
其中：w0{w_0}w0、w1{w_1}w1是任意实数。
　　为建立这直线方程就要确定w0{w_0}w0、w1{w_1}w1系数，将实测值YiY_iYi与利用计算值YjY_jYj（式1-1）的误差平方和函数：
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　∑(Yi−Yj)2{\sum {\left( {{Y_i} - {Y_j}} \right)} ^2}∑(Yi−Yj)2
的最小值为“优化判据”，有时候会看到使用L2L_2L2范数，即：S(w)=∥y−wx∥2S\left( w \right) = {\left\| {y - wx} \right\|^2}S(w)=∥y−wx∥2，其中x、y都为一组向量。

往往特征不只一个，即xix_ixi=(x1x_1x1，x2x_2x2，… xnx_nxn)，有n种特征。
考虑超定方程组（超定指未知数大于方程个数）：
　　　　　　　　　　　　　　　　　　　　　　　y^=w0+w1x1+w2x2+...+wnxn=X⋅w\hat y = {w_0} + {w_1}{x_1} + {w_2}{x_2} + ... + {w_n}{x_n} = X \cdot wy^=w0+w1x1+w2x2+...+wnxn=X⋅w （1-1）

y^\hat yy^表示线性回归模型的预测值（相对于真实观察值)
nnn表示特征的数量
xix_ixi表示第i个特征的观察值
wjw_jwj表示第j个参数的值

其中m代表有m个等式，n代表有 n 个未知数www，m>n ；将其进行向量化后为：
　　　　　　　　　　　　　　　　　　　　　　　　　　　Ｘw=y^Ｘw=\hat yＸw=y^

显然该方程组一般而言没有解，所以为了选取最合适的解www，让该等式"尽量成立"，引入残差平方和函数SSS：
　　　　　　　　　　　　　　　　　　　　　　　　S(w)=∥Xw−y∥2S\left( w \right) = {\left\| {Xw-y } \right\|^2}S(w)=∥Xw−y∥2 （1-2）

（在统计学中，残差平方和函数可以看成n倍的均方误差MSE）

当w=w^w = \hat ww=w^时，S(w)S\left( w \right)S(w)的值取最小，记作：
　　　　　　　　　　　　　　　　　　　　　　　　w^=arg⁡min⁡(S(w))\hat w = \arg \min \left( {S\left( w \right)} \right)w^=argmin(S(w)) （1-3）
通过对S(w)S\left( w \right)S(w)进行微分求最值，可以得到：
　　　　　　　　　　　　　　　　　　　　　　　　　XTXw^=XTy{X^T}X\hat w = {X^T}yXTXw^=XTy （1-4）
如果矩阵XTX{X^T}XXTX是非奇异矩阵，则www有唯一解：
　　　　　　　　　　　　　　　　　　　　　　　　w^=(XTX)−1XTy\hat w = {\left( {{X^T}X} \right)^{ - 1}}{X^T}yw^=(XTX)−1XTy （1-5）

什么是非奇异？

若n阶方阵A的行列式不为零，即 |A|≠0

也就是说如果XTX{X^T}XXTX是奇异矩阵，用最小二乘法求解不出解。也就是为什么会引入岭回归与Lasso回归对奇异矩阵XTX{X^T}XXTX进行变换使得可逆。

1.1 正规方程

正规方程也就是在最小二乘法中求解参数的方程：
　　　　　　　　　　　　　　　　　　　　　　　　w^=(XTX)−1XTy\hat w = {\left( {{X^T}X} \right)^{ - 1}}{X^T}yw^=(XTX)−1XTy （1-6）
公式推导过程并不复杂，这里有我手写的推导过程：

1.2 梯度下降法

梯度下降法是迭代法的一种，可以用于求解最小二乘问题(线性和非线性都可以)。在求解机器学习算法的模型参数，即无约束优化问题时，梯度下降（Gradient Descent）是最常采用的方法之一。在求解损失函数的最小值时，可以通过梯度下降法来一步步的迭代求解，得到最小化的损失函数和模型参数值。在机器学习中，基于基本的梯度下降法发展了两种梯度下降方法，分别为随机梯度下降法和批量梯度下降法。
　　
如果是之前讲的输入特征只有一个x时，也就是直线方程：y=w0+w1xy={w_0}+{w_1}xy=w0+w1x时，我们可以直接手动进行求解系数www，可以直接求导，也可以使用正规方程进行求解。但是往往特征不只一个，这样用矩阵的形式求解就会很困难。于是使用迭代法来逼近最优参数www。

1：梯度下降的一般步骤：
　　1：参数的初始化：通常所有参数都初始化为1；
　　2：确定学习率；
　　3：求代价函数的梯度（所有参数的偏导数）；
　　4：所有参数都沿梯度方向移动一步，步长就是学习率的大小；
　　5：重复步骤4直到参数不再发生变化（此时取到极值点，梯度为0）或达到预先设定的迭代次数。

2：学习率α：
　　学习率一般用希腊字母α表示，可能需要多尝试几次，才能找到合适的学习率。过大的学习率会导致梯度下降时越过代价函数的最小值点，随着训练步数的增加，代价函数不减反增；如果学习率太小，训练中的每一步参数的变化会非常小，这时可以看到代价函数的值在不断减小，但是需要非常大的迭代次数才能到达代价函数的最小值点。

　　图1，学习率过大会导致参数的取值越过最小值点；学习率过小会导致参数变化缓慢
3：代价函数的梯度：
在机器学习中，对代价函数中的每一个参数求偏导数，这些偏导数组成的向量就是代价函数的梯度。

首先，极小化的损失函数：
　　　　　　　　　　　　　　　　　　　　　　　J(w)=MSE=12m∑i=1m(y(i)−(wx(i)+b))2J\left(w \right)=MSE=\frac{1}{2m}{\sum\limits_{i = 1}^m {\left( {{y^{\left( i \right)}} - \left( {w{x^{\left( i \right)}} + b} \right)} \right)} ^2}J(w)=MSE=2m1i=1∑m(y(i)−(wx(i)+b))2 （1-7）

为了求导方便，添加了一个系数1/2{1}/{2}1/2，实际的MSE的定义中是没有的；

公式（1-1）向量化形式可以表示为：hw(x)=wT⋅xh_w(x)=w^T⋅xhw(x)=wT⋅x
　　1 ：J(w)J\left(w\right)J(w)对 w0w_0w0求偏导：
　　　　　　　　　　　　　　　　　∂∂w0J(w)=1m∑i=1m[(hw(x(i))−y(i))⋅x0(i)]\frac{\partial }{{\partial {w_0}}}J\left( w \right) = \frac{1}{m}\sum\limits_{i = 1}^m {\left[ {({h_w}({x^{(i)}}) - {y^{(i)}}) \cdot x_0^{(i)}} \right]}∂w0∂J(w)=m1i=1∑m[(hw(x(i))−y(i))⋅x0(i)]
　　2 ：J(w)J\left(w\right)J(w)对 w1w_1w1求偏导：
　　　　　　　　　　　　　　　　　∂∂w1J(w)=1m∑i=1m[(hw(x(i))−y(i))⋅x1(i)]\frac{\partial }{{\partial {w_1}}}J\left( w \right) = \frac{1}{m}\sum\limits_{i = 1}^m {\left[ {({h_w}({x^{(i)}}) - {y^{(i)}}) \cdot x_1^{(i)}} \right]}∂w1∂J(w)=m1i=1∑m[(hw(x(i))−y(i))⋅x1(i)]
　　3 ：J(w)J\left(w\right)J(w)对 w2w_2w2求偏导：
　　　　　　　　　　　　　　　　　∂∂w2J(w)=1m∑i=1m[(hw(x(i))−y(i))⋅x2(i)]\frac{\partial }{{\partial {w_2}}}J\left( w \right) = \frac{1}{m}\sum\limits_{i = 1}^m {\left[ {({h_w}({x^{(i)}}) - {y^{(i)}}) \cdot x_2^{(i)}} \right]}∂w2∂J(w)=m1i=1∑m[(hw(x(i))−y(i))⋅x2(i)]
．．．
对代价函数的梯度更新，学习率为ααα：
w0=w0−α1m∑i=1m[(hw(x(i))−y(i))⋅x0(i)]{w_0} = {w_0} - α\frac{1}{m}\sum\limits_{i = 1}^m {\left[ {({h_w}({x^{(i)}}) - {y^{(i)}}) \cdot x_0^{(i)}} \right]}w0=w0−αm1i=1∑m[(hw(x(i))−y(i))⋅x0(i)]
w1=w1−α1m∑i=1m[(hw(x(i))−y(i))⋅x1(i)]{w_1} = {w_1} - α\frac{1}{m}\sum\limits_{i = 1}^m {\left[ {({h_w}({x^{(i)}}) - {y^{(i)}}) \cdot x_1^{(i)}} \right]}w1=w1−αm1i=1∑m[(hw(x(i))−y(i))⋅x1(i)]
w2=w2−α1m∑i=1m[(hw(x(i))−y(i))⋅x2(i)]{w_2} = {w_2} - α\frac{1}{m}\sum\limits_{i = 1}^m {\left[ {({h_w}({x^{(i)}}) - {y^{(i)}}) \cdot x_2^{(i)}} \right]}w2=w2−αm1i=1∑m[(hw(x(i))−y(i))⋅x2(i)]
．．．

符号：
1：m为样本个数，i=1,2,3…m
2：ααα为学习率
2：yiy_iyi表示第i个样本的观察值；x1(i){x_1^{(i)}}x1(i)表示第i个样本的第1个特征的观察值；

矩阵形式： w=w−1mαXT(XTw−y)w = w - \frac{1}{m}\alpha {X^T}\left( {{X^T}w - y} \right)w=w−m1αXT(XTw−y)

1.3 代码实现

第一部份： 是使用最小二乘法的正规方程来求解参数，即w^=(XTX)−1XTy\hat w = {\left( {{X^T}X} \right)^{ - 1}}{X^T}yw^=(XTX)−1XTy ：

# 用正规方程求解theta
theta = pinv(X.T @ X) @ X.T @ y  # A@B 等于 np.dot(A, B)

第一部分完整代码：

import matplotlib.pyplot as plt
import numpy as np# 准备数据：X为(n_samples, n_features) y为(n_samples)
X = np.array([[1, 2], [3, 2], [1, 3], [2, 3], [3, 3], [3, 4]])
y = np.array([3.1, 5.1, 4.2, 5.2, 5.9, 6.8])
n_samples, n_features = X.shapeX = np.concatenate((np.ones(n_samples).reshape((n_samples, 1)), X), axis=1)  # 给X添加一列：截距
y = y.reshape((n_samples, 1))  # 将y转换成(n_samples, 1) 便于计算# 用正规方程求解theta
theta = pinv(X.T @ X) @ X.T @ y  # A@B 等于 np.dot(A, B)
intercept = theta[0, 0]  # 截距项
coef = theta[2:, 0]  # 系数print("截距项：%s" % intercept)
print("系数：%s" % coef)

结果：

截距项：0.24662162162161216
系数：[0.90675676 0.91486486]

第二部份： 是使用梯度下降法来更新参数w，矩阵形式： w=w−1mαXT(XTw−y)w = w - \frac{1}{m}\alpha {X^T}\left( {{X^T}w - y} \right)w=w−m1αXT(XTw−y)

theta_neat = theta - alpha*(X.T) @(X@theta - y)/n_samples

第二部分完整代码：

import matplotlib.pyplot as plt
import numpy as np# 准备数据
X = np.array([[1, 2], [3, 2], [1, 3], [2, 3], [3, 3], [3, 4]])
y = np.array([3.1, 5.1, 4.2, 5.2, 5.9, 6.8])
n_samples, n_features = X.shape
X = np.concatenate((np.ones(n_samples).reshape((n_samples, 1)), X), axis=1)  # X添加一列：截距
y = y.reshape((n_samples, 1))# 预定义超参数
alpha = 0.1  # 步长
max_iter = 1e6  # 最大迭代次数，防止模型发散，陷入死循环
epsilon = 1e-6  # 当前后两次迭代theta的变化小于epsilon则视为收敛，停止迭代，返回此时的theta
theta = np.zeros((n_features + 1, 1))  # 初始化theta
costs = []  # 用来存储每一次迭代时损失函数值的变化情况def cost_fuc(X, y, theta, n_samples):# 损失函数return (X @ theta - y).T @ (X @ theta - y) /2 /n_samples# 迭代求解
for iter in range(int(max_iter)):theta_neat = theta - alpha*(X.T) @(X@theta - y)/n_samplescost_new = cost_fuc(X, y, theta, n_samples)[0][0]costs.append(cost_new)if np.abs(theta - theta_neat).sum() < epsilon:theta = theta_neatprint('merge')breaktheta = theta_neat
else:print("get the max_iter, stop iter.")# 求出参数
intercept = theta[0, 0]
coef = theta[1:, 0]
print(costs)
print("截距项：%s" % intercept)
print("系数：%s" % coef)# 可视化迭代过程中损失函数值的变化情况
plt.plot(list(range(len(costs))), costs)
plt.xlabel("iter count")
plt.ylabel("cost function")
plt.show()

结果：

截距项：0.4930883696433331
系数：[0.90676633 0.91490961]

2 岭回归

2.1 岭回归原理

注：正则化是用来防止过拟合的方法。在最开始学习机器学习的课程时，只是觉得这个方法就像某种魔法一样非常神奇的改变了模型的参数。当正则化的使用对象不同时，分为岭回归与Lasso回归。

加上所有参数（不包括θ0θ_0θ0）的绝对值之和，即L1L_1L1范数，此时叫做Lasso回归
加上所有参数（不包括θ0θ_0θ0）的平方和，即L2L_2L2范数，此时叫做岭回归

岭回归与多项式回归唯一的不同在于代价函数上的差别，也就是在代价函数MSE(θ)MSE\left( \theta \right)MSE(θ)上添加了多项式系数www的L2L_2L2范数，岭回归的代价函数如下：
　　　　　　　　　　　　J(w)=1m∑i=1m(y(i)−(wx(i)+b))2+λ∥w∥22=MSE(w)+λ∑i=1nwi2J\left( w\right) = \frac{1}{m}{\sum\limits_{i = 1}^m {\left( {{y^{\left( i \right)}} - \left( {w{x^{\left( i \right)}} + b} \right)} \right)} ^2} + \lambda \left\| w \right\|_2^2 = MSE\left( w \right) + \lambda \sum\limits_{i = 1}^n {w _i^2}J(w)=m1i=1∑m(y(i)−(wx(i)+b))2+λ∥w∥22=MSE(w)+λi=1∑nwi2
(1-1)
为了方便计算导数，通常在公式前乘上1/2，也写成下面的形式：
　　　　　　　　　　　J(w)=12m∑i=1m(y(i)−(wx(i)+b))2+λ2∥w∥22=12MSE(w)+λ2∑i=1nwi2J\left(w \right) = \frac{1}{{2m}}{\sum\limits_{i = 1}^m {\left( {{y^{\left( i \right)}} - \left( {w{x^{\left( i \right)}} + b} \right)} \right)} ^2} + \frac{\lambda }{2}\left\| w \right\|_2^2 = \frac{1}{2}MSE\left( w \right) + \frac{\lambda }{2}\sum\limits_{i = 1}^n {w _i^2}J(w)=2m1i=1∑m(y(i)−(wx(i)+b))2+2λ∥w∥22=21MSE(w)+2λi=1∑nwi2
(1-2)
　　上式中的www是长度为n的多项式系数向量，不包括截距项的系数θ0θ_0θ0；θθθ是长度为n+1的向量，包括截距项的系数θ0θ_0θ0；m为样本数；n为特征数。
　　岭回归的代价函数仍然是一个凸函数，因此可以利用梯度等于0的方式求得全局最优解（正规方　程）：
　　　　　　　　　　　　　　　　　　　　　　　　w=(XTX+λI)−1(XTy)w = {\left( {{X^T}X + \lambda I} \right)^{ - 1}}\left( {{X^T}y} \right)w=(XTX+λI)−1(XTy)
推导公式过程：

　　上述正规方程与一般线性回归的正规方程相比，多了一项λI\lambda IλI，其中I表示单位矩阵。假如XTX{X^T}XXTX是一个奇异矩阵（不是满秩），添加这一项后可以保证该项可逆。由于单位矩阵的形状是对角线上为1其他地方都为0，看起来像一条山岭，因此而得名。
　　
　　除了上述正规方程之外，还可以使用梯度下降的方式来进行岭回归方程代价函数的求解。这里对式子1−2来求导：
　　　　　　　　　　　　　　　　　　　　　　∇wJ(w)=1mXT(X⋅w−y)+λw{\nabla _w }J\left(w \right) = \frac{1}{m}{X^T}\left( {X \cdot w - y} \right) + \lambda w∇wJ(w)=m1XT(X⋅w−y)+λw
(1-3)
因为式子1−2中和式第二项不包含θ0θ_0θ0，因此求导后，上式第二项中的www本来也不包含θ0θ_0θ0。为了计算方便，添加θ0θ_0θ0到www。因此在梯度下降的过程中，参数的更新可以表示成下面的公式：
　　　　　　　　　　　　　　　　　　　　　　w=w−(αmXT(X⋅w−y)+λw)w =w - \left( {\frac{\alpha }{m}{X^T}\left( {X \cdot w - y} \right) + \lambda w} \right)w=w−(mαXT(X⋅w−y)+λw)

其中α\alphaα为学习率，λ\lambdaλ为正则化项的参数。
算法求解步骤：
　　这里使用梯度下降法来求解参数系数的过程：
　　　1 首先给定参数初值
　　　2 给定步长和最大的迭代次数
　　　3 更新w
　　更新参数系数：
　　　　　　　　　　　　　w1=w1−α1m∑i=1m(((hw(x(i))−y(i))x1(i))+λmw1)j=1,2,....n{w _1} = {w _1} - \alpha \frac{1}{m}\sum\limits_{i = 1}^m {\left( {\left( {\left( {{h_w }\left( {{x^{\left( i \right)}}} \right) - {y^{\left( i \right)}}} \right)x_1^{\left( i \right)}} \right) + \frac{\lambda }{m}{w _1}} \right)} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j = 1,2,....nw1=w1−αm1i=1∑m(((hw(x(i))−y(i))x1(i))+mλw1)j=1,2,....n
　　　　　　　　　　　　　w2=w2−α1m∑i=1m(((hw(x(i))−y(i))x2(i))+λmw2)j=1,2,....n{w _2} = {w _2} - \alpha \frac{1}{m}\sum\limits_{i = 1}^m {\left( {\left( {\left( {{h_w }\left( {{x^{\left( i \right)}}} \right) - {y^{\left( i \right)}}} \right)x_2^{\left( i \right)}} \right) + \frac{\lambda }{m}{w _2}} \right)} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j = 1,2,....nw2=w2−αm1i=1∑m(((hw(x(i))−y(i))x2(i))+mλw2)j=1,2,....n
　　　　　　　　　　　　　　　　　　　　　　　　……
　　　　　　　　　　　　　wj=wj(1−αλm)−α1m∑i=1m(hw(x(i))−y(i))xj(i){w _j} = {w _j}\left( {1 - \alpha \frac{\lambda }{m}} \right) - \alpha \frac{1}{m}\sum\limits_{i = 1}^m {\left( {{h_w }\left( {{x^{\left( i \right)}}} \right) - {y^{\left( i \right)}}} \right)x_j^{\left( i \right)}}wj=wj(1−αmλ)−αm1i=1∑m(hw(x(i))−y(i))xj(i)
　　4 重复以下步骤，直到收敛。
可以看出给定的步长α\alphaα越大，收缩率越大，系数对于共线性的鲁棒性更强。最终得到收敛下的wjw_jwj，其中参数可以表示为矩阵形式：
　　　　　　　　　　　　　　　　　　　w=w(1−αλm)−α1mαXT(X⋅w−y)w = w \left( {1 - \alpha \frac{\lambda }{m}} \right) - \alpha \frac{1}{m}\alpha {X^T}\left( {X \cdot w - y} \right)w=w(1−αmλ)−αm1αXT(X⋅w−y)

2.1 代码实现

import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_errordata = np.array([[ -2.95507616,  10.94533252],[ -0.44226119,   2.96705822],[ -2.13294087,   6.57336839],[  1.84990823,   5.44244467],[  0.35139795,   2.83533936],[ -1.77443098,   5.6800407 ],[ -1.8657203 ,   6.34470814],[  1.61526823,   4.77833358],[ -2.38043687,   8.51887713],[ -1.40513866,   4.18262786]])
m = data.shape[0]  # 样本大小
X = data[:, 0].reshape(-1, 1)  # 将array转换成矩阵
y = data[:, 1].reshape(-1, 1)# 代价函数
def L_theta(theta, X_x0, y, lamb):"""lamb: lambda, the parameter of regularizationtheta: (n+1)·1 matrix, contains the parameter of x0=1X_x0: m·(n+1) matrix, plus x0"""h = np.dot(X_x0, theta)  # np.dot 表示矩阵乘法theta_without_t0 = theta[1:]L_theta = 0.5 * mean_squared_error(h, y) + 0.5 * lamb * np.sum(np.square(theta_without_t0))return L_theta# 梯度下降
def GD(lamb, X_x0, theta, y, alpha):"""lamb: lambda, the parameter of regularizationalpha: learning rateX_x0: m·(n+1), plus x0theta: (n+1)·1 matrix, contains the parameter of x0=1"""for i in range(T):h = np.dot(X_x0, theta)theta_with_t0_0 = np.r_[np.zeros([1, 1]), theta[1:]]  # set theta[0] = 0theta -= (alpha * 1/m * np.dot(X_x0.T, h - y) + lamb*(theta_with_t0_0))  # add the gradient of regularization termif i%50000==0:print(L_theta(theta, X_x0, y, lamb))return thetaT = 1200000  # 迭代次数
degree = 11
theta = np.ones((degree + 1, 1))  # 参数的初始化，degree = 11，一个12个参数
alpha = 0.0000000006 # 学习率
lamb = 0.0001
poly_features_d = PolynomialFeatures(degree=degree, include_bias=False)
X_poly_d = poly_features_d.fit_transform(X)
X_x0 = np.c_[np.ones((m, 1)), X_poly_d]  # X_poly_d添加一列：截距
theta = GD(lamb=lamb, X_x0=X_x0, theta=theta, y=y, alpha=alpha)
print(theta)

结果：

3.599476296986373
[[ 1.00078848e+00][-1.03862735e-05][ 3.85144400e-05][-3.77233288e-05][ 1.28959318e-04][-1.42449160e-04][ 4.42760996e-04][-5.11518471e-04][ 1.42533716e-03][-1.40265037e-03][ 3.13638870e-03][ 1.21862016e-03]]

3 Lasso回归

Lasso回归于岭回归非常相似，它们的差别在于使用了不同的正则化项。最终都实现了约束参数从而防止过拟合的效果。但是Lasso之所以重要，还有另一个原因是：Lasso能够将一些作用比较小的特征的参数训练为0，从而获得稀疏解。也就是说用这种方法，在训练模型的过程中实现了降维(特征筛选)的目的。
　　Lasso回归的代价函数为：
　　　　　　　　　J(w)=12m∑i=1m(y(i)−(wx(i)+b))2+λ∥w∥1=12MSE(w)+λ∑i=1nwi2J\left(w \right) = \frac{1}{{2m}}{\sum\limits_{i = 1}^m {\left( {{y^{\left( i \right)}} - \left( {w{x^{\left( i \right)}} + b} \right)} \right)} ^2} + {\lambda }{}\left\| w \right\|_1 = \frac{1}{2}MSE\left( w \right) + {\lambda }\sum\limits_{i = 1}^n {w _i^2}J(w)=2m1i=1∑m(y(i)−(wx(i)+b))2+λ∥w∥1=21MSE(w)+λi=1∑nwi2 (3-1)
　　上式中的w是长度为n的向量，不包括截距项的系数θ0θ_0θ0, θ是长度为n+1的向量，包括截距项的系数θ0θ_0θ0，m为样本数，n为特征数。∣∣w∣∣1||w||_1∣∣w∣∣1表示参数w的L1L_1L1范数，也是一种表示距离的函数。
　　
　　假如：w表示3维空间中的一个点(x,y,z)，那么∣∣w∣∣1=∣x∣+∣y∣+∣z∣||w||_1=|x|+|y|+|z|∣∣w∣∣1=∣x∣+∣y∣+∣z∣，即各个方向上的绝对值（长度）之和。式子(3-1)的梯度为：
　　　　　　　　　　　　　　　　　　　　　　∇wJ(w)+λ(sign(w1)sign(w2)...sign(wn)){\nabla _w}J\left( w \right) + \lambda \left( \begin{array}{l} sign\left( {{w_1}} \right)\\ sign\left( {{w_2}} \right)\\ ...\\ sign\left( {{w_n}} \right) \end{array} \right)∇wJ(w)+λ⎝⎜⎜⎛sign(w1)sign(w2)...sign(wn)⎠⎟⎟⎞
　　其中sign(wi)sign(w_i)sign(wi)由wiw_iwi的符号决定:
　　　　　　　　　　　　　　　　　　　　　　sign(wi)={1wi>00wi=0−1wi<0sign\left( {{w_i}} \right) = \left\{ \begin{array}{l} 1{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {w_i} > 0\\ 0{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {w_i} = 0\\ -1{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {w_i} < 0 \end{array} \right.sign(wi)=⎩⎨⎧1wi>00wi=0−1wi<0

Lasso回归直接使用scikit-learn中的函数：

from sklearn.linear_model import Lasso

可以参考官方文档：http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html

4 案例实战1-车流量预测

scikit-learn中有专门计算岭回归的函数，而且效果要比上面的方法好。使用scikit-learn中的岭回归，只需要输入以下参数：

Ridge(alpha=1.0, fit_intercept=True, normalize=False, copy_X=True, max_iter=None, tol=0.001, solver=‘auto’, random_state=None)
• alpha：用于指定λλ \lambdaλ值参数，默认为1。
• fit_intercept：bool类型，是否需要拟合截距项，默认为True。
• normalize：bool类型，建模时是否对数据集做标准化处理，默认为False。
• copy_X：bool类型，是否复制自变量X的数值，默认为True。
• max_iter：指定模型的最大迭代次数。
• solver：指定模型求解最优化问题的算法，默认为’auto’。
• random_state：指定随机生成器的种子。


import numpy as np
import pandas as pd
from sklearn.linear_model import Ridge  # 通过sklearn.linermode1加载岭回归方法
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures  # 通过sklearn.preprocessing加载PolynomialFeatures用于创建多项式特征
from sklearn.model_selection import train_test_split  # 交叉验证data_df = pd.read_csv('岭回归.csv')
data = np.array(data_df)
plt.plot(data[:, 5])  # 展示车流量信息X = data[:, 1:5]  # X用于保存1-4维数据,即属性
y = data[:, 5]  # y用于保存第5维数据,即车流量
poly = PolynomialFeatures(6)  # 用于创建最高次数6次方的的多项式特征,多次试验后决定采用6次
X = poly.fit_transform(X)  # X为创建的多项式特征train_set_X, test_set_X, train_set_y, test_set_y = train_test_split(X, y, test_size=0.3, random_state=0)
# 将所有数据划分为训练集和测试集, test_ size表示测试集的比例， random_state是随机数种子clf = Ridge(alpha=1.0, fit_intercept=True)  # 创建岭回归实例
clf.fit(train_set_X, train_set_y)  # 调用fit函数使用训练集训练回归器
score = clf.score(test_set_X, test_set_y)  # 评价拟合值的好坏(最大值：1)
# print(score)
'''
利用测试集计算回归曲线的拟合优度，clf. score返回值为0.7620拟合优度,
用于评价拟合好坏,最大为1,无最小值,当对所有输入都输出同一个值时,拟合优度为0。
'''# 绘制拟合曲线
start = 200  # 画一段200到300范围内的拟合曲线
end = 300
y_pre = clf.predict(X)  # 是调用predict函数的拟合值
time = np.arange(start, end)fig = plt.figure()  # 定义一个图片
ax = fig.add_subplot(1, 1, 1)
ax.plot(time, y[start:end], label='real')
ax.plot(time, y_pre[start:end],'r', label='predict')  # 展示真实数据(蓝色)以及拟合的曲线(红色)
plt.legend(loc = 'upper left') # 设置图例的位置
props = {'title' : 'Traffic flow forecast','xlabel' : 'Period of time[200-300]','ylabel' : 'Number of traffic'
}
ax.set(**props)
plt.show()

总结：
1：岭回归与Lasso回归的出现是为了解决线性回归出现的过拟合以及在通过正规方程方法求解θ的过程中出现的x转置乘以x不可逆这两类问题的，这两种回归均通过在损失函数中引入正则化项来达到目的。
2：正则化主要针对自变量之间存在多重共线性或者自变量个数多于样本量的情况。

参考文献：
[1] Peter Harrington.《机器学习实战》.人民邮电出版社,2013-6
[2] 刘顺祥.《从零开始学Python数据分析与挖掘》.清华大学出版社,2018
[3] http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html