线性回归模型

日常生活中，我们往往会碰到连续特征的预测问题。例如房价预测，股价预测等。这种类型的问题统称“回归”问题。解决这类问题的学习方法是有监督的。设有一个数据集$D$，其特征$\mathbf{x}$的维度为$d$，输出为$y$，$y$是连续型的。则回归问题致力于寻找一个函数$f$以建立$\mathbf{x}$到$y$的映射$y=f(\mathbf{x})$。

常见的回归模型分为==线性模型==和==非线性模型==两种。其中线性模型形式简单，易于建模。非线性模型可在线性模型的基础上通过引入层级结构或高维映射得到。本小节我们先讨论线性模型。

基本形式

设输入数据$\mathbf{x}=(x_1,x_2,\cdots,x_d)$，则线性模型是通过将输入数据的属性进行线性组合来得到的模型：

$$
f(\boldsymbol{x})=w_{1} x_{1}+w_{2} x_{2}+\ldots+w_{d} x_{d}+b
$$

我们也可以用向量表示：

$$
f(\boldsymbol{x})=\boldsymbol{w}^{\mathrm{T}} \boldsymbol{x}+b
$$

该模型的学习目的就是学得参数$\mathbf{w}$和$b$，使得$f(\boldsymbol{x})$尽可能的与$y_{i}$相等。

线性回归

刚刚我们说到线性模型的学习目的就是学得参数$\mathbf{w}$和$b$。那么如何衡量$f(\boldsymbol{x})$与$y_{i}$的差别呢？在前一篇我们已经提到过回归任务中的性能度量——均方误差。因此我们可以通过最小化均方误差来让模型参数达到最优。即

$$
\begin{aligned}
\left(w^{}, b^{}\right) &=\underset{(w, b)}{\arg \min } \sum_{i=1}^{m}\left(f\left(x_{i}\right)-y_{i}\right)^{2} \\
&=\underset{(w, b)}{\arg \min } \sum_{i=1}^{m}\left(y_{i}-w x_{i}-b\right)^{2}
\end{aligned}
$$

实际上均方误差对应于“欧氏距离”。而基于均方误差最小化的模型求解策略称为“最小二乘法”。通俗理解最小二乘法在线性模型中就是找到一条直线能够使所有样本到该直线的欧式距离之和最小。求解参数的过程称为最小二乘的“参数估计”。设$E_{(w,b)}=\sum^m_{i=1}(y_i-wx_{i}-b)$。求解即将$E_{(w,b)}$分别对$w$和$b$求导：

$$
\begin{array}{l}
\frac{\partial E_{(w, b)}}{\partial w}=2\left(w \sum_{i=1}^{m} x_{i}^{2}-\sum_{i=1}^{m}\left(y_{i}-b\right) x_{i}\right), \\
\frac{\partial E_{(w, b)}}{\partial b}=2\left(m b-\sum_{i=1}^{m}\left(y_{i}-w x_{i}\right)\right)
\end{array}
$$

然后令二者为零即可解得关于$w$和$b$的最优闭式解：

$$
w=\frac{\sum_{i=1}^{m} y_{i}\left(x_{i}-\bar{x}\right)}{\sum_{i=1}^{m} x_{i}^{2}-\frac{1}{m}\left(\sum_{i=1}^{m} x_{i}\right)^{2}}
$$

$$
b=\frac{1}{m} \sum_{i=1}^{m}\left(y_{i}-w x_{i}\right)
$$

注意：这里的$E_{(w,b)}$实际上是关于$w$，$b$的凸函数，这样求导令为零才能得到闭式解。式中$\overline{x}$为输入样本均值。

以上为一元线性回归的参数求解。那么在多元线性回归中，输入特征的维度由一维增加到了$d$维。其模型为：

$$
f(\boldsymbol{x})=\boldsymbol{w}^{\mathrm{T}} \boldsymbol{x}+b
$$

