Understanding the Under-Coverage Bias in Uncertainty Estimation

Abstract:
Estimating the data uncertainty in regression tasks is often done by learning a quantile function or a prediction interval of the true label conditioned on the input. It is frequently observed that quantile regression—a vanilla algorithm for learning quantiles with asymptotic guarantees—tends to under-cover than the desired coverage level in reality. While various fixes have been proposed, a more fundamental understanding of why this under-coverage bias happens in the first place remains elusive. In this paper, we present a rigorous theoretical study on the coverage of uncertainty estimation algorithms in learning quantiles. We prove that quantile regression suffers from an inherent under-coverage bias, in a vanilla setting where we learn a realizable linear quantile function and there is more data than parameters. More quantitatively, for α > 0.5 and small d/n, the α-quantile learned by quantile regression roughly achieves coverage α− (α− 1/2) · d/n regardless of the noise distribution, where d is the input dimension and n is the number of training data. Our theory reveals that this under-coverage bias stems from a certain highdimensional parameter estimation error that is not implied by existing theories on quantile regression. Experiments on simulated and real data verify our theory and further illustrate the effect of various factors such as sample size and model capacity on the under-coverage bias in more practical setups.

heading: 1 Introduction
text: This paper is concerned with the problem of uncertainty estimation in regression problems. Uncertainty estimation is an increasingly important task in modern machine learning applications—Models should not only make high-accuracy predictions, but also have a sense of how much the true label may deviate from the prediction. This capability is crucial for deploying machine learning in the real world, in particular in risk-sensitive domains such as medical AI [15, 29], self-driving cars [47], and so on. A common approach for uncertainty estimation in regression is to learn a quantile function or a prediction interval of the true label conditioned on the input, which provides useful distributional information about the label. Such learned quantiles are typically evaluated by their coverage, i.e., probability that it covers the true label on a new test example. For example, a learned 90% upper quantile function should be an actual upper bound of the true label at least 90% of the time.
Algorithms for learning quantiles date back to the classical quantile regression [35], which estimates the quantile function by solving an empirical risk minimization problem with a suitable loss function that depends on the desired quantile level α. Quantile regression is conceptually simple, and is
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
theoretically shown to achieve asymptotically correct coverage as the sample size goes to infinity [34] or approximately correct coverage in finite samples under specific modeling assumptions [46, 60, 56]. However, it is observed that quantile regression often under-covers than the desired coverage level in practice [53]. Various alternative approaches for constructing quantiles and confidence intervals are proposed in more recent work, for example by aggregating multiple predictions using Bayesian neural networks or ensembles [24, 37], or by building on the conformal prediction technique to construct prediction intervals with finite-sample coverage guarantees [68, 66, 39, 53]. However, despite these advances, a more fundamental understanding on why vanilla quantile regression exhibits this under-coverage bias is still lacking.
This paper revisits quantile regression and presents a first precise theoretical study on its coverage, in a new regime where the number of samples n is proportional to the dimension d, and the ratio d/n is small (so that the problem is under-parametrized). Our main result shows that quantile regression exhibits an inherent under-cover bias under this regime, even in the well-specified setting of learning a linear quantile function when the true data distribution follows a Gaussian linear model. To the best of our knowledge, this is the first rigorous theoretical justification of the under-coverage bias. Our main contributions are summarized as follows.
• We prove that linear quantile regression exhibits an inherent under-coverage bias in the wellspecified setting where the data is generated from a Gaussian linear model, and the number of samples n is proportional to the feature dimension d with a small d/n (Section 3). More quantitatively, quantile regression at nominal level α ∈ (0.5, 1) roughly achieves coverage α− (α− 1/2)d/n regardless of the noise distribution. To the best of our knowledge, this is the first rigorous characterization of the under-coverage bias in quantile regression.
• Towards understanding the source of this under-coverage bias, we disentangle the effect of estimating the bias and estimating the linear coefficient on the coverage of the learned linear quantile (Section 4). We show that the estimation error in the bias can have either an undercoverage or over-coverage effect, depending on the noise distribution. In contrast, the estimation error in the linear coefficient always drives the learned quantile to under-cover, and we show this effect is present even on broader classes of data distributions beyond the Gaussian linear model.
• We perform experiments on simulated and real data to test our theory (Section 5). Our simulations show that the coverage of quantile regression in Gaussian linear models agrees well with our precise theoretical formula as well as the α− (α− 1/2)d/n approximation. On real data, we find quantile regression using high-capacity models (such as neural networks) exhibits severe under-coverage biases, while linear quantile regression can also have a mild but non-negligible amount of under-coverage, even after we remove the potential effect of model misspecification.
• On the technical end, our analysis builds on recent understandings of empirical risk minimization problems in the high-dimensional proportional limit with a small d/n, and develops new techniques such as a novel concentration argument to deal with an additional learnable variable in learning linear models with biases, which we believe could be of further interest (Section 6).

