PaperlayerView exact revision on arXiv

Source-linked AI summary

Smoothed Quantile Regression with Large-Scale Inference

Xuming He, Xiaoou Pan, Kean Ming Tan, Wen-Xin Zhou

arXiv:2012.05187v2math.STstat.ME

TL;DR

Large-scale quantile regression requires inference and computation that remain reliable as predictor dimension grows, while conventional nonsmooth optimization and density-based variance estimation can be problematic. The paper develops conquer, a convolution-smoothed quantile regression method with scalable gradient optimization and multiplier bootstrap inference, and establishes non-asymptotic error bounds, asymptotic normality under a weaker dimension condition than standard QR, and bootstrap validity.

  • Problem

    Large-scale quantile regression faces computational bottlenecks from nonsmooth linear-programming optimization and instability in conditional-density-based variance estimation when p grows with n.

  • Method

    Conquer applies convolution-type smoothing to obtain a convex, twice-differentiable, locally strongly convex loss, then uses gradient-based optimization and multiplier bootstrap inference.

  • Results

    The analysis establishes non-asymptotic estimation and Bahadur-Kiefer bounds, asymptotic normality when p8/3/n → 0, and validity of multiplier bootstrap confidence constructions.

  • Takeaways & Limitations

    Conquer provides a computationally scalable framework for large-scale quantile regression inference in increasing-dimension settings.

  • Takeaways & Limitations

    Normal-based confidence intervals can be exceedingly wide and inaccurate at high or low quantile levels, while one warm-start idea can have higher estimation error than conventional QR and conquer.

Abstract

from arXiv · show

Quantile regression is a powerful tool for learning the relationship between a response variable and a multivariate predictor while exploring heterogeneous effects. In this paper, we consider statistical inference for quantile regression with large-scale data in the "increasing dimension" regime. We provide a comprehensive and in-depth analysis of a convolution-type smoothing approach that achieves adequate approximation to computation and inference for quantile regression. This method, which we refer to as {\it{conquer}}, turns the non-differentiable quantile loss function into a twice-differentiable, convex and locally strongly convex surrogate, which admits a fast and scalable Barzilai-Borwein gradient-based algorithm to perform optimization, and multiplier bootstrap for statistical inference. Theoretically, we establish explicit non-asymptotic bounds on both estimation and Bahadur-Kiefer linearization errors, from which we show that the asymptotic normality of the conquer estimator holds under a weaker requirement on the number of the regressors than needed for conventional quantile regression. Moreover, we prove the validity of multiplier bootstrap confidence constructions. Our numerical studies confirm the conquer estimator as a practical and reliable approach to large-scale inference for quantile regression. Software implementing the methodology is available in the \texttt{R} package \texttt{conquer}.

1 Introduction

The paper addresses computational and inferential challenges for quantile regression when both sample size and predictor dimension grow. It analyzes conquer, a convolution-smoothed approach designed to retain convexity while enabling scalable optimization and large-scale inference.

  • Motivation: Large-scale quantile regression becomes difficult because linear-programming approaches can be slow or memory-intensive when both n and p are large.The paper illustrates this challenge with a dataset exceeding 4 million observations and 920 stock-level covariates.
  • Motivation: Increasing-dimension inference is also challenged by unstable conditional-density variance estimation and uncertain finite-sample accuracy of normal approximations.The paper studies settings with p ≍ n^a for a ∈ (0,1).
  • Approach: Conquer uses convolution-type smoothing to replace the non-differentiable quantile loss with a convex, twice-differentiable, locally strongly convex objective.The method is motivated partly by the optimization difficulties of alternative smoothing that sacrifices convexity.
  • Approach: The method combines bandwidth selection adapted to n and p with a Barzilai-Borwein gradient algorithm and Huber-type initialization for scalable computation.The proposed implementation is provided in the R package conquer.
  • Theory: With probability at least 1−3e−t, the Bahadur-Kiefer linearization error is bounded at order (p + t)/(nh1/2) + h3/2 p.Under a suitable bandwidth, the linear functional is asymptotically normal when p8/3/n → 0, improving the cited condition for standard QR.
  • Empirical study: Simulation studies compare conquer with quantreg’s pfn implementation using estimation error and elapsed time across increasing sample sizes, supporting its practical feasibility for large-scale analysis.The experiment uses 100 Monte Carlo samples with p = ⌊n1/2⌋.