求解方法与一元线性回归类似。为了便于讨论，我们将$\boldsymbol{w}$和$b$合并为$\hat{\boldsymbol{w}}$。设数据集为$m\cdot(d+1)$大小的矩阵$\mathbf{X}$。即：

$$
\mathbf{X}=\left(\begin{array}{ccccc}
x_{11} & x_{12} & \ldots & x_{1 d} & 1 \
x_{21} & x_{22} & \ldots & x_{2 d} & 1 \
\vdots & \vdots & \ddots & \vdots & \vdots \
x_{m 1} & x_{m 2} & \ldots & x_{m d} & 1
\end{array}\right)=\left(\begin{array}{cc}
\boldsymbol{x}{1}^{\mathrm{T}} & 1 \
\boldsymbol{x}{2}^{\mathrm{T}} & 1 \
\vdots & \vdots \
\boldsymbol{x}_{m}^{\mathrm{T}} & 1
\end{array}\right)
$$

设标记为$\boldsymbol{y}=\left(y_{1} ; y_{2} ; \ldots ; y_{m}\right)$，则：

$$
\hat{\boldsymbol{w}}^{*}=\underset{\hat{\boldsymbol{w}}}{\arg \min }(\boldsymbol{y}-\mathbf{X} \hat{\boldsymbol{w}})^{\mathrm{T}}(\boldsymbol{y}-\mathbf{X} \hat{\boldsymbol{w}})
$$

==【推导过程】==如下：

先将$E_{\hat{\boldsymbol{w}}}=(\boldsymbol{y}-\mathbf{X} \hat{\boldsymbol{w}})^{\mathrm{T}}(\boldsymbol{y}-\mathbf{X} \hat{\boldsymbol{w}})$展开，得到

$$
E_{\hat{\boldsymbol{w}}}=\boldsymbol{y}^{\mathrm{T}} \boldsymbol{y}-\boldsymbol{y}^{\mathrm{T}} \mathbf{X} \hat{\boldsymbol{w}}-\hat{\boldsymbol{w}}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \boldsymbol{y}+\hat{\boldsymbol{w}}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \mathbf{X} \hat{\boldsymbol{w}}
$$

对$\hat{\boldsymbol{w}}$求导即得：

$$
\frac{\partial E_{\hat{\boldsymbol{w}}}}{\partial \hat{\boldsymbol{w}}}=\frac{\partial \boldsymbol{y}^{\mathrm{T}} \boldsymbol{y}}{\partial \hat{\boldsymbol{w}}}-\frac{\partial \boldsymbol{y}^{\mathrm{T}} \mathbf{X} \hat{\boldsymbol{w}}}{\partial \hat{\boldsymbol{w}}}-\frac{\partial \hat{\boldsymbol{w}}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \boldsymbol{y}}{\partial \hat{\boldsymbol{w}}}+\frac{\partial \hat{\boldsymbol{w}}^{\mathrm{T}} \mathbf{X}^{\mathrm{T}} \mathbf{X} \hat{\boldsymbol{w}}}{\partial \hat{\boldsymbol{w}}}
$$

由矩阵微分公式$\frac{\partial \boldsymbol{a}^{\mathrm{T}} \boldsymbol{x}}{\partial \boldsymbol{x}}=\frac{\partial \boldsymbol{x}^{\mathrm{T}} \boldsymbol{a}}{\partial \boldsymbol{x}}=\boldsymbol{a}, \frac{\partial \boldsymbol{x}^{\mathrm{T}} \mathbf{A} \boldsymbol{x}}{\partial \boldsymbol{x}}=\left(\mathbf{A}+\mathbf{A}^{\mathrm{T}}\right) \boldsymbol{x}$可得：

