Lasso 回归与 Ridge 回归

Reference

• Lasso回归
• Lasso—原理及最优解
• 机器学习算法系列（五）- Lasso回归算法（Lasso Regression Algorithm）
• 岭回归
• 岭回归详解 从零开始 从理论到实践
• Tikhonov regularization 吉洪诺夫正则化（L2正则化）
• 机器学习算法系列（四）- 岭回归算法（Ridge Regression Algorithm）
• Lasso (statistics))
• Tikhonov regularization

【概述】

对于线性回归来说，其损失函数的解析解为：

$$\boldsymbol{\theta}=(XX^T)^{-1}XY$$

当 $XX^T$ 不是满秩矩阵时，没有逆矩阵，此时会出现过拟合现象，其本质就是自变量 $x^{(j)}$ 之间存在多重共线性（Multicollinearity），即自变量 $x^{(j)}$ 间存在相关关系

对于线性回归中的多重共线性问题，其中一个解决方法就是剔除掉共线的自变量，可以通过计算方差膨胀因子（Variance inflation factor，VIF）来量化自变量之间的相关关系，方差扩大因子越大说明自变量的共线性越严重

另一种解决方法是通过对损失函数正则化来解决，使用 L1 范数正则化的回归模型一般称为 Lasso 回归，使用 L2 范数正则化的回归模型一般称为 Ridge 回归

本文仅介绍 Lasso 回归与 Ridge 回归，关于多重共线性与方差膨胀因子，详见：多重共线性与方差膨胀因子

【正则化】

对于给定的容量为 $n$ 的样本集 $D=\{(\mathbf{x_1},y_1),(\mathbf{x_2},y_2),…,(\mathbf{x_n},y_n)\}$，第 $i$ 组样本中的输入 $\mathbf{x_i}$ 具有 $m$ 个特征值，即：$\mathbf{x_i}=(x_i^{(1)},x_i^{(2)},…,x_i^{(m)})^T\in \mathbb{R}^m$，输出为 $y_i$

用假设函数 $f(\mathbf{x_i};\boldsymbol{\theta})$ 来表示对第 $i$ 组数据的预测结果：

$$f(\mathbf{x_i};\boldsymbol{\theta})=\theta^{(0)} + \theta^{(1)} x_i^{(1)} + \theta^{(2)} x_i^{(2)} + ... + \theta^{(m)} x_i^{(m)}$$

其中，特征参数 $\boldsymbol{\theta}$ 为 $(m+1)\times 1$ 的列向量，即：

$$\boldsymbol{\theta}=[\theta^{(0)},\theta^{(1)},...,\theta^{(m)}]^T\in \mathbb{R}^{m+1}$$

现在，希望求出相应的 $\{\theta^{(i)}\}^{m+1}_{i=0}$ 来使得 $f(\mathbf{x_i};\boldsymbol{\theta})$ 能够尽量的拟合好样本集 $D$

为了表述方便，对假设函数进行简化，定义一个额外的第 $0$ 个特征量，这个特征量对所有样本的取值全部为 $1$，这使得特征量从过去的 $m$ 个变为 $m+1$ 个，即设：$x_i^{(0)}=1$

那么假设函数就可以写为：

$$f(\mathbf{x_i};\boldsymbol{\theta})=\theta^{(0)} x_i^{(0)} + \theta^{(1)} x_i^{(1)} + \theta^{(2)} x_i^{(2)} + ... + \theta^{(m)} x_i^{(m)}$$

选用残差平方和 RSS 作为损失函数，则有：

$$\begin{align} J_0(\boldsymbol{\theta}) &= \sum_{i=1}^n (f(\mathbf{x_i})-y_i)^2 \notag \\ &= \sum_{i=1}^n(\sum_{j=0}^{m+1}\theta_jx_{i}^{(j)}-y_i)^2 \notag \\ &= \sum_{i=1}^n(\boldsymbol{\theta}^T\mathbf{x_i}-y_i)^2 \notag \end{align}$$

