Finite-Sample Distribution Theory and Efficient Large-Scale Inference for Online Quantile Regression
Authors: Ziyang Wei, Jiaqi Li, Lan Wang, Wei Biao Wu
Organizations: Department of Statistics, University of Chicago · Department of Statistics, Rice University · Department of Management Science, University of Miami
This paper studies online quantile regression for large-scale and streaming data using Stochastic SubGradient Descent (SSGD) with constant learning rates. Classical offline inference for quantile regression is computationally and memory intensive. Existing works of online inference for quantile regression provide only asymptotic guarantees and typically require sub-exponential tail conditions for distribution theory. To bridge these gaps, we introduce new techniques to prove a quenched central limit theorem (CLT) and finite-sample Gaussian approximation for SSGD under a finite-moment assumption. We further show that Ruppert-Polyak averaging with a constant learning rate has a non-vanishing bias and fails to satisfy CLT centering at the population target. Hence we propose suffix averaging to address this issue and establish its finite-sample Gaussian approximation. Based on these results, we provide an efficient online inference method for quantile regression that avoids covariance estimation. Numerical experiments show that our method achieves desirable empirical coverage rates and competitive performance compared to other inference methods. We also apply our approach to U.S. wage data to demonstrate its practical effectiveness.
Figures & tables
Figure 1 : Convergence, inference speed, and bias of RP averaging. (a) compares convergence under constant and decaying learning rates with dimension d=100 . (b) shows the computational time of three inference methods with dimension d=100 . Parallel SSGD is the one proposed in this paper, which achieves a speedup of more than 100× over the two competing methods. (c) shows that the n -scaled bias of full RP averaging remains non-negligible across different constant learning rates, whereas suffix averaging substantially reduces this scaled bias. The 1/2 -suffix estimator averages the second half of the iterates, while the 4/5 -power-suffix estimator approximately averages the final n−n4/5 iterates.
Figure 2 : Asymptotic normality: the density plot and QQ plot of standardized estimation errors. (a) last-iterate SSGD. (b) suffix averaged SSGD.
Figure 3 : Estimation and inference wall-clock time vs. sample size and dimension. The time for the first two panels is plotted on a logarithmic scale on the y-axis. The SSGD algorithm and parallel inference method are orders of magnitude faster.
Figure 4 : Parallel inference and random scaling – empirical coverage and length of confidence intervals. The first two rows: (n,d,1−α)∈{(107,1500,95%),(3×106,500,95%)} and 1000 replications. The last two rows: (n,d,1−α)∈{(107,1500,99%),(3×106,500,99%)} and 10000 replications. The shadows represent the error bars with one standard deviation.
Algorithm 1 Efficient Online Inference for Large-Scale Quantile Regression
α=0.05,d=500
Coverage (Relative error)
Length (std)
k=7.5×105
Parallel Inference
0.9490 (0.020)
0.011 (0.00333)
Random Scaling
0.9850 (0.700)
0.023 (0.01515)
k=1.5×106
Parallel Inference
0.9470 (0.060)
0.006 (0.00208)
Random Scaling
0.9840 (0.680)
0.012 (0.00766)
k=3×106
Parallel Inference
0.9470 (0.060)
0.004 (0.00137)
Random Scaling
0.9740 (0.480)
0.007 (0.00388)
Appendix
Table 1 : Comparison of Parallel Inference and Random Scaling for α=0.05 , d=500 , and n=3×106 . The value of k represents the k -th iteration or k -th sample. Coverage relative error is defined by Dα in ( 106 ). CI lengths are reported as mean with standard deviation in parentheses.
α=0.01,d=500
Coverage (Relative error)
Length (std)
k=7.5×105
Parallel Inference
0.9894 (0.060)
0.017 (0.00535)
Random Scaling
0.9975 (0.750)
0.036 (0.02287)
k=1.5×106
Parallel Inference
0.9885 (0.150)
0.010 (0.00333)
Random Scaling
0.9968 (0.680)
0.019 (0.01145)
k=3×106
Parallel Inference
0.9904 (0.040)
0.007 (0.00222)
Random Scaling
0.9955 (0.550)
0.010 (0.00574)
Appendix
Table 2 : Comparison of Parallel Inference and Random Scaling for α=0.01 , d=500 , and n=3×106 . The value of k represents the k -th iteration or k -th sample. Coverage relative error is defined by Dα in ( 106 ). CI lengths are reported as mean with standard deviation in parentheses.
α=0.05,d=1500
Coverage (Relative error)
Length (std)
k=2.5×106
Parallel Inference
0.9420 (0.160)
0.006 (0.00209)
Random Scaling
0.9890 (0.780)
0.021 (0.01486)
k=5×106
Parallel Inference
0.9350 (0.300)
0.004 (0.00126)
Random Scaling
0.9880 (0.760)
0.011 (0.00736)
k=107
Parallel Inference
0.9510 (0.020)
0.003 (0.00086)
Random Scaling
0.9850 (0.700)
0.006 (0.00369)
Appendix
Table 3 : Comparison of Parallel Inference and Random Scaling for α=0.05 , d=1500 , and n=107 . The value of k represents the k -th iteration or k -th sample. Coverage relative error is defined by Dα in ( 106 ). CI lengths are reported as mean with standard deviation in parentheses.
α=0.01,d=1500
Coverage (Relative error)
Length (std)
k=2.5×106
Parallel Inference
0.9884 (0.160)
0.010 (0.00319)
Random Scaling
0.9988 (0.880)
0.033 (0.02242)
k=5×106
Parallel Inference
0.9886 (0.140)
0.006 (0.00196)
Random Scaling
0.9991 (0.910)
0.017 (0.01118)
k=107
Parallel Inference
0.9902 (0.020)
0.004 (0.00131)
Random Scaling
0.9971 (0.710)
0.009 (0.00554)
Appendix
Table 4 : Comparison of Parallel Inference and Random Scaling for α=0.01 , d=1500 , and n=107 . The value of k represents the k -th iteration or k -th sample. Coverage relative error is defined by Dα in ( 106 ). CI lengths are reported as mean with standard deviation in parentheses.
Method
n=104
n=3×104
n=105
n=5×105
Parallel Inference
0.0009
0.0022
0.0071
0.0339
Conquer Bootstrap
0.0474
0.1116
0.3303
1.4667
Standard Quantile Regression
0.3874
0.7623
1.5802
4.6028
Appendix
Table 5 : Average inference time in seconds for different sample sizes ( d=20 ).
This paper considers the estimation of quantiles via a smoothed version of the stochastic gradient descent (SGD) algorithm. By smoothing the score function with a bandwidth tied to the learning rate, we obtain estimates that are monotone in the quantile level at every iteration, while retaining the memory and computational efficiency required for streaming data. We establish non-asymptotic tail probability bounds for the smoothed estimate with and without Polyak-Ruppert averaging, which are sub-exponential with a multi-regime structure. For the averaged estimate we further derive a Bahadur representation that is uniform in the quantile level and across coordinates, and a resulting Gaussian approximation by the maximum of Brownian bridges, with the dimension p allowed to grow exponentially in the sample size. This yields simultaneous inference across coordinates and quantile levels. As an alternative that avoids estimating the sparsity function, we propose an online multiplier bootstrap that preserves monotonicity, runs in a single pass and is asymptotically valid. Extending the theory to a localized recursion, we obtain online nonparametric conditional quantile estimates with uniform bands over design points and quantile levels. Simulations confirm accurate finite-sample coverage, and we illustrate the method on conditional value-at-risk curves.
Likai Chen, Georg Keilbar, Wei Biao Wu
Department of Mathematics and Statistics · Washington University in St.Louis · St.Louis, MO, USA +6
Stochastic gradient descent (SGD) is a foundational algorithm for large-scale statistical learning and stochastic optimization. However, statistical inference based on SGD iterates remains challenging when stochastic gradients have infinite variance, as the relevant limiting distributions depend on unknown nuisance parameters. In this paper, we develop an efficient, model-agnostic methodology for constructing confidence regions from SGD trajectories that applies in both finite- and infinite-variance regimes. The procedure is based on a joint weak convergence result for the Polyak-Ruppert averaged estimator and an empirical second-moment normalizer constructed from stochastic gradients along the SGD trajectory. This joint limit yields a self-normalized statistic in which the leading tail-dependent scaling terms cancel. We then use a subsampling calibration scheme to estimate the relevant critical values, avoiding explicit estimation of tail indices, slowly varying functions, or stable-law parameters. The resulting confidence regions are straightforward to implement and are asymptotically valid under both the finite- and infinite-second-moment regimes. Simulation studies show reliable coverage in various settings, supporting the proposed method as a practical tool for uncertainty quantification in stochastic optimization.
Jose Blanchet, Peter Glynn, Wenhao Yang
Management Science and Engineering, Stanford University
Online high-dimensional regression requires algorithms that can update sequentially while preserving structural sparsity. We propose \textit{Adaptive Iterative Hard Thresholding (AIHT)}, an online sparse-regression framework that alternates stochastic subgradient updates with adaptively scheduled hard-thresholding steps. The key idea is to separate support discovery from local refinement: early in the learning process, AIHT delays thresholding so that weak but informative coordinates have time to accumulate signal, while later it increases the projection frequency to stabilize the sparse estimator and exploit local curvature. We develop the theory for high-dimensional online quantile regression, a challenging setting in which the loss is nonsmooth and the data may exhibit heterogeneity or heavy-tailed noise. Under restricted curvature and gradient-leakage conditions, AIHT remains in an inflated sparse cone, exhibits a two-phase convergence behavior, and attains logarithmic regret for the sliding-window objective. Simulations for online quantile regression, together with threshold-scheduling ablations, support the proposed mechanism and illustrate its advantage over standard online sparse-learning baselines.