手撕统计机器学习 - LASSO 变量选择

本节所需前置知识：

• 线性回归
• 最优化方法（二次规划问题、Lagrange 函数与 KKT 条件、次梯度）

---

高维小样本问题

高维小样本问题是在 $n$ 个样本，$p$ 个维度的情况下，有 $n<<p$ 的前提的问题. 这类问题在特殊场景下非常常见，包括：

• 图片模式识别（每个像素视为一个特征）；
• 生物信息（例如 DNA 甲基化数据，单个节点均视为一个特征）；

同时，在一般的情况下，我们也想知道什么样的特征是比较重要的，对于一个背包，影响它价格的因素里面，品牌的影响更大还是背包容量的影响更大？我们因此希望剔除掉不重要的特征，提高我们的预测准确率.

在线性回归当中，我们得到的系数：
$$y={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_p{x}_p$$
实际上就代表着各个特征的重要程度，当 ${\beta }_i>0$ 时，说明第 $i$ 个特征具有正影响，绝对值越大说明特征影响越大.

---

变量选择问题

现在我们要进行变量的选择，也就是说，我们想筛除掉一些没那么重要的特征，应该如何解决这个问题？既然线性回归的系数代表着重要性，那么如果系数为 $0$ ，就说明这个变量不被选择，反之则被选择，这个时候我们就进行了变量选择.

一个简单粗暴的想法是：我们先进行 OLS，然后设定一个阈值，当系数绝对值小于这个阈值时，我们就强行压缩到 $0$ ，这样就能实现变量选择. 例如：
$$y=0.6+0.1x_1+0.8x_2-0.4x_3$$
如果设定阈值为 $0.2$ ，那么我们就选出了后两个特征.

上述的方法看似简单，但是有个问题就是：我们到底要如何控制选择的力度？阈值选择到什么程度才能尽可能地完成变量选择的同时保障 $\mathrm{MSE}$ 还不会有明显增大？

---

稀疏优化问题

为了解决刚刚提到的问题，我们不妨将原来的损失函数进行改造：
$$L:=\mathrm{MSE}+\lambda \Vert \mathbf{\beta }{\Vert }_0$$
这里 $\lambda$ 为惩罚参数，${\ell }_0$ 范数表示 $\mathbf{\beta }$ 当中非零元素的个数，这也就是说，我们既希望 $\mathbf{\beta }$ 足够稀疏，又希望 $\mathrm{MSE}$ 不会过大，那么这个时候只要 $\lambda$ 选的足够好，就能达到我们希望的效果. 这种希望解能稀疏的问题我们称为稀疏优化问题.

看起来这个问题很简单，但是以 $L$ 作为损失函数有个致命的问题在于：我们很难求解！我们当然可以预设 $\mathbf{\beta }$ 哪些元素为 $0$ ，然后进行优化，但是所有的选择方案加在一起将会是一个组合的复杂度：
$$\left(\genfrac{}{}{0ex}{}{p}{1}\right)+\left(\genfrac{}{}{0ex}{}{p}{2}\right)+\cdots +\left(\genfrac{}{}{0ex}{}{p}{p-1}\right)+\left(\genfrac{}{}{0ex}{}{p}{p}\right)$$
这个时间就会爆炸，因此这不是一个合适的计算方案. 同时，${\ell }_0$ 范数作为泛函，也不是凸的，因此常用的凸优化方法都无法找到其全局最优解.

---

LASSO

LASSO (Least Absolute Shrinkage and Selection Operator) 就是用来解决上述问题的，它最早由 Robert Tibshirani 于 1996 年提出，具有划时代的意义. 它的改动很简单，就是将刚刚所说的 ${\ell }_0$ 范数损失函数以及岭回归的损失函数改为：
$$L=\mathrm{MSE}+\lambda \Vert \mathbf{\beta }{\Vert }_1$$
也就是说，我们将其改为了 ${\ell }_1$ 范数，一个 ${\ell }_0$ 范数的凸近似. 也可以写为如下的形式：
$$\begin{array}{rl}& \min \mathrm{MSE}\\ & s.t.\sum _{i=1}^p|{\beta }_i|<t\end{array}$$

---

正交设计下的 LASSO

首先解释标题，正交设计是统计学试验设计当中的概念，表示数据 $\mathbf{X}$ 是正交矩阵，满足 ${\mathbf{X}}^{\mathrm{T}}\mathbf{X}=\mathbf{I}$ . 在这个条件下，LASSO 将会更明显地展示出其性质.

对 $\mathrm{MSE}$ （在原论文当中使用残差平方和 $\mathrm{RSS}$ 代替，由于两者仅差一个 $\frac{1}{n}$ ，因此这里不更改），可以将其展开为
$$\mathrm{MSE}=(\mathbf{X\beta }-\mathbf{y}{)}^{\mathrm{T}}(\mathbf{X\beta }-\mathbf{y})={\mathbf{\beta }}^{\mathrm{T}}{\mathbf{X}}^{\mathrm{T}}\mathbf{X}\mathbf{\beta }-{\mathbf{y}}^{\mathrm{T}}\mathbf{X\beta }-{\mathbf{\beta }}^{\mathrm{T}}{\mathbf{X}}^{\mathrm{T}}\mathbf{y}+{\mathbf{y}}^{\mathrm{T}}\mathbf{y}$$
现在这一步可以看到，第四项和 $\mathbf{\beta }$ 无关，而 ${\mathbf{X}}^{\mathrm{T}}\mathbf{X}$ 会影响优化问题的难度，最简单的时候，就是正交设计的时候. 此时不妨使用正交设计的条件继续，有
$$\mathrm{MSE}={\mathbf{\beta }}^{\mathrm{T}}\mathbf{\beta }-2{\mathbf{\beta }}^{\mathrm{T}}{\mathbf{X}}^{\mathrm{T}}\mathbf{y}+{\mathbf{y}}^{\mathrm{T}}\mathbf{y}$$
假设 OLS 已经计算出系数 ${\mathbf{\beta }}^{\circ }$ ，那么根据 OLS 系数的公式，有 ${\mathbf{\beta }}^{\circ }={\mathbf{X}}^{\mathrm{T}}\mathbf{y}$ ，因此
$$\mathrm{MSE}=\Vert \mathbf{\beta }{\Vert }_2^2-2{\mathbf{\beta }}^{\mathrm{T}}{\mathbf{\beta }}^{\circ }$$

---

Lagrange 函数的应用

此时利用 Lagrange 函数将约束优化转化为无约束优化：
$$\mathcal{L}(\mathbf{\beta },\gamma )=\Vert \mathbf{\beta }{\Vert }_2^2-2{\mathbf{\beta }}^{\mathrm{T}}{\mathbf{\beta }}^{\circ }+\gamma \left(\sum _i|{\beta }_i|-t\right)$$
对分量求导，我们可以发现有 $|{\beta }_i|$ 项不好求导，此时应该分情况求导得到 $\mathrm{sign}({\beta }_i)$ （尽管这里 $0$ 处是尖点，但是考虑到系数为 $0$ 就说明变量被筛出，因此可以直接视为 $0$），随后有：
$$\frac{\partial \mathcal{L}}{\partial {\beta }_i}=2{\beta }_i-2{\beta }_i^{\circ }+\gamma \mathrm{sign}({\beta }_i)$$
令上式为 $0$ ，有
$${\beta }_i={\beta }_i^{\circ }-\frac{\gamma }{2}\mathrm{sign}({\beta }_i)$$
其中 $\gamma ⩾0$，此时要分情况讨论，当 ${\beta }_i^{\circ }>0$ 时，要保证 ${\beta }_i>0$ 必须 ${\beta }_i^{\circ }>\frac{\gamma }{2}$ ，对 ${\beta }_i^{\circ }<0$ 时讨论也是类似的，最终得到
$${\beta }_i=\mathrm{sign}({\beta }_i^{\circ }){\left({\beta }_i^{\circ }-\frac{\gamma }{2}\right)}^{+}$$
这就是正交设计下的系数. 其中 $\gamma$ 利用 KKT 条件进行确定.

---

软阈值解

对于最终的系数：
$${\beta }_i=\mathrm{sign}({\beta }_i^{\circ }){\left({\beta }_i^{\circ }-\frac{\gamma }{2}\right)}^{+}$$
可以看到，当 ${\beta }_i\ne 0$ 时，第 $i$ 个特征被选择出来，此时 ${\beta }_i$ 和 ${\beta }_i^{\circ }$ 仅差一个常数，这个解我们称为软阈值解. 对于特征，软阈值解只有两种行为：选择和收缩，并且由于收缩的力度 $\gamma$ 是一致的，所以软阈值解可以看为下图所示的函数：

---

几何意义

在二维情况下，我们可以研究 LASSO 的几何意义，有
$$\mathrm{MSE}=(\mathbf{\beta }-{\mathbf{\beta }}^{\circ }{)}^{\mathrm{T}}{\mathbf{X}}^{\mathrm{T}}\mathbf{X}(\mathbf{\beta }-{\mathbf{\beta }}^{\circ })+C$$
其中 $C$ 为与 $\mathbf{\beta }$ 无关的常数，上述二次型可看作二维下的椭圆等高线，由于 $p=2$ ，可以以 ${\beta }_1,{\beta }_2$ 为两个轴绘制图像：

这个就是 LASSO 几乎最经典的图像了，${\ell }_1$ 范数的单位圆顶点是尖点，在特征选择之后，如果有特征被筛除，就会类似 (a) 当中与尖点相切，此时 ${\beta }_1=0$ ，符合特征筛选的几何直觉.

---

非正交设计下的 LASSO

然而，${\mathbf{X}}^{\mathrm{T}}\mathbf{X}=\mathbf{I}$ 肯定是一个过于苛刻的条件，绝大多数时候都不会是正交的，那么此时又该如何求解？事实上我们就要回到
$$L=\mathrm{MSE}+\lambda \Vert \mathbf{\beta }{\Vert }_1$$
损失函数的最小化上面了，这个函数最小化的解没有解析解，因此需要用数值计算的方法来求解其最终选择结果.

我们接下来的思路为：

1. 证明该问题为凸优化问题，而非之前所说的组合优化问题.
2. 引入优化方法：坐标轴下降法和最小角回归法.

---

前置准备：凸函数证明

即证明
$$L(\alpha {\mathbf{\beta }}_1+(1-\alpha ){\mathbf{\beta }}_2)⩽\alpha L({\mathbf{\beta }}_1)+(1-\alpha )L({\mathbf{\beta }}_2)$$
左侧展开有
$$L(\alpha {\mathbf{\beta }}_1+(1-\alpha ){\mathbf{\beta }}_2)=\frac{1}{n}{\left(\sum _{i=1}^n{\left(\alpha {\mathbf{\beta }}_1+(1-\alpha ){\mathbf{\beta }}_2\right)}^{\mathrm{T}}{\mathbf{x}}_i-y_i\right)}^2+\lambda \Vert \alpha {\mathbf{\beta }}_1+(1-\alpha ){\mathbf{\beta }}_2{\Vert }_1$$
右侧有
$$\alpha L({\mathbf{\beta }}_1)+(1-\alpha )L({\mathbf{\beta }}_2)=\frac{\alpha }{n}{\left(\sum _{i=1}^n{\mathbf{\beta }}_1{\mathbf{x}}_i-y_i\right)}^2+\frac{1-\alpha }{n}{\left(\sum _{i=1}^n{\mathbf{\beta }}_2{\mathbf{x}}_i-y_i\right)}^2+\alpha \lambda \Vert {\mathbf{\beta }}_1{\Vert }_1+(1-\alpha )\lambda \Vert {\mathbf{\beta }}_2{\Vert }_1$$
左侧减右侧有
$$\lambda (\Vert \alpha {\mathbf{\beta }}_1+(1-\alpha ){\mathbf{\beta }}_2{\Vert }_1-\alpha \Vert {\mathbf{\beta }}_1\Vert -(1-\alpha )\Vert {\mathbf{\beta }}_2\Vert )$$
${\ell }_1$ 范数是凸范数，因此上面的式子 $⩽0$ ，从而我们证明了凸函数.

---

坐标轴下降法

简单思路

坐标轴下降法的思路非常简单：

• 初始解：$\hat{\mathbf{\beta }}=0$ .
• 在未收敛的情况下，循环执行：
  • 选取坐标位置 $j$ .
  • 解优化子问题：${\hat{\beta }}_j\leftarrow \arg {\min }_{{\beta }_j}L({\hat{\beta }}_1,{\hat{\beta }}_2,\cdots ,{\hat{\beta }}_{j-1},{\beta }_j,{\hat{\beta }}_{j+1},\cdots ,{\hat{\beta }}_p)$ .

也就是每次都挑选一个坐标，下降到最小值，然后挑选出来.

---

下降方法

可以看到，算法当中只有子问题的求解是关键的，我们按照分量形式展开：
$$\begin{array}{rl}L(\mathbf{\beta })& =\sum _{i=1}^n{\left(y_i-\sum _{j=1}^p{x}_{ij}{\beta }_j\right)}^2+\lambda \sum _{i=1}^p|{\beta }_i|\end{array}$$
此时对 ${\beta }_k$ 求导
$$\frac{\partial L}{\partial {\beta }_k}=-2\sum _{i=1}^n\left(y_i-\sum _{j=1}^p{x}_{ij}{\beta }_j\right)x_{ik}+\lambda \frac{\partial |{\beta }_k|}{\partial {\beta }_k}$$
现在有个不符合直觉的地方，$|{\beta }_k|$ 是有不可导点的，但是在优化上却是可行的，因此这里引入最优化当中次梯度的概念，详见附录内容. 求其次梯度有
$$\lambda \frac{\partial |{\beta }_k|}{\partial {\beta }_k}=\{\begin{array}{ll}\lambda ,& {\beta }_k>0\\ [-\lambda ,\lambda ],& {\beta }_k=0\\ -\lambda ,& {\beta }_k<0\end{array}$$
然后令偏导为 $0$ 即可达到想要的最小值.

---

${\beta }_k$ 的最优值

下面进行分类讨论，为方便解方程，将 ${\beta }_k$ 分离出来：
$$\begin{array}{rl}\frac{\partial L}{\partial {\beta }_k}& =-2\sum _{i=1}^n\left(y_i-\sum _{j=1}^p{x}_{ij}{\beta }_j\right)x_{ik}+\lambda \frac{\partial |{\beta }_k|}{\partial {\beta }_k}\\ & =-2\sum _{i=1}^n\left(y_i-\sum _{j\ne k}{x}_{ij}{\beta }_j\right)x_{ik}+2{\beta }_k\sum _{i=1}^n{x}_{ik}^2+\lambda \frac{\partial |{\beta }_k|}{\partial {\beta }_k}\\ & =r_k+z_k{\beta }_k+\frac{\partial (\lambda |{\beta }_k|)}{\partial {\beta }_k}\end{array}$$
这里 $r_k$ 对应第二个等号后的第一项，$z_k$ 对应 $2\sum _{i=1}^n{x}_{ik}^2$ .

当 ${\beta }_k>0$ 时，令上式为 $0$ ，解得 ${\beta }_k=\frac{r_k-\lambda }{z_k}$ ，当 ${\beta }_k<0$ 时，解得 ${\beta }_k=\frac{r_k+\lambda }{z_k}$ ，此时对 $r_k$ 也有相应要求. 最终有
$${\beta }_k=\{\begin{array}{ll}\frac{r_k-\lambda }{z_k},& r_k>\lambda \\ 0,& -\lambda ⩽r_k⩽\lambda \\ \frac{r_k+\lambda }{z_k},& r_k<-\lambda \end{array}$$

---

代码实现

这里还是先构造数据：

import numpy as np
import pandas as pd
from sklearn.datasets import make_regression

# 生成包含 100 个样本、4 个特征，噪声为 10 的数据集
X, y = make_regression(n_samples=50, n_features=4, noise=10, random_state=42)

然后就实现 LASSO ：

class Lasso:
    def __init__(self, alpha=0.1, fit_intercept=True, max_iter=1000, tol=1e-4):
        """
        Lasso回归实现（坐标下降法）
        
        参数:
        alpha: float, 正则化强度
        fit_intercept: bool, 是否拟合截距项
        max_iter: int, 最大迭代次数
        tol: float, 收敛阈值
        """

        self.alpha = alpha
        self.fit_intercept = fit_intercept
        self.max_iter = max_iter
        self.tol = tol
        self.coef_ = None
        self.intercept_ = 0.0 if fit_intercept else None
        self._is_fit = False
        self.n_iter_ = 0
    
    def soft_threshold(self, rho, lamda):
        """软阈值函数"""
        if rho < -lamda:
            return rho + lamda
        elif rho > lamda:
            return rho - lamda
        else:
            return 0.0
    
    def loss_(self, y_pred, y_true):
        RSS = np.sum((y_pred - y_true) ** 2)
        penalty = np.sum(np.abs(self.coef_))

        # RSS 乘以 0.5 是为了符合 sklearn 的结果.
        return 0.5 * RSS + self.alpha * penalty
    
    def fit(self, X, y):

        n_samples, n_features = X.shape
        
        # 初始化系数
        if self.coef_ is None:
            self.coef_ = np.zeros(n_features)
        
        # 如果需要截距项，对数据进行中心化
        if self.fit_intercept:
            X_mean = np.mean(X, axis=0)
            y_mean = np.mean(y)
            X_centered = X - X_mean
            y_centered = y - y_mean
        else:
            X_centered = X.copy()
            y_centered = y.copy()
        
        # 坐标下降
        prev_coef = self.coef_.copy()
        for iteration in range(self.max_iter):
            for j in range(n_features):

                r_j = y_centered - np.dot(X_centered, self.coef_) + self.coef_[j] * X_centered[:, j]
                rho_j = np.dot(X_centered[:, j], r_j)
                self.coef_[j] = self.soft_threshold(rho_j, self.alpha * n_samples) / np.sum(X_centered[:, j] ** 2)
            # 检查收敛
            coef_change = np.max(np.abs(self.coef_ - prev_coef))
            if coef_change < self.tol:
                break
                
            prev_coef = self.coef_.copy()
        
        # 计算截距
        if self.fit_intercept:
            self.intercept_ = y_mean - np.dot(X_mean, self.coef_)
        
        self.n_iter_ = iteration + 1
        self._is_fit = True
        
        return self
    
    def predict(self, X):
        """使用拟合的模型进行预测"""
        if not self._is_fit:
            raise ValueError("[ERROR] 模型尚未训练，请先调用fit方法！")
        
        if self.fit_intercept:
            return np.dot(X, self.coef_) + self.intercept_
        else:
            return np.dot(X, self.coef_)

---

然后分别用 scikit-learn 和自己的对比：

lasso = Lasso(alpha=50)
lasso.fit(X,y)
lasso.coef_

from sklearn import linear_model
model = linear_model.Lasso(alpha=50)
model.fit(X,y)
model.coef_

出来的结果会有些许小数点上的区别，但是筛选变量会是一致的. 结果的不一样其实也是因为 sklearn 官方使用的是最小角回归法进行优化的，这里读者可以尝试阅读别的资料来进行实现。