对于原始的损失函数 $J_0(\boldsymbol{\theta})$，在加入正则化项 $R(\boldsymbol{\theta})$ 后，新的损失函数为：

$$J(\boldsymbol{\theta}) = J_0(\boldsymbol{\theta})+\lambda R(\boldsymbol{\theta}),\quad \lambda\geq 0$$

其中，$R(\boldsymbol{\theta})$ 是正则化项，$\lambda\geq0$ 为正则化参数，用于控制两者之间的平衡关系

如果 $\lambda$ 设置的过小，那么 $\boldsymbol{\theta}$ 参数的惩罚程度为起到应有的效果，仍会出现过拟合问题

如果 $\lambda$ 设置的过大，那么 $\boldsymbol{\theta}$ 参数的惩罚程度太大，所有的 $\boldsymbol{\theta}$ 参数几乎全部等于 $0$，这就使得把假设函数的除去 $\theta^{(0)}$ 的项全都忽略掉了，相当于直接用 $f(x)=\theta^{(0)}$ 来拟合数据，这就出现了欠拟合问题

因此，为了让正则化起到应用的效果，应当根据实际情况去选择一个合适的正则化参数 $\lambda$

【Lasso 回归】

一种对损失函数正则化的方式是将正则化项 $R(\boldsymbol{\theta})$ 取全部权重向量 $\boldsymbol{\theta}$ 的 L1 范数 $||\boldsymbol{\theta}||_1$，这个过程被称为 L1 正则化，经过这种损失函数正则化后的线性回归，被称为 Lasso 回归，其完整名称是最小绝对值收敛和选择算子（Least Absolute Shrinkage and Selection Operator）回归

$$\begin{align*} J(\boldsymbol{\theta}) &= \sum_{i=1}^n(\boldsymbol{\theta}^T\mathbf{x_i}-y_i)^2 + \lambda ||\boldsymbol{\theta}||_1\\ &= \sum_{i=1}^n(\boldsymbol{\theta}^T\mathbf{x_i}-y_i)^2 + \lambda\sum_{j=0}^m |\theta^{(j)}| \end{align*}$$

此时最优化目标是求损失函数 $J(\boldsymbol{\theta})$ 最小时 $\boldsymbol{\theta}$ 的大小，即：

$$\boldsymbol{\theta}^* = \arg \min_{\boldsymbol{\theta}} \Big( \sum_{i=1}^n(\boldsymbol{\theta}^T\mathbf{x_i}-y_i)^2 + \lambda\sum_{j=0}^m |\theta^{(j)}| \Big)$$

由于加入了 L1 范数，存在绝对值，因此损失函数 $J(\boldsymbol{\theta})$ 并不是处处可导的，无法通过直接求导的方式来直接得到 $\boldsymbol{\theta}$ 的解析解，通常使用坐标下降法（Coordinate Descent）或最小角回归法（Least Angle Regression，LARS），来求 $\boldsymbol{\theta}$ 的最优解

如图所示，展示了采用最小角回归法时惩罚系数 $\lambda$ 对各个自变量的权重系数的影响，图中横轴为惩罚系数 $\lambda$，纵轴为权重系数，每一个颜色表示一个自变量的权重系数

可以看到当 $\lambda$ 逐渐增大时（向左），某些特征的权重系数会快速变成零，通过这个性质说明 Lasso 回归可以用来进行特征选择，即通过控制 $\lambda$ 的大小来选择关键特征

• 关于坐标下降法的详细介绍，详见：坐标下降法
• 关于最小角回归法，详见：最小角回归法

【Ridge 回归】

最优化过程