2 Smoothed quantile regression

Convolution-type smoothing produces the conquer estimator, a twice-differentiable and globally convex alternative to the nonsmooth quantile-regression objective. The method supports efficient gradient-based computation and theoretical inference in increasing dimensions, including bootstrap validity and explicit error bounds.

  • Convolution-type smoothing: The conquer loss is globally convex, whereas Horowitz’s smoothed loss is not.This distinction is central to the computational advantages of convolution smoothing.
  • Bandwidth selection: The bandwidth should adapt to sample size n and dimension p to balance statistical accuracy and computational stability.The paper treats bandwidth selection as depending on the increasing-dimension regime.
  • Convolution-type smoothing: Convolution smoothing replaces the nonsmooth check loss with an objective that is twice continuously differentiable and convex when the kernel is non-negative.The resulting estimator satisfies a first-order condition, enabling gradient-based optimization.
  • Statistical properties: Conquer has lower asymptotic mean squared error than Horowitz’s smoothed estimator and smaller Bahadur linearization error than standard QR when p is fixed.These comparisons are reported for earlier fixed-dimensional analyses of the conquer estimator.
  • Statistical properties: The paper derives finite-sample bounds for estimation, Bahadur-Kiefer linearization, and multiplier-bootstrap counterparts in increasing dimensions.The analysis also establishes concentration inequalities and non-asymptotic Bahadur representations.
  • Inference: Theoretical analysis establishes asymptotic normality and valid multiplier-bootstrap inference under increasing dimension, with a variance expression involving J^-1ΣJ^-1.The bootstrap can remain valid after variance adjustment even when the random weights do not have unit variance.
  • Inference: Normal-based confidence intervals can be excessively wide and inaccurate at high or low quantile levels when estimated density matrices are unstable.The issue is linked to sensitivity in estimating and inverting the conditional-density matrix when p is large.
  • Computation: Conquer’s gradient-based algorithm improves computational efficiency for smoothed quantile regression with many regressors and massive sample sizes.The method uses non-negative kernels and local strong convexity to support tailored computation.

3 Computational methods for conquer

The computational methods smooth the quantile loss into a differentiable objective and solve it with gradient-based procedures designed for large-scale optimization. The section combines Barzilai–Borwein updates with asymmetric Huber warm starts, while addressing conditioning and implementation stability.

  • The smoothed objective is twice differentiable, enabling gradient descent while avoiding computationally intensive line searches.Its Hessian is explicit, but standard stepsize schedules may be inefficient in large-scale settings.
  • The convolution-smoothed check function has explicit forms for Gaussian, logistic, uniform, Epanechnikov, and triangular kernels.The uniform-kernel form is proportional to a shifted Huber loss, while the other kernels yield analogous smooth approximations.
  • The GD-BB algorithm uses Barzilai–Borwein stepsizes to reduce the cost of approximating inverse-Hessian updates in high dimensions.The method derives stepsizes from successive iterates and gradients rather than explicitly inverting the Hessian.
  • Newton optimization can converge slowly when the bandwidth is small or the quantile level is near 0 or 1 because the Hessian becomes ill-conditioned.Backtracking line searches also require repeated loss evaluations, motivating explicit smoothed-loss formulas for common kernels.
  • 3.2 Warm start via asymmetric Huber regression: An asymmetric Huber M-estimator provides a warm start for GD-BB because its loss is twice continuously differentiable, convex, and locally strongly convex.The asymmetric Huber loss approximates the check function as its shape parameter γ approaches 0.
  • 3.2 Warm start via asymmetric Huber regression: The shrinking-γ alternative is not numerically stable across simulated datasets without careful control of the minimum γ and has higher estimation error than conventional QR and conquer.This result limits the practical appeal of using asymmetric Huber regression with a progressively shrinking shape parameter.

4 Statistical analysis

