3.线性回归

线性回归

最小二乘法

前提

假设数据集$$\mathcal{D}={(x_1, y_1),(x_2, y_2),\cdots,(x_N, y_N)}$$，$$x_i\in\mathbb{R},y_i\in\mathbb{R},i=1,2,\cdots,N$$,记为： $$ \begin{align} X&=(x_1,x_2,\cdots,x_N)^T \ &= \begin{bmatrix} x_{1}^T \ x_{2}^T \ \vdots \ x_{N}^T \ \end{bmatrix}{N\times p} \ & = \begin{bmatrix} x{11} & x_{12} & \cdots & x_{1p} \ x_{21} & x_{22} & \cdots & x_{2p} \ \vdots & \vdots & \ddots & \vdots \ x_{N1} & x_{N2} & \cdots & x_{Np} \ \end{bmatrix}_{N\times p} \ \end{align} $$

$$ \begin{align} Y&=(y_1,y_2,\cdots,y_N)^T \ &= \begin{bmatrix} y_{1}^T \ y_{2}^T \ \vdots \ y_{N}^T \ \end{bmatrix}{N\times p} \ & = \begin{bmatrix} y{11} & y_{12} & \cdots & y_{1p} \ y_{21} & y_{22} & \cdots & y_{2p} \ \vdots & \vdots & \ddots & \vdots \ y_{N1} & y_{N2} & \cdots & y_{Np} \ \end{bmatrix}_{N\times p} \ \end{align} $$

线性回归假设： $$ f(w)=w^Tx $$

构造数据

$$X=[x^T,1]^T$$,利用等式$$y=3x+2$$造些伪数据，并给$$x$$添加一些噪声数据，其中$$w=[3,2]^T$$，即$$X.dot(w)$$。

import numpy as np

"""
	构造数据
"""
X = np.linspace(0,100,100) # 生成100个样本数据
X = np.c_[X,np.ones(100)]
w = np.asarray([3,2])
Y = X.dot(w) # 矩阵点乘
X = X.astype('float')
Y = Y.astype('float')
X[:,0] += np.random.normal(size=(X[:,0].shape))*3 #添加噪声
Y = Y.reshape(Y.size,1)

$$x$$如图所示

$$X$$如图所示

绘制图形

import matplotlib.pyplot as plt

# 绘制(x,y)
plt.scatter(X[:,0],Y)
plt.plot(np.arange(0,100).reshape((100,1)),Y,'y.')
plt.xlabel("X")
plt.ylabel("Y")

如图：

Text(0, 0.5, 'Y')

定义

直接求闭式解

采用二范数定义的平方误差来定义损失函数： $$ \begin{align} loss function： L(w)&=\sum\limits_{i=1}^N||w^Tx_i-y_i||^2_2 \ & = \sum\limits_{i=1}^N(w^Tx_i-y_i)^2 \ &=(w^Tx_1-y_1,\cdots,w^Tx_N-y_N)\cdot (w^Tx_1-y_1,\cdots,w^Tx_N-y_N)^T\nonumber \ &=(w^TX^T-Y^T)\cdot (Xw-Y) \ &=w^TX^TXw-Y^TXw-w^TX^TY+Y^TY\nonumber \ &=w^TX^TXw-2w^TX^TY+Y^TY \end{align} $$

$$ (w^Tx_i-y_i)^2=(Y−Xw)^T(Y−Xw) $$

最小化值$$ \hat{w}$$ ,，进行求导： $$ \begin{align} \hat{w}=\mathop{argmin}\limits_wL(w)&\longrightarrow\frac{\partial}{\partial w}L(w)=0\nonumber\ &\longrightarrow2X^TX\hat{w}-2X^TY=0\nonumber\ &\longrightarrow \hat{w}=(X^TX)^{-1}X^TY=X^+Y \end{align} $$

$$(X^TX)^{-1}X^T$$记为$$X^+$$(伪逆)

这个式子中 $$(X^TX)^{-1}X^T$$ 又被称为伪逆。对于行满秩或者列满秩的 $$X$$，可以直接求解，但是对于非满秩的样本集合，需要使用奇异值分解（SVD）的方法，对 $$X$$ 求奇异值分解，得到： $$ X=U\Sigma V^T $$ 于是： $$ X^+=V\Sigma^{-1}U^T $$

