使用基于 L1 的模型进行稀疏信号回归

简介

在本实验中，我们将演示如何使用基于 L1 的回归模型来处理高维稀疏信号。具体来说，我们将比较三种流行的基于 L1 的模型：套索（Lasso）、自动相关性确定（ARD）和弹性网络（ElasticNet）。我们将使用一个合成数据集来说明这些模型在拟合时间、R2 分数和估计系数的稀疏性方面的性能。

虚拟机使用提示

虚拟机启动完成后，点击左上角切换到“笔记本”标签页，以访问 Jupyter Notebook 进行练习。

有时，你可能需要等待几秒钟让 Jupyter Notebook 完成加载。由于 Jupyter Notebook 的限制，操作验证无法自动化。

如果你在学习过程中遇到问题，随时向 Labby 提问。课程结束后提供反馈，我们会及时为你解决问题。

Skills Graph

生成合成数据集

首先，我们生成一个样本数量少于特征总数的数据集。这会导致一个欠定系统，即解不唯一，我们不能单独应用普通最小二乘法。正则化会在目标函数中引入一个惩罚项，它会修改优化问题，并有助于缓解系统的欠定性质。我们将生成一个目标 y，它是正弦信号的线性组合，其系数正负交替。在 X 的 100 个频率中，只有最低的 10 个频率用于生成 y，而其余特征则没有信息价值。这就导致了一个高维稀疏特征空间，在这种情况下，某种程度的 L1 惩罚是必要的。

import numpy as np

rng = np.random.RandomState(0)
n_samples, n_features, n_informative = 50, 100, 10
time_step = np.linspace(-2, 2, n_samples)
freqs = 2 * np.pi * np.sort(rng.rand(n_features)) / 0.01
X = np.zeros((n_samples, n_features))

for i in range(n_features):
    X[:, i] = np.sin(freqs[i] * time_step)

idx = np.arange(n_features)
true_coef = (-1) ** idx * np.exp(-idx / 10)
true_coef[n_informative:] = 0  ## 使系数稀疏化
y = np.dot(X, true_coef)

## 使用 numpy.random.random_sample 引入随机相位
## 使用 numpy.random.normal 添加一些高斯噪声
for i in range(n_features):
    X[:, i] = np.sin(freqs[i] * time_step + 2 * (rng.random_sample() - 0.5))
    X[:, i] += 0.2 * rng.normal(0, 1, n_samples)

y += 0.2 * rng.normal(0, 1, n_samples)

## 使用 sklearn 中的 train_test_split 将数据拆分为训练集和测试集
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, shuffle=False)

套索（Lasso）

在这一步中，我们将演示如何使用套索回归模型来估计数据集的稀疏系数。我们将使用正则化参数 alpha 的一个固定值。在实际应用中，应该通过将 TimeSeriesSplit 交叉验证策略传递给 LassoCV 来选择最优参数 alpha。为了使示例简单且执行速度快，我们在此直接设置 alpha 的最优值。

from sklearn.linear_model import Lasso
from sklearn.metrics import r2_score
from time import time

t0 = time()
lasso = Lasso(alpha=0.14).fit(X_train, y_train)
print(f"套索（Lasso）拟合完成，耗时{(time() - t0):.3f}秒")

y_pred_lasso = lasso.predict(X_test)
r2_score_lasso = r2_score(y_test, y_pred_lasso)
print(f"套索（Lasso）在测试数据上的 R 平方值 : {r2_score_lasso:.3f}")

自动相关性确定（ARD）

自动相关性确定（ARD）回归是套索（Lasso）的贝叶斯版本。如果需要，它可以为所有参数（包括误差方差）生成区间估计。当信号具有高斯噪声时，它是一个合适的选择。

from sklearn.linear_model import ARDRegression

t0 = time()
ard = ARDRegression().fit(X_train, y_train)
print(f"自动相关性确定（ARD）拟合完成，耗时{(time() - t0):.3f}秒")

y_pred_ard = ard.predict(X_test)
r2_score_ard = r2_score(y_test, y_pred_ard)
print(f"自动相关性确定（ARD）在测试数据上的 R 平方值 : {r2_score_ard:.3f}")

弹性网络（ElasticNet）

弹性网络（ElasticNet）是套索（Lasso）回归和岭回归之间的一种折中方法，因为它结合了 L1 和 L2 惩罚。正则化的程度由两个超参数 l1_ratio 和 alpha 控制。对于 l1_ratio = 0，惩罚是纯粹的 L2 惩罚，模型等同于岭回归。类似地，l1_ratio = 1 是纯粹的 L1 惩罚，模型等同于套索回归。对于 0 < l1_ratio < 1，惩罚是 L1 和 L2 的组合。

from sklearn.linear_model import ElasticNet

t0 = time()
enet = ElasticNet(alpha=0.08, l1_ratio=0.5).fit(X_train, y_train)
print(f"弹性网络（ElasticNet）拟合完成，耗时{(time() - t0):.3f}秒")

y_pred_enet = enet.predict(X_test)
r2_score_enet = r2_score(y_test, y_pred_enet)
print(f"弹性网络（ElasticNet）在测试数据上的 R 平方值 : {r2_score_enet:.3f}")

结果的绘制与分析

在这一步中，我们使用热图来可视化各个线性模型的真实系数和估计系数的稀疏性。

import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
from matplotlib.colors import SymLogNorm

df = pd.DataFrame(
    {
        "真实系数": true_coef,
        "套索（Lasso）": lasso.coef_,
        "自动相关性确定（ARD）回归": ard.coef_,
        "弹性网络（ElasticNet）": enet.coef_,
    }
)

plt.figure(figsize=(10, 6))
ax = sns.heatmap(
    df.T,
    norm=SymLogNorm(linthresh=10e-4, vmin=-1, vmax=1),
    cbar_kws={"label": "系数值"},
    cmap="seismic_r"
)
plt.ylabel("线性模型")
plt.xlabel("系数")
plt.title(
    f"模型的系数\n套索（Lasso） $R^2$: {r2_score_lasso:.3f}, "
    f"自动相关性确定（ARD） $R^2$: {r2_score_ard:.3f}, "
    f"弹性网络（ElasticNet） $R^2$: {r2_score_enet:.3f}"
)
plt.tight_layout()

总结

众所周知，套索（Lasso）能有效地恢复稀疏数据，但在处理高度相关的特征时表现不佳。实际上，如果几个相关特征对目标有贡献，套索（Lasso）最终只会选择其中一个。在特征稀疏但不相关的情况下，套索（Lasso）模型会更合适。

弹性网络（ElasticNet）在系数上引入了一定的稀疏性，并将其值收缩至零。因此，在存在对目标有贡献的相关特征时，该模型仍能够降低它们的权重，而不会将其精确设置为零。这导致模型的稀疏性低于纯套索（Lasso）模型，并且可能还会捕捉到非预测性特征。

自动相关性确定（ARD）回归在处理高斯噪声时表现更好，但仍然无法处理相关特征，并且由于拟合先验而需要更多时间。