PaperlayerView exact revision on arXiv

Source-linked AI summary

Adaptive Huber Regression

Qiang Sun, Wenxin Zhou, Jianqing Fan

arXiv:1706.06991v2math.STstat.ME

TL;DR

Heavy-tailed and asymmetric data make conventional regression methods unreliable, motivating robust estimation without sub-Gaussian assumptions. The paper proposes adaptive Huber regression, tuning robustification to sample size, dimension, and moments. It establishes a sharp phase transition: sub-Gaussian-type deviations when δ ≥1 and slower rates when 0 < δ < 1, with near-optimal guarantees across low- and high-dimensional settings.

  • Problem

    Heavy-tailed data and asymmetric errors challenge conventional regression, while fixed Huber robustification can prevent consistent estimation under asymmetric distributions.

  • Method

    The paper uses Huber regression with a robustification parameter adapted to sample size, dimension, and available moments for robust estimation and inference.

  • Results

    The estimator shows a sharp, near-optimal phase transition: sub-Gaussian-type deviations when δ ≥1 and slower concentration when 0 < δ < 1.

  • Takeaways & Limitations

    Adaptive robustification combines robustness with asymptotic unbiasedness and supports regression with heavy-tailed errors, including heteroscedastic models.

  • Takeaways & Limitations

    For heavy-tailed observation noise, whether the sharper high-dimensional bound can be achieved by a Huber-type regularized estimator remains future work.

Abstract

from arXiv · show

Big data can easily be contaminated by outliers or contain variables with heavy-tailed distributions, which makes many conventional methods inadequate. To address this challenge, we propose the adaptive Huber regression for robust estimation and inference. The key observation is that the robustification parameter should adapt to the sample size, dimension and moments for optimal tradeoff between bias and robustness. Our theoretical framework deals with heavy-tailed distributions with bounded $(1+δ)$-th moment for any $δ> 0$. We establish a sharp phase transition for robust estimation of regression parameters in both low and high dimensions: when $δ\geq 1$, the estimator admits a sub-Gaussian-type deviation bound without sub-Gaussian assumptions on the data, while only a slower rate is available in the regime $0<δ< 1$. Furthermore, this transition is smooth and optimal. In addition, we extend the methodology to allow both heavy-tailed predictors and observation noise. Simulation studies lend further support to the theory. In a genetic study of cancer cell lines that exhibit heavy-tailedness, the proposed methods are shown to be more robust and predictive.

1 Introduction

Heavy-tailed data challenge methods built on sub-Gaussian assumptions, while fixed-parameter Huber regression can be inconsistent under asymmetric errors. The paper develops adaptive Huber regression and establishes nonasymptotic, near-optimal guarantees across moment and dimensional regimes.

  • Motivation: Heavy-tailed distributions occur in applications including fMRI, gene expression, and finance, making sub-Gaussian assumptions and least-squares procedures potentially unrealistic or unreliable.The cited discussion gives examples of non-Gaussian fMRI data, heavy-tailed gene expression levels, power-law financial returns, and poor least-squares behavior under heavy tails.
  • Related work: Fixed robustification based on a 95% asymptotic-efficiency rule can prevent consistent estimation of model-generating parameters when the sample distribution is asymmetric.The issue arises because the classical robustification parameter is held fixed rather than adapted to the estimation setting.
  • Contribution: Adaptive Huber regression uses a robustification parameter that adapts to sample size, dimension, and moments for robust estimation and inference.The proposed estimator has exponential-type concentration under low-order moments and is asymptotically unbiased without symmetry or homoscedasticity assumptions on errors.
  • Contribution: The paper establishes matching nonasymptotic upper and lower bounds, a sharp phase transition, high-dimensional regularized results, and a Bahadur representation when δ ≥1.The framework unifies low- and high-dimensional analyses through effective dimension and effective sample size, while localized analysis removes an artificial bounded-parameter constraint.
  • Related work: The paper distinguishes its conditional-mean target from quantile regression, which is generally biased for mean-regression coefficients under asymmetric errors.This distinction follows from the different target and the role of error-distribution symmetry.

2 Methodology

