稀疏信号的基于 L1 的模型#

本示例比较了三种基于 L1 的回归模型在从稀疏和相关特征获得的合成信号上的性能，这些特征进一步被加性高斯噪声破坏

一个 Lasso；

一个自动相关性确定 - ARD；

一个弹性网络。

众所周知，当数据维度增加时，Lasso 估计值接近模型选择估计值，前提是不相关变量与相关变量的相关性不太强。在存在相关特征的情况下，Lasso 本身无法选择正确的稀疏模式 [1]。

在这里，我们比较了三个模型在 $R^2$ 分数、拟合时间和估计系数的稀疏性（与真实值相比）方面的性能。

# Author: Arturo Amor <[email protected]>

生成合成数据集#

我们生成一个样本数量少于特征总数的数据集。这会导致欠定系统，即解不是唯一的，因此我们不能单独应用普通最小二乘法。正则化在目标函数中引入了一个惩罚项，它修改了优化问题，并且可以帮助缓解系统的欠定性质。

目标 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  # sparsify coef
y = np.dot(X, true_coef)

一些信息特征具有接近的频率以引起（反）相关性。

freqs[:n_informative]

array([ 2.9502547 , 11.8059798 , 12.63394388, 12.70359377, 24.62241605,
       37.84077985, 40.30506066, 44.63327171, 54.74495357, 59.02456369])

使用 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)

这种稀疏、嘈杂和相关的特征可以从监测某些环境变量的传感器节点获得，因为它们通常根据其位置记录相似的值（空间相关性）。我们可以可视化目标。

import matplotlib.pyplot as plt

plt.plot(time_step, y)
plt.ylabel("target signal")
plt.xlabel("time")
_ = plt.title("Superposition of sinusoidal signals")

为了简单起见，我们将数据分成训练集和测试集。在实践中，应该使用 TimeSeriesSplit 交叉验证来估计测试分数的方差。在这里，我们设置 shuffle="False"，因为在处理具有时间关系的数据时，我们不能使用测试数据之后的训练数据。

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)

在下文中，我们根据拟合优度 $R^2$ 分数和拟合时间来计算三个基于 L1 的模型的性能。然后，我们绘制一个图表，比较估计系数相对于真实系数的稀疏性，最后分析前面的结果。

Lasso#

在本例中，我们演示了一个具有固定正则化参数 alpha 值的 Lasso。在实践中，最佳参数 alpha 应该通过将 TimeSeriesSplit 交叉验证策略传递给 LassoCV 来选择。为了保持示例简单快速执行，我们在此处直接设置 alpha 的最佳值。

from time import time

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

t0 = time()
lasso = Lasso(alpha=0.14).fit(X_train, y_train)
print(f"Lasso fit done in {(time() - t0):.3f}s")

y_pred_lasso = lasso.predict(X_test)
r2_score_lasso = r2_score(y_test, y_pred_lasso)
print(f"Lasso r^2 on test data : {r2_score_lasso:.3f}")

Lasso fit done in 0.002s
Lasso r^2 on test data : 0.480

自动相关性确定 (ARD)#

ARD 回归是 Lasso 的贝叶斯版本。如果需要，它可以为所有参数（包括误差方差）生成区间估计。当信号具有高斯噪声时，它是一个合适的选择。有关 ARDRegression 和 BayesianRidge 回归器的比较，请参见示例比较线性贝叶斯回归器。

from sklearn.linear_model import ARDRegression

t0 = time()
ard = ARDRegression().fit(X_train, y_train)
print(f"ARD fit done in {(time() - t0):.3f}s")

y_pred_ard = ard.predict(X_test)
r2_score_ard = r2_score(y_test, y_pred_ard)
print(f"ARD r^2 on test data : {r2_score_ard:.3f}")

ARD fit done in 0.015s
ARD r^2 on test data : 0.543

弹性网络#

ElasticNet 是 Lasso 和 Ridge 之间的折衷方案，因为它结合了 L1 和 L2 正则化。正则化的数量由两个超参数 l1_ratio 和 alpha 控制。对于 l1_ratio = 0，惩罚是纯 L2，模型等效于 Ridge 。类似地，l1_ratio = 1 是纯 L1 惩罚，模型等效于 Lasso 。对于 0 < l1_ratio < 1，惩罚是 L1 和 L2 的组合。

如前所述，我们使用 alpha 和 l1_ratio 的固定值训练模型。为了选择它们的最佳值，我们使用了 ElasticNetCV ，为了保持示例简单，此处未显示。

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 fit done in {(time() - t0):.3f}s")

y_pred_enet = enet.predict(X_test)
r2_score_enet = r2_score(y_test, y_pred_enet)
print(f"ElasticNet r^2 on test data : {r2_score_enet:.3f}")

ElasticNet fit done in 0.001s
ElasticNet r^2 on test data : 0.636

结果的绘图和分析#

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

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

df = pd.DataFrame(
    {
        "True coefficients": true_coef,
        "Lasso": lasso.coef_,
        "ARDRegression": 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": "coefficients' values"},
    cmap="seismic_r",
)
plt.ylabel("linear model")
plt.xlabel("coefficients")
plt.title(
    f"Models' coefficients\nLasso $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()

Models' coefficients Lasso $R^2$: 0.480, ARD $R^2$: 0.543, ElasticNet $R^2$: 0.636

在本例中，ElasticNet 产生了最佳分数并捕获了大部分预测特征，但仍然无法找到所有真实成分。请注意，ElasticNet 和 ARDRegression 产生的模型都比 Lasso 的稀疏性低。

结论#

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

ElasticNet 对系数引入了一些稀疏性，并将它们的值缩小到零。因此，在存在对目标有贡献的相关特征的情况下，模型仍然能够在不将它们的权重设置为零的情况下降低它们的权重。这导致模型的稀疏性低于纯 Lasso ，并且也可能捕获非预测特征。

ARDRegression 在处理高斯噪声时表现更好，但仍然无法处理相关的特征，并且由于需要拟合先验，因此需要更多的时间。

参考文献#

[1]
“高维数据的稀疏表示的 Lasso 类型恢复” N. Meinshausen, B. Yu - The Annals of Statistics 2009, Vol. 37, No. 1, 246-270

脚本总运行时间：（0 分 0.506 秒）