含污染数据最小二乘插值的双重下降现象：一项模拟研究

# 受污染数据上最小二乘插值的双重下降：一项模拟研究 来源：https://arxiv.org/html/2605.21494 \(2026年4月15日\)

###### 摘要

过参数化模型可以展现出卓越的泛化性能，尽管根据经典统计理论，它们应该容易出现过拟合。 "双重下降" 现象的发现（表明泛化误差在达到某个模型复杂度后会下降）开启了一条新的研究路线。 稳健统计研究的是受污染数据上的统计估计问题。由于真实数据不满足假设，导致数据点相对于假设的"理想"分布表现为异常值，并可能严重扭曲任何经典估计量。 我们探讨的问题是：在训练数据受污染的线性回归设定中，是否能够观察到双重下降现象。 我们将高度非稳健的最小二乘插值估计器与几种稳健替代方法的性能进行了比较。 结果表明，较大的过参数化确实允许出现双重下降现象，使得最小二乘插值器的泛化性能非常好，超越了那些稳健替代方法。

**关键词：** 稳健统计，双重下降，插值区间，受污染数据

## 1 引言

经典统计理论根据所考虑模型类的复杂度来量化机器学习算法的泛化能力，复杂度通过 Rademacher 复杂度、Vapnik-Chervonenkis 维数 (vapnik71) 或伪维数 (pol90) 来衡量（例如，kolt01; kolt02; kolt06; bart02b; bart05b）。 这些结果促使通过偏差-方差权衡来寻找模型类的最优复杂度，即既要避免因过于稀疏的模型导致高偏差，也要避免因过于复杂的模型导致高方差。 由于神经网络类的复杂度随其参数数量增长 (bart03,maass)，大规模神经网络的泛化界限非常宽松，因此无法解释它们通常出色的泛化能力。 在最近一条研究路线中，从实证工作 belkin19 和理论工作 belkin 开始，发现了"双重下降"现象，表明模型类的泛化损失随着复杂度的增加而上升，但随着模型复杂度的进一步增加而下降。 泛化损失的增加被称为过拟合，而下降行为被称为"插值区间"。 虽然该主题的早期工作考虑了最小二乘插值 (belkin;neal;muthu20)，但随后出现了关于过参数化分类 (muthu21;dar21) 和神经网络分类器 (frei;frei23;george;zhu23) 的结果。

稳健统计 (huber;hampel;rieder;maronna) 考虑的是受污染数据的分析问题。 这种污染源于模型误设，使得来自未知真实数据分布的实现相对于假设的"理想"分布表现为异常值。 由于经典的非稳健估计器在应用于受污染数据时容易失真，导致泛化能力差，稳健统计提供了在受污染数据上进行估计同时保持良好泛化能力的方法，代价是效率降低。

在本工作中，我们纯粹从实证角度研究在受污染数据上训练并在干净数据上评估的过参数化模型的泛化性能。 我们考虑最小二乘插值器和不同的稳健对应方法，即使用 Huber 损失、Tukey 损失的回归、稀疏最小截断平方 (SLTS; alfons13) 和稳健 Boosting (RRBoost; ju21)。 我们还提供了一种思路，将最小 \(l_2\) 范数插值与干净子集选择相结合，方法是先分别应用 SLTS 或 RRBoost，然后仅对识别出的干净子集训练一个最小 \(l_2\) 范数插值器。