在几何上，最小二乘法相当于模型（这里就是直线）和试验值的距离的平方求和，假设我们的试验样本张成一个 $p$ 维空间（满秩的情况）：$$X=Span(x_1,\cdots,x_N)$$，而模型可以写成 $$f(w)=w^Tx=X^T\beta$$，反过来看，也就是 $$x_1,\cdots,x_N$$ 的某种组合，而最小二乘法就是说希望 $$Y$$ 和这个模型距离越小越好，于是它们的差应该与这个张成的空间垂直：

$$ X^T\cdot(Y-X\beta)=0\longrightarrow\beta=(X^TX)^{-1}X^TY $$

提示：$$ \vec a \perp\vec b => a^T \cdot b = 0 $$

通过伪逆求解到的结果有如下优点：

（1）当$$ w $$有解时，$$w=X^+Y $$是所有解中欧几里得距离$$||w||^2$$最小的一个；

（2）当$$w$$无解时，通过伪逆得到的$$w$$是使得$$Xw$$与$$Y $$的欧几里得距离$$ (w^Tx_i-y_i)^2 $$最小

• 求逆

# np.linalg.inv(X) <==> np.matrix(X).I

• 求伪逆

np.linalg.pinv(X)

$$ X $$数据如图所示，

$$ Y $$数据如图所示，

求出参数$$ w $$结果为：

np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y)

求出$$Y1$$​为：

w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y)
Y1 = X.dot(w)
Y1[:15]

绘图：

# 绘制(x,y)
plt.scatter(X[:,0],Y)
plt.plot(X[:,0],Y1,'y')
plt.xlabel("X")
plt.ylabel("Y")

梯度下降求解

随机梯度下降法

但对于数据量很大的情况，求闭式解的方式会让内存很吃力，我们可以通过随机梯度下降法（SGD）对$$ w $$进行更新，首先随机初始化$$ w $$，然后使用如下的迭代公式对$$w$$进行迭代更新： $$ w:=w-\eta \frac{dL}{dw}，其中\eta 参数取值一般很小，比如0.0001 $$ 其中$$\frac{dL}{dw}=2X^TX\hat{w}-2X^TY=2X^T(X\hat{w}-Y) （1）$$，$$ (w^Tx_i-y_i)^2=(Y−Xw)^T(Y−Xw) （2）$$

# 初始化
w = np.random.random(size=(2,1))
# 迭代次数
epoches = 100
# eta参数
eta = 0.0000001

#记录loss变化(损失)
losses = []
for e in range(epoches):
    nw = 2*X.T.dot(X.dot(w)-Y) #（1）
    w = w - eta*nw
    # reshape(-1)：变成一维数组
    losses.append((Y-X.dot(w)).T.dot(Y-X.dot(w)).reshape(-1)) # （2）
w

参数$$\eta$$以及迭代次数，不要设置太低或者太高，可能会使loss变化达不到预期目标

求出$$w$$结果为：

可视化：

loss变化：

plt.plot(losses)
plt.xlabel("iterations")
plt.ylabel("loss")

小批量梯度下降法

问题

在上面的梯度下降的例子中存在一个问题，$$w1$$基本能收敛到3附近，而$$w2$$却基本在0附近，很难收敛到2，说明$$w1比w2$$更容易收敛($$ w=[w1,w2]^T $$)，这很容易理解，模型可以写作：$$f(x)=x∗w1+1⋅w2$$，如果$$x$$量纲比1大很多，为了使$$ f(x)$$变化，只需更新少量的$$w1$$就能达到目的，而$$w2$$的更新动力略显不足；可以考虑用如下方式：

（1）对输入$$X $$进行归一化，使得x无量纲，$$w1,w2 $$的更新动力一样（后面封装代码时添加上），如下图；

（2）梯度更新时，w1,w2使用了一样的学习率，可以让w1,w2使用不一样的学习率进行更新，比如对w2使用更大的学习率进行更新（可以利用学习率自适应一类的梯度下降法，比如adam），如下图：

使用sklearn

from sklearn.linear_model import LinearRegression