The methodology replaces least squares or fixed-parameter Huber fitting with an adaptively tuned Huber loss that balances bias and robustness. Its guarantees exhibit a phase transition at δ = 1 across low- and high-dimensional settings.

  • Huber loss: Adaptive Huber loss is quadratic for small residuals and linear beyond τ, down-weighting outliers while τ balances approximation bias against robustness.A larger τ reduces bias but weakens robustness; the adaptive choice depends on sample size, dimension, and moments.
  • Phase transition: For fixed effective dimension, the ℓ2-error rate is n^-δ/(1+δ), with the phase transition at δ = 1 yielding sub-Gaussian-type behavior when δ ≥1.The unified bounds use effective dimension d_eff and effective sample size n_eff, with n_eff = n in low dimensions and n_eff = n/log d in high dimensions.
  • High dimensions: The regularized high-dimensional estimator achieves ℓ2-error of order s^1/2{(log d)/n}^min{δ/(1+δ),1/2} with high probability.Here s is the size of the true support, and the high-dimensional analysis uses regularized adaptive Huber regression.

3 Nonasymptotic Theory

The theory tunes Huber’s robustification parameter to balance approximation bias and robustness under finite (1+δ)-th moments. Matching bounds establish near-optimal rates and a sharp phase transition across low- and high-dimensional settings.

  • Low-dimensional theory: Adaptive Huber regression tunes τ to balance estimation error and approximation bias, whose order is τ^-δ.Increasing τ reduces bias but compromises robustness.
  • Low-dimensional theory: Theorem 1 uses τ = τ0(n/t)^max{1/(1+δ),1/2} and yields exponential-type deviation bounds under finite (1+δ)-th moments.The result assumes bounded predictor coordinates and a fixed-design scaling condition.
  • Extensions and assumptions: The method accommodates heteroscedastic regression without requiring E(|ε_i|^(1+δ)|x_i) to be constant.The random-design analysis is provided in supplementary material.
  • Phase transition: The adaptive estimator converges at rate n^-1/2 when δ ≥ 1, but only n^-δ/(1+δ) when the second moment does not exist.Matching lower bounds show the transition is near-optimal, and root-n exponential concentration is impossible for 0 < δ < 1.
  • Inference and efficiency: The estimator has a nonasymptotic Bahadur representation and asymptotically the same efficiency as ordinary least squares.The representation provides a linear approximation for future inference and confidence-set construction.
  • High-dimensional theory: In high dimensions, ℓ1-regularized adaptive Huber regression achieves rate s^1/2{(log d)/n}^min{δ/(1+δ),1/2} with overwhelming probability.With finite-variance observation noise, the stated rate matches the Lasso rate for sub-Gaussian errors.
  • High-dimensional theory: The regularized estimator attains minimax-optimal ℓ1-, ℓ2-, and prediction-error bounds under random designs with n ≳ s log d.Under fixed designs, the stated scaling is n ≳ s^2 log d.

4 Extension to Heavy-tailed Designs

The paper extends adaptive Huber regression to high-dimensional settings with heavy-tailed covariates and errors by truncating covariates and regularizing the estimator. Under finite fourth-moment covariate conditions, the modified estimator has exponential concentration, with near-optimal rates under additional boundedness assumptions.

  • Scope: The extension targets sparse high-dimensional regression where both covariates and regression errors may be heavy-tailed.The method focuses on d ≫ n with sparsity s = ∥β∗∥0 ≪ n.
  • Method: Covariates are robustified coordinatewise by truncating each component to the interval [−ϖ, ϖ].The truncated vector is combined with a modified adaptive Huber estimator and regularization.
  • Theory: The modified estimator admits exponentially fast concentration when covariates have only finite fourth moments, subject to stronger scaling conditions.
  • Theory: Under bounded v3, M4, and ∥β∗∥2, the regularized estimator achieves a near-optimal convergence rate.The stated result uses t ≍ log d in the high-dimensional bound.
  • Theory: The theoretically optimal robustification parameter differs from the one used for sub-Gaussian designs.

5 Algorithm and Implementation