The statistical analysis establishes non-asymptotic guarantees for the smoothed conquer estimator and its bootstrap inference in increasing dimensions. It derives concentration, Bahadur-Kiefer, normal-approximation, and bootstrap results under explicit bandwidth and dimensionality conditions.

  • Smoothing bias and estimation: Under random-design moment and density conditions, the analysis bounds smoothing bias and estimation error for the conquer estimator.The assumptions include bounded third moments for standardized covariates, a positive conditional error density at zero, and Lipschitz density variation.
  • Bahadur-Kiefer representation: The paper establishes concentration and Bahadur-Kiefer representations for conquer, yielding explicit bounds for the estimator and its linearization remainder.The Bahadur representation requires a stronger sub-Gaussian covariate condition and supports subsequent linear-functional inference.
  • Smoothing bias and estimation: Theorem 4.1 provides a high-probability estimation bound combining a bandwidth-dependent bias term with a statistical rate involving (p + t)/n.The confidence level is controlled by t through 1 − 2e^-t, while an additional log^2(1/h) factor arises from the proof’s peeling argument.
  • Normal approximation: With h = {(p + log n)/n}^2/5, linear functionals are asymptotically normal when p^8/3/n → 0, improving the reported growth condition for standard quantile regression.The corresponding Gaussian approximation error is of order (p + log n)^4/5 n^-3/10.
  • Approximation accuracy: The inference approximation has error faster than n^-1/4 but slower than n^-1/2, while the paper notes that alternative standard-quantile expansions can attain n^-1/2 up to logarithmic factors.This comparison qualifies the accuracy rate of the conquer-based Bahadur and normal approximations.
  • Bootstrap inference: The bootstrap analysis develops a Bahadur representation for the bootstrap estimator and establishes multiplier-bootstrap validity under explicit sample-size, bandwidth, and dimensionality conditions.Rademacher multiplier bootstrap avoids estimating potentially unstable variance-covariance matrices, especially near extreme quantile levels; RMB can also operate on a subsample of about n/2 observations.

5 Numerical studies

Numerical studies compare conquer with standard and smoothed quantile regression for estimation, computation, bandwidth sensitivity, and bootstrap inference. Conquer generally improves estimation accuracy, remains computationally efficient, and provides confidence intervals balancing coverage, width, and runtime.

  • Estimation: Conquer has lower estimation error than classical QR across all scenarios and than Horowitz’s smoothing in most settings.Advantages are especially evident at low and high quantile levels.
  • Computational efficiency: Conquer is computationally efficient and stable across scenarios, while classical QR becomes slower as sample size and dimension increase.Its runtime is insensitive to the choice of kernel functions.
  • Bandwidth sensitivity: Conquer’s estimation error remains uniformly below QR across a range of bandwidths, including the default hde = {(p + log n)/n}2/5.The default bandwidth performs comparably to an AMSE-based bandwidth using oracle knowledge.
  • Bootstrap inference: The multiplier-bootstrap methods mb-per, mb-piv, and mb-norm balance high coverage probability with narrow confidence intervals and achieve the lowest runtime.At τ = 0.9, pwy is most conservative with the widest intervals, whereas wild gives the narrowest intervals at the cost of coverage.
  • Bootstrap inference: At (n, p) = (5000, 250), rank inversion takes approximately 80 minutes, compared with 41 seconds for conquer with multiplier bootstrap.Rank inversion is therefore omitted from the larger-scale comparison.
  • Bootstrap inference: Bootstrap calibration produces more reliable and accurate confidence intervals than normal calibration for extreme quantile regression under the reported t1.5-error setting.The comparison concerns empirical coverage and confidence-interval width.

6 Discussion

The discussion presents conquer as a convex, twice-differentiable approach for increasing-dimensional quantile regression, with gradient-based computation and multiplier-bootstrap inference. Higher-order one-step conquer further relaxes the dimension constraint under additional smoothness, while regularized extensions are motivated by the same computational structure.

  • Conquer properties: With a non-negative kernel, the smoothed objective is convex, twice continuously differentiable, and locally strongly convex near β∗ with high probability.These properties support efficient gradient-based optimization.
  • Conquer properties: Convolution smoothing balances statistical accuracy and computational complexity for large-scale quantile regression.The discussion contrasts this balance with generic optimization that can be slow or memory-intensive.
  • Higher-order smoothing: Higher-order one-step conquer is asymptotically normal under a dimension constraint of roughly p^2/n → 0 when increasing smoothness is available.It uses higher-order kernels without solving a large-scale non-convex optimization problem.
  • Regularized extensions: Because conquer is convex and twice differentiable, gradient-based algorithms are expected to improve computational efficiency for regularized conquer without sacrificing statistical accuracy.The discussion specifically mentions coordinate gradient descent and proximal gradient descent.

