3.1线性回归：最小二乘法+岭回归

简介

介绍了线性回归模型，包括最小二乘法的原理、最小二乘法中的注意点以及代码实现方法。

什么是回归

回归分析是一种统计学方法，当给定样本的时候，用于估测多个变量之间的关系。这种技术通常用于预测分析，时间序列模型以及发现变量之间的因果关系。举个简单的例子，二元回归就是学习自变量x和因变量y之间的关系，当给定x之后能够大致估测出y的数值。当数据的维数比较高的时候，学习各个变量之间的关系的过程为多元回归。其中，自变量即相当于机器学习中的特征，因变量相当于机器学习中的标签。

线性回归模型

假定有$m$个样本，其特征包含了$d$个维度，第$i$样本可以表示为 x i = { x i 1 , x i 2 , x i 3 , x i 4 . . . . . . . x i d } \bm{x_i}=\{x_{i1},x_{i2},x_{i3},x_{i4}.......x_{id}\} xi​={xi1​,xi2​,xi3​,xi4​.......xid​}，对应的因变量为$y_i$。

所谓线性回归，即为特征与标签之间的关系为线性关系。这里假设多元特征${\mathbf{x}}_{\mathbf{i}}$和$y_i$之间的关系式$f({\mathbf{x}}_{\mathbf{i}})$
$$f({\mathbf{x}}_{\mathbf{i}})={\mathbf{\omega }}^{\mathbf{T}}{\mathbf{x}}_{\mathbf{i}}+b,使得f({\mathbf{x}}_{\mathbf{i}})近似等于{\mathbf{y}}_{\mathbf{i}},称为“\mathbf{多元线性回归}"。$$上面的表达式中$\mathbf{\omega }$和$b$就是在回归分析中需要学习的参数。回归的目的就是使得$f({\mathbf{x}}_{\mathbf{i}})$近似等于${\mathbf{y}}_{\mathbf{i}}$，这一点很重要，对于参数的学习也就是从这里来的。

最小二乘法估计参数

所以如何获取这两个参数呢？具体地，我们可以从使得预测值和真实值相等着手，也就是误差。

定义损失函数

误差越小，一般就意味着我们学习得到的模型越好 。误差可以用预测值和真实值之间的差距表示，即误差可以表示为$\sum _{i=1}^m|f({\mathbf{x}}_{\mathbf{i}})-{\mathbf{y}}_{\mathbf{i}}|$。就如同图中预测值和真实值之差一样，使其最小。图源知乎

为了运算方便，进一步地我们可以表示为差值的平方，这样就得到了新的误差函数$$E=\sum _{i=1}^m(f({\mathbf{x}}_{\mathbf{i}})-{\mathbf{y}}_{\mathbf{i}}{)}^2$$
这个表达式，称为“均方误差”。

矩阵形式的均方误差

把矩阵$\mathbf{\omega }$和$b$合在一起，得到一个$d+1$行的列向量$\hat{\mathbf{\omega }}=(\mathbf{\omega };b)$。同样把数据集的矩阵表示为$m\times (d+1)$的矩阵$X$，最后一列补充为1。具体如下

标签也用一个列向量表示$\mathbf{y}=(y_i;y_2;......;y_m)$。这样根据上面所说的，我们的均方误差为$E=(\mathbf{y}-\mathbf{X}\hat{\mathbf{\omega }}{)}^T(\mathbf{y}-\mathbf{X}\hat{\mathbf{\omega }})$，只需要求得$\hat{\mathbf{\omega }}$使得E最小即可。

求解