除上述的 L1 正则化外，另一种对损失函数正则化的方式是将正则化项 $R(\boldsymbol{\theta})$ 取全部权重向量 $\boldsymbol{\theta}$ 的 L2 范数 $||\boldsymbol{\theta}||_2$ 的平方，这种方式被吉洪诺夫正则化（Tikhonov Regularization），又称为 L2 正则化，经过这种损失函数正则化后的线性回归，被称为 Ridge 回归，即岭回归

$$\begin{align*} J(\boldsymbol{\theta}) &= \sum_{i=1}^n(\boldsymbol{\theta}^T\mathbf{x_i}-y_i)^2 + \lambda ||\boldsymbol{\theta}||_2^2, \\ &= \sum_{i=1}^n(\boldsymbol{\theta}^T\mathbf{x_i}-y_i)^2 + \lambda \sum_{j=0}^{m}(\theta^{(j)})^2 \\ &= \sum_{i=1}^n(\boldsymbol{\theta}^T\mathbf{x_i}-y_i)^2 + \lambda \boldsymbol{\theta}^T\boldsymbol{\theta} \\ \end{align*}$$

此时最优化目标是求损失函数 $J(\boldsymbol{\theta})$ 最小时 $\boldsymbol{\theta}$ 的大小，即：

$$\boldsymbol{\theta}^* = \arg \min_{\boldsymbol{\theta}} \Big( \sum_{i=1}^n(\boldsymbol{\theta}^T\mathbf{x_i}-y_i)^2 + \lambda \boldsymbol{\theta}^T\boldsymbol{\theta} \Big)$$

对于损失函数 $J(\boldsymbol{\theta})$ 可以通过求导来直接得到 $\boldsymbol{\theta}$ 的解析解，即：

$$\boldsymbol{\theta}=(XX^T+\lambda I)^{-1} XY$$

可以看到岭回归的损失函数的解析解，相较于线性回归来说多了一个可以人为控制的对角矩阵 $I$，这时可以通过调整不同的 $\lambda$ 来使得矩阵 $XX^T+\lambda I$ 整体可逆

如图所示，展示了惩罚系数 $\lambda$ 对各个自变量的权重系数的影响，图中横轴为惩罚系数 $\lambda$，纵轴为权重系数，每一个颜色表示一个自变量的权重系数

可以看到当 $\lambda$ 越大时（向左），惩罚项的影响越来越大，逐渐占据主导地位，会使得每个自变量的权重系数趋近于零；当 $\lambda$ 越小时（向右），惩罚项的影响越来越小，会导致每个自变量的权重系数震荡的幅度变大

推导过程

将数据集 $D$ 写为 $(m+1)\times n$ 的矩阵，即：

$$X=\begin{bmatrix} x_{1}^{(0)} & x_{2}^{(0)} & ... & x_{n}^{(0)} \\ x_{1}^{(1)} & x_{2}^{(1)} & ... & x_{n}^{(1)} \\ \vdots & \vdots & \ddots & \vdots \\ x_{1}^{(m)} & x_{2}^{(m)} & ... & x_{n}^{(m)} \end{bmatrix}$$

同时，将样本中的 $y_i$ 也写为矩阵形式，即输出变量 $Y$ 为 $n\times 1$ 的列向量：

$$Y=[y_1,y_2,...,y_n]^T \in \mathbb{R}^{n}$$

那么损失函数 $J(\boldsymbol{\theta})$ 可写为：

$$\begin{align*} J(\boldsymbol{\theta}) &= \sum_{i=1}^n(\boldsymbol{\theta}^T\mathbf{x_i}-y_i)^2 + \lambda \boldsymbol{\theta}^T\boldsymbol{\theta} \\ &= \boldsymbol{\theta}^TXX^T\boldsymbol{\theta} -\boldsymbol{\theta}^TXY -Y^TX^T\boldsymbol{\theta} +Y^TY + \lambda \boldsymbol{\theta}^T\boldsymbol{\theta} \\ &= \boldsymbol{\theta}^T(XX^T+\lambda I)\boldsymbol{\theta} -\boldsymbol{\theta}^TXY -Y^TX^T\boldsymbol{\theta} +Y^TY \end{align*}$$