A Asymptotic mean squared error

The analysis characterizes the asymptotic mean squared error of the smoothed estimator through bandwidth-dependent bias and variance terms. It also connects kernel choice and local convexity to estimator accuracy and computational feasibility.

  • Conditions and derivation: Under the stated density, covariate, and bandwidth conditions, the estimator’s asymptotic behavior is obtained by Taylor expansions, concentration bounds, and explicit remainder control.The derivation uses conditional-density smoothness and moment assumptions to control bias and variance terms.
  • AMSE: The asymptotic mean squared error is determined by the leading term of the smoothed estimator’s mean squared error under the prescribed bandwidth constraint.The analysis derives this leading term after expanding the relevant variance expression.
  • Kernel smoothing: Second-order kernels yield smoothing bias of order h^2, while the smoothed empirical loss remains twice-differentiable and convex.The loss is also strongly convex locally with high probability.
  • Higher-order kernels: Higher-order kernels can reduce smoothing bias, but the resulting convolution-smoothed objective is non-convex and may be computationally intractable to optimize globally.The paper instead motivates a one-step estimator to avoid solving this large-scale non-convex problem.
  • One-step estimator: The one-step conquer estimator starts from a second-order-kernel estimator and applies a higher-order-kernel Newton-type correction.The correction is defined through a gradient equation and uses a positive-definite Hessian-like matrix.
  • Local analysis: The analysis establishes local convergence control through convexity, first-order optimality, and restricted strong convexity of the empirical smoothed loss.Uniform empirical-process bounds control the loss over local ellipsoidal neighborhoods.

C.3 An alternative proof to Theorem 4.1

The alternative proof removes an extra logarithmic factor by refining the local restricted strong convexity argument. It first obtains a coarse radius and then sharpens the estimator bound within that neighborhood.

  • Alternative proof: The peeling proof introduces an additional log^2(1/h) term, which is bounded by the slowly growing log(log n) under the prescribed bandwidth constraint.The alternative analysis treats this term as a proof artifact rather than an intrinsic rate component.
  • Alternative proof: The refined argument conditions on the estimator lying in a local ellipsoid and uses the first-order optimality condition to improve its convergence rate.The key condition is ∇b Qh(bβh) = 0.
  • Local curvature: A symmetrized Bregman divergence is used to quantify local curvature of the convex empirical loss.The divergence is nonnegative and compares gradients at two parameter values.
  • Assumptions: The local curvature bounds rely on upper and lower conditional-density assumptions and standardized-covariate moment control.The quantity ιδ summarizes a tail condition for standardized covariate projections.
  • Result: The resulting probability bounds hold uniformly over local neighborhoods and yield the refined convergence bound without the additional logarithmic factor.The lemmas are combined on a high-probability good event, with bandwidth restrictions controlling the radius.

C.3.1 Proof of Lemma C.2

The proof of Lemma C.2 controls the stochastic gradient discrepancy using sphere coverings, bounded increments, conditional variance bounds, and Bernstein’s inequality.

  • Gradient representation: The normalized gradient discrepancy is represented as an average of centered vectors involving smoothed residual scores and standardized covariates.The representation defines ξi from the smoothed indicator and wi from the standardized covariates.
  • Sphere covering: An ε-net of the unit sphere reduces uniform control over vector directions to finitely many scalar projections.The net cardinality grows exponentially with the dimension p.
  • Moment bounds: The projected variables are bounded, and their conditional second moments are controlled by a constant involving τ, the density smoothness, and h^2.These bounds provide the inputs for a Bernstein concentration argument.
  • Concentration result: Bernstein’s inequality combined with the covering argument yields a high-probability uniform bound for the gradient discrepancy.The final probability bound is obtained by selecting a specific net resolution and setting the deviation parameter relative to p and t.

C.3.2 Proof of Lemma C.3