heading: 1.1 Related work
text: Algorithms for uncertainty estimation in regression The earliest methods for uncertainty estimation in regression adopted subsampling methods (bootstrap) or leave-one-out methods (Jackknife) for assessing or calibrating prediction uncertainty [52, 63, 58, 25]. More recently, a growing line of work builds on the idea of conformal prediction [55] to design uncertainty estimation algorithms for regression. These algorithms provide confidence bounds or prediction intervals by post-processing any predictor, and can achieve distribution-free finite-sample marginal coverage guarantees utilizing exchangeability of the data [50, 66, 39, 53, 33, 67, 69, 68, 11]. Further modifications of the conformal prediction technique can yield stronger guarantees such as group coverage [10] or coverage under distribution shift [9] under additional assumptions. Our under-coverage results advocate the necessity of such post-processing techniques, and are complementary in the sense that we provide understandings on the more vanilla quantile regression algorithm. Quantiles and prediction intervals can also be obtained by aggregating multiple predictors, such as using Bayesian neural networks [41, 24, 32, 44, 42] or ensembles [37, 49, 28, 45]. These methods offer an alternative approach for uncertainty estimation, but do not typically come with coverage guarantees.
Theoretical analysis of quantile regression Linear quantile regression with the pinball loss dates back to the late 1970s [34]. The same work proved the asymptotic normality of the regression coefficients in the n → ∞, fixed d limit. Takeuchi et al. [60] studied non-parametric quantile regression using kernel methods, and provided generalization bounds (with the pinball loss) based on the Rademacher complexity. Meinshausen [46] studied non-parametric quantile regression using random forest and showed its consistency under proper assumptions. Christmann and Steinwart [19], Steinwart et al. [56] established a “self-calibration” inequality for the quantile loss, which, when combined with standard generalization bounds, can be translated to an estimation error bound for quantile regression. These works all focus on bounding the parameter or function estimation error, which can be translated to bounds on the coverage bias, but does not tell the sign of this coverage bias as we do in this paper. We also remark that conformalization can be used in conjunction with quantile regression to correct its coverage bias [53].
Uncertainty quantification for classification For classification problems, two main types of uncertainty quantification methods have been considered: outputting discrete prediction sets with guarantees of covering the true (discrete) label [70, 71, 38, 7, 12, 18, 17], or calibrating the predicted probabilities [51, 72, 73, 37, 26]. The connection between prediction sets and calibration was discussed in [27]. The sample complexity of calibration has been studied in a number of theoretical works [36, 27, 54, 30, 40, 8]. Our work is inspired by the recent work of Bai et al. [8], which showed that logistic regression is over-confident even if the model is correctly specified and the sample size is larger than the dimension.
High-dimensional behaviors of empirical risk minimization There is a rapidly growing literature on limiting characterizations of convex optimization-based estimators in the n ∝ d regime [21, 13, 23, 31, 57, 61, 20, 62, 43, 59, 16]. Our analysis builds on results for unregularized Mestimator derived in [62] and generalizes theirs in certain aspects (see also [23, 20, 31]).