要令目标函数最小，显然要令 $\frac{\partial}{\partial\boldsymbol{\theta}}J(\boldsymbol{\theta})=0$

对于 $\frac{\partial}{\partial\boldsymbol{\theta}}J(\boldsymbol{\theta})$，有：

$$\begin{align*} \frac{\partial}{\partial \boldsymbol{\theta}}J(\boldsymbol{\theta}) &= \frac{\partial}{\partial \boldsymbol{\theta}} \Big ( \boldsymbol{\theta}^T(XX^T+\lambda I)\boldsymbol{\theta} -\boldsymbol{\theta}^TXY -Y^TX^T\boldsymbol{\theta} +Y^TY \Big) \\ &= \frac{\partial}{\partial \boldsymbol{\theta}} \boldsymbol{\theta}^T(XX^T+\lambda I)\boldsymbol{\theta} - \frac{\partial}{\partial \boldsymbol{\theta}} \boldsymbol{\theta}^TXY - \frac{\partial}{\partial \boldsymbol{\theta}}Y^TX^T\boldsymbol{\theta} + \frac{\partial}{\partial \boldsymbol{\theta}}Y^TY \end{align*}$$

对于上式的第一部分，由于 $\boldsymbol{\theta}$ 为 $(m+1)\times 1$ 的列向量，$XX^T+\lambda I$ 为 $(m+1)\times (m+1)$ 的矩阵，故 $\theta^T(XX^T+\lambda I)\boldsymbol{\theta}$ 为标量，求导为分母布局下的标量/向量的形式，故有：

$$\begin{align} \frac{\partial}{\partial \boldsymbol{\theta}} \theta^T(XX^T+\lambda I)\boldsymbol{\theta} &= \big[(XX^T+\lambda I)+(XX^T+\lambda I)^T\big] \boldsymbol{\theta} \notag \\ &= \big[ (XX^T+\lambda I)+(XX^T+\lambda I) \big] \boldsymbol{\theta} \notag \\ &= 2(XX^T+\lambda I) \boldsymbol{\theta} \notag \end{align}$$

对于剩余的三部分（证明详见：最小二乘法），有：

$$- \frac{\partial}{\partial \boldsymbol{\theta}} \boldsymbol{\theta}^TXY - \frac{\partial}{\partial \boldsymbol{\theta}}Y^TX^T\boldsymbol{\theta} + \frac{\partial}{\partial \boldsymbol{\theta}}Y^TY = -2XY$$

因此，对于 $\frac{\partial}{\partial\boldsymbol{\theta}}J(\boldsymbol{\theta})$，有：

$$\frac{\partial}{\partial\boldsymbol{\theta}}J(\boldsymbol{\theta})=2(XX^T+\lambda I) \boldsymbol{\theta} - 2XY$$

令 $\frac{\partial}{\partial\boldsymbol{\theta}}J(\boldsymbol{\theta})=0$，则有：

$$\begin{gather} && 2(XX^T+\lambda I) \boldsymbol{\theta} -2XY = 0 \notag \\ &\Rightarrow& (XX^T+\lambda I) \boldsymbol{\theta}=XY \notag \\ &\Rightarrow& (XX^T+\lambda I)^{-1} (XX^T+\lambda I)\boldsymbol{\theta} = (XX^T+\lambda I)^{-1} XY \notag \\ &\Rightarrow& \boldsymbol{\theta}=(XX^T+\lambda I)^{-1} XY \notag \end{gather}$$

【sklearn 实现】

Lasso 回归

以 sklearn 中的波士顿房价数据集为例，实现 Lasso 回归

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import Lasso
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import r2_score

# 特征提取
def deal_data():
    boston = load_boston()  # sklearn的波士顿房价数据集
    df = pd.DataFrame(boston.data, columns=boston.feature_names)
    df['result'] = boston.target
    data = np.array(df)
    return data[:, :-1], data[:, -1]

