This paper studies M-estimators with gradient-Lipschitz loss function regularized with convex penalty in linear models with Gaussian design matrix and arbitrary noise distribution. A practical example is the robust M-estimator constructed with the Huber loss and the Elastic-Net penalty and the noise distribution has heavy-tails. Our main contributions are three-fold. (i) We provide general formulae for the derivatives of regularized M-estimators $\hatβ(y,X)$ where differentiation is taken with respect to both $y$ and $X$; this reveals a simple differentiability structure shared by all convex regularized M-estimators. (ii) Using these derivatives, we characterize the distribution of the residual $r_i = y_i-x_i^\top\hatβ$ in the intermediate high-dimensional regime where dimension and sample size are of the same order. (iii) Motivated by the distribution of the residuals, we propose a novel adaptive criterion to select tuning parameters of regularized M-estimators. The criterion approximates the out-of-sample error up to an additive constant independent of the estimator, so that minimizing the criterion provides a proxy for minimizing the out-of-sample error. The proposed adaptive criterion does not require the knowledge of the noise distribution or of the covariance of the design. Simulated data confirms the theoretical findings, regarding both the distribution of the residuals and the success of the criterion as a proxy of the out-of-sample error. Finally our results reveal new relationships between the derivatives of $\hatβ(y,X)$ and the effective degrees of freedom of the M-estimator, which are of independent interest.
Derivatives and residual distribution of regularized M-estimators with application to adaptive tuning 论文ID : 2107.05143标题 : Derivatives and residual distribution of regularized M-estimators with application to adaptive tuning作者 : Pierre C. Bellec (Rutgers University), Yiwei Shen (Rutgers University)分类 : math.ST stat.ML stat.TH发表会议 : Proceedings of Machine Learning Research vol 178:1–36, 2022论文链接 : https://arxiv.org/abs/2107.05143 本文研究在高斯设计矩阵和任意噪声分布的线性模型中,具有梯度Lipschitz损失函数和凸惩罚项的M-估计器。主要贡献包括:(1) 提供了正则化M-估计器 β ^ ( y , X ) \hat{\beta}(y,X) β ^ ( y , X ) 关于 y y y 和 X X X 的导数的一般公式,揭示了所有凸正则化M-估计器共享的简单可微结构;(2) 利用这些导数,在维数和样本量同阶的中等高维regime下刻画了残差 r i = y i − x i ⊤ β ^ r_i = y_i-x_i^\top\hat{\beta} r i = y i − x i ⊤ β ^ 的分布;(3) 基于残差分布提出了一个新的自适应准则来选择正则化M-估计器的调优参数,该准则能够逼近样本外误差且不需要知道噪声分布或设计协方差。
在高维统计中,M-估计器是处理异常值和重尾噪声的重要工具。典型的M-估计器形式为:
β ^ ( y , X ) = arg min b ∈ R p 1 n ∑ i = 1 n ρ ( y i − x i ⊤ b ) + g ( b ) \hat{\beta}(y,X) = \arg\min_{b\in\mathbb{R}^p} \frac{1}{n}\sum_{i=1}^n \rho(y_i - x_i^\top b) + g(b) β ^ ( y , X ) = arg min b ∈ R p n 1 ∑ i = 1 n ρ ( y i − x i ⊤ b ) + g ( b )
其中 ρ \rho ρ 是凸损失函数(如Huber损失),g g g 是凸惩罚项(如Elastic-Net)。
参数调优的困难性 :现有的调优方法通常需要知道噪声分布或设计协方差矩阵,在实际应用中往往不可获得。理论理解的不足 :对于一般的M-估计器,其可微性结构和残差分布的理论理解还不够深入。实用性需求 :需要一个完全自适应的调优准则,既不依赖于未知参数,又能有效选择最优的损失-惩罚对。大多数现有工作局限于平方损失 需要知道设计协方差矩阵 Σ \Sigma Σ 对非光滑惩罚函数缺乏理论保证 导数公式的统一框架 :为任意凸正则化M-估计器提供了关于 ( y , X ) (y,X) ( y , X ) 的导数的一般公式,揭示了统一的可微结构。残差分布的随机表示 :在中等高维regime下,给出了个体残差的精确随机表示和渐近正态性结果。自适应调优准则 :提出了完全自适应的参数选择准则,不需要知道噪声分布或设计协方差。有效自由度的新关系 :建立了M-估计器导数与有效自由度之间的新联系。考虑线性模型 y = X β ∗ + ε y = X\beta^* + \varepsilon y = X β ∗ + ε ,其中:
X ∈ R n × p X \in \mathbb{R}^{n \times p} X ∈ R n × p 的行向量独立同分布于 N ( 0 , Σ ) N(0,\Sigma) N ( 0 , Σ ) ε \varepsilon ε 独立于 X X X ,具有连续分布维数 p p p 与样本量 n n n 同阶 对于几乎所有 ( y , X ) (y,X) ( y , X ) ,存在矩阵 A ^ ∈ R p × p \hat{A} \in \mathbb{R}^{p \times p} A ^ ∈ R p × p 使得:
∂ ∂ y i β ^ ( y , X ) = A ^ X ⊤ e i ψ ′ ( r i ) \frac{\partial}{\partial y_i}\hat{\beta}(y,X) = \hat{A}X^\top e_i \psi'(r_i) ∂ y i ∂ β ^ ( y , X ) = A ^ X ⊤ e i ψ ′ ( r i )
∂ ∂ x i j β ^ ( y , X ) = A ^ e j ψ ( r i ) − A ^ X ⊤ e i ψ ′ ( r i ) β ^ j \frac{\partial}{\partial x_{ij}}\hat{\beta}(y,X) = \hat{A}e_j\psi(r_i) - \hat{A}X^\top e_i \psi'(r_i)\hat{\beta}_j ∂ x ij ∂ β ^ ( y , X ) = A ^ e j ψ ( r i ) − A ^ X ⊤ e i ψ ′ ( r i ) β ^ j
其中 r i = y i − x i ⊤ β ^ r_i = y_i - x_i^\top\hat{\beta} r i = y i − x i ⊤ β ^ ,ψ = ρ ′ \psi = \rho' ψ = ρ ′ ,∥ Σ 1 / 2 A ^ Σ 1 / 2 ∥ o p ≤ ( n μ ) − 1 \|\Sigma^{1/2}\hat{A}\Sigma^{1/2}\|_{op} \leq (n\mu)^{-1} ∥ Σ 1/2 A ^ Σ 1/2 ∥ o p ≤ ( n μ ) − 1 。
对于每个 i = 1 , … , n i = 1,\ldots,n i = 1 , … , n ,存在 Z i ∼ N ( 0 , 1 ) Z_i \sim N(0,1) Z i ∼ N ( 0 , 1 ) 独立于 ε i \varepsilon_i ε i 使得:
∣ r i + tr [ Σ A ^ ] ψ ( r i ) − ( ε i + ∥ Σ 1 / 2 ( β ^ − β ∗ ) ∥ Z i ) ∣ ≤ O P ( n − 1 / 4 ) ( 误差项 ) \left|r_i + \text{tr}[\Sigma\hat{A}]\psi(r_i) - (\varepsilon_i + \|\Sigma^{1/2}(\hat{\beta}-\beta^*)\|Z_i)\right| \leq O_P(n^{-1/4})(\text{误差项}) r i + tr [ Σ A ^ ] ψ ( r i ) − ( ε i + ∥ Σ 1/2 ( β ^ − β ∗ ) ∥ Z i ) ≤ O P ( n − 1/4 ) ( 误差项 )
这给出了残差的随机表示:
r i + tr [ Σ A ^ ] ψ ( r i ) ≈ ε i + ∥ Σ 1 / 2 ( β ^ − β ∗ ) ∥ Z i r_i + \text{tr}[\Sigma\hat{A}]\psi(r_i) \approx \varepsilon_i + \|\Sigma^{1/2}(\hat{\beta}-\beta^*)\|Z_i r i + tr [ Σ A ^ ] ψ ( r i ) ≈ ε i + ∥ Σ 1/2 ( β ^ − β ∗ ) ∥ Z i
基于残差分布,提出调优准则:
Crit ( ρ , g ) = ∥ r + d f ^ tr [ V ] ψ ( r ) ∥ 2 \text{Crit}(\rho, g) = \left\|r + \frac{\hat{df}}{\text{tr}[V]}\psi(r)\right\|^2 Crit ( ρ , g ) = r + tr [ V ] df ^ ψ ( r ) 2
其中:
r = y − X β ^ ρ , g r = y - X\hat{\beta}_{\rho,g} r = y − X β ^ ρ , g d f ^ = tr [ X ( ∂ / ∂ y ) β ^ ρ , g ] \hat{df} = \text{tr}[X(\partial/\partial y)\hat{\beta}_{\rho,g}] df ^ = tr [ X ( ∂ / ∂ y ) β ^ ρ , g ] V = diag { ψ ′ ( r ) } ( I n − X ( ∂ / ∂ y ) β ^ ρ , g ) V = \text{diag}\{\psi'(r)\}(I_n - X(\partial/\partial y)\hat{\beta}_{\rho,g}) V = diag { ψ ′ ( r )} ( I n − X ( ∂ / ∂ y ) β ^ ρ , g ) 统一的可微结构 :首次为一般凸M-估计器建立了统一的导数公式,包括非光滑惩罚。有效自由度估计 :提出 d f ^ / tr [ V ] \hat{df}/\text{tr}[V] df ^ / tr [ V ] 作为 tr [ Σ A ^ ] \text{tr}[\Sigma\hat{A}] tr [ Σ A ^ ] 的估计,避免了对 Σ \Sigma Σ 的依赖。概率工具的创新使用 :巧妙运用Stein公式和Gaussian积分技巧处理高维M-估计器。样本量 :n = 1001 n = 1001 n = 1001 ,维数 :p = 1000 p = 1000 p = 1000 设计矩阵 :X X X 的行独立同分布于 N ( 0 , Σ ) N(0,\Sigma) N ( 0 , Σ ) ,其中 Σ = R ⊤ R / ( 2 p ) \Sigma = R^\top R/(2p) Σ = R ⊤ R / ( 2 p ) ,R R R 为Rademacher矩阵真实参数 :β ∗ \beta^* β ∗ 的前100个分量为 10 / 10 \sqrt{10}/10 10 /10 ,其余为0噪声 :ε i \varepsilon_i ε i 独立同分布于自由度为2的t分布(重尾)使用Huber-Elastic-Net估计器:
损失函数 :ρ ( u ; Λ ) = Λ 2 H ( Λ − 1 u ) \rho(u;\Lambda) = \Lambda^2 H(\Lambda^{-1}u) ρ ( u ; Λ ) = Λ 2 H ( Λ − 1 u ) ,其中 H H H 为Huber损失惩罚项 :g ( b ; λ , τ ) = λ ∥ b ∥ 1 + ( τ / 2 ) ∥ b ∥ 2 2 g(b;\lambda,\tau) = \lambda\|b\|_1 + (\tau/2)\|b\|_2^2 g ( b ; λ , τ ) = λ ∥ b ∥ 1 + ( τ /2 ) ∥ b ∥ 2 2 样本外误差:∥ Σ 1 / 2 ( β ^ − β ∗ ) ∥ 2 \|\Sigma^{1/2}(\hat{\beta}-\beta^*)\|^2 ∥ Σ 1/2 ( β ^ − β ∗ ) ∥ 2 调优准则的逼近误差 残差正态性检验 图1展示了在 ( λ , τ ) (\lambda,\tau) ( λ , τ ) 网格上:
真实样本外误差 ∥ Σ 1 / 2 ( β ^ − β ∗ ) ∥ 2 \|\Sigma^{1/2}(\hat{\beta}-\beta^*)\|^2 ∥ Σ 1/2 ( β ^ − β ∗ ) ∥ 2 调优准则的逼近 ∥ r + ( d f ^ / tr [ V ] ) ψ ( r ) ∥ 2 / n − ∥ ε ∥ 2 / n \|r + (\hat{df}/\text{tr}[V])\psi(r)\|^2/n - \|\varepsilon\|^2/n ∥ r + ( df ^ / tr [ V ]) ψ ( r ) ∥ 2 / n − ∥ ε ∥ 2 / n 逼近误差 结果显示调优准则能够准确逼近样本外误差的相对大小。
图2展示了标准化残差 ζ 1 \zeta_1 ζ 1 的直方图和QQ图,在不同参数组合下都很好地符合标准正态分布,验证了理论预测。
表1显示 ∣ tr [ Σ A ^ ] − d f ^ / tr [ V ] ∣ |\text{tr}[\Sigma\hat{A}] - \hat{df}/\text{tr}[V]| ∣ tr [ Σ A ^ ] − df ^ / tr [ V ] ∣ 的值很小(约0.002),证实了 d f ^ / tr [ V ] \hat{df}/\text{tr}[V] df ^ / tr [ V ] 是 tr [ Σ A ^ ] \text{tr}[\Sigma\hat{A}] tr [ Σ A ^ ] 的良好估计。
定理7-8 :证明了基于调优准则选择的估计器以高概率达到最优样本外误差定理9 :E [ ∣ tr [ Σ A ^ ] tr [ V ] / n − d f ^ / n ∣ ] ≤ C ( γ , μ ) n − 1 / 2 E[|\text{tr}[\Sigma\hat{A}]\text{tr}[V]/n - \hat{df}/n|] \leq C(γ,μ)n^{-1/2} E [ ∣ tr [ Σ A ^ ] tr [ V ] / n − df ^ / n ∣ ] ≤ C ( γ , μ ) n − 1/2 定理6 :∥ Σ 1 / 2 ( β ^ − β ∗ ) ∥ 2 + ∥ ε ∥ 2 / n = ( 1 + O P ( n − 1 / 2 ) ) ∥ r + tr [ Σ A ^ ] ψ ( r ) ∥ 2 / n \|\Sigma^{1/2}(\hat{\beta}-\beta^*)\|^2 + \|\varepsilon\|^2/n = (1+O_P(n^{-1/2}))\|r + \text{tr}[\Sigma\hat{A}]\psi(r)\|^2/n ∥ Σ 1/2 ( β ^ − β ∗ ) ∥ 2 + ∥ ε ∥ 2 / n = ( 1 + O P ( n − 1/2 )) ∥ r + tr [ Σ A ^ ] ψ ( r ) ∥ 2 / n 本文建立在以下工作基础上:
Bayati & Montanari (2012) :LASSO的风险分析El Karoui et al. (2013) :无惩罚M-估计器的研究Thrampoulidis et al. (2018) :一般损失-惩罚对的精确误差分析与现有方法的比较:
ALO准则 (Rad et al., 2020) :需要二阶连续可微性假设基于Σ的准则 (Bellec, 2020) :需要知道设计协方差本文方法 :完全自适应,适用于非光滑函数本文首次将可观测量(仅依赖于数据)用于描述M-估计器行为,而非依赖不可观测的先验分布或协方差矩阵。
统一理论框架 :为凸正则化M-估计器建立了统一的可微性理论。实用调优工具 :提供了不需要先验知识的自适应参数选择方法。理论保证 :在合理假设下证明了方法的有效性。高斯设计假设 :主要理论结果需要高斯设计矩阵,尽管模拟显示对Rademacher设计也有效。强凸性要求 :部分结果需要惩罚项的强凸性,虽然第7节提供了放松方法。计算复杂性 :对于某些非光滑惩罚,矩阵 A ^ \hat{A} A ^ 没有闭式表达式。扩展到非高斯设计 处理更一般的损失函数类 开发计算高效的实现算法 理论贡献显著 :首次为一般M-估计器提供了统一的导数理论,填补了重要理论空白。实用价值高 :提出的调优准则完全自适应,在实际应用中具有重要价值。技术创新性强 :巧妙结合了凸分析、随机矩阵理论和Stein方法。实验验证充分 :通过多种设置验证了理论预测的准确性。假设限制性 :高斯设计假设限制了方法的普适性。计算考虑不足 :对于实际计算中的数值稳定性和效率讨论较少。比较不够全面 :与其他自适应方法的经验比较有限。理论影响 :为高维M-估计器理论提供了新的分析工具。实践价值 :为鲁棒回归中的参数选择提供了实用方法。方法论贡献 :展示了如何将高维概率论与统计推断相结合。高维鲁棒回归问题 存在异常值或重尾噪声的数据分析 需要自适应参数选择的机器学习应用 金融、生物信息学等对鲁棒性要求高的领域 主要参考文献包括:
Bayati, M. and Montanari, A. (2012). The lasso risk for gaussian matrices. El Karoui, N. et al. (2013). On robust regression with high-dimensional predictors. Thrampoulidis, C. et al. (2018). Precise error analysis of regularized m-estimators in high dimensions. Bellec, P.C. (2020). Out-of-sample error estimate for robust m-estimators with convex penalty.