The proof of Lemma C.3 establishes a uniform lower bound on the local Bregman divergence. It separates the expected curvature from empirical fluctuations and controls the latter with empirical-process inequalities.

  • Curvature quantity: The empirical loss is convex and twice continuously differentiable, and its local curvature is represented through a symmetrized Bregman divergence.The proof uses this divergence as the quantity requiring a uniform lower bound.
  • Pointwise lower bound: For bounded parameter neighborhoods, events controlling residuals and covariate projections provide a pointwise lower bound on the curvature contribution.The kernel’s positive lower bound on a central interval is used on these events.
  • Truncation device: The truncated quadratic function ϕR is introduced to express the curvature contribution while preserving Lipschitz and scaling properties.These properties enable subsequent empirical-process control.
  • Empirical-process control: The expected curvature and random fluctuation are analyzed separately.The fluctuation term is controlled using a refined Talagrand inequality, symmetrization, and contraction.
  • Concentration result: The contraction argument reduces the relevant process to a class indexed by standardized covariate projections and applies a 1-Lipschitz transformation.The resulting concentration bound is combined with earlier inequalities to prove the lemma.

C.4 Proof of Theorem 4.2

The proof controls local fluctuations of the smoothed score process using moment, exponential-tail, and empirical-process bounds. These controls yield a high-probability bound for the linearization remainder.

  • Local fluctuation control: The local remainder Δ(δ) is decomposed into deterministic and centered stochastic fluctuations around the population Hessian J_h.The centered component is represented through the gradient process G(β), while the deterministic component is handled separately.
  • Moment bounds: Sub-Gaussian design assumptions provide exponential-moment and fourth-moment bounds needed to control products of projected covariates.The proof uses bounds on E exp{⟨w,u⟩^2/(2υ_1)^2} and related weighted fourth moments.
  • Concentration argument: The empirical-process argument applies a concentration theorem after verifying a moment condition with g = min{h/(4κ_u), 1/ f̄}(n/2)^1/2 and ν_0 = C_0h^−1/2.The resulting control is combined with earlier bounds on the estimator's local neighborhood.
  • Result: With probability at least 1 − 3e^−t, the remainder satisfies ∥Δ(bδ)∥2 ≲ (p + t)/(h^1/2n) + h^3/2 p.This follows after combining the local fluctuation bound with the estimator localization result.

C.5 Proof of Theorem 4.3

The proof derives variance and third-moment controls for a linear functional of the smoothed estimator, then applies Berry–Esseen arguments to establish normal approximation under bandwidth conditions.

  • Variance control: For a linear functional a, the smoothed score variance is {τ(1 − τ) + O(h)}∥J_h^−1a∥2_Σ for sufficiently small h.The variance expression is obtained from the kernel-smoothed score representation.
  • Normal approximation: The centered score is expressed through ξ_i = τ − K(−ε_i/h) and γ_i = ⟨J_h^−1a, x_i⟩, enabling moment-based Berry–Esseen control.The proof separately bounds the conditional bias and the third moment of γ_iξ_i.
  • Conclusion: The resulting bounds hold uniformly over x and a, and the proof obtains the refined Berry–Esseen bound (4.9).The uniformity follows because the preceding bounds are independent of x and a.

C.6 Proof of Theorem 4.4

The proof establishes bootstrap control for the smoothed loss and its local error process using Rademacher multipliers, covering arguments, and concentration inequalities.

  • Bootstrap construction: The weighted bootstrap loss remains convex, and its bootstrap counterpart is defined through non-negative multiplier weights.Convexity supports the local analysis of the bootstrap estimator.
  • Bootstrap error bounds: Lemma C.4 provides high-probability bounds for the bootstrap error process when n ≳ p + t.The bounds hold over prescribed annuli in the standardized parameter norm.
  • Conclusion: Combining the concentration bounds with the good-event construction yields the claimed bootstrap bound under the scaling n ≳ p + t.The argument concludes after applying the bounds across the relevant radial regions.
  • Concentration argument: The proof controls the bootstrap loss difference with independent Rademacher multipliers using conditional Chernoff bounds, contraction, discretization, and peeling.The resulting probability statements are conditional on the observed data and cover local parameter regions.

C.7 Proof of Theorem 4.5