$$
\begin{array}{c}
\frac{\partial E_{\hat{\boldsymbol{w}}}}{\partial \hat{\boldsymbol{w}}}=0-\mathbf{X}^{\mathrm{T}} \boldsymbol{y}-\mathbf{X}^{\mathrm{T}} \boldsymbol{y}+\left(\mathbf{X}^{\mathrm{T}} \mathbf{X}+\mathbf{X}^{\mathrm{T}} \mathbf{X}\right) \hat{\boldsymbol{w}} \
\frac{\partial E_{\hat{\boldsymbol{w}}}}{\partial \hat{\boldsymbol{w}}}=2 \mathbf{X}^{\mathrm{T}}(\mathbf{X} \hat{\boldsymbol{w}}-\boldsymbol{y})
\end{array}
$$

令上式为零，可以得到最优参数的闭式解：

$$
\hat{\boldsymbol{w}}^{*}=\left(\mathbf{X}^{\mathrm{T}} \mathbf{X}\right)^{-1} \mathbf{X}^{\mathrm{T}} \boldsymbol{y}
$$

当特征数量$d$大于训练样本数$n$时，以上闭式解很容易产生过拟合的问题。但是我们可以通过正则化来解决该问题。下面我们来讨论一下两种常见的线性回归正则化方法——岭回归和LASSO回归。

线性回归正则化

所谓正则化，就是通过在模型中加入一些惩罚项或约束条件来控制模型的复杂度，以达到减轻过拟合的目的。岭回归（$ridge\ regression$）和LASSO是目前最为流行的两种正则化方法。

岭回归&LASSO

根据定义，岭回归和LASSO分别对应$l_2$和$l_1$正则化，其对$\boldsymbol{w}$所做出的先验限制分别为$||\boldsymbol{w}||_2\le C$和$||\boldsymbol{w}||_1\le C$，$C$为需要我们预先设定的常数。由此问题转为约束优化问题，对于岭回归有如下：

$$
\min {\boldsymbol{w}}|\boldsymbol{y}-\boldsymbol{X} \boldsymbol{w}|{2}^{2}, \quad \text { s.t. } \quad|\boldsymbol{w}|_{2} \leqslant C
$$

对于LASSO，则有如下：

$$
\min {w} \quad|\boldsymbol{y}-\boldsymbol{X} \boldsymbol{w}|{2}^{2}, \quad \text { s.t. } \quad|\boldsymbol{w}|_{1} \leqslant C
$$

针对此类约束问题我们可以利用拉格朗日乘子法来解决。

关于拉格朗日乘子法可以参考我的这篇博客：
链接: https://ylxy123.github.io/2021/05/06/优化算法（一）拉格朗日乘子法/
来源: Welcom to YLXY’s blog.
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

通过该方法引入拉格朗日乘子来构造拉格朗日函数$L_1(\boldsymbol{X},\lambda)=||\boldsymbol{y}-\boldsymbol{X} ||^2_2+\lambda||\boldsymbol{w}||^2_2$和$L_2(\boldsymbol{X},\lambda)=||\boldsymbol{y}-\boldsymbol{X} ||^2_2+\lambda||\boldsymbol{w}||_1$。

原问题即转化为如下无约束带惩罚函数的优化问题：

$$
\min {\boldsymbol{w}}|\boldsymbol{y}-\boldsymbol{X} \boldsymbol{w}|{2}^{2}+\lambda|\boldsymbol{w}|_{2}^{2}
$$

和

$$
\min {\boldsymbol{w}} \frac{1}{2}|\boldsymbol{y}-\boldsymbol{X} \boldsymbol{w}|{2}^{2}+\lambda|\boldsymbol{w}|_{1}
$$

其中系数$\lambda>0$是依赖$C$的常数。

针对岭回归问题，对如上目标函数直接求导并令梯度为零可得：

$$
\begin{aligned}
\boldsymbol{w}^{\text {ridge }} &=\underset{\boldsymbol{w}}{\operatorname{argmin}}\left(|\boldsymbol{y}-\boldsymbol{X} \boldsymbol{w}|{2}^{2}+\lambda|\boldsymbol{w}|{2}^{2}\right) \
&=\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{X}+\lambda \boldsymbol{I}\right)^{-1} \boldsymbol{X}^{\mathrm{T}} \boldsymbol{y}
\end{aligned}
$$