heading: 2 Preliminaries
text: In this paper we focus on the problem of learning quantiles. Suppose we observe a training dataset {(xi, yi)}ni=1 drawn i.i.d. from some joint distribution P on Rd × R, where xi ∈ Rd is the input features and y ∈ R is the real-valued response (label). Let F (t|x) := P(Y ≤ t|X = x) denote the conditional CDF of Y |X. Our goal is to learn the α-(conditional) quantile of Y |X:
q?α(x) := inf {t ∈ R : F (t|x) ≥ α}. For example, q?0.95(x) is the ground truth 95% quantile of the true conditional distribution Y |X, and can be seen as the “ideal” 95% upper confidence bound for the label y given the features x. Throughout this paper we work with upper quantiles, that is, α ∈ (0.5, 1) (some typical choices are α ∈ {0.8, 0.9, 0.95}); by symmetry our results hold for learning lower quantiles as well.
Coverage For any learned quantile function f̂ : Rd → R, the marginal coverage (henceforth “coverage”) of f̂ is the probability of y ≤ f̂(x) on a new test example (x, y):
Coverage(f̂) := P(x,y) ( y ≤ f̂(x) ) = Ex [ P ( y ≤ f̂(x)|x )] . (1)
For learning the α-quantile (α > 0.5), we usually expect Coverage(f̂) ≈ α, i.e. f̂(x) covers the label y on approximately α proportion of the data, under the ground truth data distribution.
We say that f̂ has under-coverage if Coverage(f̂) < α and over-coverage if Coverage(f̂) > α. Note that these two notions are not symmetric: Over-coverage means that the learned upper quantile f̂(x) is overly conservative (higher than enough), and is typically tolerable; In contrast, under-coverage means that f̂(x) fails to cover y with α probability, and is typically considered as a failure. We remark that while there exist more fine-grained notions of coverage such as conditional coverage [10], the (marginal) coverage is still a basic requirement for any quantile learning algorithm.
Quantile regression We consider quantile regression, a standard method for learning quantiles from data [35]. Quantile regression estimates the true quantile function qα(·) via the pinball loss [34, 56]
`α(t) = −(1− α)t1 {t ≤ 0}+ αt1 {t > 0} . (2)
Note that in the special case of α = 0.5, we have `0.5(t) = |t|/2, and thus the pinball loss strictly generalizes the absolute loss (for learning medians) to learning any quantile. Given the training dataset and any function class {fθ : θ ∈ Θ} (e.g. linear models or neural networks), quantile regression solves the (unregularized) empirical risk minimization (ERM) problem
θ̂ = arg min θ∈Θ
R̂n(θ) := 1
n n∑ i=1 `α(yi − fθ(xi)). (3)
(We take θ̂ as any minimizer of R̂n when the minimizer is non-unique.) Let R(θ) := E[R̂n(θ)] denote the corresponding population risk. It is known that the population risk over all (measurable) functions is minimized at the true quantile q?α = arg minf R(f) under minimal regularity conditions (for completeness, we provide a proof in Appendix B.1).

heading: 3 Quantile regression exhibits under-coverage
text: We analyze quantile regression in the vanilla setting where the input distribution is a standard Gaussian and y follows a linear model of x:
y = w>? x + z, where x ∼ N(0, Id), z ∼ Pz. (4)
Above, w? ∈ Rd is the ground truth coefficient vector, and the noise z ∼ Pz is independent of x. The Gaussian input assumption is required only for technical convenience in the high-dimensional limiting analysis, and we believe it is not strictly required for the same result to hold1 (an extension to more general input distributions can also be found in Section 4). The noise distribution Pz is required to satisfy the following smoothness assumption, but can otherwise be arbitrary: Assumption A (Smooth density). The noise distribution Pz has a smooth density φz ∈ C∞(R) (with corresponding CDF Φz), with bounded derivatives: supt∈R |φ (k) z (t)| < ∞ for any k ≥ 0. We further assume that φz(zα) > 0, where α ∈ (0.5, 1) is our pre-specified quantile level, and zα := inf {t ∈ R : Φz(t) ≥ α} is the α-quantile of Pz .
Under the above model, it is straightforward to see that the true α-conditional quantile of y|x is also a linear model (with bias):
q?α(x) = w > ? x + zα. (5)
Given the training data {(xi, yi)}ni=1, we learn a linear quantile function f̂(x) = ŵ>x + b̂ via quantile regression:
(ŵ, b̂) = arg min w,b
R̂n(w, b) := 1
n n∑ i=1 `α(yi − (w>xi + b)), (6)
where `α is the pinball loss in (2). As our linear function class realizes the true quantile function (5), the population risk is minimized at the true quantile: arg minw,bR(w, b) = (w?, zα).
We are now ready to state our main result, which shows that quantile regression exhibits an inherent under-coverage bias even in this vanilla realizable setting. Theorem 1 (Quantile regression exhibits under-coverage bias). Suppose the data is generated from the linear model (4) and the noise satisfies Assumption A. Let f̂(x) = ŵ>x + b̂ be the output of quantile regression (6) at level α ∈ (0.5, 1). Then, in the limit of n, d → ∞ and d/n → κ where κ ∈ (0, κ0] for some small κ0 > 0, for the coverage (1), we have (
p→ denotes convergence in probability)
Coverage(f̂) p→ α− Cα,κ for some Cα,κ > 0.
That is, the limiting coverage of the learned quantile function is less than α. Further, for small enough κ we have the local linear expansion
Cα,κ = (α− 1/2)κ+ o(κ). (7) 1Our results can be extended directly to any correlated Gaussian input x ∼ N(0,Σ) by the transform x̃ = Σ−1/2x and w̃? = Σ1/2w?. We believe our results also hold for i.i.d. sub-Gaussian inputs by the universality principle (e.g. [14]).
Theorem 1 builds on the precise characterization of ERM problems in the high-dimensional proportional limit [62], along with new techniques over existing work for dealing with the unique challenges in quantile regression (such as analyzing the additional learnable bias b). An overview of the main technical steps is provided in Section 6, and the full proof is deferred to Appendix C.
Implications Theorem 1 can be illustrated by the following numeric example. Suppose we perform quantile regression at α = 0.9, where the data follows the linear model (4), and our κ = d/n = 0.1 (so that the sample size is 10x number of parameters). Then Theorem 1 shows that, even in this realizable, under-parametrized setting, the coverage of the learned quantile f̂ is going to be roughly 0.9− Cα,κ when n, d are large, and further Cα,κ ≈ (α− 1/2)κ = 0.04. Thus the actual coverage is around 0.9− 0.04 = 0.86, and such a 4% under-coverage bias can be rather non-negligible in reality. We remark that a symmetric conclusion of Theorem 1 also holds for lower quantiles, and we further expect similar results also hold for learning prediction intervals, where the coverage is defined as the two-sided coverage of the prediction interval formed by the learned {lower quantile, upper quantile}. To the best of our knowledge, this offers a first precise theoretical understanding of why practically trained quantiles or prediction intervals often under-cover than the desired coverage level [53].
Comparison against existing theories An important feature of the under-coverage bias shown in Theorem 1 is that it only shows up in the n, d proportional regime, and is not implied by existing theories on quantile regression. Classical asymptotic theory only shows asymptotic normality √ n([ŵ, b̂] − [w?, zα]) → N(0,V) in the n → ∞, fixed d limit [34, 64]. Under this limit, Coverage(f̂) is consistent at α with O(1/ √ n) deviation. Christmann and Steinwart [19], Steinwart et al. [56] consider the finite n, d setting and establish self-calibration inequalities (similar to strong
convexity) that bounds the quantile estimation error by the square root excess loss √ R(f̂)−R(q?α). Combined with standard generalization theories (e.g. via Rademacher complexities) and Lipschitzness, this can be turned into a bound on |Coverage(f̂) − α|, but does not tell the sign (positive or negative) of the coverage bias.
Large κ; extension to over-parametrized learning While Theorem 1 requires a small κ = d/n, the approximation formula (7) suggests that the under-coverage should get more severe as κ—the measure of over-parametrization in this problem—gets larger. We confirm this trend experimentally in our simulations in Section 5.1.
As an extension to Theorem 1, we also show theoretically that the under-coverage bias indeed becomes even more severe in over-parametrized learning, under the same linear model (4): When d/n > Õ(1), and the noise Pz is sub-Gaussian and symmetrically distributed about 0, the convergence point of the gradient descent path on the quantile regression risk R̂n is the minimum-norm interpolator of the data, which has coverage 0.5± Õ(1/ √ d) with high probability (see Appendix D for the formal statement and the proof). Notably, this 0.5 coverage does not depend on α and exhibits a severe under-coverage.

heading: 4 Understanding the source of the under-coverage bias
text: In this section, we take steps towards a deeper understanding of how the under-coverage bias shown in Theorem 1 happens. Recall that the quantile regression returns f̂(x) = ŵ>x+ b̂ where (ŵ, b̂) is a solution to the ERM problem (6) and estimates the true parameters (w?, zα). Our main approach in this section is to disentangle the effect of the two sources—the estimation error in b̂ and the estimation error in w—on the coverage of f̂ .
We show that the estimation error in b̂ can have either an under-coverage or an over-coverage effect, depending on the noise distribution (Section 4.1). In contrast, the estimation error in ŵ always has an under-coverage effect; this holds not only for the linear model assumed in Theorem 1, but also on more general data distributions (Section 4.2). In the setting of Theorem 1, this under-coverage effect of ŵ is always strong enough to dominate the effect of b̂, leading to the overall under-coverage.

heading: 4.1 Effect of estimation error in b̂
text: To study the effect of b̂, we use the quantity b̂ − zα as a measure for its effect on the coverage— Recall that the true quantile is q?α(x) = w > ? x + zα, thus having b̂ < zα means that b̂ contributes to under-coverage, whereas b̂ > zα means b̂ contributes to over-coverage. (This can be seen more straightforwardly in the easier case where we know w? and only output b̂ to estimate zα.)
The following corollary shows that, under the same settings of Theorem 1, the error b̂− zα can be understood precisely. The proof can be found in Appendix E.1.
Corollary 2 (Effect of b̂ on coverage depends on noise distribution). Under the same settings as Theorem 1, for any α ∈ (0.5, 1), as n, d→∞ with d/n→ κ ∈ (0, κ0], we have
(a) The learned bias b̂ from quantile regression (6) converges to the following limit:
b̂− zα p→ Cbα,κ = b̄0κ+ o(κ),
where b̄0 has a closed-form expression:
b̄0 := −α(1− α)φ′z(zα)− (2α− 1)φ2z(zα)
2φ3z(zα) . (8)
(b) For any α ∈ (0.5, 1), when Pz is the Gaussian distribution (with arbitrary scale), we have b̄0 < 0 in which case Cbα,κ < 0 for small enough κ. Conversely, for any α ∈ (0.5, 1), there exists some noise distribution Pz for which b̄0 > 0, in which case Cbα,κ > 0 for small enough κ.
Corollary 2 shows that the sign ofCbα,κ in the limiting regime (and thus the effect of b̂ on the coverage) depends on b̄0, which in turn depends on the noise distribution Pz . For common noise distributions such as Gaussian we have Cbα,κ < 0 at small κ, but there also exists Pz such that C b α,κ > 0. Note that the second claim in part (b) follows directly from (8): we can always design the density φz by varying φz(zα) and φ′z(zα) so that b̄0 > 0. Overall, this result shows that the under-coverage bias in Theorem 1 cannot be simply explained by the under-estimation error in b̂.

heading: 4.2 Effect of estimation error in ŵ; relaxed data distributions
text: We now show that the primary source of the under-coverage is the estimation error in ŵ, which happens not only on the linear data distribution assumed in Theorem 1, but also on a broader class of data distributions. We consider the following relaxed data distribution assumption
y = µ?(x) + σ?(x)z, (9)
where the noise z ∼ Pz . We do not put structural assumptions on (µ?, σ?), except that we assume the true α-quantile is still a linear function of x, that is, there exists (w?, b?) for which
q?α(x) = µ?(x) + σ?(x)zα = w > ? x + b?. (10)
Since here we are interested in the effect of estimating w?, for simplicity, we assume that we know b? and only estimate w? via some estimator ŵ. We now collect our assumptions and state the result.
Assumption B (Relaxed data distribution). The data is distributed as model (9) with a linear αquantile function (10). Further, the data distribution satisfies the following regularity conditions:
(a) The distribution of x ∈ Rd is symmetric about 0, has a lower bounded covariance E[xx>] γId, and is K-sub-Gaussian, for constants γ,K > 0.
(b) The variance function σ?(·) is bounded and symmetric: For all x ∈ Rd we have σ ≤ σ?(x) ≤ σ for some constants σ, σ > 0, and σ?(x) = σ?(−x).
(c) The noise density φz is continuously differentiable and symmetric about 0, i.e. φz(t) = φz(−t) for all t ∈ R. Further, φz is uni-modal, i.e. φ′z(t)|t<0 > 0 and φ′z(t)|t>0 < 0.
Theorem 3 (Estimation error in ŵ leads to under-coverage on a family of data distributions). Under the relaxed data distribution assumption (Assumption B), for any α > 3/4, there exists constants c, r0 > 0 such that for any learned quantile estimate f̂(x) = ŵ>x + b? with small estimation error ‖ŵ −w?‖2 ≤ r0, we have
Coverage(f̂) ≤ α− cγ/σ2 · ‖ŵ −w?‖22 ,
that is, the learned quantile under-covers by at least Ω(‖ŵ −w?‖22). Above, c > 0 is an absolute constant, and r0 > 0 depends on (γ, σ,K,Φz, α) but not (n, d).
Implications; proof intuition Theorem 3 shows that, for a broad class of data distributions, any estimator ŵ will under-cover by at least Ω(‖ŵ −w?‖22). If particular, any estimator satisfying ‖ŵ −w?‖2 Õ( √ d/n) (e.g. from standard generalization theory) will under-cover by Õ(d/n). This confirms that the estimation error in the (bulk) regression coefficient ŵ is the primary source of the under-coverage bias, under assumptions that are more general than Theorem 1 in certain aspects (such as the distribution of x and y|x). We remark that as opposed to Theorem 1, Therorem 3 does not give an end-to-end characterization of any specific algorithm, but assumes we have an estimator ŵ with a small error.
At a high-level, Theorem 3 follows from the fact that any estimator ŵ>x+ b? must be lower than the true quantile w>? x + b? for some x and higher for some other x. Averaging the coverage indicator over x, such under-coverage and over-coverage cancel out on the first-order if x has a symmetric distribution, but aggregate to yield a under-coverage effect on the second-order as long as Φz(t)|t>0 is concave (which holds if Pz is unimodal). The proof of Theorem 3 can be found in Appendix E.2.

heading: 5 Experiments

heading: 5.1 Simulations
text: Setup We first test our Theorem 1 via simulations. We generate data from the linear model (4) in d = 100 dimensions with ‖w?‖2 = 1 and noise distribution Pz = N(0, 0.25). We vary κ = d/n ∈ {0.02, 0.04, . . . , 0.5} where κ determines a sample size n, and vary α ∈ {0.5, 0.52, ..., 0.98}. For each combination of (α, κ), we generate 8 random problem instances, and solve the quantile regression ERM problem (6) on each instance via (sub)-gradient descent. We evaluate the coverage of the learned quantile f̂ (thanks to the linear model (4), the coverage can be computed exactly without needing to introduce a test set). Additional details about the setup can be found in Appendix F.1.
Results Figure 1 plots the coverage of the learned quantiles. Observe that quantile regression exhibits under-coverage consistently across different values of (α, κ). Figure 1a shows that at fixed κ, the amount of under-coverage gets more severe at a higher α, which is qualitatively consistent with our approximation formula (α− 1/2)κ. Figure 1b further compares simulations with our analytical formula α − Cα,κ (found numerically through solving the system of equations 12), as well as the
local linear approximation α− (α− 1/2)κ claimed in Theorem 1. Note that the simulations agree extremely well with the analytical formula. The approximation α− (α− 1/2)κ is also very accurate for almost all κ at α = 0.8, and accurate for small κ at α = 0.9, 0.95. These verify our Theorem 1 and suggests it holds at rather realistic values of the dimension (d = 100).

heading: 5.2 Real data experiments
text: Datasets and models We take six real-world regression datasets: community and crimes (Community) [2], bike sharing (Bike) [1], Tennessee’s student teacher achievement ratio (STAR) [6], as well as the medical expenditure survey number 19 (MEPS_19) [3], number 20 (MEPS_20) [4], and number 21 (MEPS_21) [5]. All datasets are pre-processed to have standarized features and randomly split into a 80% train set and 20% test set.
To go beyond linear quantile functions, we perform quantile regression with one of the following four models as our fθ: linear model (Linear), a 3-layer MLP (two non-linear layers) with width 64 (MLP-3-64), 512 (MLP-3-512), and a variant of the width-512 MLP where all representation layers are frozen and only the last linear layer is trained (MLP-freeze-3-512). All linear layers include a trainable bias. We minimize the α-quantile loss (3) via momentum SGD with batch size 64. For each setting, we average over 8 random seeds where each seed determines the train-validation split, model initialization, and SGD batching. In our real experiments we fix α = 0.9. (Results at α ∈ {0.8, 0.95} as well as additional experimental setups can be found in Appendix F.2).
Results Table 1 reports the coverage of the learned quantile functions (evaluated on the test sets). Observe that all MLPs exhibit under-coverage compared with the nominal level 90%. Additionally, the amount of under-coverage correlates well with model capacity—the two vanilla MLPs undercovers more severely than the MLP-freeze and the linear model. Notice that the linear model does not have a notable under-coverage on most datasets—we believe this is a consequence of d/n being small on these datasets. The only exception is the Community dataset with the highest d/n ≈ 1/16, on which the linear model does under-cover mildly by roughly 1%.

heading: 5.3 Linear quantile regression on pseudo-labels
text: To further test the coverage of linear quantile regression on real data distributions, we make two modifications: (1) We subset the training data by fixing d and reducing n, so as to test the coverage across differerent values of κ = d/n; (2) We compare linear quantile regression on both true labels yi, and pseudo-labels y pseudo i generated from estimated linear models. These pseudo-labels are generated by first fitting a linear model ŵ ∈ Rd (with square loss) on the training data, and then generating a new label using the fitted linear model ŵ:
ypseudoi = ŵ >xi + σ̂zi,
where σ̂ is estimated as √
Ê(x,y)[(y − ŵ>x)2] on a separate hold-out split, and zi ∼ N(0, 1). The motivation for the pseudo-labels is to make sure that the data comes from a true linear model, removing the potential effect of model misspecification.
Table 2 shows that on the MEPS_20 dataset, linear quantile regression exhibits under-coverage at relatively large values of κ (0.1, 0.2, 0.5) for both kinds of labels. Also, there is no notable difference
between pseudo-labels and true labels. This provides evidence that our theory on linear quantile regression may hold broadly on real-world data distributions.

heading: 6 Proof overview of Theorem 1
text: Closed-form expression for coverage Our first step is to obtain a closed-form expression for the coverage. Recall that
Coverage(f̂) := P(x,y)(y ≤ f̂(x)) = P(x,z)(〈w?,x〉+ z ≤ 〈ŵ,x〉+ b̂).
As x is standard Gaussian, and the random variable z has cumulative distribution function Φz , standard calculation then yields the closed form expression (Lemma B.1)
Coverage(f̂) = EG∼N(0,1)[Φz(‖ŵ −w?‖2G+ b̂)].
Concentration of ‖ŵ − w?‖2 and b̂ We generalize results from recent advances in highdimensional M-estimator in linear models [23, 20, 31, 62] to show that ‖ŵ −w?‖2 and b̂ obtained by quantile regression 6 concentrates around fixed values in the high-dimensional limit. We show that, in the limit of d, n→∞ and d/n→ κ, the following concentration happens:
‖ŵ −w?‖2 p→ τ?(κ), and b̂ p→ b?(κ). (11)
Above, τ? and b? are determined by the solutions of a system of nonlinear equations with three variables (τ, λ, b):  τ2κ = λ2 · E(G,Z)∼N(0,1)×Pz [e ′ `αb (τG+ Z;λ)2], τκ = λ · E(G,Z)∼N(0,1)×Pz [e ′ `αb (τG+ Z;λ)G],
0 = E(G,Z)∼N(0,1)×Pz [e ′ `αb (τG+ Z;λ)],
(12)
where e`(x; τ) = minv 12τ (x − v) 2 + `(v) and `αb = ` α(t − b) is the shifted pinball loss (2). (See Theorem C.1 for the formal statement.) This is established via two main steps: We first build on the results of Thrampoulidis et al. [62] to show that a variant of the risk minimization problem with a fixed bias b concentrates around the solution to the first two equations in (12). We then develop a novel concentration argument to deal with the additional learnable bias b in the minimization problem (6), which introduces the third equation in (12) that will be used in characterizing the limiting value of the minimizer b̂.
The concentration (11) implies that Coverage(f̂) also converges to the following limiting coverage value (Lemma C.3):
Coverage(f̂) p→ EG∼N(0,1)[Φz(τ?(κ)G+ b?(κ))] =: α− Cα,κ. (13)
Calculating the limiting coverage via local linear analysis In this final step, as another technical crux of the proof, we further evaluate the small κ approximation of coverage value (13), and determine the sign of Cα,κ. This is achieved by a local linear analysis on the solutions of the aforementioned system of equations at small κ (Lemma C.2) in a similar fashion as in [8], and a precise analysis on the interplay between the concentration values τ?, b?, and the noise density φz . Combining these calculations yields that Cα,κ/κ = −(α−1/2) +o(1) for small enough κ (Lemma C.4). As α > 1/2, this establishes Theorem 1. All details on these analyses can be found in our proofs in Appendix C.