The implementation solves regularized adaptive Huber regression with a local adaptive majorize-minimization algorithm. Each iteration uses a locally valid isotropic quadratic surrogate, soft-thresholding, and parameter inflation until the surrogate majorizes the objective.

  • Motivation: Standard cutting-plane and interior-point methods are not scalable for the large-scale convex optimization problem.The implementation therefore focuses on a local adaptive majorize-minimization approach.
  • LAMM algorithm: LAMM iteratively minimizes a local quadratic majorizer of the adaptive Huber objective.The majorization requirement needs to hold locally at the next iterate.
  • LAMM algorithm: An isotropic quadratic surrogate yields a simple analytic update for the regularized problem.The update is expressed through the soft-thresholding operator.
  • LAMM algorithm: The algorithm starts with a small quadratic parameter and repeatedly multiplies it by γu > 1 until the majorization condition is satisfied.The paper gives γu = 2 as an example.
  • Stopping: Iterations stop when successive coefficient vectors differ by at most ϵ, with ϵ = 10^-4 suggested as a simple criterion.The resulting sequence is described as convergent.

6 Numerical Studies

Numerical studies compare adaptive Huber regression with least squares across light- and heavy-tailed settings, examine phase-transition and effective-sample-size scaling, and evaluate genomic prediction. The results support robustness under heavy tails with little loss under normal errors.

  • Finite-sample performance: Adaptive Huber regression performs comparably to least squares under normal noise but significantly better under Student’s t and log-normal errors.The comparison uses ℓ2-error averaged over 100 simulations.
  • Phase transition: The empirical negative log ℓ2-error curve closely follows theory in both low and high dimensions, supporting the predicted phase transition.Figure 2 uses 200 repetitions for each (n, d) combination.
  • Phase transition: Adaptive Huber regression has a significant advantage over OLS when δ is small, while OLS gradually catches up as δ increases.
  • Effective sample size: In regularized experiments, ℓ2-error decreases with sample size, and plotting against n/log d aligns curves across dimensions.The alignment is reported for d ∈ {100, 500, 5000}.
  • NCI-60 data: NCI-60 protein and gene-expression variables show substantial empirical heavy-tailedness relative to the normal kurtosis benchmark of 3.More than 89.5% of protein variables have heavier tails than the normal distribution, and about 36.5% of gene variables exceed the kurtosis of t5.
  • NCI-60 data: For KRT19 prediction, TAHuber has the smallest MAE, followed by AHuber and Lasso, while Lasso selects a comparatively large model.The methods are evaluated by leave-one-out cross-validation.

A A Lepski-type method

The method adapts the Huber robustification parameter using a Lepski-type procedure, without requiring the variance to be known in advance. It computes estimators across a geometrically spaced parameter grid and selects a data-driven estimator with an explicit probability bound.

  • A A Lepski-type method: The robustification parameter is adapted through Lepski’s method without knowing the variance in advance.The construction requires preliminary upper and lower variance bounds.
  • A A Lepski-type method: The procedure assumes crude bounds σmin ≤ v1/2 ≤ σmax and forms a geometric grid σj = σmin a^j.The grid includes indices satisfying σmin ≤ σj < aσmax.
  • A A Lepski-type method: For each grid value, it computes a Huber estimator with τj = σj(n/t)1/2 before selecting the final data-driven estimator.The candidate set has cardinality at most 1 + log_a(σmax/σmin).
  • A A Lepski-type method: The selected estimator satisfies an explicit deviation bound with probability at least 1 − (2d + 1) log_a(aσmax/σmin)e^−t.The result requires n ≥ 8 max(4eL^2d, eL^4d^2)t.
  • A A Lepski-type method: In practice, σmin and σmax can be constructed from the least-squares residual variance estimate using a multiplicative factor K > 1.The effectiveness depends on the sharpness of theoretical constants, which may not be sharp.

B.1 Huber regression in low dimensions

In low dimensions, Huber regression uses a calibrated robustification parameter to obtain exponential-type concentration under finite (1 + δ)-th moments. With finite variance, it also provides a finite-sample linear approximation and high efficiency as the parameter diverges.

  • B.1 Huber regression in low dimensions: The low-dimensional Huber estimator is defined by minimizing empirical Huber loss with a tunable robustification parameter τ.The estimator is analyzed under random-design moment conditions.
  • B.1 Huber regression in low dimensions: For finite (1 + δ)-th moments, properly calibrated τ yields exponential-type concentration inequalities for the estimator.The theorem applies when n ≥ C2(d + t), with constants depending only on A0.
  • B.1 Huber regression in low dimensions: Under finite variance, the estimator admits a nonasymptotic Bahadur representation with a higher-order remainder having sub-exponential tails.Its distribution is governed mainly by a linear stochastic term.
  • B.1 Huber regression in low dimensions: The adaptive estimator balances nonasymptotic robustness against heavy-tailed errors with high efficiency when τ diverges to infinity.This conclusion follows from the behavior of truncated moments and the estimator’s asymptotic representation.
  • B.1 Huber regression in low dimensions: With t = log n and n ≳ d, the recommended scaling is τ ≍ [n/(d + log n)] raised to the exponent determined by δ.The resulting bound holds with probability at least 1 − O(n^−1).