lr = LinearRegression()
lr.fit(X[:,:-1],Y) #训练数据
predict = lr.predict(X[:,:-1])
#查看w系数,b常量
print('w:',lr.coef_,'b:',lr.intercept_)
#查看标准差
np.std(Y-predict)

假设函数为$$ y=wx+\epsilon $$，

• lr.coef_:为$$ w $$，即系数

• lr.intercept_:为$$ \epsilon $$，即常量

噪声为高斯分布的 MLE（极大似然估计）

对于一维的情况，记 $$y=w^Tx+\epsilon,\epsilon\sim\mathcal{N}(0,\sigma^2)$$，那么 $$y\sim\mathcal{N}(w^Tx,\sigma^2) $$。代入极大似然估计中：

其中正态分布函数的概率密度为$$p(y|x_iw)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y_i-w^Tx_i)^2}{2\sigma^2}}$$ $$ \begin{align} L(w)=\log p(Y|X,w)&=\log\prod\limits_{i=1}^Np(y_i|x_i,w)\nonumber\ &=\sum\limits_{i=1}^N\log(\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y_i-w^Tx_i)^2}{2\sigma^2}})\ &=\sum\limits_{i=1}^N\big ( \log\frac{1}{\sqrt{2\pi}\sigma}+\log e^{-\frac{(y_i-w^Tx_i)^2}{2\sigma^2}} \big )\ &=\sum\limits_{i=1}^N \big ( \log\frac{1}{\sqrt{2\pi}\sigma}- \frac{(y_i-w^Tx_i)^2}{2\sigma^2} \big )\

\mathop{argmax}\limits_wL(w) &=\mathop{argmax}\limits_w\sum \limits_{i=1}^N -\frac 1{2\sigma^2} (y_i-w^Tx_i)^2\ &=\mathop{argmin}\limits_w\sum\limits_{i=1}^N (y_i-w^Tx_i)^2 (进而转化为求最小) \end{align} $$

LSE（最小二乘估计）<==> MLE（极大似然估计），其中噪声为高斯分布

（正则化）Regularized LSE <==> MAP，其中噪声为高斯分布

权重先验为高斯分布的 MAP

取先验分布 $$w\sim\mathcal{N}(0,\sigma_0^2)$$。于是： 

概率密度为$$p(w)=\frac{1}{\sqrt{2\pi\sigma}}e^{-\frac{w^2}{2\sigma_0^2}}$$，$$p(y|w)=\frac{1}{\sqrt{2\pi\sigma}}e^{-\frac{(y-w^Tx)^2}{2\sigma^2}}$$。

有$$p(y|w) \cdot p(w) =\frac{1}{\sqrt{2\pi}\sigma_0}\cdot \frac{1}{\sqrt{2\pi}\sigma} \cdot e^{-\frac{(y-w^Tx)^2}{2\sigma^2}} \cdot e^{-\frac{||w||^2}{2\sigma^2}}$$ $$ \begin{align} \hat{w}&=\mathop{argmax}\limits_wp(w|y) \ &=\mathop{argmax}\limits_wp(y|w)p(w)\nonumber\ &=\mathop{argmax}\limits_w\log p(y|w)p(w)\nonumber\ &=\mathop{argmax}\limits_w(\log p(y|w)+\log p(w))\nonumber\ &=\mathop{argmin}\limits_w[\frac {(y-w^Tx)^2} {2\sigma^2} +\frac {||w||^2} {2\sigma_0^2}] \ &=\mathop{argmin}\limits_w[(y-w^Tx)^2+\frac{\sigma^2}{\sigma_0^2}w^Tw]\ \end{align} $$ 这里省略了 $X$，$p(Y)$和 $w$ 没有关系，同时也利用了上面高斯分布的 MLE的结果。

将会看到，超参数 $\sigma_0$的存在和下面会介绍的 Ridge 正则项可以对应，同样的如果将先验分布取为 Laplace 分布，那么就会得到和 L1 正则类似的结果。

正则化

过拟合

演示

• 伪造数据

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

#造伪样本
X=np.linspace(0,100,100)
X=np.c_[X,np.ones(100)]
w=np.asarray([3,2])
Y=X.dot(w)
X=X.astype('float')
Y=Y.astype('float')
X[:,0]+=np.random.normal(size=(X[:,0].shape))*3#添加噪声
Y=Y.reshape(100,1) # 一维数组转变为二维矩阵