heading: 7 Conclusion
text: This paper presents a first theoretical justification of the under-coverage bias in quantile regression. We prove that quantile regression suffers from an inherent under-coverage bias even in well-specified linear settings, and provide a precise quantitative characterization of the amount of the under-coverage bias on Gaussian linear models. Our theory further identifies the high-dimensional estimation error in the regression coefficient as the main source of this under-coverage bias, which holds more generally on a broad class of data distributions. We believe our work opens up several interesting directions for future work, such as analyzing non-linear quantile regression, as well as analyzing other notions of uncertainty in regression problems.
Funding transparency statement
YB, HW, CX are funded through employment with Salesforce.

Review summaries:
In this work, the authors consider linear quantile regression, paying particular attention to the "coverage" properties of the resulting predictor. The traditional goal of quantile regression is to learn a quantile of the conditional distribution of output $Y$, conditioned on the inputs $X$. Denoting the predictor returned by any quantile regression algorithm by $\hat{f}(\cdot)$, one can use this to construct a predictive interval for a freshly-drawn $Y$. The authors consider the one-sided interval $C(X) = (-\infty, \hat{f}(X)]$, and the event $\\{ Y \in C(X) \\}$. The authors call the probability of this event (over the draw of a pair $(X,Y)$) the coverage of $\hat{f}$. Their main results elucidate conditions under which the traditional pinball-loss based ERM quantile regression solution is biased below the desired quantile level, that is, the predictive interval tends to be "too small," due to $\hat{f}$ under-estimating the desired quantile.  In particular, they show that in the limit where the sample size and dimensionality are proportional, the coverage (the aforementioned probability, which depends on $\hat{f}$ and thus the training data) converges below the desired quantile, with this error depending directly on the ratio of dimension to sample size $d/n$. Furthermore, they show (under fixed $n$ and $d$) how for a certain family of data distributions and a linear quantile model, given any estimator of the weights, the downward bias term scales with the weight estimation error. Their empirical tests also effectively highlight the phenomena described formally.