B.2 Huber regression in high dimensions

In high dimensions, the paper studies an ℓ1-regularized Huber estimator for sparse regression under heavy-tailed errors. Under finite variance and suitable tuning, it achieves minimax-order estimation rates, while the analysis relies on restricted strong convexity.

  • B.2 Huber regression in high dimensions: The high-dimensional method uses an ℓ1-regularized Huber estimator when d ≫ n and the true coefficient vector is sparse.τ and λ serve as the robustification and regularization parameters.
  • B.2 Huber regression in high dimensions: Theorem 8 establishes the estimator’s statistical consistency for suitable robustification and regularization parameters.The result assumes a sparse β∗ and Condition 5.
  • B.2 Huber regression in high dimensions: Under finite variance, ℓ1- and ℓ2-errors scale as s√(log d/n) and √(s log d/n), respectively, when n ≳ s log d.These are the minimax rates associated with standard Lasso under Gaussian or sub-Gaussian errors.
  • B.2 Huber regression in high dimensions: Under random designs, the required sample size is O(s log d), compared with O(s^2 log d) under fixed designs.The random-design scaling is attributed to restricted strong convexity of Huber loss near β∗.
  • B.2 Huber regression in high dimensions: The fixed-design sparsity factor in one theorem is identified as a proof artifact, while achieving an oracle excess-risk rate requires sparsity of order O(√(n/log n)).Equivalently, the required sample size scales as s^2 log n.
  • B.2 Huber regression in high dimensions: The restricted strong convexity analysis controls the adaptive Huber loss locally over an ℓ1-cone around the true parameter.This analysis uses moment conditions on errors and covariates.

C.3 Proof of Theorem 1

The proof combines concentration for truncated score terms with local curvature arguments to establish the low-dimensional upper bound. It also uses a construction showing that the rate cannot generally be improved under only finite (1 + δ)-th moments.

  • C.3 Proof of Theorem 1: An intermediate estimator is introduced along the segment from β∗ to the Huber estimator to localize the proof within a radius-r ball.The argument then shows the intermediate estimator lies in the ball’s interior, forcing it to equal the Huber estimator.
  • C.3 Proof of Theorem 1: The proof uses the stationarity condition ∇Lτ(β̂) = 0 and the vector-valued mean value theorem to relate estimation error to local curvature and score terms.This connects the estimator’s optimality equation to the deviation analysis.
  • C.3 Proof of Theorem 1: Truncated score coordinates are controlled using Markov-type tail bounds, yielding exponential control when τ is sufficiently large relative to the moment scale.The resulting coordinate bound is combined across dimensions.
  • C.3 Proof of Theorem 1: The proof establishes the required event with probability at least 1 − 2e^−t under n ≥ C1(d + t).The constant C1 depends only on A0.
  • C.3 Proof of Theorem 1: A matching lower-bound construction shows that root-n consistency with exponential concentration is impossible when 0 < δ < 1.The construction uses distributions with bounded (1 + δ)-th moments and a non-vanishing design-angle condition.
  • C.3 Proof of Theorem 1: The lower-bound argument specializes naturally to the mean model with X = (1, …, 1)^T and uses a sign vector satisfying a minimum-angle condition.This links the regression lower bound to the corresponding univariate mean-estimation phenomenon.

C.5 Proof of Theorem 3

