Technical Summary: Robust Regularised M-Estimators for High-Dimensional Regression with Heavy-Tailed and Skewed Errors

1. Problem Statement

The paper addresses the challenge of high-dimensional linear regression where the number of predictors ($p$) substantially exceeds the sample size ($n$). While penalised regression methods like the Lasso are standard in this regime, their theoretical guarantees and practical performance rely heavily on the assumption that model errors are sub-Gaussian or possess at least four finite moments.

In many empirical domains—such as financial log-returns, gene expression data, and environmental monitoring—error distributions are often heavy-tailed (e.g., Student's $t$ with low degrees of freedom), skewed, or contaminated. Under these conditions, standard squared-loss estimators (like the Lasso) can suffer from infinite variance, instability, or complete divergence. Existing robust high-dimensional regression literature has largely been constrained by requiring at least three or four finite moments, leaving a gap in theoretical and practical guidance for distributions with only $\bar{\alpha} > 1$ finite absolute moments (e.g., near-Cauchy tails).

2. Methodology

2.1 Model and Assumptions

The authors consider the linear model $y_i = x_i^T \beta^* + \varepsilon_i$, where $\beta^*$ is a sparse vector. The framework relies on four key assumptions:

1. Sparsity: The true support size $s^*$ satisfies $s^* \ll p$ and $s^* \log p = o(n)$.
2. Moment Condition: Errors have finite absolute moments of order $\bar{\alpha} \in (1, 2]$. Symmetry is assumed for Huber loss results, while asymmetric extensions are handled by Catoni and rank-based losses.
3. Design Normalization: Predictors are normalized to unit $\ell_2$ norm.
4. Restricted Eigenvalue: The design matrix satisfies a restricted eigenvalue condition to ensure identifiability.

2.2 Estimators and Loss Functions

The paper proposes a unified framework combining convex M-estimation losses with sparsity-inducing penalties ($\ell_1$ and adaptive $\ell_1$). Three specific loss functions are analyzed:

• Huber Loss: Piecewise quadratic, differentiable, and minimax optimal for contaminated normal models. Requires a threshold parameter $c$.
• Catoni-type Loss: A smooth, bounded influence function that does not require a growing threshold, making it suitable when $\bar{\alpha}$ is unknown and potentially close to 1.
• Convoluted Rank-Based Loss: A smoothed rank-based objective that naturally accommodates asymmetric errors without bias correction.

These losses are paired with the Lasso penalty ($\lambda \|\beta\|_1$) and the Adaptive Lasso penalty ($\lambda \sum \hat{w}_j |\beta_j|$), where weights are derived from an initial Huber-Lasso estimate.

2.3 Parameter Selection (RGCV)

A critical methodological contribution is the Robust Generalised Cross-Validation (RGCV) criterion. Standard cross-validation fails under heavy tails because the squared prediction error on held-out folds has infinite variance when $\bar{\alpha} < 2$. The RGCV replaces the squared loss in the numerator with the Huber loss and uses a stochastic estimator (Hutchinson's method) to approximate the trace of the hat matrix, allowing for the joint selection of the robustification parameter ($c$) and the regularisation parameter ($\lambda$).

2.4 Computational Implementation

The estimators are solved using Iteratively Reweighted Least Squares (IRLS) combined with coordinate descent. The implementation is provided in R and Python.

3. Key Theoretical Contributions

3.1 Oracle Inequalities under Weak Moments

The paper establishes non-asymptotic oracle inequalities for the proposed estimators under the condition $\bar{\alpha} > 1$. This relaxes the standard requirement of $\bar{\alpha} \ge 4$ found in prior work (e.g., Loh, 2015).

• Convergence Rates: The prediction error rate adapts continuously to the tail index $\bar{\alpha}$. For $\bar{\alpha}=2$ (Gaussian), the rate is $O(s^* \log p / n)$. As $\bar{\alpha} \to 1^+$, the rate degrades to $O(s^* (\log p / n)^{2\bar{\alpha}/(\bar{\alpha}+1)})$, reflecting the fundamental difficulty of estimation under near-Cauchy tails.
• Generality: The results hold for Huber, Catoni, and Rank-based losses, with Catoni and Rank-based losses explicitly handling asymmetric errors.

3.2 Negative Result on Variable Selection

A significant theoretical finding is the characterization of the breakdown of variable selection under near-Cauchy errors.

• Sign Consistency: The paper proves that for sign consistency (recovery of the true support) to hold, the condition $s^* \log p = o(n^{\bar{\alpha}-1})$ must be met.
• Implication: For $\bar{\alpha} \in (1, 2)$, this condition becomes practically unachievable for growing $n$ and $s^*$. Specifically, under near-Cauchy errors ($\bar{\alpha} \approx 1.2$), reliable variable selection is theoretically impossible regardless of the loss function used. The role of robust methods in this regime is limited to stable coefficient estimation, not support recovery.

4. Empirical Results

4.1 Simulation Study

Simulations were conducted across $p/n \in [1, 25]$ and five error distributions (Gaussian, Student's $t(3)$, Cauchy, Contaminated Normal, Log-Normal).

• Prediction Error: Under Gaussian errors, robust methods incur a minor efficiency loss (4–14%) compared to standard Lasso. Under heavy-tailed errors ($t(3)$), the Adaptive Huber-Lasso achieved a 50% reduction in Mean Squared Prediction Error (MSPE) compared to standard Lasso. Under Cauchy errors, standard Lasso diverged in 34% of replications, whereas robust methods remained stable.
• Variable Selection: Under $t(3)$ errors, the Adaptive Huber-Lasso achieved high F1 scores (0.90), approaching the Gaussian benchmark. Under Cauchy errors, no method achieved an F1 score above 0.65, confirming the theoretical negative result regarding support recovery.
• Parameter Selection: The RGCV criterion was shown to track oracle selection with a correlation above 0.77.

4.2 Real-Data Application

The methods were applied to TCGA Breast Cancer genomic data ($n=312, p=8,942$).

• Diagnostics: Residual analysis indicated heavy tails (excess kurtosis 4.7, consistent with $t(4.3)$).
• Performance: The Adaptive Huber-Lasso reduced cross-validated prediction error by 30% (0.203 vs. 0.291) compared to standard Lasso.
• Biological Relevance: The selected genes showed a 46% improvement in enrichment for BRCA1 regulatory pathways (Enrichment Ratio 0.704 vs. 0.477 for Lasso), with the result being statistically significant ($p_{perm} < 0.01$).

5. Significance and Claims

The paper claims to provide a unified theoretical and practical framework for high-dimensional regression that extends robustness to distributions with only $\bar{\alpha} > 1$ finite moments, a regime previously excluded by standard theory.

• Theoretical Advancement: It bridges the gap between classical robust statistics and modern high-dimensional theory by proving oracle inequalities under weaker moment conditions and explicitly characterizing the limits of variable selection under heavy tails.
• Practical Utility: It introduces the RGCV criterion to solve the joint tuning problem, which is critical for heavy-tailed data where standard cross-validation fails.
• Empirical Validation: Through simulations and genomic application, the paper demonstrates that robust regularised estimators (particularly Adaptive Huber-Lasso and Catoni-Lasso) offer substantial gains in prediction accuracy and biological interpretability when data deviates from Gaussian assumptions, while maintaining near-optimal efficiency under Gaussian conditions.

The authors maintain a modest stance, acknowledging that variable selection is not reliable under near-Cauchy errors and that the computational cost of RGCV (approx. 37x that of standard Lasso CV) may be prohibitive for extremely high-dimensional datasets ($p > 50,000$) without specialized hardware.