The proof controls bootstrap Hessian and variance fluctuations over local parameter sets using isotropic sub-Gaussian design bounds, nets, and Bernstein-type inequalities. These controls establish the stated result under the prescribed bandwidth constraint.

  • Standardization: The local fluctuation is analyzed after transforming δ to v = Σ^1/2δ, so the parameter region becomes an Euclidean ball.The standardized covariates w_i = Σ^−1/2x_i are isotropic random vectors.
  • Hessian control: The empirical Hessian fluctuation is decomposed into net approximation, pointwise stochastic error, and population smoothness terms.The decomposition is controlled over sphere nets and local parameter nets.
  • Design concentration: Sub-Gaussian design assumptions yield a uniform covariate-norm bound ∥w_i∥2^2 ≤ C_1υ_1^2(p + log n + t) with probability at least 1 − e^−t.This bound is obtained by applying a concentration inequality and taking a union bound over observations.
  • Uniform control: The proof combines net cardinality bounds, higher-order moment estimates, Bernstein inequalities, and union bounds to control the local processes.These tools produce the required high-probability bounds under the sample-size scaling n ≳ p + t.
  • Variance control: The same argument controls the bootstrap variance process because the squared kernel score is Lipschitz and bounded.This yields a uniform bound for the standardized variance fluctuation over the local parameter set.

D Theoretical properties of one-step conquer

The one-step conquer estimator uses higher-order convolution smoothing to reduce bias and support asymptotic normality under weaker dimensional growth. Its theory establishes Hessian positivity, Bahadur-error control, and bandwidth conditions under smoothness assumptions.

  • Higher-order smoothing: The one-step estimator uses a kernel of order ν ≥ 4 and requires a sufficiently smooth conditional error density to reduce smoothing bias.The density must be (ν − 1)-times differentiable under the stated conditions.
  • Higher-order smoothing: Higher-order kernels reduce the population-score bias to O(b^(ν−1)), while the population Hessian remains positive definite for sufficiently small bandwidth.These properties support the one-step construction despite higher-order kernels having negative parts.
  • Local curvature: With a consistent initial estimate, the empirical smoothed-loss Hessian is positive definite with high probability in a local neighborhood of the target.The result follows from uniform convergence of the empirical Hessian to its population counterpart.
  • Bahadur representation: Theorem D.1 provides a Bahadur linearization bound when n ≳ p(log n)^2+t and (p+t)/n ≲ b ≲ {(p+t)/n}^(1/(2ν)).The theorem controls the one-step estimator under the stated sample-size and bandwidth restrictions.
  • Bahadur representation: Higher-order kernels permit larger bandwidths, reducing the variance and total Bahadur linearization error.This is the stated consequence of Theorem D.1.
  • Asymptotic normality: For ν = 4, asymptotic normality holds when p(log n)^(9/16) = o(n^(7/16)) using h ≍ {(p + log n)/n}^(2/5) and b = {(p + log n)/n}^(2/9).The result applies to n^(1/2)⟨a, β̂ − β⟩ for arbitrary a ∈ R^p.

E Additional simulation results

The appendix extends the paper’s simulation settings across additional quantile levels, noise distributions, estimation-error and runtime measures, bandwidth choices, and confidence-interval diagnostics.

  • Simulation extensions: Figure E.1 extends simulation results under models (5.1)–(5.3) to τ ∈ {0.1, 0.3, 0.5, 0.7} with N(0, 4) noise.Results are averaged over 500 simulations.
  • Runtime: Figures E.2 and E.4 extend elapsed-time comparisons across τ ∈ {0.1, 0.3, 0.5, 0.7} under N(0, 4) and t2 noise.The comparisons include standard QR, Horowitz’s smoothing, and conquer.
  • Estimation results: Figures E.3 and E.5 extend estimation-result comparisons across additional quantile levels under t2 noise and multivariate normal random design.Figure E.5 includes τ = 0.9 in addition to the lower quantile levels.
  • Bandwidth sensitivity: Figure E.7 evaluates conquer’s sensitivity to the Gaussian-kernel bandwidth for n = 2000 and p ∈ {25, 50, 75, 100, 125, 150}.It compares the default bandwidth, an oracle-informed bandwidth, and standard QR’s estimation error.
  • Inference diagnostics: Figures E.8–E.10 report empirical coverage, confidence-interval width, and elapsed time for six methods under τ = 0.5 or 0.9 and different error distributions.The appendix includes both N(0, 4) and t2 error settings.
Loading 2012.05187v2…