通过与多元线性回归的最小二乘闭式解相比，我们可以发现岭回归得到的解只多了一个正则项$\lambda\boldsymbol{I}$，该项的存在使得$\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{X}+\lambda \boldsymbol{I}\right)^{-1}$在数值计算上的表现更加稳定。特别是当多重共线性的情况发生时，即使$\boldsymbol{X}^{\mathrm{T}}\boldsymbol{X}$接近奇异，岭回归还是能产生较稳定的结果。我们甚至还能通过观察岭迹表来剔除具有多重共线性的特征。

与岭回归不同的是，LASSO无法得到解析解，但$\boldsymbol{w}^{LASSO}$是稀疏的，即有很多分量接近零。因此LASSO模型的解可以帮我们做特征选择。当某个$x_i$和$y$的相关性不高或与其他$x_j$存在多重共线性时，$w^{LASSO}$有很大的可能性为零。

正则化的本质实际上是当我们在对参数进行估计时做一个偏差和方差的平衡。不带正则项的线性回归模型得到的参数往往是最小方差的无偏估计。通俗地说，正则化是通过牺牲参数估计的无偏性来换取模型的效果的。岭回归和LASSO作为有偏估计，期望以无偏性为代价来获得更符合实际效果的回归系数。

说了这么多，下面我们来讲讲LASSO的求解方法。既然LASSO没办法得到解析解，那么我们可以通过数值迭代的方法来逼近真实解。在今天已经有很多方法可以用来求解LASSO，例如坐标下降、LARS算法和基于近似梯度方法的ISTA及FISTA算法（ISTA的加速版本）。在此我们简单讲讲ISTA算法：

我们已经知道对于光滑函数解的近似逼近最常用的方法是梯度下降法：
$$
\boldsymbol{w}^{(t+1)}=\boldsymbol{w}^{(t)}-\eta \nabla f\left(\boldsymbol{w}^{(t)}\right)
$$
其中$\eta$为学习率。为了方便理解，在这我先介绍一下近似点算法PPA的定义：设一闭凸函数$f$，有如下迭代格式：

$$
\begin{aligned}
x_{k+1} &=\operatorname{prox}{t_k f}\left(x{k}\right) \
&=\arg \min {u}\left(f(u)+\frac{1}{2 t{k}}\left|u-x_{k}\right|_{2}^{2}\right)
\end{aligned}
$$

这个式子的意思是寻找一个距离$x_k$不太远的一个点$u$，使得$f(u)$尽可能小。即$f(u)\le f(x_k)$。

接下来介绍近似梯度优化方法。设我们待优化的目标函数为$F(x)=g(x)+h(x)$，其中$g(x)$为连续可微函数，$h(x)$不连续。LASSO目标函数恰好符合该条件，因此可以用近似梯度优化求解：

近似优化算法（Proximal Gradient Algorithm）
1
2
3
4
5
for t in range(n):
    v_t = x_t - lr*grad(F(x_t))  # grad(·)为梯度函数，v_t为沿梯度方向找到的一点，lr为学习率，F(x)为目标函数
    x_t = prox_h(lambda, v_t)  # 使用prox式来优化h(x)
    if 收敛或达到最大迭代次数:
        break
注意算法中prox_h(x)函数中参数$\lambda$的选择。一般取$\lambda\in(0,\frac{1}{L})$。若$L$未知，可通过线搜索寻找：
1
2
3
4
5
while True:
    z = prox_h(lambda, v_t)
    if 终止条件:
        break
    lambda *= 1/2
其中终止条件为$F(z)<=F\left(v^{t}\right)+\nabla F^\text{T}\left(v^{t}\right)\left(v^{2}-z\right)+\frac{1}{2 \lambda}\left|v^{t}-z\right|^{2}$

算法最终输出$x^{(t+1)}=z$