本文的组织结构如下。 第 2 节 (https://arxiv.org/html/2605.21494#S2) 列出了关于双重下降现象的相关工作。 第 3 节 (https://arxiv.org/html/2605.21494#S3) 汇集了稳健统计、过参数化回归以及我们方法的基本概念，我们首先使用稳健算法识别干净子集，然后对该子集应用最小 \(l_2\) 范数插值器。 第 4 节 (https://arxiv.org/html/2605.21494#S4) 描述了我们的模拟设置、所应用的算法、R 包和评估方法。 关于测试误差、训练误差、系数和迭代次数的结果分别以图形方式呈现在第 5 节 (https://arxiv.org/html/2605.21494#S5)、第 6 节 (https://arxiv.org/html/2605.21494#S6)、第 7 节 (https://arxiv.org/html/2605.21494#S7) 和第 8 节 (https://arxiv.org/html/2605.21494#S8) 中。 我们在第 9 节 (https://arxiv.org/html/2605.21494#S9) 中讨论结果并得出结论。

## 2 相关工作

发现双重下降现象的实证工作 belkin19 已成为一些关于最小二乘插值的理论工作的起点。 在 belkin 中，计算了最小二乘估计量的预测风险，确认了在信噪比 (SNR) 足够高的情况下存在理论上的双重下降。 在神经网络训练的 epoch 数量方面也可能出现双重下降 (nakkiran)。 在一些工作中（例如，belkin;neal;muthu20），计算了高斯设定下最小 \(l_2\) 范数插值器的测试 MSE，表明当 \(p > n\) 时，MSE 随着 \(p\) 的增长而下降。 此外，测试 MSE 可以分解为一个依赖于噪声的分量（取决于噪声方差，因此在无噪声数据上消失）和一个依赖于信号的分量（取决于系数向量的范数）。 在 bartlett20; hastie22; muthu; dar21; tsigler 中计算了依赖于噪声和依赖于信号的误差界。 在 dar21; bartlett21; tsigler 中已表明，双重下降需要低维信号，即必须与 \(s < n\) 的特征向量对齐，然后才能观察到 \(p > n\) 时的插值区间。 vilucchi 研究了 Huber 损失的回归，并在展望中提及可以研究污染与良性过拟合之间的相互作用。 kausik 考虑了有噪声的输入，即只能访问 \(X+A\)，其中噪声矩阵 \(A\) 来自旋转双不变分布，其特征值仍然允许 Marchenko-Pastur 律。 他们表明，最小 \(l_2\) 范数插值器也会出现双重下降。 tripuraneni 考虑了以幂律位移形式出现的协变量偏移，该律可以抑制最大的特征方向。 他们证实了过参数化线性模型对此类偏移具有稳健性。

## 3 预备知识

### 3.1 稳健回归

令 \(n\) 和 \(p\) 分别表示实例数和预测变量数。 在线性回归设定中，令 \(D = (X, Y) \in \mathbb{R}^{n \times (p+1)}\) 为一个数据集，其中 \(Y \in \mathbb{R}^n\) 是响应向量，\(X \in \mathbb{R}^{n \times p}\) 是预测变量矩阵。 我们假设线性结构
\[
Y_i = X_i \beta + \epsilon_i \quad (3.1)
\]
其中 \(\beta \in \mathbb{R}^p\) 是系数向量，\(\epsilon_i\) 是噪声项。 在最小二乘设定中，假设 \(\epsilon_i \sim \mathcal{N}(0, \sigma^2)\) 且独立同分布。 最小二乘估计量在经典设定中，即当 \(p < n\) 时，由下式给出：
\[
\hat{\beta}^{LS} = \arg\min_{\beta \in \mathbb{R}^p} \frac{1}{n} \sum_{i=1}^n (Y_i - X_i \beta)^2 = (X^T X)^{-1} X^T Y \quad (3.2)
\]
稳健 M-估计量 (huber) 通过将平方损失替换为增长较慢的损失函数 \(L(\cdot)\) 来定义：
\[
\hat{\beta}^{M} = \arg\min_{\beta \in \mathbb{R}^p} \frac{1}{n} \sum_{i=1}^n L(Y_i - X_i \beta) \quad (3.3)
\]
这种损失函数的选择影响了 M-估计量的稳健性，即其崩溃点 (BDP)。 一个凸且单调的损失函数是 Huber 损失，\(\delta\) 是调节参数：
\[
L_{\delta}^{\text{Huber}}(r) = \begin{cases} r^2/2, & |r| \leq \delta \\ \delta |r| - \delta^2/2, & |r| > \delta \end{cases}
\]
其导数是有界的但仍然是单调的。 其 BDP 最多为 25%，但根据预测变量的分布可以任意接近 0（参见 hampel, Sec. 6.4）。 一个具有甚至重新下降导数（即 \(\lim_{|r| \rightarrow \infty} (\psi(r)) = 0\)）的损失函数是 Tukey 双权损失函数：
\[
L_k^{\text{Tukey}}(r) = \begin{cases} 1 - [1 - (r/k)^2]^3, & |r| \leq k \\ 1, & |r| > k \end{cases}
\]
另一种策略是最小二乘截断 (LTS) 估计量，由 rous84 提出。 不同于优化所有 \(n\) 个实例的残差平方和，rous84 建议用截断和替换完全和，即
\[
\hat{\beta}^{\text{LTS}} = \operatorname*{argmin}_{\beta \in \mathbb{R}^p} \left( \frac{1}{h} \sum_{i=1}^h r(\beta)_{i:n}^2 \right)
\]
其中 \(r(\beta) = Y - X\beta\)，\(z_{i:n}\) 表示向量 \(z\) 中第 \(i\) 个最小的元素。 换句话说，优化是在一个基数为 \(h\) 的"干净"子集上的平方残差。 rous84 表明，如果选择干净子集的大小 \(h = \lfloor n/2 \rfloor + 1\)，则其 BDP 为 \(\lfloor n/2 \rfloor + \lfloor (p+1)/2 \rfloor\)。

### 3.2 \(p > n\) 的情况

对于过参数化设定，即 \(p > n\)，由于 \(X^T X\) 的奇异性，最小二乘解 (3.2) 中的逆被替换为 Moore-Penrose 逆，从而得到闭式解 \(X^+ Y\)。 至于优化 Huber 或 Tukey 损失的稳健对应方法，由于解缺乏闭式表示，通常应用迭代重加权最小二乘算法 (IRWLS, huber)。 加权最小二乘估计量由下式给出：
\[
\hat{\beta}^{\{WLS\}} = \operatorname*{argmin}_{\beta \in \mathbb{R}^p} \left( \frac{1}{n} \sum_{i=1}^n w_i (Y_i - X_i \beta)^2 \right) = (X^T W X)^{-1} X^T W Y \quad (3.4)
\]
其中 \(W\) 是对角矩阵，其第 \(i\) 个对角元素为 \(w_i\)。 权重 \(w_i\) 满足 \(w_i \geq 0\) 且 \(\sum_i w_i = 1\)。 对于 \(p > n\)，加权最小二乘解被替换为 \((\tilde{X})^+ Y\)，其中 \(\tilde{X} = W^{1/2} X\)。 因此，用最小二乘估计量初始化 IRWLS 方案会将任何 IRWLS 过程简化为最小 \(l_2\) 范数插值，因为最小二乘估计量已经插值了数据。 因此，在我们的模拟中，无法使用 R 包 `robustreg` 中提供的分别最小化 Huber 损失和 Tukey 损失的 IRWLS 估计量的实现。 因此，我们使用 R 包 `MTE` (qin17,MTE)，其中通过梯度下降优化给定的损失函数。 我们将最大迭代次数增加到了 100。 LTS 估计量的计算通过迭代更新干净子集并在更新后的子集上计算最小二乘估计量来完成。 对于 \(p > h < n\)，再次面临干净子集上的拟合是完美的问题。 在干净子集上优化其他损失函数（如 Huber 或 Tukey 损失）是没有意义的，因为这将对应于双重稳健化，既用有界导数损失替换平方损失又进行截断，这至少会显著降低估计量的效率并增加计算成本。 此外，我们在实验中观察到，至少当 \(p\) 足够大时，Huber M-估计量会插值训练数据，这使得任何干净子集的识别都变得毫无意义。 我们提出应用一种允许处理 \(p > n\) 数据但仍优化最小二乘损失的稳健回归技术，并决定考虑两个候选方法：稀疏 LTS (alfons13)，在 R 包 `robustHD` (alfons16) 中提供；以及 ju21 的稳健 Boosting 算法 RRBoost，在 R 包 `RRBoost` (rrboost) 中提供。 为了在干净子集上执行最小 \(l_2\) 范数插值，我们在整个训练集上应用稀疏 LTS 或稳健 Boosting 算法，并基于平方残差识别干净子集。 在该干净子集上，我们计算最小 \(l_2\) 范数插值器。

## 4 模拟

### 4.1 数据生成

本文中，我们仅考虑线性回归。 受 kobak 和 muthu20 的实验启发，我们使用稀疏的真实系数向量，并根据高斯分布生成预测变量，要么是独立特征，要么是尖峰协方差矩阵。 预测变量 \(X_i\) 是来自多元正态分布 \(N_p(\mu \cdot 1_p, \Sigma)\) 的独立同分布实现，其中 \(1_p\) 表示长度为 \(p\) 的全 1 向量，\(\mu \in \mathbb{R}\)，\(\Sigma\) 是某个协方差矩阵。 在独立情况下，我们设 \(\Sigma = I_{p \times p}\)，表示 \(p \times p\) 维的单位矩阵。 在尖峰协方差情况下，我们设 \(\Sigma = I_{p \times p} + \rho 1_p 1_p^T\)。 我们始终使用 \(\rho = 0.25\)。 真实系数向量 \(\beta\) 是随机生成的，要么是高斯的，即 \(\beta \sim N_p(0_p, I_{p \times p})\)，要么是均匀的，即 \(\beta_j \sim U([1,2])\) 且独立同分布。 在我们的模拟中，为了生成关于 \(p\) 的测试损失曲线，我们使用不同的 \(p\)，但始终让真实预测变量数为 \(s = 20\)。 对于 \(p > s\)，\(\beta\) 的 \(p-s\) 个分量被随机选择并设为零。 对于 \(p < s\)，为了保持可识别性，我们使用 \(\beta_j \sim N(0,1)\) 生成所有系数，因为这仍然适用于非常稀疏的情况。 对于 \(s=20\) 且 \(p < 20\)，我们使用所有 \(p\) 个系数为非零，但服从 \(N(0,1)\) 分布，这仍产生稀疏性吗？很好。 然后我们生成 \(\epsilon_i \sim N(0, \sigma^2)\)，其中 \(\sigma^2\) 使得 SNR（定义为 \(\text{Var}(\mathbb{E}[Y|X]) / \mathbb{E}[\text{Var}(Y|X)] = \beta^T \Sigma \beta / \sigma^2\)）等于某个值。 我们使用三个不同的 SNR 值：0.5, 1, 2。 我们生成了 \(n_0 = 5000\) 个干净的测试点。 对于训练数据，我们生成了 \(n\) 个点，然后污染这些点。 我们考虑三种数据污染情况：无污染、\(Y\) 污染和 \(X\) 污染。 对于 \(Y\) 污染，我们随机选择 \(n_{\text{cont}} = \max( r \cdot n, c_{\text{out}} )\) 个实例，其中 \(r \in \{0.1, 0.25, 0.5, 0.75\}\) 是污染比例，\(c_{\text{out}} \in \{100, 10000\}\) 是最小污染绝对数，并将相应的响应值设置为随机值 \(Y_i \sim N(\mu_Y, \sigma_Y^2)\)，其中 \(\mu_Y = 100\)，\(\sigma_Y = 50\)。 对于 \(X\) 污染，我们随机选择 \(n_{\text{cont}}\) 个实例，并将相应的预测变量设置为随机值 \(X_i \sim N(\mu_X \cdot 1_p, I_{p \times p})\)，其中 \(\mu_X = 100\)，\(\sigma_X = 50\)。

### 4.2 评估与表示

我们使用三种不同的稳健回归方法：最小化 Huber 损失的估计量（简称为 Huber）、最小化 Tukey 损失的估计量（简称为 Tukey）、以及稀疏 LTS 估计量（简称为 SLTS）和稳健 Boosting 算法 RRBoost（简称为 RRB）。 通过调用 R 包 `MTE` 中的 `MTE` 函数来训练 Huber 和 Tukey 估计量，而 SLTS 使用 `robustHD` 包中的 `sparseLTS` 函数，RRBoost 使用 `RRBoost` 包中的 `RRBoost` 函数。 作为性能指标，我们计算测试 MSE，即真实值 \(Y_i\) 与预测值 \(\hat{Y}_i\) 之间的均方误差：\(\frac{1}{n_0} \sum_{i=1}^{n_0} (Y_i - \hat{Y}_i)^2\)。 我们也在训练数据上计算训练 MSE。 此外，我们计算估计系数与真实系数之间的 \(l_2\) 距离，以及——对于需要迭代次数的方法——收敛所需的迭代次数。

### 4.3 模拟参数

在所有模拟中，我们使用 \(n = 100\) 或 \(n = 200\) 个训练实例。 对于预测变量数 \(p\)，我们使用以下网格：\(p \in \{10, 20, 50, 100, 200, 500, 1000, 2000, 5000\}\)。 我们为 \((n, p)\) 的每种组合生成 \(N_{\text{sim}} = 100\) 次独立重复。 对于每种重复，我们生成新的系数向量 \(\beta\)、新的训练数据（可能受污染）和新的测试数据。 训练用 \(n\) 个实例生成，测试用 \(n_0 = 5000\) 个实例生成。 我们报告所有重复的平均 MSE 以及系数距离，并针对污染比例 \(r\)、协方差结构（独立或尖峰）和 SNR（0.5, 1, 2）的各种组合进行。 为了简洁，我们仅呈现 \(n = 100, \mu = 0\) 的独立特征情况的结果以及 \(n = 200, c_{\text{out}} = 10000\) 的 \(Y\) 污染结果。

## 5 测试误差

### 5.1 独立特征，\(\mu = 0\)

#### 5.1.1 最小 \(l_2\) 范数插值

图 1：在干净训练数据上训练的最小 \(l_2\) 范数插值的测试 MSE。 当在干净数据上训练时（图 1），MSE 在 \(p = s\) 处达到最小值，然后随着 \(p\) 增长而下降，最终稳定在某个水平。 在 \(Y\) 污染下（图 2 顶部一行），对于 \(c_{\text{out}} = 100\)，MSE 行为类似于干净情况，但数值更高。 对于 \(c_{\text{out}} = 10000\)，曲线在 \(p = 20\) 后急剧上升，然后随着 \(p\) 增长而下降。 在 \(X\) 污染下（图 2 底部一行），MSE 在所有情况下都单调下降。

#### 5.1.2 Huber 损失插值

图 3：在干净训练数据上训练的 Huber 损失插值的测试 MSE。 在干净数据上（图 3），MSE 在 \(p\) 增长时首先上升，达到峰值，然后下降。 对于低污染（\(r = 0.1, 0.25\)，图 4 顶部两行），Huber 插值的 MSE 曲线遵循类似的模式，但峰值更高且出现在更大的 \(p\) 值。 对于高污染（\(r = 0.5, 0.75\)），MSE 在 \(p = 200\) 前单调上升，然后似乎下降但波动很大。

### 5.2 尖峰协方差，\(\mu = 0\)

#### 5.2.1 最小 \(l_2\) 范数插值

图 5：在干净训练数据上训练的最小 \(l_2\) 范数插值的测试 MSE。 在干净数据上（图 5），MSE 在 \(p = s\) 后下降，但下降速度比独立情况慢。 在污染下（图 6），行为与独立情况定性相似，但 MSE 值更高，且峰值更不明显。

### 5.3 \(n = 100, c_{\text{out}} = 100\)

#### 5.3.1 最小 \(l_2\) 范数插值

图 7：在 \(Y\) 污染训练数据上训练的最小 \(l_2\) 范数插值的测试 MSE。 对于 \(c_{\text{out}} = 100\)，MSE 曲线（图 7）在低污染（\(r = 0.1\)）时类似于干净情况，但在 \(r \ge 0.25\) 时，MSE 在 \(p = 20\) 后上升，然后在 \(p = 1000\) 后下降。

### 5.4 \(n = 100, c_{\text{out}} = 10000\)

#### 5.4.1 最小 \(l_2\) 范数插值

图 38：在 \(Y\) 污染训练数据上训练的最小 \(l_2\) 范数插值的测试 MSE。 对于 \(c_{\text{out}} = 10000\)（图 38），MSE 曲线在 \(p = 20\) 后急剧上升，达到远高于干净情况的峰值，然后在 \(p = 5000\) 时下降但仍高于干净情况。

### 5.5 \(n = 200, c_{\text{out}} = 10000\)

#### 5.5.1 最小 \(l_2\) 范数插值

图 39：在 \(Y\) 污染训练数据上训练的最小 \(l_2\) 范数插值的测试 MSE。 图 39 中的 MSE 曲线与 \(c_{\text{out}} = 100\) 时的图 37 相似，尽管 MSE 值明显更大。

#### 5.5.2 Huber 损失插值

图 40：在 \(Y\) 污染训练数据上训练的 Huber 损失插值的测试 MSE。 对于 \(r \in \{0.1, 0.25\}\)，图 40 中的 MSE 曲线在 \(p = s\) 处达到最小值，之后单调增加。 与 \(c_{\text{out}} = 100\) 时图 38 所示的情况相反，即使考虑的最大的 \(p\) 值，曲线也没有下降。 对于 \(r \ge 0.5\)，MSE 单调增加，并在 \(p = 2000\) 处达到峰值，然后似乎随着 \(p\) 的进一步增长而下降。

## 6 训练误差

### 6.1 独立特征，\(\mu = 0\)

#### 6.1.1 最小 \(l_2\) 范数插值

图 41：在干净训练数据上训练的最小 \(l_2\) 范数插值的训练 MSE。 通过插值，一旦 \(p > n\)，训练误差消失，如图 41 和图 42 所示。 正如预期，在 \(Y\) 污染和 \(X\) 污染数据上，低 \(p\) 时的训练误差高于干净数据。

图 42：在受污染训练数据上训练的最小 \(l_2\) 范数插值的训练 MSE。

#### 6.1.2 Huber 损失插值

图 43：Huber 损失