上面的关于$\hat{\mathbf{\omega }}$的函数是一个凸函数，具有最小值（可以自己求海森矩阵试试，这里不再证明）。我们对$\hat{\mathbf{\omega }}$求导，得到$$\frac{\partial E}{\partial \hat{\mathbf{\omega }}}=2{\mathbf{X}}^T(\mathbf{X}\hat{\mathbf{\omega }}-\mathbf{y})$$
令上面的式子为0，就可以得到最优取值点${\hat{\mathbf{\omega }}}^{\ast }$。当 ${\mathbf{X}}^{\mathbf{T}}\mathbf{X}$ 为满秩矩阵或正定矩阵时，其结果为$${\hat{\mathbf{\omega }}}^{\ast }=({\mathbf{X}}^{\mathbf{T}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathbf{T}}\mathbf{y}$$
令${\hat{\mathbf{x}}}_{\mathbf{i}}=({\mathbf{x}}_{\mathbf{i}};1)$因此我们最终学习得到的线性回归模型为$$f({\hat{\mathbf{x}}}_{\mathbf{i}})={\hat{\mathbf{x}}}_{\mathbf{i}}^{\mathbf{T}}({\mathbf{X}}^{\mathbf{T}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathbf{T}}\mathbf{y}$$

$({\mathbf{X}}^{\mathbf{T}}\mathbf{X}{)}^{-1}$的存在问题

不满秩

$\mathbf{X}$是一个$m\times (d+1)$的矩阵，所以${\mathbf{X}}^{\mathbf{T}}\mathbf{X}$应该是一个$d+1$阶的方阵。如果$({\mathbf{X}}^{\mathbf{T}}\mathbf{X}{)}^{-1}$是存在的，则意味着矩阵${\mathbf{X}}^{\mathbf{T}}\mathbf{X}$满秩。那么${\mathbf{X}}^{\mathbf{T}}\mathbf{X}$是否满秩呢，即秩为$d+1$？显然是可能不满秩的，有两种情况。

• $d+1>m$:当特征值维度大于样本数目的时候，对于矩阵$\mathbf{X}$来说，它满秩即为$m$，而两者相乘的情况下,秩最多也就是m，所以不可能满秩。（原因： r a n k A B ≤ m i n { r a n k A , r a n k B } rankAB\leq min\{rankA,rankB\} rankAB≤min{rankA,rankB}）
• 矩阵$\mathbf{X}$的列秩小于$d+1$:在这种情况下，乘积得到的矩阵${\mathbf{X}}^{\mathbf{T}}\mathbf{X}$的秩小于d+1，也是不满秩的。
  上述两种情况，第一种情况发生的可能性不大，因为数据一般比较多，第二种发生的可能性偏大一些。

取值无法计算出来

一般在训练的时候样本数目是比较多的，这就意味着矩阵${\mathbf{X}}^{\mathbf{T}}\mathbf{X}$的行列式会非常大，随着数据数量越多，行列式越大。而矩阵的逆为矩阵的伴随除以矩阵的行列式，这种情况下矩阵的逆就会变得很小，如果小于计算机浮点数所能表示的情况，那该矩阵就相当于也是不存在了。

代码实现

# @Date 2023/8/11 16:53
# @Discription 最小二乘法线性回归,二元数据回归
import numpy as np
import matplotlib.pyplot as plt
# 梯度下降法
def gradientDescent(iter, x, y, w, alpha):
    # w_trans = w.reshape(-1,1) # w转置,变成列向量,对于向量来说好像用.T和.Transpose都不太好用
    x_trans = x.transpose()
    for i in range(0, iter):
        pre = np.dot(x, w)  #m*2与2*1相乘得到为一个m行的列向量
        loss = (pre - y)
        gradient = 2 * np.dot(x_trans, loss)/m
        w = w - alpha * gradient
        print("第{}次得到的w为{}".format(i,w))
    return w


# 输入数据
# delimiter参数的作用是指定分隔符,dtype参数的作用是指定数据类型
data = np.genfromtxt('train.csv', delimiter=',', dtype=float, skip_header=1,)  # 将csv数据读取为numpy
print('df的数据类型为:{}'.format(type(data)))
x = data.copy()
m, n = np.shape(x)
x_data = np.ones((m, n))
x_data[:, :-1] = x[:, :-1]  # 扩展为特征矩阵加一列1
y_data = np.ones((m, 1))
y_data[:,0] = x[:, -1]  # 取出最后一列,即标签
m, n = np.shape(x_data)
w = np.ones((n,1))
# 这个下降的学习率可能和数据量有关，数据越多，则上面求出来的梯度越大，可能会导致w变化过大，无法学习得到结果
result = gradientDescent(1000, x_data, y_data, w, 0.0001)
y_pre = np.dot(x_data, result)  # 699×2乘以2×1得到的是m行的列向量
print("线性回归模型 w: ", result)

plt.rc('font', family='FangSong', size=14)
plt.scatter(x[:, 0], x[:, 1], color='b', s=6, label='tran data')
plt.plot(x[:, 0], y_pre, color='r', label='forecast data')
plt.xlabel('x')
plt.ylabel('y')
plt.title('线性回归预测(梯度下降)')
plt.legend()
plt.savefig('result.png', dpi=300)
plt.show()

正由于上面所述问题，所以在实际的操作中很少有人直接通过求极值点的方法获得拟合参数，更多的方法是使用梯度下降法。梯度下降法的核心就是选择按照当前的导数方向搜索更优点，直至得到一个稳定的结果。在上面的程序段里，因为问题比较简单，所以我直接选择了迭代1000次。

岭回归

上面说到，在矩阵$\mathbf{X}$的秩小于$d+1$的时候，$({\mathbf{X}}^{\mathbf{T}}\mathbf{X}{)}^{-1}$肯定是不存在的。换句话说，就是当矩阵$\mathbf{X}$中存在两个特征具有线性关系的时候，用最小二乘法是无法得到最优的${\hat{\mathbf{\omega }}}^{\ast }$的，往往得到的结果会偏离正常现象。那么，针对于特征存在共线的情况，岭回归往往会得到一更优结果。

岭回归原理

由于共线问题导致的是$|{\mathbf{X}}^{\mathbf{T}}\mathbf{X}|=0$，同时根据${\hat{\mathbf{\omega }}}^{\ast }=({\mathbf{X}}^{\mathbf{T}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathbf{T}}\mathbf{y}$，我们可以考虑通过放弃一点最小二乘法的无偏性，以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归模型。具体操作就是给${\mathbf{X}}^{\mathbf{T}}\mathbf{X}$加上一个$kI$,这样子得到的新的模型参数为$({\mathbf{X}}^{\mathbf{T}}\mathbf{X}+\mathbf{kI}{)}^{-1}{\mathbf{X}}^{\mathbf{T}}\mathbf{y}$,可以记为$${\mathbf{\beta }}^{\ast }=({\mathbf{X}}^{\mathbf{T}}\mathbf{X}+\mathbf{kI}{)}^{-1}{\mathbf{X}}^{\mathbf{T}}\mathbf{y}$$
由于仅仅加了一个$kI$，一方面对精度损失不太多，另一方面求解的时候也能够得到一个较为稳定的模型，这种方法在实际情况下应用很多。
根据上面的式子反推我们相当于在一开始的损失函数中引入了一个惩罚项，称为L2正则化。新的损失函数为$$\hat{E}=\sum _{i=1}^m(f({\mathbf{x}}_{\mathbf{i}})-{\mathbf{y}}_{\mathbf{i}}{)}^2+k\sum _{j=1}^{d+1}{\beta }_j^2$$$$=||{\mathbf{\beta }}^{\mathbf{T}}\mathbf{X}|{|}^2+k||\mathbf{\beta }|{|}^2$$
这里的$k$称之为超参数。