Expected Pinball Loss For Quantile Regression And Inverse CDF Estimation
Authors: Taman Narayan
Serena Lutong Wang
Kevin Robert Canini
Maya Gupta
Venue: Accepted by TMLR ()

Abstract: We analyze and improve a recent strategy to train a quantile regression model by minimizing an expected pinball loss over all quantiles.  Through an asymptotic convergence analysis, we show that minimizing the expected pinball loss can be more efficient at estimating single quantiles than training with the standard pinball loss for that quantile, an insight that generalizes the known deficiencies of the sample quantile in the unconditioned setting. Then, to guarantee a legitimate inverse CDF, we propose using flexible deep lattice networks with a monotonicity constraint on the quantile input to guarantee non-crossing quantiles, and show lattice models can be regularized to the same location-scale family. Our analysis and experiments on simulated and real datasets show that the proposed method produces state-of-the-art legitimate inverse CDF estimates that are likely to be as good or better for specific target quantiles.

Published in Transactions on Machine Learning Research (01/2024)
Expected Pinball Loss For Quantile Regression And Inverse CDF Estimation
Taman Narayan
tamann@google.com
Google Research
Screna Wang
serenawaangkggoogle.ccmm
Google Research
University of California, Berkeley
Kevin Canini
conéni@dgoogle.com
Google Research
Maya R. Gupta University of Washington
wlatknventropyy@mmall.com
Reviewed on OpenReview: https: // openreview. net/ forum? id Eg8RnbOHad
Abstract
We analyze and improve recent strategy train a quantile regression model by minimizing an expected pinball loss over all quantiles. Through an asymptotic convergence analysis, show that minimizing the expected pinball loss can be more efficient. at estimating single quantiles than training with the standard pinball loss for that quantile, an insight that generalizes the known deficiencies of the sample quantile in the unconditioned setting. Then, to guarantee a legitimate inverse CDF, we propose using flexible deep lattice networks with a monotonicity constraint on the quantile input to guarantee non-crossing quantiles, and show lattice models can be regularized the same location-scale family. Our analysis and experiments simulated and real datasets show that the proposed method produces state- of-the-art legitimate inverse CDF estimates that are likely to be good or better for specific target quantiles.
1 Introduction
Real world applications often soek estimates of the quantiles of a random variable. For example, an airline would like to tell passengers that 90% of the time their flight will arrive within hours, so they know they can make their connecting flights. In addition, one might condition that estimate on features. For example, conditioned on a forecast of inches of snow at the destination airport, the conditional estimate might worsen to 90% of the time their flight arriving within 6 hours and minutes.
Quantile regression learns model from training examples to produce such estimates for quantiles. If the model can estimate all quantiles, it is an inverse CDF estimator.
1.1 Formal Definitions
Formally, let random variable Y E IR have cumulative distribution function (CDF) F, and for E 10, the r-quantile of Y is defined as = F-1(T) inf(g : P(Y < g) T). In the conditional setting where one also has random feature vector X E RDD the conditional r-quantile of for feature vector X is defined as g. (z(z)= = inf(g P(Y < q|X = x) T).
Quontile regression takes as training data set of n pairs (z, y) E R from joint distribution over (X, Y), and forms estimator for one more of the conditional quantiles of Y for any value of X. A

---

Published in Transactions on Machine Learning Research (01/2024)
standard objective to train model to estimate the quantile for is to minimize the pinball loss (Koenker & Bassett 1978), Thr (v. 0) = max(r(y - 0), (r 1) 0)) for V. E IR. In the unconditioned case with no feature vector X, the training data is simply a set of n scalar values Hy and minimizing the pinball loss has the satisfying property that it selects the empirical quantile of the training set (Chernozhukov 2005
1.2 Expected Pinball Loss
Recent work by Tagasovska & Lopez-Paz 2019) fI,, T) that takes 7 as an imput, and is trained to minimize an expected pinball loss over random T, drawn from a uniform distribution. That work is important because it provides simultaneous quantile regression of the entire inverse CDF. However, it left open a couple of important theoretical and practical issues that we address in this paper.
proposed training one deep neural network (DNN) model
1.3 Is Expected Pinball Loss Worse At Specific Quantiles?
The first. issue we address is whether one pays a price in estimation accuracy for a model f(x; T) that can predict any quantile, compared to a model trained specifically for one target quantile. Surprisingly, we show theoretically with an analysis of asymptotic convergence properties that learning the full inverse CDF model (Ir,T) r) can produce more efficient. estimator for single quantile than minimizing the pinball loss for just single quantile, depending
the true distribution, quantile, and function class. Our simulations and real-world experiments confirm that one does often win single quantile estimates by training with the expected pinball loss, though the full inverse CDF model f(x,T) may take a bit longer to train and a bit more memory to store.
We also derensstrate a novel use case of the expected pinball loss for single quantile regression - training with a Beta distribution over 7 and show how that can also outperform single-quantile regression (without the added bonus of accurately estimating the full inverse CDF).
1.4 The Problem Of Non-crossing Quantiles
The second issue we address is that training a DNN model with the expected pinball loss as proposed by Tagasovska & Lopez-Pazz Paz (2019) may fail to produce a legitimate inverse CDF because it does not guarantee non-crossing quantiles: that any two quantile estimates satisfy the common sense requirement that 9- (x())= 2 (z) for r' at every
For example, model that does not guarantee non-crossing quantiles might tell passenger there is a 90% chance their flight will arrive in 4 hours, but that there is a 80% chance their flight will arrive in 4 hours and 12 minutes. Such non-crossing quantile estimates have long been thought objectionable (Koenker & Bassett 1978).
To lest whether non-experts people would actually notice or mind non-crossing estimates, we emailed 100 frequent customers of Artifact Puzzles, who are estimated to be 70% college-educated and 20% with post- graduate education, and told them we were considering giving them estimates of how long it would take for their package to arrive, and gave them example arrival time estimates for the 50%, 80%, 90% quantiles, and asked for their feedback on how useful they would find such estimates. However, the example estimates were all made-up: we gave 50 of the customers estimates with crossing 80% and 90% quantiles (3 days, 8 days and 7 days), and the other 50 customers got estimates with non-crossing quantile predictions (3 days, 7 days, 8 days). From the 50 emails with non-crossing estimates, we received 9 emails back to say they would apprcciate such estimates, and no concerns. From the 50 emails with crossing estimates, we received 16 emails back all expressing enthusiasm for having such predictions, but 11 of the 16 emails pointed out the crossing quantiles with negative adjectives including 'wrong,' "broken," "buggy," 1 'didn't make sense to me."
Embarrassing quantile-crossing mistakes are known to happen often in quantile regression He 1997 Koenker & Bassett 1978); see also Table 1. Tagasovska & Paz (2019) hypothesized that training with an expected pinball loss induces smoothing across т that would reduce non-crossing. However, we show exper- imentally (Table 3) that a flexible model like a DNN optimized to minimize the expected pinball loss casily suffers a problematic amount of quantile crossing. To address non-crossing without losing model flexibility, we propose using deep lattice networks (DLNss) You et al. 2017] with a monotonicity shape constraint on
2

---