The proof of Theorem 3 constructs an intermediate estimator within an ℓ1-radius and uses cone and gradient bounds to establish the stated result. It also derives high-probability bounds under moment-dependent choices of the Huber parameter.

  • Curvature control: The proof also establishes uniform control of the Hessian quadratic form over the cone with probability at least 1 − e^-t, under a lower bound on τ and n.The upper curvature bound ⟨u, Hτu⟩ ≤ κu completes the lemma used in the theorem proof.
  • Proof strategy: Lemma 8 places the estimator in an ℓ1-cone whenever the gradient at β∗ is bounded by λ/2.For any E containing the support S, the off-support ℓ1-error is at most three times the on-set error.
  • Proof strategy: The proof restricts attention to δ ∈ (0, 1] and constructs bβη between β∗ and bβ so that its ℓ1-error is at most r.The interpolation parameter is η = 1 when the target radius is already met and otherwise is chosen to reach the radius exactly.
  • Error control: The intermediate estimator satisfies ∥bβη − β∗∥1 ≤ 4s^1/2∥bβη − β∗∥2 ≤ 12κ_l^-1 sλ < r, forcing bβη = bβ.This closes the error argument for the theorem once the stated inequalities hold.
  • Probability bound: Choosing τ = τ0(n/t)^(1/(1+δ)) with τ0 ≥ νδ gives P{∥∇Lτ(β∗)_S∥∞ ≥ 2Lτ n^-1 t} ≤ 2s e^-t.The bound is combined with the preceding inequality to prove the stated result.
  • Probability bound: Taking t = (1 + c) log d yields probability at least 1 − (2s + 1)d^(-1-c) for the corresponding gradient bound.Condition 3 with k = 2s implies 2s + 1 ≤ d, giving the displayed conclusion.

C.7 Proof of Theorem 5

The proof of Theorem 5 follows the argument for Theorem 3 and derives a coordinatewise probability bound for the relevant gradient event. Bernstein’s inequality and a union bound over the support yield the stated result.

  • Proof strategy: The proof is almost identical to that of Theorem 3 and focuses on deriving a probability bound for the event involving ξ∗.The argument proceeds under the restriction 0 < δ ≤ 1.
  • Proof strategy: The proof decomposes heavy-tailed predictors into truncated components and residuals, incorporating the residual contribution into ϵi = εi + ⟨zi, β∗⟩.The resulting gradient quantity is represented using ψϖ(xij) and the residual-adjusted errors.
  • Concentration: For each fixed coordinate, Bernstein’s inequality provides the required deviation bound.The coordinatewise argument is then aggregated across the support.
  • Concentration: A union bound over j ∈ S gives probability at least 1 − 2s e^-t for the resulting simultaneous bound.This is identified as the stated result of the theorem.

C.8 Proof of Theorem 7

The proof of Theorem 7 combines deviation inequalities for the normalized Huber gradient with restricted strong convexity and an intermediate estimator localized in a covariance norm. It then calibrates τ, λ, and r to obtain high-probability error bounds.

  • Proof of (19): The proof derives deviation inequalities for ∥Σ^-1/2∇Lτ(β∗)∥2 and establishes restricted strong convexity for the Huber loss under vδ < ∞.The parameter set Θ0(r) localizes β around β∗ in the Σ,2 norm.
  • Proof of (19): With probability at least 1 − 2e^-t, the intermediate estimator satisfies ∥bβτ,η − β∗∥Σ,2 ≤ 4r0 < r when n ≥ C1(d + t).This conclusion follows by combining the preceding bounds and the localized construction.
  • Proof of (20): The empirical-process argument bounds B(bβτ) by separately controlling its centered process and its expectation over Θ0(r).The mean-value theorem is used for the expectation component, while Theorem A.3 controls the empirical process.
  • Proof of (20): The Huber-score bias and second-moment terms are bounded using tail integrals, including a τ^-κ dependence for the higher-moment contribution.These calculations use E(ε) = 0 and moment bounds on ε.
  • Sparse error analysis: The sparse estimator is placed in an ℓ1-cone under λ ≥ 2∥∇Lτ(β∗)∥∞, enabling the subsequent error analysis.The cone follows from convexity, optimality, and the ℓ1 penalty structure.
  • Sparse error analysis: Choosing λ proportional to max_j σjj^1/2 τ0{(log d)/n}^(δ/(1+δ)) controls the gradient with probability at least 1 − 2d^-1.The restricted strong convexity condition is then combined with this calibration to prove the stated sparse bound.
Loading 1706.06991v2…