• 拟合数据并可视化

#拟合数据并可视化
w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y)
Y1 = X.dot(w)
# 绘制(x,y)
plt.scatter(X[:,0],Y)
plt.plot(X[:,0],Y1,'y')
plt.xlabel("X")
plt.ylabel("Y")

# 初始化
w = np.random.random(size=(2,1))
# 迭代次数
epoches = 100
# eta参数
eta = 0.0000001

#记录loss变化(损失)
losses = []
for e in range(epoches):
    nw = 2*X.T.dot(X.dot(w)-Y) #（1）
    w = w - eta*nw
    # reshape(-1)：变成一维数组
    losses.append((Y-X.dot(w)).T.dot(Y-X.dot(w)).reshape(-1)) # （2）

plt.plot(losses)

• 加入几个异常点

# np.concatenate(x,y)将x和y两个矩阵垂直拼接在一起
X=np.concatenate([X,np.array([[100,1],[101,1],[102,1],[103,1],[104,1]])])
Y=np.concatenate([Y,np.array([[3000],[3300],[3600],[3800],[3900]])])

#拟合数据并可视化
w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y)
Y1 = X.dot(w)
# 绘制(x,y)
plt.scatter(X[:,0],Y)
plt.plot(X[:,0],Y1,'y')
plt.xlabel("X")
plt.ylabel("Y")

因此在实际应用时，如果样本容量不远远大于样本的特征维度，很可能造成过拟合，对这种情况，有下面三个解决方式：

1. 加数据
2. 特征选择/特征提取（降低特征维度）如 PCA 算法。
3. 正则化

正则化一般是在损失函数（如上面介绍的最小二乘损失）上加入正则化项（表示模型的复杂度对模型的惩罚），下面我们介绍一般情况下的两种正则化框架。 $$ \begin{align} (lasso \quad regression) \quad L1&:\mathop{argmin}\limits_wL(w)+\lambda||w||_1,\lambda\gt0\ (ridge \quad regression/岭回归/权值衰减) \quad L2&:\mathop{argmin}\limits_wL(w)+\lambda||w||^2_2,\lambda \gt 0 \end{align} $$

$$||w||^2 = w^T \cdot w$$

---

$\text{当}||w|| = |w_1| + |w_2| + ··· + |w_n| \text{，称之为L1正则化}$。

$$\text{当}||w|| = w_1^2 + w_2^2 + ··· + w_n^2$$。

这里$$||w||$$也被称为惩罚因子，那么，为什么加上惩罚因子之后，就能避免过拟合呢？

---

L1 Lasso

L1正则化可以引起稀疏解。

从最小化损失的角度看，由于 L1 项求导在0附近的左右导数都不是0，因此更容易取到0解。

从另一个方面看，L1 正则化相当于： $$ \mathop{argmin}\limits_wL(w)\ s.t. ||w||_1\lt C\ \hat{w}=(X^TX)^{-1}(X^TY-\frac{\lambda I}2)（后面会用这个） $$ 我们已经看到平方误差损失函数在 $w$ 空间是一个椭球，因此上式求解就是椭球和 $||w||_1=C$的切点，因此更容易相切在坐标轴上。

下面分别对LASSO的cost function的两部分求解:

1. $$RSS$$部分

$$ \begin{align} RSS(w) &= \sum_{i=1}^m(w^Tx_i-y_i)^2 \ &= \sum_{i=1}^m(\sum_{j=1}^n{x_{ij}w_j}-y_i)^2 \end{align} $$

求导：

$$ (对w_k进行求导，w_{i≠k}为常数) $$ $$ \begin{align} \frac{\partial}{\partial w_k}RSS(w) &= 2\sum_{i=1}^m{x_{ik}}(\sum_{j=1}^nx_{ij}w_j-y_i) \ &= 2\sum_{i=1}^m(x_{ik}\sum_{j=1}^nx_{ij}w_j-y_ix_{ik}) \ &=2\sum_{i=1}^m(x_{ik}\sum_{j=1,j≠k}^nx_{ij}w_j+x_{ik}^2w_k-y_ix_{ik}) \ &=2\sum_{i=1}^mx_{ik}(\sum_{j=1,j≠k}^nx_{ij}w_j-y_i)+2w_k\sum_{i=1}^m x_{ik}^2 \ \end{align} $$