Published in Transactions on Machine Learning Research (01/2024)
to guarantee non-crossing quantiles, thus providing legitimate inverse CDF estimates. Additionally, we show that the DIN function class is amenable to two additional kinds of useful regularization for quantile regression: restricting the learned inverse CDF to location-scale family, and restricting other features to have only monotonic effect on the predictions.
1.5 Other Uncertainty Estimation Problems
In this paper, we focus estimating multiple quantiles and the extreme case of all quantiles: inverse CDF estimation. Such problems are only part of the landscape of estimating uncertainty. A closely related problem is estimating prediction intervals such that an a-interval contains the true label alpha is and sharp around a. Prediction intervals can be generated by optimizing for the calibration and sharpness of the interval (Pearce et al. 2018 Chung et al., 2021). Tagasovska & Lopez-Pas Paz (2019) instead formed them by choosing a pair of quantiles from multiple quantile regression model, which works well because the standard pinball lose used for training quantile regression models balances calibration and sharpness. A different type of uncertainty estimation problem is estimating the uncertainty of specific statistic, such the mean in a regression model; such problems are often handled with Bayesian methods. All of these estimation problems become more challenging in the out-of-domain (OOD) setting.
1.6 Paper Organization
Next, in Section 2 we formalize the problem of producing legitimate inverse CDF estimates with quantile regression. Then in Section 3 we dig into what is known about beating the standard pinball loss by smooth- ing over multiple quantile estimates, and give a asymptotic convergence result for minimizing the expected pinball loss. In Section 4 we propose using DLNs as the model architecture. Experiments in Sec- tion 5 on simulations and real-world data show that the proposed non-crossing DLNs provide competitive, trust-worthy estimates.
2 Estimating A Legitimate Inverse CDF
We give constrained optimization problem train quantile regression estimator without crossing quan- tiles, which can produce legitimate inverse CDF. Let (.),,,)) be a training set where each X2 E RD is feature vector and E R is corresponding label. Denote the model f(z, T:f), where E (0, 1) specifies the quantile of interest, and the model f :R RRD+1 - IR is parameterized by 0 E IR" Note that we sometimes stick to f(x,T) r) for notational simplicity.
Recall the standard pinball loss given is (y. () = max(r(y 0). 1) 9)). Let T~ Py denote random so can minimize the expected pinball loss over T to train f(x, F; as in Tagasovska & Lopez-Paz 2019), though here generalize to Pr, rather than assuming uniform distribution for T. The standard pinball loss can be written as the special case where Py is the Dirac delta distribution on the target quantile.
Also, we constrain the empirical loss minimization with non-crossing constraints, producing the following constrained training objective:
min E Σ Ly (ys, LT f(r,, T:0))
(1)
s.t. f(x, + 2 ( 25.
(2)
2.1 Related Work In Minimizing Sum Of Pinball Losses
The idea of training one model that can predict multiple quantiles by minimizing a sum of pinball losses with non-crossing constraints dates back to at least the 2006 work of Takeuchi et al. (2006) for a pre- specified, discrete set of quantiles. Other prior work used similar mechanism for a pre-specified discrete set of quantiles, and for more-restriclive function classes that are monotonic on v, such as linear models
3

---

Published in Transactions on Machine Learning Research (01/2024)
which are known to have limited flexibility (Daniels & Velikova 2010). Monotonic DNN models with more layers (Minin et al. 2010: Lang 2005) or min-max layers (Sill] 1998) can theoretically provide universal approximations. but have not
[Bondell et al. 2010). and two-layer monotonic neural networks (Cannon 2018
performed as well experimentally as DLNs You et al. 2017). Monotonic neural nets have also been proposed
for estimating a CDP (Chilinski vil 2018))
Training only for (I) without (2) Tagasovska & Lopez-Paz 2019: does not guarantee a legitimate inverse CDF because the estimated quantiles can cross. An in-between strategy is to instead add to [1 penalty violations of the monotonicity constraint training samples (and one can add virtual examples as well for more coverage) (Sill & Abu-Mostafa 1997), but applying that strategy to quantile regression cannot guarantee non-crossing quantiles everywhere.
2.2 Related Work In Estimating The Inverse CDF
There are two main prior approaches to estimating a complete legitimate inverse CDF. The first is nonpara- metric. For example, k-NN can be extended to predict quantiles by taking the quantiles rather than the mean from within a neighborhood (Bhattacharya & Gangopadhyay 1990). Similarly, quantile regression forests (Meinshausen 2006) use random forests leaf nodes to generate local empirical quantiles.
The second strategy is to predict parametric inverse CDF for any x. Traditionally these methods have been fairly rigid. He 1997) developed a method to fit a shared but learned location-scale family across x. Yan et al. (0018) proposed modified 4-parameter Gaussian whose skew and variance was dependent Recently, Gasthaus et al. 2201 proposed the spline quontile function DNN (SQF-DNN), which is DNN model that takes an input x, and outputs the paramcters for a monotonic piccewise-linca quantile function on т. SQF-DNN can fit any continuous bounded distribution, given sufficient output parameters, and because for any x the output is a monotonic function on т, guarantees no quantile crossing. Though Gasthaus et al. (2019) focused on RNNs, their framework is easy to apply to the generic quantile regression setting as well, which we do in our experiments.
3 Comparison To The Single Quantile Pinball Loss
We show through asymptotic analysis and simulations that minimizing the expected pinball loss can be better at estimating the r-quantile than minimizing the pinball loss just for т, and that these positive results are foreshadowed by classic results for unconditional quantile estimates. This section focuses largely on the unconditoneed case (no feature vector æ);; will return to the conditional problem given features in Section 4
3.1 Deficiency Of The Sample Quantile
Given only samples from Y, minimizing a pinball loss produces the corresponding sample quantile (Cher- nozhukov 2005). However, it has long been known that simply taking the sample quantile is not necessarily the best quantile estimator, and in particular that it can be beaten by strategies that smooth multiple sample
quantiles (Harrell & Davis 1982 Kaigh & Lachenbruch, 1982 David & Steinberg 1986 Reiss 1980 Falk
1984). For example, one important classic result that follows from the Lehnann-Schaffe theorem, says that if the data are known to be from a uniform (a, b] distribution, then the minimum variance unbiased estimator for the median not the sample median, rather it is the average of the sample min and sample max.
Even if we do not know anything about the true distribution, general-purpose methods that smooth multiple quantiles of the data when estimating particular quantile have proven effective. In fact, the widely used Harrell-Davis estimator is a weighted average of all the sample order statistics where the weights depend on T, and the estimator for the median is asymptotically equivalent to computing the mean of bootstrapped medians Harrell & Davis 1982)-
The prior work above suggests that there can be a better training loss than the r-specific pinball loss for the general problem of quantile regression given features one that increases efficiency by better smoothing across the quantiles. Indeed, by minimizing the expected pinball loss as in ( the estimate of any quantile
4

---

Published in Transactions on Machine Learning Research (01/2024)
is affected by the other quantiles (given some smoothness in f). We provide some further intuition about how it smooths in Appendix A.7
3.2 Asymptotic Properties
We give new asymptotic convergence result for estimating the full inverse CDF by minimizing the expected
pinball loss D optionally with monotonicity constraint (2), and compare the resulting asymptotic variance to that of single-quantile estimation learned by minimizing individual pinball loss.
We show that depending on the underlying data distribution and the function class chosen to estimate the inverse CDF, the simultancous estimation can produce a more efficient estimator for a single target quantile than minimizing the pinball loss for only that quantile.
For tractability, we work in the setting where there are no features æ condition on, only set of training samples {yi}, and we are interested in estimating the quantile F-) for a random variable Y E R given IID samples {y. Such analysis is relevant to the conditional case because for function classes that are flexible over X, they intuitively will have similar dynamics at each value of X as the unconditioned case we work with here. We formalize this in basic corollary in Section 3.4
We also focus our theory on the case where the quantile model is well-specified, i.e. the used function class fxx 7) contains the true data generating model. In practice, this often not the case, but we further explore in Section 3.5 how and why the use of (lunably) flexible function class can still improve quantile estimation in the general case. In particular, moves us from global smoothing across all of the data towards more local smoothing in the vicinity of the target quantiles.
explore how adding more flexibility the function class
General theory for each of these issues - conditional CDF estimation and misspecified function classes - remain open problems. Our main experiments in Section 5.5 though, do consider the general case of conditional quantile regression with unknown output distributions.
3.2.1 Review of Single Pinball Loss Asymptotic Properties
We start with the well-understood case of single-quantile estimation. Let 0 arg min, Ey (Y, (Y,0) be the true minimizer for the corresponding pinball loss and 0.²)) arg ming a/ Σn (T>((, 0)] be the empirical minimizer over our data. (In context, (T) the true quantile of Y and ofm) is the 7 sample quantile of the data (Chernozhukov 2005 .) The asymptotic distribution for this single quantile estimate is given by,
vū(((()) - 07) 4 N (o, -(1--))
p(07)²
where p(y) is the density for Y (Koenker 2005).
3.2.2 Expected Pinball Loss Asymptotic Properties
Next, consider learning the full inverse CDF f(T:0) parameterized by 0 E θ C R" by minimizing the expected pinball loss:
arg min
8E0
where T E (0, 1) is random variable drawn from some distribution independently of Y. Note that we do not separately analyze the monotonicity constraint in Equation (2) except to note that said constraints can be wrapped into the function class ff and feasible parameter space θ It may be that explicitly modeling such constraints would improve the asymptotic variance we derive; converscly, imposing the monotonicity constraint on certain function classes may lead violations of some of the sufficient conditions set forth in Lemma 1
Now let ((n) be the corresponding minimizer under empirical distribution over Y with samples (y In Theorem 1 below, we present general asymptotic distribution for g(x) that depends the properties of
5

---

Published in Transactions on Machine Learning Research (01/2024)
the parametric function f. feasible parameter space θ, and the distribution of T. A natural choice for the distribution of T is Unif(O, 1), but in fact, depending on the distribution of Y and properties of the function class 0(T;0), other distributions could lead to lower asymptotic variance. In choosing the distribution of T, note that choosing distributions with wider support will make it casier to satisfy the requirements for asymplotic normality in Theorem [ below (see note in Appendix A.4).
Theorem 1. Suppose f(T;0) is well specified with a unique 0" E C ORC Th such that [(T; F-)(T) for all E (0, 1), and that 0 is also the unique minimizer for the risk R(0) = (Y,,LL(Y,,,(TT 0))] Suppose the estimator Φ(n) is weakly consistent, g(n) Ly f Suppose further that the function ++ [(T; 0) is continuous, locally Lipschitz, and twice differentiable
Hr for all (0, Then,
vn(f(r;B(m) - f(T; 3 N (0, (of (7 -- vo[(T;
(4)
where Q -=-E(U((T:0))NT
with r(r) -
Lemma ! in general depends on properties of the function f and parameter space O
provides a simple example of sufficent conditions for the consistency condition in Theorem which
Lemma 1. If f(T, ()) is convex in f for all r E (0, 1): θ is convex, and 0" is in the interior of θ then ((x) is weakly consistent, ((m) 2 0".
Lemma applics when the optimization problem in [1 is convex. Aside from convexity, other example conditions could include compactness θ or further Lipschitz continuity or bounded moment conditions on f (Powell 2010).
Given a specific distribution of T, function class of the inverse CDF f. and data distribution Y, directly compare the asymplotic variance of the direct quantile estimates (Equation 3 to that of quantile estimates from the learned inverse CDF (Equation We illustrate this below with several examples.
3.2.3 Example 1: Single Parameter Uniform
To build intuition for the implications of Theorem 1 we consider centered uniform distribution Y Unif(- -.), whose inverse CDF F (r) - + Let f(T;0) = 0(T f). Expanding the asymptotic variance term in Theorem 1,
Q IV(--))+= a a²Erl-T' + 27 5T2 + +1
ET - -T+T 1
This depends on the first four moments of the distribution of T, the choice which is up to the algorithm designer. For the natural choice of T Unif(00 1), the asymptotic variance for estimating a specific To- quantile using the function f(ru; ((0)) = ((m) (To - }1) is 1.2o 1). This variance is lowest when To = 0.5, and increases to 0.3a² the most extreme quantiles when 7 = or r 1.
We next compare this to estimating the To-quantile using the single pinball loss Top Equation (3) shows that the estimate g(*) has asymptotic variance a²ro(1 - mp). This shows that depending on the exact value of To, the single quantile estimate could be less or more efficient than the estimate resulting from evaluating the function [T;; ((m)). The single quantile estimate more efficient for the extreme quantiles, whercas the inverse CDF estimate is more efficient for estimating the median. Intuitively, the inverse CDF median estimate benefits from fitting the inverse CDF function to the surrounding data points, which is reminiscent of classic quantile smoothing.
3.2.4 Example 2: Two-Parameter Uniform
We next consider a two-parameter uniform distribution example: for Y Unil(a, B), p(y) aーa and (-+(7) + (8 a)r. Let e IR with f(r;0) 0o + 0T.. That is, we learn two-parameter line as our inverse CDF. We then have Val(T;0)= 1 and can therefore write the Q and V matrices from Theorem

---

Published in Transactions on Machine Learning Research (01/2024)
as
T T T
T T-T² T
T T3 73
3-a ET
= ET
-2
Evaluating each of these matrices for T ~ Unif(O, 1) we arrive at
Q v'Q'"=") = (8 - a)² 15
and therefore an asymptotic variance for particular 70-quantile estimate f(T;; 0(1)) 0 01 of (B () + To (o)). This is quadratic function that finds its minimum at 0.5 and is highest (on the unit interval) at 0 and 1.
We can compare it to the asymptotic variance for single-pinball regression of (3 - a) T), which is a quadratic function that attains its marimum at r 0.5.
As we see in Figure 1 the variance is smaller for the expected pinball regression for quantiles between roughly 0.3 and 7 0.7. The opposite is true for the extreme quantiles.
We validate these results empirically through simulation in Table In each simulation, draw 1,000 samples of Y from Uniform(0,1) and seek to best estimate the truc 10th, 50th, and 90th quantiles (i.c. 0.1, 0.5, and 0.9). We train a linear model over T using the expected pinball loss (equivalent to a 2-keypoint 1-dimensional DLN as explained Section [ and compare it to taking the sample quantile or using the Hlarrell-avvis estimator. We find the expected pinball loss model performs best at estimating the median but worse at more extreme quantiles. Note that this is specific the uniform distribution simulation; we would not necessarily expect. the same outcome for any quantile estimation problem.
Asymptotic Convergence Rates
Figure 1; Plots of the asymptotic variances for the expected pinball regression (orange) and individual single-pinball re- gressions (blue) for different quantiles, when assume that the true distribution is Uniform(a, 3), the expected pin- ball loss regression learns linear model, and we assume B - a = 1. Relaxing the last assumption would scale the y-axis but not change the intersection points.
0.35
0.30
0.15
0.10
Individual OR Simultaneous OR
01
02
03
04
0.5
06
0.7
0.9
Quantile
Table Simulation results for unconditional quantile estimation given 1,000 samples of Uniform (0.1) random variable. Results are shown for monotonic DLNs trained with the expected pinball loss over uniformly sampled T, sample quantiles, and the Harrell- estimator. DLNs are learned with keypoints in order to embed the correct linear assumption on the inverse CDF.
Model
T= 0.1
T 0.5
== 0.9
DLN, Expected Pinball
0.000109 + 0.000005
0.000077 + 0.000004
0.000109 + 0.000005
Sample Quantile
0.000097 0.000004
0.000245 + 0.000010
0.000086 0.000004
Harrell-Davis
0.000089 0.000004
0.000247 + 0.000011
0.000083 0.000004
Overall, these examples illustrate the sometimes counterintuitive fact that it can be more asymptotically efficient estimate a quantile through learning full inverse CDF than it would be to just estimate the single quantile alone.

---

Published in Transactions on Machine Learning Research (01/2024)
3.3 Expected Pinball Loss with Beta Distribution
So far, our examples have analyzed the case where the expected pinball loss is taken over the uniform distribution; that is, equal weight is put on the pinball loss for each quantile. Theorem is more general, though, allowing for any distribution over T.
We now explore the case where we put more weight on pinball losses around a desired quantile using the Beta distribution. This makes less likely that the overall inverse CDF will be high quality, but offers the possibility that an individual quantile estimate may improve by focusing the Beta distribution on it. We parameterize the Beta distribution by its mode and concentration. This allows us to fix the mode parameter M (the peak of the beta distribution) interpolating between the uniform distribution (C 2) and the point-mass distribution with all weight at the mode (C + 00).
desired quantile and the concentration parameter C as smoothly
Consider again Figure which found that taking the sample quantile achioved similar performance to training with an expected pinball loss when comes to estimating the 0.7 quantile. What if we only cared about that quantile; could we do better?
We analyze this with slight change to setup in Section 3.2.4
where replace T Unif(O, 1) with T Beta(M [M 0.7, C' æ) for different values of z. We reuse the Q and V matrices and only change the moments we plug in. The Beta(a, (3) distribution has E(x²) = I atr and we can convert from a+8+r1 {M, C} to {a, 3} via the formulae a = M(C - 2) + and B = (1 - M) - -
Asymptotic Convergence Rates, 0.7 quantile
Figure 2: Plots of the asymptotic variances for the 0.7 quan- tile for linear models trained with a Beta distribution ex- pected pinball loss. The Beta distributions have modes of 0.7 and varied concentrations. The dotted lines show the asymptotic variances for uniform expected pinball loss es- timator and the sample quantile; these values are equivalent to the values of the curves in Figure | for the 0.7 quan- tile. As in Figure [ we assume the true distribution is Uniform(a, 3) with 3 - ==
02150
02125
0.100
0.2500
0.2025
00000
01975
Beta Pinball Loss Sample Quantile
01150
10
10a
10a
Beta Concentration
Figure shows the results. C=2 equivalent the uniform pinball loss, which is slightly worse than the sample quantile. As the concentration parameter increases, we improve to a point before getting worse and asymptoting to the performance of the sample quantile. These results show the influence of local quantile smoothing; we can sometimes do better by primarily considering nearby data points instead of all data points, even when our chosen function matches the true distribution. We will see in Section 3.5 that local smoothing becomes particularly important when the function class is misspecified.
We further explore the performance of the Beta distribution in Section
3.4 Extension to Conditional Quantile Estimation
The asymptotic theory so far shows that expected pinball loss estimators can sometimes outperform sin- gle quantile estimators in the absence of conditioning on any additional features X. Conditional quantile estimation in general adds considerable amount of theoretical complexity which do not fully address here, and to the best of our knowledge, characterization of asymptotic properties for full conditional in- verse CDF estimation remains an open problem. However, the provided unconditional asymptotic theory can yield intuition and extreme examples for cases when estimating the full conditional inverse CDF would outperform single conditional quantile estimation. For instance, when X is a finite categorical feature, and the function f and parameter space are expressive enough in X, Theorem I can be dircctly applicd the individual quantile estimation problems conditioned on each finite value of X. We formalize this in the following Corollary.
8

---

Published in Transactions on Machine Learning Research (01/2024)
Corollary 1. Let X be distributed uniformly over finite set of categorical values X, letting (2). Suppose that a dataset {(r,, u:) is created by sampling {uş { values IID from the conditional distribution _m
Y|X- = x for each value of x (assuming m divides n),, and laking the union of these sels.
Let f: X, (0, 1), + R be fully separable over x; that is, let e given parameter space Gu for each value of x X. Let f take the form f(x, T;0) EzJEx 1 =T;)g(T; Oxo) for some function 9 (0, (), Oo + R.
mixx...xmm where O2 represents copy of
Suppose f(x,T;0) is well specified with unique 0" E for all T E (0,1) and all x E X, and that H* is also a unique minimizer for the rish R(0) = L Ezex (Y, T;0)|X = x|. Suppose the estimator ((m) is weakly consistent, ((n) _x 0-. Sup- pose further that the function 8 + f(z, T: continuous, locally Lipschitz, and twice differentiable at f o* for all (0,1) and all X. Then for each E X,
θ S Rn such that f(x,T;0" = gr (z)
/n/m(f(z, f(x, 3 N (0. Vof(z,T;0") Q (of(a,r;;))
where Q = (T = Er(T(1- T)r(T)), with pz (v) being the density of Y = x and [(r) = V (1, Vo_
More generally, function classes with high degree of fiexibility such as neural networks may approach this extreme categorical separability described above. We describe one such flexible function class below which, given enough parameters, can be specified to meet this separability condition for finite categorical features. Of course, computational constraints and generalizability must also be considered.
3.5 Model Misspecification And Local Quantile Smoothing
We now expand our simulations to a case where the model's functional form f(r:0) is not well-specified,
specifically estimating the true quantiles of the exponential distribution (with X = 1) given a sample of data. While the true quantile function is -In(1 -T), we learn a linear or piecewise-lineat f (i.e. 1-dimensional DIN; see Section 4 trained with the proposed expected pinball loss and monotonicity constraint т. We compare to the sample quantile and the Harrell Davis estimator (discussed in Section 3.1. which also does not rely on knowing the underlying distribution, but does not extend lo quantile regression with features.
Results in Table show the average MSE when we use 10-fold cross-validation over the expected pinball loss to choose the number of knots K in the piccewiso-lincar function f(r). The DLN is better than the sample quantile al predicting the median, but not as good as the Harrell-Davis estimator, but outperforms both for the more extreme quantile 7 0.99.
Table 2: Estimating the median and 99th percentile of exponential distribution via cross-validated DLN training. Here, the number of knots in the DLN (a piecewise linear function on 7) is tuned via 10-fold cross-validation over the N = 51 or N = 505 training points, respectively, and the model is then retrained on the full data with the selected number of keypoints.
Model
0.5
0.99
CV-DLN, Expected Pinball
0.019263 + 0.00086
0.17566 0.0080
Sample Quantile
0.019968 + 0.00046
0.19552 + 0.0043
Harrell-Davis
0.018081 + 0.00038
0.18887 = 0.0042
We can better understand how the chosen ficxibility of f(T) affects the quality of our quantile estimator by considering the performance for different numbers of keypoints K in the piecewise linear function, as shown in the top row of Figure 3 At a high-level, we see a consistent and intuitive pattern: low performance for very low number of keypoints which oversmooths across т, great performance for an intermediate number performance for a high number of keypoints
of keypoints which regularizes, and dlse-to-ssamll-quuee the effect of very flexible model ((r) is similar to the individual pinball loss.
9

---

Published in Transactions on Machine Learning Research (01/2024)
0.021
Avg MSE to true q59
0.020
0.20
0.019
0.18
MSE
0.018
0.16
HD
0.017
Sample
Sample
DLN
DLN
100
150
250
Number of keypoints
Number of keypoints
0.021
median
0.020
Avg MSE to true q99
0.20
0.019
0.018
0.19
0.017
0.18
HD
0.016
Sample DLN
Sample
0.17
DLN
0.015
10'
10*
10a
Beta concentration
Beta concentration
Figure 3: Estimating the median (left plots) and 99th percentile (right plots) of an exponential distribution (no features) by the sample median, by Harrell-Davis, and by solving ① and (2) with DLN function class. The median estimates were from N = 51 training samples. The т .99 estimates were from N 505 training samples. Error bars are computed over 1,000 repeats. Top: The Pr is the uniform distribution, and the DLN is a piece-wise linear function (PLF) with K keypoints. Larger
makes the model more flexible, and imposcs less smoothing across quantiles. Results for different choices of K are shown across the x-axis. Bottom: The DLN just a linear model on t E [0, 1] with 0 E R and Pr a beta distribution. Results are shown for six choices of the beta concentration across the x-axis, where bigger concentration makes Py more peaked around the target quantile, imposing less smoothing across quantiles.
With a low number of keypoints (note: results are even worse for very-low K which we cut off for ease of visualization), the model misspecification is the dominant effect. Limited to only few keypoints, the function class cannot capture the exponential distribution's quantile function, and the quantile estimates are poor. With high number of keypoints, meanwhile, we mechanistically approach the performance of the sample quantile; note that at the extreme when we have the same number of keypoints as data points, our expected pinball loss model will exactly recover every sample quantile.
In the middle, we find that the disadvantages of misspecification are outweighted by the benefit of smoothing across ncarby data samples when estimating any one quantile. For example, consider the high-performing model with K 100 when predicting the 0.99 quantile on an N 505 dataset; that suggests that each keypoint is influenced by around data points each, and locally linear assumption on that scale improves performance.
We confirm the key role of local smoothing across nearby quantiles by running another experiment. Again, we seek to predict the exponential distribution; but this time, we do so with a linear model f rather than piccewise-lincar model. We focus this lincar model around the neighborhood of our desired quantile by training with Beta distribution over pinball losses instead of Uniform distribution over pinball losses. This setup is very similar to the theoretical analysis conducted in Figure 2
10

---

Published in Transactions on Machine Learning Research (01/2024)
We observe, in the bottom row of Figure + qualitatively similar pattern in the top row. A Beta distribution with low concentration around the desired quantile performs poorly, and one with a very high concentration performs similarly to the sample quantile. Intermediate concentrations, ie. concentrations that fit a linear approximation over nearby sample quantiles, see improved performance.
In practice, whether we train with a flexible f over all pinball losses equally or with potentially less-flexible f over neighboring pinball losses depends on whether we are interested in the full distribution, as well as the particulars of our problem. Section 5.5 contains comprehensive set of experiments that compare various strategies.
To recap, we have shown that we can replicate or even beat the performance of not only the sample quantile but unconditional quantile estimators such as the Harrell-Davis estimalor as well. We do this, though, with a strategy that can seamlessly scale to the problems of conditional quantile regression and full inverse CDF estimation.
4 DLNs for Legitimate but Flexible Inverse CDFs
We propose using DLNs You et al. 2017 as the function class in (1] because of their flexibility and because they efficiently enable three key regularization stratogics: non-crossing constraints, restricting the Icarned distribution to location-scale families, and monotonic features, as explained below.
4.1 Review of Deep Lattice Networks
DLNs are multi-layer models where some layers are lattices or ensembles of latticos. Lattices build upon a well-established strategy for approximating functions over a bounded domain (Campbell- el. al.] 2003): record the function's values on regular grid of imputs, and then one can linearly interpolate between the recorded values. Generalizing to multiple dimensions, lattice function is multi-dinensional look-up table whose parameters are the look-up table values, and the function is produced by linearly interpolating the look-up table values surrounding an input point. See Fig. 4 for example 2D lattices.
In the special case of a one-dimensional domain, a lattice is just a piecewise-linees function. We recommend reading Gupta et al.] (2016) for the basics of lattices with monotonicity constraints. See also Garcia et al. (202) for basics of lattice models and their relationship to other splines.
Mathematically, lattice function f(x) R + R can be expressed f(x) = 0 (b(x), where the parameters 0 E R² are the look-up table values, and the nonlinear transformation v(xx) R2 maps an input x E R to the ultra-sparse vector of linear interpolation weights over the look-up table values.
Garcia & Gupta (2099) showed that lattice models can be trained using a standard empirical risk minimiza- tion framework. Lattices can approximate any continuous bounded function if they are sampled on regular grid with enough knots.
A key reason to use lattice functions in an application like quantile regression is that the regular structure of their parameters makes them amenable to shape constraints. Here, the constraint 2 requires f(x, r; 0) to be monotonically increasing in which is imposed efficiently with sparse linear inequality constraints on the lattice parameters corresponding to forcing any neighboring gridpoints in the direction of r to be increasing Gupta et al. 2016). To more efficiently represent. functions, lattice models are usually architected with a first layer that separately calibrates cach input with a learned ono-dimensional piecewise linear function (PLF) whose knots need-not be regularly-spacee referred to as calibrator, then fuses the calibrated inputs together with coarser multi-dimmensioall lattice (Garcia et al.) 2012 Gupta et al. 2016).
The main drawback to latticos is that the number of lattice parameters needed to approximate D- dimensional functions blows up exponentially as 0(22). This problem is handled by making an ensemble of calibrated lattices (Canini et al. 2016) much as one makes a random forest, and these calibration and lattice layers can be mixed with other layers like lincar layers, and cascaded ad infinitum, forming arbitrarily deep lattice networks (DLNs) YYou et al., 2017 Cotter et al. 2019). Importantly, as long as each layer of DLN is trained to respect a monotonic feature like T, the entire DI.N will be monotonic in 7.
11

---

Published in Transactions on Machine Learning Research (01/2024)
θ
0
e
0
1.00
0.75
0.50
0.25
0.00
0
0/ 0
02 0
0
0 = (0.0, 0.5, 0.5, 1.0)
0 = (0.2, 0.7, 1.0, 0.4)
(0.9, 0.9, 0.5, 0.0)
Figure Example 2D lattices sized 2 x 2 have four parameters. Each parameter represents the value the lattice function takes at one of the four corners of the input domain [0, The lattice function values in- between the vertices computed by bilinear interpolation of the four parameter values, these functions can be re-parameterized and expressed as bilinear polynomial functions over [0, 1] of the form f(x) a + bæı cr2 + dæ1x2. A 2D lattice with more knots, e.g. 3 x 3, would be locally bilinear across cach cell of the lattice.
More generally, one could unrestricted ReLU embedding layers for the first few model layers ø; then illustrates how arbitrary neural net architecture
fuse in later with 7-monotonic DLN layers. Figure Bigusteeillsstrttss be used to reduce the dimensionality of the input features of before fusing with allowing our approach to be used for essentially any type and size of imput.
Figure 5: Block diagram showing how an arbi- trary neural net archilecture can be used re- duce the dimensionality of to a more manage- able size, before combining in a DLN architecture with T, which is constrained to be monotonic.
Neural Net
Predicted quantile
X
DLN
Unconstrained
Monotonic
In our experiments, our DLNs are built with the standard two-layer calibrated lattice (Gupta al. 2016) and
(Canini et al. 2016) architectures offered by the open-source TensorFlow Lattice library (TensorFlow Blog Post 2020). The TensorFlow Lattice library enables solving the constrained optimization problem ! and (2) with Dykstra's projection algorithm. In our experiments, we found our DLN models trained in a similar amount of time as the DNN and SQF-DNN models.
4.2 Location-Scale Family Models
We show that by different architecture choices for the DLN model used can regularize the estimated inverse CDF in caaaantcalll-mmannnen example, one constrains the DLN architocture to not have any interactions between т and any of the x features, then one learns a regression model with homoskedastic errors. As another example, the basic two- layer DLN called calibrated lottice model (Gupta al.] 2016) described Section 4.1 can be constrained to learn distributions across & that come from shared, learned, location-scale family:
ways, similar to prior more rigid distributional approaches. As a simple
Let f(x, be a calibrated lattice model with piece-wise linear calibrator c(r) [0, 1] - [0, for т, and let there be two knots in the regular grid in the direction of т, and the lattice parameters are interpolated with standard multilincar interpolation. Then f(x, T) represents an inverse CDF function py1(y)>) where the estimated distribution for every from the same location-scale family as the calibrator c(r).
Lemma 2.
12

---

Published in Transactions on Machine Learning Research (01/2024)
The proof is in the Appendix. Note the number of keypoints in c(r) controls the complexity of the learned base distribution, allowing the model to approximate location-scale families like the Gaussian, gamma, or Pareto distributions. Then in the second layer of f, the number of knots in the lattice over the c(T) input controls how much the distribution should be allowed to vary across z. Two knots in the lattice limits us to a shared location-scale family, as noted above, while three knots in the lattice gives an extra degree of freedom to shrink or stretch one side of the distribution differently across x. More lattice knots across T, or ensembling of lattices, or more layers can all be used to move the model towards full generality.
4.3 Other Monotonic Features
Using DLNS also enables imposing monotonicity on the features if domain knowledge says they should have monotonic impact the estimates (Gupta et al., 2016). For example, a pizza delivery time estimator might treat constrain an input quantifying how bad traffic is to only increase estimated quantiles (that is, worse traffic never makes the estimated delivery time shorter).
Our experience with real-world applications of quantile estimation is that real-world quantile estimates do often use features that are past measurements (or otherwise strong correlates) to predict the future distribution of measurements. For example, in the real-world Puzzles experiment in Section 5 we predict the quantilos of how long Hoefnagel Puzzle Club member will keep a rented puzzle, and one of the features is how long they kept their last puzzle. We constrain that input to have monolonically positive effect on the estimate. (However, earlier hold-times are not necessarily monotonically indicative of the next future hold-time, because they may indicate a strongly docreasing trend over time). Domain knowledge about the relationship between the model inputs-and-ouput that can be captured with monotonicity constraints is a tuning-free, semnantiaall-
regularizer that can improve test performance and increase model explainability (Gupta et al. 2016), and where applicable, even help capture ethical rules (Wang & Gupta 2020).
DLNs also enable multi-dimensional shape constraints to capture domain knowledge like complementary features (Gupta et al. 2020), and functions defined on sets (Cotter et al., 2019
5
Experiments
Using simulations and real data, we first show that training a DIN that is monotonic in 7 to satisfy both (1) and 2 works well compared to training a DNN for [I as done in Tagasovska & Lopez-Pat (2019) or the SQF-DNN parametric approach of Gasthaus et al. (2019); the results are in Tables 3 and Then we use the proposed r-monotonic DLN models to compare training to minimize expected pinball loss versus training to minimize the pinball loss for specific t; the results are in Tables 7 Fadd g
Bolded results in tables indicate that the metric is not statistically siguificantly different from the best metric among the models being compared, using an unpaired t-test.
5.1 Model Training Details
All hyperparameters were optimized on validation sets. We used Keras models in TensorFlow 2.2 for the unre- stricted DNN comparisons that optimize (1) (Tagasovska & Lopez-Pas 2019 and the SQF-DNN (Gasthaus
et al. 2019), which optimizes the same objective a different manner while also guarantceing non-crossing
quantile estimates. For DLNs, we used the TensorFlow Lattice library TensorFlow Blog Post 2020 For all DNN and DLN experiments, we use the Adam optimizer (Kingma & Baj 2015) with its default learning rate of 0.001, except where noted. For DNN models, validated the number of hidden layers and the hidden dimension. For the SQF-DNN, we also validated the number of distribution keypoints. For the smaller DLN models, we used the common two-layer calibrated lattice architect (Gupta et al.] (2016) and validated over its number of calibration keypoints and lattice vertices. For the larger datasets, we used an ensemble of calibrated lattices for the DLN model, and then also validated over the the number and dimensionality of the base models in the ensemble Canini et al. 2016). For both DLNs and DNNs, we additionally validated over the number of training epochs.
13

---

Published in Transactions on Machine Learning Research (01/2024)
5.2 Benchmark and Real Datasets Used
We tested on three publicly available datasets and one proprietary dataset for real-world problems. Code is available
wetpeeaarrc/goolleeeeeee researe trre/ maatea//uuuuaan reiillerrrrre.....
Air Quality: The Beijing Multi-Site Air-Quality dataset from UCI (Zhang al. 2017 contains hourly air quality data from 12 monitoring regions around Beijing. We predict the quantiles of the PM2.5 concentration from D = 7 features: temperature, pressure, dew point, rain, wind speed, region, and wind direction. The DLN model is an ensemble of two-layer calibrated lattice models. We split the data by time (not IID) with earlier examples forming a training set of size 252,481, later examples a validation set of size 84,145, and most recent examples a test set of size 84,145.
Puzzles: This dataset is from the Hoefnagel Puzzle Club, private library that rents wooden jigsaw puzzles. We predict how long library member will hold borrowed puzzle given their five previous hold- times and sixth categorical feature that indicates a user's past activity (one of factive users, high-variance users, new users)) The DLN model is an ensemble of two-layer calibrated lattice models. Based on their domain knowledge, we also constrain the DLN output to be monotonically increasing in the most-recent-pest hold-time feature. The 984 train and 247 validation examples are IID from past data, while the 211 test samples are the most recent samples (so the test samples are not IID with the anonymized dataset is publicly available at www.mayagupia.org /atag//ata///a maata//uzzleCubb Holdoooes.css
set). The
Traffic: This is a proprietary dataset from Google for estimating the time to drive route. The DLN is a two-layer calibrated lattice whose D = 4 inputs are 1 categorical and 3 continuous features. We used 1,000 examples each for training, validation, and testing, with the training examples occurring earlier in time than the validation and test examples (not IID).
We predict the quantiles of wine quality on 100-point scale. We used D 42 features, including price, 1-d learned embedding for the country of origin (aka categorical calibration (Gupta et al. 2016 and 40 Boolean features describing cach wine. The DLN model is an ensemble of two-layer calibrated lattice models. The DLN outpul is constrained to be monotonically increasing in the price feature, as suggested in prior work (Gupta et al.] 2018). The data was split IID with 84,641 examples for training, 12,091 for validation, and 24,184 for testing.
Wine: We used the Wine Reviews dataset from Kaggle (Bahri 2018
5.3 Model Architecture Experiments
We start by demonstrating the efficacy of using monotonic DLNs to predict the inverse CDF on simulations, and then on the real data.
5.3.1 Simulations
We used simulations from a recent quantile regression survey paper (Torossian al. 2020) based on the sine, Griewank, Michalewiez, and Ackley functions with carefully designed noise distributions to represent range of variances and skews across the respective input domains. We used 250 training examples for the 1-D sine and Michalewicz functions, 1,000 examples for the 2-D Griewank function, and 10,000 examples for the 9-D Ackley function.
Our metric the average squared difference between the estimated and true quantile curves sampled at 99 uniform quantiles 7 e 0.01,..., 0.99, averaged over 100 repeats, and averaged over values of across the domain (we use fine grids of 1,000 and 10,000 & values for the 1-D and 2-D experiments, respectively, and random sample of 10,000 v-values for the 9-D experiment). We also report the fraction of test points for which at least two of their 99 quantiles crossed.
Results in Table 3 demonstrate that the proposed monotonic-in-r DLNs are the best or statistically tied for the best across all four disparate simulations. The DNNs were trained with the same sampled expected pinball loss, and are sometimes close in performance, but suffer substantially from crossing quantiles. For example, on the sine-skew distribution, quantile-crossing was observed for the DNN model between at least two quantiles on 38% of test x values!
14

---

Published in Transactions on Machine Learning Research (01/2024)
Spline quantile functions (Gasthaus al. 2019 avoid crossing by construction, but their accuracy was much worse than the DIN accuracy on three of the four datasets, and not better on the fourth. In particular, the DLNs do much better at the Griewank and Ackley simulations, which have the lowest signal-to-noise ratios Torossian et al. 2020
Table 3: Simulations: Quantile MSE and percent crossing violations for E {0.01, 0.02, 0.99).
Sine-skew (1,7)
Griewank
Michalewiez
Ackley
Model
MSE, Crossing
MSE, Crossing
MSE, Crossing
MSE, Crossing
DNN
3.53 = 0.09, 38%
1.29 + 0.02, 5%
0.311 + 0.011, 12%
237 = 6, 0.9%
SQF-DNN
5.68 + 0.11, 0%
1.05 + 0.01, 0%
0.232 + 0.006, 0%
667 + 81, 0%
DLN
3.51 = 0.13, 0%
0.55 + 0.01, 0%
0.219 + 0.006, 0%
206 = 2, 0%
5.4 Training Time Comparison
We also provide Liming results in Table 4 for the Ackley simulation, which complex 9-dimensional problem required non-trivial training times. We compared timings of the different models with their validated hyperparameters. The validated DNN chose half the number training batches as the DLN, and this caused it to be twice as fast. Each training batch took roughly the same amount of time though, and more gramular hyperparameter tuning the number of training batches could have led to different results. The SQF-DNN, meanwhile, noticeably slower to train. However, while these implementations were comparable in their use of software and processors, all of these approaches could certainly be better optimized for speed. For example, compiler-optimiae DL.N models can execute ten times faster than an optimized C++ implementation, which was 100 times faster than TensorFlow implementation (Zhang et al. 2021
Table 4: Simulations: Training time, in seconds, for the validated model of each type the Ackley simulation from Table 3 We also report the number of training steps used by each model (i.e. the number of batches seen by the model) and the time taken per training step.
Model
Train Time (s) Training Steps
Time Per Batch (ms)
DNN
3446
500000
6.9
SQF-DNN
28042
100000
280
DLN
8444
1000000
8.4
5.4.1 Real Data
Table compares these models on the four real datasets. We report the pinball loss averaged over T E 0.02,, .., 0.99), which is a standard metric for judging conditional quantile estimates, and one which approximates the continuous ranked probability score metric used for measuring the quality of probabilistic forecasts (Gasthaus et al. 2019 Yan ct al. 2018).
The 7-monotonic DLNs performed the best on-average for all four problems. Unlike the simulation results in Table 3 the SQF-DNN (monotonic by construction) performs more similarly to the DNN on these real datasots.
We also report the absolute devialion calibration error CChung et al. 2021 Table t (some authors use a squared variant [Kuleshov al. 2018)), the formula is given in A.8 By that metric the DLN is better or stat. sig. the same SQF-DNN, with substantial wins over SQF-DNN on Traffic and Wine.
5.5 Expected Pinball Loss Experiments
We next empirically show that training a full inverse CDF by minimizing the expected pinball loss is at least good at estimating specific quantiles training with the pinball loss just. for those specific quantiles.
15

---

Published in Transactions on Machine Learning Research (01/2024)
Table
Real data experiments: Pinball loss on the test set, averaged over 7
0.02, (0.01,0.02.....0.9999)) 0.999..
Model
Air Quality
Puzzles
Traffic
Wine
DNN
18.041 + 0.088
3.229 + 0.017
0.0494 + 0.00037
0.6603 + 0.0088
SQF-DNN
17.356 + 0.106
3.634
= 0.142
0.0491 + 0.00022
0.6552 = 0.0037
DLN
17.280 + 0.157
3.188 + 0.041
0.0479 + 0.00003
0.6381 + 0.0005
Table 6:
Real data experiments:
Calibration on the test set, averaged over {0.01, 0.02, 0.99).
Model
Air Quality
Puzzles
Traffic
Wine
DNN
0.0761 + 0.0038
0.0821 + 0.0032
0.0562 0.0014
0.0253 + 0.0010
SQF-DNN
0.0798 0.0042
0.0614 = 0.0038
0.0600 0.0062
0.0671 0.0083
DLN
0.0588 0.0057
0.0739 0.0078
0.0210 0.0027
0.0125 + 0.0031
5.5.1 Simulations
We ran simulations on the 1D sine-skew function as per Torossian al. (0020) with three noise scenarios illustrated in Figure [ noise parameters (a, (1,1) for symmetric low-noise, (7,7) for symmetric high- noise, and (1, for asymmetric high-noise. These simulations are the simplest extension of the unconditioned context studied empirically and theoretically in Section B to the conditioned case, in that there a single æ feature, and enough data for a flexible quantile regression model to approximate the corresponding sample quantile y within a small neighborhood of any value x.
0-
-20
-20
-20
-30
-30
-30
-1.00
-0.5
-0.00
-0.5
0.5
050
0.5
100
-1.00
-0.75
-0.00
-0.2
600
025
030
1.00
-1.00
-035
-0.0
-0.5
0.5
050
075
100
(a) Sine-skew (1,1)
(b) Sine-skew (7,7)
(c) Sine-skew (1,7)
Figure 6: The colored lines show the true quantiles for = 0.1, 0.5, 0.6, 0.9 for the sine-skew distribution simulations Torossian et al. (2020).
The results shown in Table 7 are the error in predicting the median given either 100 or = 1, 000 training samples. Overall, the single pinball loss model struggles relative to the expected pinball loss models. even we increase the number of data points from 100 to 1000, at which point there are about 50 data points within each 0.1 region of E [-1,1].
The uniform pinball loss model which trains equally across all quantiles performs quite well in both symmetric-noise cases 1) and (7,7) However, for the simulation with highly skewed noise (1, it strug- gles a bit relative to the models trained with a Beta pinball loss, perhaps because the extreme quantiles start becoming very uninformative about the median. Still, even in that case, performs similarly to the single pinball loss models.
5.5.2 Real Data
In Table g we give the accuracy at predicting three specific target quantiles on the four real datascts by either training single model estimating the full inverse CDF by minimizing Ey where is drawn uniformly from (0, 1); training single model to minimize the discrete of the three pinball losses for the
16

---

Published in Transactions on Machine Learning Research (01/2024)
Table 7: Sine-skew experiment with different sine-skew noise choices (1,1), (7,7) and (1,7). We present results for DLNs trained with four different distributions over т: the uniform distribution, the beta distribution with a mode of 0.5 and concentration of 10, the beta distribution with a mode of 0.5 and concentration of 100, and just the 0.5 pinball loss. The metric is the MSE between the true median and estimated median, averaged over x 1ni
(a, b)
N
min ET, Unif
min ET, Beta(10)
ET, (eta
min single-7 loss
(1,1)
100
0.421 + 0.025
0.348 + 0.015
0.301 + 0.016
0.588 + 0.036
(1,1)
1,000
0.050 + 0.003
0.057 = 0.002
0.056 = 0.002
0.067 0.003
(7.7)
100
7.334 + 0.367
6.654 + 0.330
7.402 + 0.397
9.860 + 0.560
(7.7)
1,000
0.967 + 0.058
0.861 = 0.036
1.057 + 0.037
1.470 + 0.086
(1, 7) 100
4.458 0.259
3.515 + 0.178
3.031 0.179
4.416 0.386
(1,7) 1,000
0.451 + 0.037
0.336 + 0.022
0.338 + 0.020
0.542 0.067
considered; or training three separate models that minimize ET where T is concentrated around the target quantiles; or training three separate models that minimize each spocific r's pinball loss.
The test metric is computed on finite test sets, unlike the simulations where we could compare to the true quantiles.
When measuring the accuracy of the predicted median, via the == 0.5 pinball loss, there is no statistically significant difference between the model trained for the single quantile individually and the model trained to estimate the full inverse CDF using a distribution of T ~ Unif(0, 1) for the loss. In other words, the estimate of the median does not suffer from trying to estimate the full inverse CDF.
At the more extreme quantiles, the single-r estimator sometimes outperforms the full inverse CDF estimator, e.g. the 99th percentile for Air Quality, and the 10th and 90th percentile for Wine. This is possibly due to the fact that the smoothing from the ET loss is more one-sided. Additionally, for Wine, we hypothesize this is because it is a fairly easy regression problem with a large number of training samples and one continuous feature (wine prico) that correlates highly with the label (wine quality), so the extra regularization is not helping. That said, the single r estimator is still not always the best even on extreme quantiles; for the large Traffic dataset, the 99th percentile is statistically tied, while on Puzzles the single 7 estimator for the 90th percentile worse
In all cases, smoothing from 7 Beta is still beneficial: training with sampled from Beta distribution is statistically similar or better than training with a single pinball loss on all datasets. This tracks with the simulated results and remains true even on extreme quantiles.
6 Conclusions And Some Open Questions
We gave theoretical and empirical evidence that training quantile regression models by minimizing the expected pinball loss can perform as well or better on any individual quantile than models trained for that quantile alone. We showed that minimizing the expected pinball loss is not sufficient to provide legitimate inverse CDF estimates, one must also produce a model whose quantiles are guarantced to not cross. We showed that DLN models can be effectively constrained be monotonic in T and produce stale-of-the-art quantile estimates. The non-crossing issue is important not just for theoretical reasons, but because it is confusing to many non-experts if model's quantile estimates cross, eroding trust in that model, and by association, all Al.
We provide a new asymptotic convergence analysis for estimation of the full inverse CDF, assuming an oracle that minimizes the empirical loss. We leave open the question of non-asymptotic variances, and whether some of our assumptions can be relaxed, or whether lower variances can be had with some more stringent assumptions. Theorem is for unconditioned quantile estimation, with simple extension to separable
17

---

Published in Transactions on Machine Learning Research (01/2024)
Table 8: Accuracy of single quantile estimates, comparing one model trained minimize the expected pinball loss with respect to uniform Py; one model trained to minimize the sum of the three pinball losses for the three quantiles shown (Discrete); three separate models trained with individual beta distributions PT with the mode set at the target quantile and the beta concentration chosen on the validation set; and three separate models trained with individual r pinball losses (Single). All models were DLNs with the same architecture as in Table
5 and monotonic on r and any input features noted in the text.
Air Quality:
Model
Pinball loss (T = 0.5)
Pinball loss (T = 0.9)
Pinball loss (r = 0.99)
Unif(O. 1)
23.576 + 0.047
14.850 0.075
2.947 + 0.025
Discrete
24.083 + 0.073
15.439 = 0.062
3.042 + 0.012
T ~ Beta
23.586 + 0.183
13.942 + 0.054
2.709 = 0.029
Single E T
23.634 + 0.068
14.908 = 0.092
2.700 = 0.013
Traffic:
Model
Pinball loss (r 0.5)
Pinball loss (r 0.9)
Pinball loss (r 0.99)
T Unif(O, 1)
0.064386 + 0.00014
0.040159 0.00018
0.010645 0.00018
T Discrete
0. 0.064578 + 0.00028
0. 0.039804 + 0.00014
0.011290 0.00014
Beta
0.064988 + 0.00030
0.039549 + 0.00009
0.01064 + 0.00011
Single T
0.064791 + 0.00012
0.039900 + 0.00010
0.01070 + 0.00018
Wine:
Model
Pinball loss (7 0.1)
Pinball loss (r 0.5)
Pinball loss (T 0.9)
T Unif(O, 1)
0.4099 + 0.0006
0.8889 + 0.0017
0.3773 = 0.0019
5 Discrete
0.4094 = 0.0005
0.8891 + 0.0010
0.3706 = 0.0003
t ~ Beta
0.4047 = 0.0005
0.8871 0.0008
0.3665 + 0.0004
Single T
0.4049 + 0.0004
0.8867 0.0006
0.3663 +'2 0.0003
Puzzles: Model
Pinball loss (r = 0.5)
Pinball loss (r = 0.7)
Pinball loss (T 0.9)
r Unif(O, 1)
4.437 + 0.148
4.427 + 0.202
2.696 + 0.091
r Discrete
4.265 + 0.052
4.192 + 0.044
2.686 + 0.035
5 ~ Beta
4.188 + 0.049
4.259 0.054
2.680 + 0.020
Single e
4.258 + 0.056
4.255 + 0.057
2.793 + 0.049
18

---

Published in Transactions on Machine Learning Research (01/2024)
conditional quantile estimation in Corollary 1. Important open work would be providing stronger results for conditional quantile estimation where there is smoothing across the features.
A key question is the choice of distribution of T when minimizing the expected pinball loss. We note one can control the regularization across quantiles both by the shape of Pr and by the model flexibility with respect to T. Thus we showed a uniform distribution Pr generally provides good results one is validating the model flexibility over T, and a uniform distribution has no hyperparameters and gives you good estimate of the entire inverse CDF for free. However, if one is only interested in a single т, also showed that one may do better with beta distribution for PT, and that a concentration of 10 likely makes good default choice, but for best results one would want to validate the concentration.
One open question is whether the proposed methodology improves the accuracy of prediction intervals [Pearee et al. 2018), as measured by traditional metrics such as interval width and calibration. Tagasovskal & Lopez-Paz (2019) have shown that minimizing the expected pinball loss can be useful strategy in that regard, and our results in this paper about improved quantile accuracy and calibration suggest there may be improved performance in prediction intervals as well.
Our experiments, like most published quantile regression experiments, were limited to datascts with D = 42 features and less. This is partly because one needs a large amount of data to make decent estimates of tail quantiles. In theory the proposed methodology can be applied to any size problem by only bringing T in for the last layers as shown in Fig. 4, but an open question is how well this or any quantile regression works in practice with a large number of features, and whether there are novel scaling issues.
In particular, we hypothesize that for larger models there may be increase in quantile collapse (Lopez- Martin et al. 2021). Quantile collapse is the edge case of non-crossing, defined as flat f(T) over some range of T. For example, the airline tells passengers there is an 80% chance the flight will arrive in 4 hours and 6 minutes, and that there is 90% chance the flight will arrive in 4 hours and 6 minutes. Our monotonicity constraint (2) allows quantile collapse because it only requires monotonicity in т, not strict monotonicity.
Quantile collapse can happen in regions of sparse data where there are no training samples to fit to such that the model only has to satisfy (2), which it can do with flat f(z, r). We belive this can partly be solved by initializing with f(z, T) that is strictly monotonically incrcasing (as we do in experiments), as then one only gets quantile collapse if there is enough relevant training data to pull the model to the edge of the feasibility of .
Quantile collapse is also a risk if the model f(x, r) has too much fiexibility in x such that the data is sparse with respect that flexibility. For example, assume x is a categorical feature with I million possible values (e.g. a language model), and that (xx, is flexible enough learn an independent estimate for each value of x. Suppose for one of those categorical values there is only one training sample w:). Then since the pinball loss (expected or standard) will be minimized when f(z, r) Yr, there is no dependence, and for that category you will have quantile collapse. More generally, if you have a very flexible model and sparse data in some regions, we expect you may get quantile collapsc if no other regularization to avoid it is used.
Broader Impact Statement
We hope the broader impact of this work be positive in increasing the deserved trustworthiness of AI by promoting quantile regression models with no crossing-quantile. The alternative of quantile estimates that cross will undoubtedly erode public confidence in AI. For example, prior work would launch a quantile regression model that predicts there is a 90 percent chance your pizza will be delivered in 30 minutes, and that there is an 80 percent chance it will be delivered in 34 minutes. Such models are unnecessarily confusing and untrustworthy for non-experts, and using them in practice may then reduce the ability of worthy AI to have a positive impact on the world.
References
D. Bahri. Wine reviews. Kaggle, 2018 URL
19

---

Published in Transactions on Machine Learning Research (01/2024)
P. K. Bhattacharya and A. K. Gangopadhyay. Kernel and nearest-neighbor estimation of conditional quantile. The Annals of Statistics, 1990.
H. D. Bondell, B. J. Reich, and H. Wang. Non-crossing quantile regression curve estimation. Biometrika, 2010.
M. Campbell- Kelly, M. Croaken, and E. Robson. The History of Mathematical Tables: From Sumer to Spreadsheets. Oxford University Press, Oxford, England, 2003.
K. Canini, A. Cotter, M. M. Fard, M. R. Gupta, and J. Pfeifer. Fast and flexible monotonic functions with ensembles of lattices. NeurIPS, 2016.
A. J. Cannon. Non-crossing nonlinear regression quantiles. Stochastic Environmental Research and Risk Assessment, 32:3207-3225, 2018.
V. Chernozhukov, Extremal quantile regression. Annals of Statistics, 33(2), 2005.
P. Chilinski and R. Silva. Neural likelihoods via cumulative distribution functions. arXiv, 2018.
Y. Chung, W. Neiswanger, Char, and J. Schneider. Beyond pinball loss: Quantile methods for calibrated uncertainty quantification. NeurIPS, 2021.
A. Cotter, M. R. Gupta, H. Jiang, E, Louidor, J. Muller, T. Narayan, S. Wang, and Zhu. Shape constraints for set functions. ICML, 2019.
Daniels and M. Velikova. Monotone and partially monotone neural networks. IEEE Trans. Neural Networks, 21(6):906-917, 2010.
C. E. David and S. M. Steinberg. Quantile estimation. In Encyclopedia of Statistical Sciences, volume 7. Wiley, New York, 1986.
M. Palk. Relative deficiency of kernel type estimators of quantiles. Annols of Statistics, 12, 1984.
E. K. Garcia and M. Gupta. Lattice regression. NewrIPS, 2009.
E. K. Garcia, R. Arora, and M. R. Gupta. Optimized regression for efficient function evaluation. 1EEE Trans. Image Processing, (199):122844444 September 2012.
J. Gasthaus, K. Benidis, Y. Wang, S. Rangapuram, D. Salinas, V. Flunkert, and T. Jamuschowski. Proba- bilistic forecasting with spline quantile function RNNs. AIStats, 2019.
M. R. Gupta, A. Cotter, J. Pfeifer, K. Voevodski, K. Canini, A. Mangylov, W. Moczydlowski, and A. Van Esbroeck. Monotonic calibrated interpolated look-up tables. Journal of Machine Learning Research (JMLR), 17(109):1-47, 2016. URL http://jmlr http://jml..or/ppprrret orggppppers//1715-223 html
M. R. Gupta, D. Bahri, A. Cotter, and K. Canini. Diminishing returns shape constraints for interpretability and regularization. NearIPS, 2018.
M. R. Gupta, E. Louidor, O. Mangylov, N. Morioka, Narayan, and S. Zhao. Multidimennional shape constraints. ICML, 2020.
F. E. Harrell and C. E. Davis. A new distribution-free quantile estimator. Biometrika, 1982.
Fumio Hayashi. Econometrics. Princeton University Press, 2011.
X. He. Quantile curves without crossing. The American Statistician, 51(2):186 192, 1997.
W. D. Kaigh and P. A. Lachenbruch. A generalized quantile estimator. Communications in Statistics Theory and Melhods, 1982.
D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015.
20

---

Published in Transactions on Machine Learning Research (01/2024)
R. Koenker. Quantile Regression. Cambridge University Press, 2005.
R. Koenker and G. Basselt. Regression quantiles. In Econometrica, 1978.
V. Kuleshov, N. Fenner, and S. Ermon. Accurate uncertainties for deep learning using calibrated regression. In ICML, 2018.
B. Lang. Monotonic multi-layer perceptron networks as universal approximators. Artificial neural networks: formal models and their applications-CAANN 2005.
M. Lopez-Martin, A. Banchez-Esgueville
Hernandez Calleej, J. I. Arribas, and B. Carro. Additive ensemble neural network with constrained weighted quantile loss for probabilistic electric-load forecasting. Sensors, 2021.
N. Meinshausen. Quantile regression forests. Journal of Machine Learning Research (JMLR), 2006.
A. Minin, M. Velikova, B. Lang, and H. Daniels. Comparison of universal approximators incorporating partial monotonicity by structure. Neural Networks, 23(4):471-475, 2010.
T. Pearce, A. Brintrup, M. Zaki, and A. Neely. High-quality prediction intervals for deep learning: distribution-fres ensembled approach. ICML, 2018.
James L. Powell. Notes on median and quantile regression, 2010.
R.-D. Reiss. Estimation of quantiles in certain nonparametric models. Annals of Statistics, 8, 1980.
J. Sill. Monotonic networks, NeurIPS, 1998.
J. Sill and Y. S. Abu-Mostafa. Monotonicity hints. NewrIPS, 1997.
N. Tagasovska and D. Lopez-Paz. Single-model uncertainties for deep learning, NeurIPS, 2019.
1. Takeuchi, Q V. Le, T. D. Sears, and A. J. Smola. Nonparametric quantile estimation. Journal of Machine Learning Research (JMLR), 2006.
TensorFlow Blog Post.
TensorFlow Lattice: flexible, controlled, and interpretable ML,
2020.
http:///ww..orrooo000000o00oooa0 controlled-ance interpettable-M.........
URL
L. Torossian, V. Picheny, R. Faivre, and A. Garivier. A review on quantile regression for stochastic computer experiments. Retiable Engineering & System Safely, 2020.
A. W. van der Vaart. Asymplotic Statislics. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 1998.
S. Wang and M. R. Gupta. Deontological ethics by monotonicity shape constraints. In ATStats, 2020.
X. Yan, W. Zhang, L. Ma, W. Liu, and Q Wu. Parsimonious quantile regression of financial asset tail dynamics via sequential learning. NewrIPS, 2018.
S. You, K. Canini, D. Ding, J. Pfeifer, and M. R. Gupta. Deep lattice networks and partial monotonic functions. NewrIPS, 2017.
N. Zhang, K. Canini, S. Silva, and M. Gupta. Fast linear interpolation. ACM Journal on Emerging Technologies in Computing Systems, 2021.
S. Zhang, B. Guo, A. Dong, J. He, Z. Xu, and S.X. Chen. Cautionary tales air-quality improvement in Beijing. Proceedings of the Royal Society A, 473, 2017.
21

---

Published in Transactions on Machine Learning Research (01/2024)
A Appendix
A.1 Proof of Lemma 2
Lemma 2 Let T(,,T) be calibrated lottice model with piece-wise linear calibrator c(T) [0, + [0,1] for and let there be two knols in the regular grid in the direction of T, and the lallice porameters are interpolated with standard multilinear interpolation. Then f(z, T) represents an inverse CDF function (ylx) where the estimated distribution for every x is from the same location-scale family as the calibrator c(T).
Proof. I random variable Y conditioned on belongs to the location-scale family then for T E (0, 1) and and 6 0, it must hold that the conditional inverse CDF satisfies Frjx= a (r) Note that for r E (0,1), interpolating a lattice with two lattice vertices in the 7 dimension yields the estimate Prix -rix=;(T) (7) = f(2,T) = f(z, 0) + (f,) f(z,0)) Thus mapping to the location-scale property, a = fzz, 0), = 1), f(=,)), and (T) = c(7). Thus every estimated conditional inverse CDF A- (T) is a translation and scaling of the piecewise linear function c(T).
A.2 Proof of Theorem 1
Theorem 1 Suppose (T: is well specified with unique 0r E O R" such that f(r; F for all T E 1), and let 'T be distributed such that f* is also a unique minimizer for the rish R(0) F ( ((T 0))] Suppose the estimator ((m) is weakly consistent, ((n) P, 0" Suppose further that the function + f (T; is continuous, localty Lipschitz, and twice differentiable at 0 = O for all T E (0, 1).
Then,
- N (O, )
where,
Q Q= Q=F = -
with I'(7) = Vof(::"")
Consider the loss function Lyy,0) = Er (u. ((; Let g(n) minimize the risk with an empirical distribution over Y with IID samples {,1)_11 ((m) arg min, 1 Σx= L(v,, 0). We assume that Y is a continuous random variable with density p(y), and T is independent of Y.
Proof.
We apply Theorem 5.23 from van der Vaarl 1998
which says that under some regularity conditions,
(n((()) - -0) 0-) -
VoL(y,,0)
where o1(1) refers to sequence of random variables that converges in probability to 0: X, = or(1) if and only if X, P, 0. Specifically, Theorem 5.23 from van der Vaart (998)) requires the following conditions:
(i) ((m) is weakly consistent for 8": ((n) Py rr This is assumed, but in general this consistency depends on properties of the function f and the feasible parameter space e See Lemma T for an example.
(ii) En= L(us, 0in) < info Σn L(,, 0) + op(1/n): true by optimality of ((m),
(iii) The loss L(y. 8) is locally Lipschitz near 0 = 0". This holds as long as f(r;0) is locally Lipschitz near = 8r.
(iv) For P-almost all y, the function 0 + L(y, 0) is differentiable at 0=0-
(v) R(O) is twice differentiable at 8* with positive definite Hessian V²R(0*) > 0.
A. W. van der Vaart. Asymptotic Statistics. Cambridge Series in Statistical and Probabilistie Mathematios. Cambridge University Press, 1998.
22

---

Published in Transactions on Machine Learning Research (01/2024)
To verify condition (iv), note that
-
Since is a continuous random variable with a density near /(T;0"), VoL, S(, ((T,0') exists for P-almost every y as long as f(T,0) is also differentiable at = = This because for continuous Y, the T-quantile is unique, and the derivative V,I+ exists with probability The full gradient is:
VoL,(w,
(1=))](((T;;)) v) -1([(T;;))) u)) VoS(T;0").
Therefore, VoL(y, ErlVo (0)] is also well defined and exists for P-almost every y.
To verify condition (v), we explicitly compute the Hessian:
VR(0") = (y(VoL(Y,00)
-BT11:(V since T TL L
- SY) VoJ(T:0)
-
= -EI((EUTT:"")
VR(0') ( +(F((T:)))) -T)V3((:88))
= -CF10(11000))/((:77))01111. =Q>0. 0.
The third equality in the computation of the Hessian above follows by assumption that f is well specified,
since f(r;0") (T), F(f(T;8")) - 0.
With all conditions verified, applying Theorem 5.23 from van der Vaari 1998 we have,
a
vn(((() 0*) Q-
VoL (l)()..
(5)
Focusing on the sum term,
ET (VH
ET 1(11((11))))) - [((1))]((((T8)) -T) s v.))
Vef(T:0))
Vn
Fr S -)
Since 1 E ((T;0)) are IID samples from Bernoulli(7) and T independent of Y, the random variable Ey < ((T:0))) T)Vof (T: has mean and variance
= FT(T(1 - 7
By the central limit theorem, the converges in distribution to (0 V). Combining this with Equation yields,
Vū((ln) 8r) S NN0.0---QQ-117
Applying the delta method with ( as a function for fixed desired quantile т, we have the final result,
(m)) $ N (o. v1
23

---

Published in Transactions on Machine Learning Research (01/2024)
A.3 Sufficient conditions for consistency
Lemma [ IJ () is convex in 0 for all E (0, 1)), O is convex, and O is in the interior of O then ((n)
is weakly consistent, ((m) Eg g".
The result follows from Proposition 7.4 from Hayashi (2011), provided we show that L 0) is also convex.
Proof.
Ty, 0) (T,
= (ay- ((T:0),, (T 1 f(T:0)))
For any fixed scalar value 7, the term -yy - [(T;0) is convex in 0 since it is a linear mapping of f(T; 0). The same applies the term (r- f(; Therefore, max(r(y- f(:: ((: (T- 1) - ff:::;(r:00 0))) as the pointwise maximun of two convex functions, is also convex in 0. Finally, the linearity of the ET operator gives that L(y, 0) is convex in 0.
A.4 A note uniqueness of 0"
Theorem I and Lemma| assume that f(;; 0) well specified with a unique 0* E 2 E 1MTh such that f(r; 0*) (r) for all 10, and that 0* is also the unique minimizer for the risk R(0)) = Ey,T(LT(Y, ((T:())) The property that f(r;0" (*) (T) directly implies that 0* is minimizer of R(0); however, whether 0" is the unique minimizer of R(O) depends on the distribution of T and function class defined by and O
Lemma Suppose (T; is well specified with unique 0" e S RT such that ((T;0") Fr) for all (0, 1). For 0 E 6, define the set Tg = f(; * (T; If for all 0 θ where 0 & 0", the set To has measure greater than O under random variable T, then 0 is also the unique minimizer for the risk R(0) (XT
Proof. Suppose there exists O' where 0' + 0* but 0' also minimizes of the risk R(O). Since O' + 0", the set To, nonempty, and for all E To:
f(T:0) + (T) = 0' mig Ey (Y,
0 (Y; (Y;
Then since Tor has measure greater than 0 under random variable T, this implies
which contradicts that 0' minimizes the risk R(0).
The conditions of uniqueness can be satisfied by a combination of wide enough support T, and a limited enough function class given by f, For example, if T supported on the full interval (0, 1), and f(r; En 0,T is a polynomial function over θ C RTm then any pair + 0' would yield set { f(r; 0) + f(r:0)) with measure greater than 0.
Uniqueness breaks down of the function class is too expressive. For example, if f(T; 0) is infinitely expressive over θ (e.g. θ = IR x (0, 1), and f(T;0) = 0, 1 =)), then it's possible to find a pair 0 + 0' that yield a set (r: (T;;0) # O(T:0')} that contains only a single value of 7 with measure 0.
Uniqueness can also break down if the support of T is too limited. For example, single value To: and f(T: = Em=1, 0,T is a polynomial function, then it's possible to find pair 0 # 8' where [(r:: 0) [To; 0'), and thus (r: f(; 0) + [(r;0)) has measure 0 for T.
is only supported on
24

---

Published in Transactions on Machine Learning Research (01/2024)
A.5 Basic extension to categorical features
Corollary 1 Let X be distribuled uniformly over fivile set of calegorical values X. For nolational convenience, let X' = fl...., m) C N. Suppose that a dataset {(Fs, y)) created by sampling {yz} values IID from the conditional distribulion Y|X = x for each value of x (assuming m divides n), and laking the union of these sets.
Let f X. (0, - R be fully separable over x; that is, let θ O1 x Om, where O₂ represents a copy of a given parameler space On for each value of x E X. Let f lake the form
[(x,T;0)=
Σ ((==xo)g(T: FOEX
for some function g: (0, Oo R.
Suppose f1),T;0) is well specified with unique 0" E θ S R. such that ((x)):00) (x(2) for all T E (0,1) and all x E X, and that H* is also a unique minimizer for the rish R(0) 1 Ezex LY,, f f(z, T:0)] = x]. Suppose the estimator giv) is weakly consistent, ((n) L 0r. Sup- pose further that the function 8 + f(z, r: continuous, locally Lipschitz, and twice differentiable at rr for all 10,1) and all X.
Then for each E X,
- N wo
where,
Q = (2: -B
with Pa (y) being the density of YIX and [(()) O"
Proof- For notational convenience, let X (1....m) m) c N. Define loss function Mx y, 0) Er[Lr(y f(x, T;0)] Let ((m) minimize the empirical risk over samples {ax ys Wim.. ((m) arg ming E Ty,, z4, 0). By assumption of the separability of f and θ by x,
n Σ
((m) e arg min dee n
L(yis xz. 8)
n
arg min
x, =x)
m
n/m
sEN
arg min
Σ m
Er (xx, (t,,T;0) I(x; = x)
n/m
lEe
TEX
arg min lie0
Σ TEX
1
(1IT-(0.9(T:0.))) 1(x; x))
n/m
arg min #.E0.
Er
1(x; x)
m
n/m
x
((n) minimizes an empirical loss over points with x; = x:
Therefore, each component
((m) arg min 1 0,e0, n/m Σ Er(LT(M.
=x).
By assumption, components of the optimum 0 also uniquely minimize each conditional risk
0% = arg min (vx, x]. 0,
25

---

Published in Transactions on Machine Learning Research (01/2024)
Therefore, we may apply Theorem 1 to describe the convergence of each component 0% to its respective optimum 04 for each X, yielding the result.
]
A.6 Additional Examples Illustrating Theorem
Single parameter uniform anchored at O. Suppose the true distribution of interest is Y Unif(O, a). Let f(T; 8) = AT, which correctly parameterizes the true lincar inverse CDF. By Theorem 1, the asymptotic variance of vn((()) - 0*) is
- -
Thus, the asymptotic variance depends the distribution of T through its third and fourth moments, and depends on the distribution of Y through For T Unif(O, 1), the exact asymptotic variance is 0.45a². To estimate a specific To-quantile, we would evaluate the function f(To; ((mi) = ((m) ro, which would have an asymptotic variance of 0.45a27o.
A.7 Connection To Regression On Sample Quantiles
We provide some more intuition here about how the training procedure we describe in this paper learning function f(r) with an expected pinball loss smooths sample quantiles in the case that there are no conditioning x-features. In particular, we can show that our proposed method (except with uniform discrete rather than uniform continuous distribution for r) is cquivalent to regression on the observed sample quantiles with a particular cost. function.
Proposition 1. Given a dataset fu of n elements, construct the dataset (T(i), ()))) pr where y(u repre- i-1
sents the ith order statistic of the data and T(i) = n-1 is the sample quantile that order statistic represents. Then training a model with an expected pinball loss tribution {r is equivalent learning single-variable regression over the data points {T(i), ) The cost function e for the regression at the specific point (r(i), is when f(T(i)) = g0() and otherwise piecewise-linear with slope
! where T is drawn uniformly from the discrete dis-
c(F(F(()) = { (i - - (U(T())
for f(T)) < <(
(1 (1- -()) + Ej=4 1(f(T(1)) > ylrrfor f(T)) > y
Proof. The r-pinball loss can be expressed as T(V ((v-f(T)) f(r))
if f(T) < (1--)(((+) - v) if f())
Lr [(T))
When training with an expected pinball loss over the discrete values of T corresponding to the sample quantiles of the data, this means that the expected pinball loss a function only of [((()), or the quantile predictions of f at those sample quantiles. This is precisely the context of ordinary regression; the loss we face is function only of the values our function takes at the r-coordinates of our data points, which here are set up to be exactly the {r(1))=
Let us consider the pinball loss applied to T(i). To optimize 1 we average the losses across our data points P for this proof, we (equivalently) work with the sums instead. Let our regression-cquivalle cost (((()())) be equal to the relative increase in the T(1) pinball loss we face, summed over the data points, compared to we had predicted f(r(i)) (+).
By the definition of the T(i) pinball loss, decreasing f(r(i) by will mean that we face an additional cost of r(i)e for every data point we are below and moving away from and recover (1-7(i))e for every data point we are above and moving towards. Say we start at f(r(1)) y(0) and choose an c such that f(r(i))-c (i-1) The point ( has (()) 1) data points below it and (1 (1-7(i)((x))))) data points above it by definition. So our marginal loss from decreasing f(r(i)) by e is 7(i)(1 (00(1-7(i))( - 1)e due to the points above y()) while we recover (1 (1-(()))(i)(zz
1)c due to the points below it. These are equal! But we also lose r(1)c because we are now moving away from the point y( itself, making the slope of c(f(((1))) equal to -7(1) in the region between (f(-1) and ( Fach time we pass neighboring data point, we start losing r(i)e from that point
26

---

Published in Transactions on Machine Learning Research (01/2024)
instead of gaining (1 11-7(6) 7())e, increasing the steepness of the slope by 1. This is exactly what we want to show: <(((T(i))) is equal to (i( minus the number of data points smaller than ylu that ((((i)) is smaller than.
The opposite case analyzing what happens when we increase f(r(i)) by is essentially analagous. Increasing [()(i)) by e leads to an additional cost of (1 (1-77 - for every data point we are above and moving away from and decreased cost of T(i)E for every data point we are below and moving towards. This implies the the cost accrued and recovered by the points above and below y( cancel out and we only have the net cost due to moving away from y() itself.
slope of the cost function between yl) and y(t+ is (1- 7(1),
(1-7(t)),
And for the same reason described above, passing data points adds to the slope by 1.
See Figure for some example cost functions.
This cost function is nonstandard in couple ways. First, it varies per point. And second, its steepness is a function of how many neighboring sample data points the prediction f(T) is off by rather than just how far away it is from the sample quantile in an absolute sense.
Still, it otherwise follows the normal properties of a cost function, such as taking the value when the prediction exactly correct and (weakly) rising in convex fashion as the prediction is too high too low. This makes it clear why sufficiently flexible function f(T) will exactly pass through all of the sample quantiles and also why more limited ((() will serve smoother, passing below some sample quantiles and above others.
Worked Example: Consider the dataset {1, 2, 5, 6, 8}. These correspond to the sample quantiles {0, 0.25, 0.5, 0.75, 1}. Figure 7 plots the sample quantiles on a graph, the least squares regression fit,, and the expected pinball loss fit with linear f(r). We can see that both are quite similar. For example, both estimate the median as being lower than the observed sample median. The figure also shows the regressonn-oquivalenn cost function for each of the quantile predictions, if we were to think of our method regression sample quantiles.
A.8 Calibration Error Metric Definition
A predicted conditional inverse CDF f(, T; is calibrated if
P(Y s f(), =T.
A calibration error metric measuring the difference between P(Y S f(X,T:0)) and can be approximated over finite dataset with data points (~i, W: and quantiles Ti.-. Tm- We consider the absolute deviation:
Σ 1 S Ty; (mi;
Tj
(6)
m
Others such as Kuleshov et al. (0018) have considered the squared deviation.
A.9 Hyperparameter Tuning Details
For each real data experiment, tune hyperparameters over validation dataset, where the validation dataset is selected according to Section 5.2 The tuning criterion was the validation pinball loss averaged over E (0.01,..., 0.99). For DLNs, the search spaces for real data experiments were the following:
The number of calibration keypoints for the seccewise-lineat calibration function over were tuned between {10, 20, 50, 100). Note that as this number goes up, you get more model flexibility with respect to - (only), which makes it effectively like just training separate models for each 7.
The number of lattice keypoints for was tuned between 32, 70). That controls the flexibility of the interactions between and the other features .c.
27

---

Published in Transactions on Machine Learning Research (01/2024)
25)
(0.5)
8-
Type
150
19.200
10.50
Cost
0.25
3.50
0.75
1.50
quantie
Figure Given the dataset {1, 2, 5, 6, 8), we show visually how learning f(r) via the expected pinball loss with a feature can be seen as a type of regression the sample quantiles. For simplicity, we take f(r) to be linear, which would be the correct functional form if the data are drawn from uniform distribution. On the left, we plot the sample quantiles and compare the lines learned by direct squared-loss regression on those data points (i.c. (r(1), (² with the line Icarned by our method (i.c. learning f(r) with the expected pinball loss). On the right, we show the equivalent costs imposed on predictions for each sample quantile we reframed training f()) with an expected pinball loss over T = {0, 0.25, 0.5, 0.75, 1} as a regression on sample quantiles. The cost when the prediction exactly the corresponding sample quantile and weakly rising it varies.
Other feature calibration keypoints were tuned between {5, 10, 15, 20}, which are common values for DLNs.
Step sizes were tuned between (0.001, 0.01, 0.05, 0.1}
minibatch sizes were tuned between {1000, 10000}.
Number of steps was tuned between {100, 1000, 10000).
To tune the DLN architecture for experiments with real data, we searched over the following hyperparameters:
Air Quality: number of lattices {4, (..8).nnmma number of features per lattice 6,
Puxzles number of lattices E (8, 16, 32, number of features per lattice =
Traffic: single lattice which includes all features (not ensemble)
Wine: number of lattices (100, 200, 400, 800}, number of fcatures per lattice (2, 428..).
28


Reviews:
[{"review_id": "41sCCHWdlu", "reviewer": "Reviewer_f79v", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "\n**Strengths:** Suggestion of deep lattice networks for quantile estimators to solve crossing quantile issue.\n\n\n\n\n\n**Weakness:**  Writing and organization of the paper", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper provides strong evidence, using both theory and real-world data analysis, that performing quantile regression at different quantiles is more beneficial than doing it at just one specific quantile, even if we are primarily interested in that particular quantile. Moreover, the paper suggests using Deep Lattice networks [1] for estimating quantiles, as these networks address significant problems observed in previous methods [2], such as the issue of quantiles crossing each other. The paper also highlights the practical implications of these problems.\n\n\n**Refs**\n\n[1] S. You, K. Canini, D. Ding, J. Pfeifer, and M. R. Gupta. Deep lattice networks and partial monotonic functions. NeurIPS, 2017.\n\n[2] N. Tagasovska and D. Lopez-Paz. Single-model uncertainties for deep learning. NeurIPS, 2019.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "**Writing:**\n\nWriting feels like not upto the mark. I will just highlight few of them. \n\n\nPage1-2:  \"*That* work is important because it provides simultaneous quantile regression of the entire inverse CDF. However, *that* work left open a couple of important theoretical and practical issues that we address in this paper.\"\n\nPage 11 : \"Lattices build upon a truly *ancient* strategy for approximating functions over a bounded domain.. \" \n\n**Relation with Calibration:**\n\nIt would be beneficial to have a discussion and explore the practical applications of quantile estimators for calibration purposes.\n\n\n\n\n", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "za5GqGYf8E", "reviewer": "Reviewer_ko41", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. The research addresses an interesting problem. \n2. The paper provides several theoretical results. \n3. The research includes a thorough analysis of numerical studies.\n\nWeakness:\nThe revised version effectively addresses the previously mentioned comments.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper investigates the training of a quantile regression model by minimizing the expected pinball loss across all quantiles. The asymptotic convergence rate indicates that minimizing the expected pinball loss can more effectively estimate single quantiles than training with the standard pinball loss for that quantile. Additionally, the authors utilize deep lattice networks to address the issue of crossing quantiles.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The revised version effectively addresses the previously mentioned comments.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "fM4Dd3rYPs", "reviewer": "Reviewer_Nbdr", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "no additional comments.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This is a resubmission of a manuscript I reviewed before. The authors have\naddressed my previous concerns, and I don't have additional major\ncomments. However, the writing, especially some notations, needs some\nclarification.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Both $f(x, \\tau)$ and $f(\\tau; x)$ are used in Section 1.3, which\nmay be a typo.\n\nLater $f(x, \\tau; \\theta)$ is used. I guess $f(x, \\tau)$ means for one sample\nproblem and $f(x, \\tau; \\theta)$ is for the conditional quantile, but it will\nsame readers a lot of efforts if the authors explicitly state the their\nmeanings.\n\nIn the paper, $f(x, \\tau; \\theta)$ is a model, so it is not appropriate to call\nit an estimator in Sections 2, 3.2, and probably some other places. I believe by\n\"i.e. our function class contains the true distribution.\" you mean the working\nmodel contains the true data generating model. Please go other the paper and try\nto improve this type of wage statements.\n\nAre the two models correctly specified in Examples 1 and 2? It would be helpful\nto write the true inverse cdf in the discussion.\n\nThe title for Section 3.3 can be more accurate, so as some statements in the\nabstract. It is about the sampling properties (or asymptotic properties) of the\nestimator. The rate of convergence for the proposed estimator and that of the\nestimator based on a single quantile are the same. It is their asymptotic\nvariances that are different.\n\n", "broader_impact_concerns": "No additional comments", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Reviewers agree that this resubmission addresses the major concerns raised in the previous round, especially around empirical evaluation, assumptions in theory, and discussion of related work. The work makes a valuable contribution to quantile regression, that will be of broad interest to the community. The theory has been improved with clearer assumptions (in particular, a well specified assumption). Some of these assumptions are strong and so empirical evaluations are important. The empirical work is now extensive enough to paint a clear picture that the method offers a valuable tool.
Decision: Accept (poster)