This paper presents a theoretical study on the coverage of quantile regression algorithms. More specifically, the authors prove theoretically that linear quantile regression exhibits an inherent under-coverage bias, with an order d/n difference in the marginal coverage regardless of the noise distribution (but Gaussian for the covariates). They also show that the primary source of this under-coverage bias is the estimation error in the linear coefficients, while the estimation error in the bias term can have either an under- or over-coverage effect. Numerical experiments further support their findings.

Quantifying predictive uncertainty is a crucial task for high-stakes prediction problems. The paper focus on quantile regression, a common approach for quantifying predictive uncertainty in regression problems. Authors consider realizable linear setting when the true model belongs to the considered functional class. With a focus on the under-parameterized setting, they theoretically justify the under-coverage of the learnt quantile regressors. Established theoretical results are backed up with supporting empirical evidence.

This paper gives theoretical insights into the under-coverage phenomenon often observed in quantile regression, i.e., for a continuous distribution, the probability of observing  a true label below a learned  $\alpha$-quantile is often below $\alpha$, which is consider a failure in estimating the uncertainty. Two important results are shown. The first one provides a characterization of the under-coverage in the setting of a Gaussian linear model in the under-parameterized high-dimensional proportional limit (i.e. $d,n \rightarrow \infty$ , $d/n\rightarrow\kappa$ with small $\kappa$) independently of the noise distribution. The experimental results show that the formulas match pretty well the observed coverage in realistic situations of dimension around 100, on synthetic and real data. The second important result shows that the source of the under-coverage stems essentially from the error in estimating the high-dimensional linear coefficient, under some assumptions on the underlying source.