$$ 令p_k=\sum_{i=1}^mx_{ik}(y_i-\sum_{j=1,j≠k}^nx_{ij}w_j),z_k=\sum_{i=1}^m x_{ik}^2 $$

有， $$ \begin{align} \frac{\partial}{\partial w_k}RSS(w) &= -2p_k+2z_kw_k \end{align} $$

1. 惩罚项的求导

$$ \begin{align} \lambda \frac{\partial \sum_{i=1}^n|w_i| }{\partial w_k} &= \left{
\begin{array}{lr}
-\lambda, & w_k<0\ [-\lambda,\lambda], & w_k=0\
\lambda, & w_k>0\
\end{array}
\right. \end{align} $$

进行整体的求导可得： $$ \begin{align} \frac{\partial}{\partial w_k}L(w) &=\frac{\partial}{\partial w_k}RSS(w) +\lambda \frac{\partial \sum_{i=1}^n|w_i| }{\partial w_k}\ &= 2z_kw_k-2p_k+ \left{
\begin{array}{lr}
-\lambda, & w_k<0\ [-\lambda,\lambda], & w_k=0\
\lambda, & w_k>0\
\end{array}
\right. \ &= \left{
\begin{array}{lr}
2z_kw_k-2p_k-\lambda, & w_k<0\ [-2p_k-\lambda,-2p_k+\lambda], & w_k=0\
2z_kw_k-2p_k+\lambda, & w_k>0\
\end{array}
\right. \ & = 0 \end{align} $$ 得到： $$ \begin{align} \hat{w_k} &= \left{
\begin{array}{lr}
(\lambda/2+p_k)/z_k, & p_k<-\lambda/2,\ 0, & -\lambda/2 \leq p_k \leq \lambda/2\
(p_k-\lambda/2)/z_k, & p_k>\lambda/2\
\end{array}
\right. \ 其中，p_k&=\sum_{i=1}^mx_{ik}(y_i-\sum_{j=1,j≠k}^nx_{ij}w_j),z_k=\sum_{i=1}^m x_{ik}^2，\& (w^Tx_i-y_i)^2=(Y−Xw)^T(Y−Xw) \end{align} $$

手动实现代码

• 利用公式（1）直接求解w

import numpy as np
import matplotlib.pyplot as plt

# 计算残差
def RSS(X,y,w):
     return (Y-X.dot(w)).T.dot(Y-X.dot(w))

#造伪样本
X=np.linspace(0,100,100)
X=np.c_[X,np.ones(100)]
w=np.asarray([3,2])
print(w)
Y=X.dot(w)
X=X.astype('float')
Y=Y.astype('float')
X[:,0]+=np.random.normal(size=(X[:,0].shape))*3#添加噪声
Y=Y.reshape(Y.shape[0],1) # 一维数组转变为二维矩阵

# np.concatenate(x,y)将x和y两个矩阵垂直拼接在一起
X=np.concatenate([X,np.array([[100,1],[101,1],[102,1],[103,1],[104,1]])])
Y=np.concatenate([Y,np.array([[3000],[3300],[3600],[3800],[3900]])])
#拟合数据并可视化

lam = 1000000
i = np.eye(X.shape[1])
w = np.linalg.pinv(X.T.dot(X)).dot((X.T.dot(Y))-(lam*i)/2)
Y1=X.dot(w[:,0])
plt.scatter(X[:,:-1],Y)
plt.plot(X[:,:-1],Y1)

• 利用公式（2）直接求解w

import numpy as np
import matplotlib.pyplot as plt

# 计算残差
def RSS(X,y,w):
     return (Y-X.dot(w)).T.dot(Y-X.dot(w))

#造伪样本
X=np.linspace(0,100,100)
X=np.c_[X,np.ones(100)]
w=np.asarray([3,2])
print(w)
Y=X.dot(w)
X=X.astype('float')
Y=Y.astype('float')
X[:,0]+=np.random.normal(size=(X[:,0].shape))*3#添加噪声
Y=Y.reshape(Y.shape[0],1) # 一维数组转变为二维矩阵

# np.concatenate(x,y)将x和y两个矩阵垂直拼接在一起
X=np.concatenate([X,np.array([[100,1],[101,1],[102,1],[103,1],[104,1]])])
Y=np.concatenate([Y,np.array([[3000],[3300],[3600],[3800],[3900]])])
#拟合数据并可视化
m,n = X.shape
lambda_params = 100

z = [0 for i in range(n)]
p = [0 for i in range(n)]
w = np.eye(n,1)

for k in range(n):
    z[k] = np.sum(X[:,k]*X[:,k],axis=0)
    p[k] = np.sum([X[i,k]* (Y[i]-np.sum([X[i,j]*w[j] for j in range(n)])) for i in range(m)],axis=0)
    if p[k] > lambda_params/2:
        w_k = (p[k]-lambda_params/2) / z[k]
    elif p[k] < -lambda_params/2:
        w_k = (lambda_params/2+p[k]) / z[k]
    else:
        w_k = 0
    w[k] = w_k
y = X.dot(w)
plt.scatter(X[:,:-1],Y)
plt.plot(X[:,:-1],y)
w

使用sklearn

from sklearn.linear_model import Lasso
from sklearn.preprocessing import StandardScaler

# 特征值和目标值是都必须进行标准化处理, 实例化两个标准化API
std_x = StandardScaler()
std_y = StandardScaler()
X = std_x.fit_transform(X)
Y = std_y.fit_transform(Y.reshape(-1,1))

lasso = Lasso()
lasso.fit(X[:,:-1],Y)
predict = lasso.predict(X[:,:-1])
plt.plot(X[:,:-1],predict)
plt.scatter(X[:,:-1],Y)

L2 Ridge

$$ s.t. ||w||_2^2\lt C\ $$

$$ \begin{align} \hat{w}=\mathop{argmin}\limits_wL(w)+\lambda w^Tw&\longrightarrow\frac{\partial}{\partial w}L(w)+2\lambda w=0\nonumber\ &\longrightarrow2X^TX\hat{w}-2X^TY+2\lambda \hat w=0\nonumber\ &\longrightarrow \hat{w}=(X^TX+\lambda \mathbb{I})^{-1}X^TY \end{align} $$

可以看到，这个正则化参数和前面的 MAP 结果不谋而合。利用2范数进行正则化不仅可以是模型选择 $$w$$ 较小的参数，同时也避免 $$ X^TX$$不可逆的问题。

$$\mathbb{I}是一个单位矩阵$$

定义符号函数并进行向量化，用于对L1正则化项的梯度计算:

手动实现代码

alpha = 100.0

w = np.linalg.pinv(
  X.T.dot(X)+alpha*np.eye(X.T.dot(X).shape[0])
).dot(X.T.dot(Y))
y = X.dot(w)
plt.scatter(X[:,:-1],Y)
plt.plot(X[:,:-1],y)

使用sklearn

from sklearn.linear_model import Ridge
from sklearn.preprocessing import StandardScaler

std_x = StandardScaler()
std_y = StandardScaler()
X2 = std_x.fit_transform(X)
Y2 = std_y.fit_transform(Y.reshape(-1,1))

lasso = Ridge(alpha=100)
lasso.fit(X2[:,:-1],Y2)
predict = lasso.predict(X2[:,:-1])
plt.plot(X2[:,:-1],predict)
plt.scatter(X2[:,:-1],Y2)

L1与L2结合--ElasticNet回归

from sklearn.linear_model import ElasticNet

std_x = StandardScaler()
std_y = StandardScaler()
X3 = std_x.fit_transform(X)
Y3 = std_y.fit_transform(Y.reshape(-1,1))

lasso = ElasticNet()
lasso.fit(X3[:,:-1],Y3)
predict = lasso.predict(X3[:,:-1])
plt.plot(X3[:,:-1],predict)
plt.scatter(X3[:,:-1],Y3)

1. 当岭参数$$\lambda=0$$时，得到的解是最小二乘解。
2. 当岭参数 $$\lambda$$ 趋向更大时，岭回归系数 $$w_i$$ 趋向于0，约束项 $$t$$ 很小。