# 模型训练
def train_model(features, labels):  
    # 建立Lasso回归模型
    model = Lasso()
    
    # 训练
    model.fit(features, labels)
    return model

# 模型评估
def estimate_model(y_true, y_pred):
    MSE = mean_squared_error(y_true, y_pred)
    RMSE = np.sqrt(MSE)
    MAE = mean_absolute_error(y_true, y_pred)
    R2 = r2_score(y_true, y_pred)
    indicators = {"MSE": MSE, "RMSE":RMSE, "MAE":MAE, "R2":R2}
    return indicators

# 可视化
def visualization(y_true, y_pred, model):
    # 绘图
    plt.plot(range(y_true.shape[0]), y_true, "b-") 
    plt.plot(range(y_true.shape[0]), y_pred, "r-.")
    plt.legend(["original value", "predicted value"])
    plt.xlabel("samples", fontsize="15")
    plt.ylabel("y", fontsize="15")
        
    plt.show()

if __name__ == "__main__":
    # 特征提取
    x, y = deal_data()
    
    # 简单交叉验证
    x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3,random_state=0)
    
    # 模型训练
    model = train_model(x_train, y_train)
    
    # 预测结果
    y_pred = model.predict(x_test) # predict()输入输出均为二维
    print("y test:", y_test[:10]) # 测试集y值
    print("y pred:", y_pred[:10]) # 预测y值
    
    # 模型评估
    indicators = estimate_model(y_test, y_pred)
    print("MSE:", indicators["MSE"])
    print("RMSE:", indicators["RMSE"])
    print("MAE:", indicators["MAE"])
    print("R2:", indicators["R2"])
    
    # 可视化
    visualization(y_test, y_pred, model)

Ridge 回归

以 sklearn 中的波士顿房价数据集为例，实现 Ridge 回归

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import r2_score

# 特征提取
def deal_data():
    boston = load_boston()  # sklearn的波士顿房价数据集
    df = pd.DataFrame(boston.data, columns=boston.feature_names)
    df['result'] = boston.target
    data = np.array(df)
    return data[:, :-1], data[:, -1]

# 模型训练
def train_model(features, labels):  
    # 建立Ridge回归模型
    model = Ridge()
    
    # 训练
    model.fit(features, labels)
    return model

# 模型评估
def estimate_model(y_true, y_pred):
    MSE = mean_squared_error(y_true, y_pred)
    RMSE = np.sqrt(MSE)
    MAE = mean_absolute_error(y_true, y_pred)
    R2 = r2_score(y_true, y_pred)
    indicators = {"MSE": MSE, "RMSE":RMSE, "MAE":MAE, "R2":R2}
    return indicators

# 可视化
def visualization(y_true, y_pred, model):
    # 绘图
    plt.plot(range(y_true.shape[0]), y_true, "b-") 
    plt.plot(range(y_true.shape[0]), y_pred, "r-.")
    plt.legend(["original value", "predicted value"])
    plt.xlabel("samples", fontsize="15")
    plt.ylabel("y", fontsize="15")
    
    plt.show()

if __name__ == "__main__":
    # 特征提取
    x, y = deal_data()
    
    # 简单交叉验证
    x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3,random_state=0)
    
    # 模型训练
    model = train_model(x_train, y_train)
    
    # 预测结果
    y_pred = model.predict(x_test) # predict()输入输出均为二维
    print("y test:", y_test[:10]) # 测试集y值
    print("y pred:", y_pred[:10]) # 预测y值
    
    # 模型评估
    indicators = estimate_model(y_test, y_pred)
    print("MSE:", indicators["MSE"])
    print("RMSE:", indicators["RMSE"])
    print("MAE:", indicators["MAE"])
    print("R2:", indicators["R2"])
    
    # 可视化
    visualization(y_test, y_pred, model)