Nonparametric Quantile Regression: Non-Crossing Constraints and Conformal Prediction

1 Introduction

How to assess uncertainty in prediction is a fundamental problem in statistics. Conformal prediction is a general distribution-free methodology for constructing prediction intervals with a guaranteed coverage probability in finite samples

(papadopoulos2002inductive; vovk2005algorithmic). We propose a method for nonparametric estimation of quantile regression functions with the constraint that two quantile functions for different quantile levels do not cross. We then use estimated non-crossing quantile regression functions for constructing conformal prediction intervals.

Since vovk2005algorithmic formally introduced the basic framework of conformal prediction, there has been a number of important advancements on conformal prediction (vovk2012conditional; lei2013distribution; lei2014distribution; lei2018distribution). lei2014distribution and vovk2012conditional showed that the conditional validity for prediction interval with finite length is impossible without any regularity and consistency assumptions on the model and the estimator. zeni2020conformal established that the marginal validity, a conventional coverage guarantee, can be achieved under the assumption that the observations are independent and identically distributed. Recently, several papers have studied the coverage probability, the prediction accuracy in terms of the length of prediction interval, and the computational complexities of conformal prediction using neural networks (Barber2021limits; lei2018distribution; papadopoulos2008inductive).

Earlier works on conformal prediction were based on estimating a conditional mean function and constructing intervals of constant width, assuming homoscedastic errors. Recently, romano2019conformalized

proposed a conformal prediction method based on quantile regression, called conformalized quantile regression. This method is adaptive to data heteroscedasticity and can have varying length across the input space. A similar construction of adaptive and distribution-free prediction intervals using deep neural networks have been considered by

kivaranovic2020adaptive. A comparison study of conformal prediction based on quantile regression with two choices of the conformity scores is given in sesia2020comparison.

Nonetheless, associated with the great flexibility of regression quantiles is the quantile-crossing phenomenon. The quantile crossing problem, due to separate estimation of regression quantile curves at individual quantile levels, has been observed in simple linear quantile regression, and can happen more frequently in multiple regression. Several papers have attempted to deal with the crossing problem. In linear quantile regression, koenker1984note studied a parallel quantile plane approach to avoid the crossing problem. he1997quantile proposed a restricted regression quantile method that avoids quantile crossing while maintaining modeling flexibility. neocleous2008monotonicity

established the asymptotic guarantees that the probability of crossing will tend to zero for the linear interpolation of the Koenker-Bassett linear quantile regression estimator. These papers focused on the linear quantile regression setting.

bondell2010noncrossing proposed a constrained quantile regression to avoid the crossing problem, and considered nonparametric non-crossing quantile regression using smoothing splines with a one-dimensional predictor. However, this approach may not work well with a multi-dimensional predictor. Recently, interesting findings on simultaneous quantile regression that alleviates the crossing quantile problem were reported. tagasovska2019single proposed simultaneous quantile regression to estimate the quantiles by minimizing the pinball loss where the target quantile is randomly sampled in every training iteration. brando2022deep proposed an algorithm for predicting an arbitrary number of quantiles, which ensures the quantile monotonicity by imposing a restriction on the partial derivative of the quantile functions.

In this paper, we make the following methodological and theoretical contributions.

We propose a penalized deep quantile regression approach, in which a novel penalty function based on the rectified linear unit (ReLU) function is proposed to encourage the non-crossing of the estimated quantile regression curves.
Based on the estimated non-crossing quantile regression curves, we study a conformalized quantile regression approach to construct non-crossing conformal prediction intervals, which are fully adaptive to heteroscedasticity and have locally varying length.
We study the properties of the ReLU-penalized nonparametric quantile regression using deep feedforward neural networks. We derive non-asymptotic upper bounds for the excess risk of the non-crossing empirical risk minimizers. Our error bounds achieve optimal minimax rate of convergence, and the prefactor of the error bounds depends polynomially on the dimension of the predictor, instead of exponentially.
We establish theoretical guarantees of valid coverage of the proposed approach to constructing conformal prediction intervals. We also give a non-asymptotic upper bound of the difference between the non-crossing conformal prediction interval and the oracle prediction interval. Extensive numerical studies are conducted to support the theory.

2 Non-crossing nonparametric quantile regression

In this section, we describe the proposed method for nonparametric estimation of quantile curves with non-crossing constraints.

For any given $τ \in (0, 1)$ , the conditional quantile function of $Y$ given $X = x$ is defined by

fτ(x)\coloneqqQτ(Y|X=x)=inf{y∈R:F(y|X=x)≥τ},x∈X,

(1)

where $Y \in R$ is a response, $X \in X \subseteq R^{d}$ is a $d$

-dimensional vector of covariates, and

F (\cdot | X = x)

is the conditional distribution function (c.d.f) of

Y

given

X = x

. It holds that

P (Y \leq f_{τ} (X) ∣ X = x) = τ

, where

P (\cdot | X = x)

is the conditional probability measure of

Y

given

X = x

. By definition (1), an inherent constraint of the conditional quantile curves is the monotonicity property: for any

0 < τ_{1} < τ_{2} < 1

, it holds that

f_{τ_{1}} (x) \leq f_{τ_{2}} (x), x \in X .

(2)

Therefore, estimated conditional quantile curves should also satisfy this property, otherwise, it would be difficult to interpret the estimated quantiles.

Quantile regression is based on the check loss function defined by

ρ_{τ} (u) = u {τ - 1 (u \leq 0)}, u \in R,

(3)

where $1 (\cdot)$ is the indicator function. The target conditional quantile function $f_{τ}$ is the minimizer of $E ρ_{τ} (Y - f (X))$ over all measurable function $f$ (koenker_2005). In applications, only a finite random sample ${(X_{i}, Y_{i})}_{i = 1}^{n}$ is available. The quantile regression estimator for a given $τ$ is

{^f}_{τ} (x) \in arg min f \in F_{n} \frac{1}{n} n \sum i = 1 ρ_{τ} (Y_{i} - f (X_{i})),

(4)

where $F_{n}$ is a class of functions which may depend on the sample size $n$ .

For two quantile levels $0 < τ_{1} < τ_{2} < 1$ , we can obtain the estimated quantile curves ${^f}_{τ_{1}}$ and ${^f}_{τ_{2}}$ by using (4) separately for $τ_{1}$ and $τ_{2}$ . However, such estimated quantile curves may not satisfy the monotonicity constraint (2), that is, there may exist $x \in X$ for which ${^f}_{τ_{1}} (x) > {^f}_{τ_{2}} (x) .$ Below, we propose a penalized method to mitigate this problem.

2.1 Non-crossing quantile regression via ReLU penalty

In this subsection, we propose a penalized quantile regression framework to estimate quantile curves that can avoid the crossing problem. We first introduce a ReLU-based penalty function to enforce the non-crossing constraint $f_{2} (x) \geq f_{1} (x)$ in quantile regression. The ReLU penalty is defined as

V(f1,f2;x)\coloneqqmax{f1(x)−f2(x),0}.

This penalty function encourages $f_{2} (x) \geq f_{1} (x)$ when combined with the quantile loss function. At the population level, for $0 < τ_{1} < τ_{2} < 1$ , the expected penalized quantile loss function for $f_{1}, f_{2}$ is

R^{τ_{1}, τ_{2}} (f_{1}, f_{2}) = E [ρ_{τ_{1}} {Y - f_{1} (X)} + ρ_{τ_{2}} {Y - f_{2} (X)}] + λ E V (f_{1}, f_{2}; X),

(5)

where $λ \geq 0$ is a tuning parameter. Note that, with $f_{τ_{1}}, f_{τ_{2}}$ defined in (1), $E V (f_{τ_{1}}, f_{τ_{2}}; X) = 0$ for $τ_{1} < τ_{2}$ . Moreover, as mentioned earlier, $f_{τ_{i}}$ is the minimizer of the expected check loss $E ρ_{τ_{i}} (Y - f (X))$ over measurable function $f$ for $i = 1, 2$ respectively. The following lemma establishes the identifiability of the quantile functions through the loss function (5). This is the basis of the proposed method.

Lemma 2.1

Suppose $Z$

is a random variable with c.d.f.

F (\cdot)

. For a given

τ \in (0, 1)

, any element in

{z : F (z) = τ}

is a minimizer of the expected check loss function

E ρ_{τ} (Z - t)

. Moreover, the pair of true conditional quantile functions

(f_{τ_{1}}, f_{τ_{2}})

is the minimizer of the loss function (5), that is,

(f_{τ_{1}}, f_{τ_{2}}) = arg min f_{1}, f_{2} \in F_{0} R^{τ_{1}, τ_{2}} (f_{1}, f_{2}),

where $F_{0}$ is a class of measurable functions that contains the true conditional quantile functions.

Lemma 2.1 shows that the ReLU penalty does not introduce any bias at the population level. This is a special property of this penalty function, which differs from the usual penalty function such as the ridge penalty used in penalized regression.

When only a random sample ${X_{i}, Y_{i}}_{i = 1}^{n}$ is available, we use the empirical loss function

R_{n}^{τ_{1}, τ_{2}} (f_{1}, f_{2}) := n \sum i = 1 [ρ_{τ_{1}} {Y_{i} - f_{1} (X_{i})} + ρ_{τ_{2}} {Y_{i} - f_{2} (x_{i})}] + λ n \sum i = 1 V (f_{1}, f_{2}; X_{i}),

(6)

and define

({^f}_{1}, {^f}_{2}) \in arg min f_{1}, f_{2} \in F R_{n}^{τ_{1}, τ_{2}} (f_{1}, f_{2}),

(7)

where $F$ is a class of functions that approximate $F_{0}$ . In (6), a positive value of $f_{1} (X_{i}) - f_{2} (X_{i})$ (i.e., the quantile curves cross at $X_{i}$ ) will be penalized with the penalty parameter $λ$ . On the other hand, a negative value of $f_{1} (X_{i}) - f_{2} (X_{i})$ (i.e., the quantile curves do not cross at $X_{i}$ ) will not incur any penalty. Therefore, with a sufficiently large penalty parameter $λ$ , quantile crossing will be prevented.

Throughout the paper, we choose the function class $F$ to be a function class of feedforward neural networks. We note that other approximation classes can also be used. However, an important advantage of neural network functions is that they are effective in approximating smooth functions in $R^{d}$ , see, for example, jiao2021deep and the references therein. A detailed description of neural network functions will be given in Section 2.2. Figure 1 previews non-crossing curve estimation via a toy example, which shows the comparison between the proposed method in (7) (with penalty) and the deep quantile estimation studied in shen2021deep (without penalty); see Section B.1.2 in Appendix for more details.

Figure 1: A toy example shows that NC-CQR can avoid quantile crossing. The quantile curves on the left panel are estimated by NC-CQR, and those on the right are estimated by the deep quantile regressionshen2021deep. The yellow color indicates quantile crossing.

Remark 2.1

The rectified linear unit (ReLU) function $σ (x) = max (0, x)$ is a piecewise-linear function. An important advantage of the ReLU penalty $max {f_{1} - f_{2}, 0}$ is that it is convex with respect to $f_{1}$ and $f_{2}$ , thus the loss function in (6) is also convex with respect to $f_{1}$ and $f_{2}$ .

Remark 2.2

The tuning parameter $λ$ controls the amount of penalty on crossing quantile curves. As shown in Section 4, a bigger $λ$ leads to a larger estimation bias. Therefore, there is a trade-off between the crossing and the accuracy of prediction intervals. Both can be achieved with proper choice of $λ$ . In the appendix, we propose a cross-validation method to select $λ$ .

Remark 2.3

bondell2010noncrossing proposed a non-crossing nonparametric quantile regression using smoothing splines with non-crossing constraints. For two given quantile levels $0 < τ_{1} < τ_{2} < 1$ , they proposed to estimate the quantile curves $f_{τ_{1}}$ and $f_{τ_{2}}$ by minimizing the following constrained loss function

	$min f_{1}, f_{2}$	$2 \sum t = 1 n \sum i = 1 ρ_{τ_{t}} (Y_{i} - f_{t} (X_{i})) + 2 \sum t = 1 λ_{τ_{t}} V^{*} (f_{t}^{'})$
		$subject to f_{2} (x) \geq f_{1} (x) for all% x \in X,$		(8)

where $f^{'}$ denotes the derivative of a function $f$ , and $V^{*} (f^{'}) = \int_{0}^{1} | f^{''} (t) | d t$ is the total variation penalty to guarantee smoothness. Such a spline-based method works well in the one-dimensional setting, however, it is difficult to apply this approach to multi-dimensional problems.

2.2 ReLU Feedforward neural networks

For the estimation of conditional quantile functions, we choose the function class $F$ in (7) to be $F_{D, W, U, S, B}$ , a class of feedforward neural networks $f_{ϕ} :$ $R^{d} \to R$ with parameter $ϕ$ , depth $D$ , width $W$ , size $S$

, number of neurons

U

and

f_{ϕ}

satisfying

{∥ ∥ f_{ϕ} ∥ ∥}_{\infty} \leq B

for some positive constant

B

, where

∥ f ∥_{\infty}

is the supreme norm of a function

f : R^{d} \to R

. Note that the network parameters may depend on the sample size

n

, but this dependence is omitted for notational simplicity. Such a network

f_{ϕ}

has

D

hidden layers and

(D + 2)

layers in total. We use a

(D + 2)

vector

(w_{0}, w_{1}, \dots, w_{D}, w_{D + 1})

to describe the width of each layer; particularly in nonparametric regression problems,

w_{0} = d

is the dimension of the input and

w_{D + 1} = 1

is the dimension of the response. The width

W

is defined as the maximum width of hidden layers, i.e.,

W = max {w_{1}, \dots, w_{D}};

the size

S

is defined as the total number of parameters in the network

f_{ϕ}

, i.e.,

S = \sum_{i = 0}^{D} w_{i + 1} (w_{i} + 1);

the number of neurons

U

is defined as the number of computational units in hidden layers, i.e.,

U = \sum_{i = 1}^{D} w_{i} .

For an MLP

F_{D, U, W, S, B}

, its size

S

satisfies

max {W, D} \leq S \leq W (D + 1) + (W^{2} + W) (D - 1) + W + 1.

From now on, we write $F_{D, W, U, S, B}$ as $F_{ϕ}$ for short. Then, at the population level, the non-crossing nonparametric quantile estimation is to find a pair of measurable functions $(f_{ϕ 1}^{*}, f_{ϕ 2}^{*}) \in F_{ϕ}$ satisfying

(f_{ϕ 1}^{*}, f_{ϕ 2}^{*}) := arg min f_{1}, f_{2} \in F_{ϕ} R^{τ_{1}, τ_{2}} (f_{1}, f_{2}) .

(9)

3 Non-crossing quantile regression for conformal prediction

Now suppose we have a new observation $X_{n + 1}$ . We are interested in predicting the corresponding unknown value of $Y_{n + 1}$ . Our goal is to construct a distribution-free prediction interval $C (X_{n + 1}) \subseteq R$ with a coverage probability satisfying

P {Y_{n + 1} \in C (X_{n + 1})} \geq 1 - α

(10)

for any joint distribution

P_{X Y}

and any sample size

n

, where

0 < α < 1

is often called a miscoverage rate. We refer to such a coverage stated in (10) as marginal coverage.

First, we define an oracle prediction band based on the conditional quantile functions. For a pre-specified miscoverage rate $α$ , we consider the lower and upper quantiles levels such as $τ_{1} = α / 2$ and $τ_{2} = 1 - α / 2$ . Then, a conditional prediction interval for $Y$ given $X = x$ with a nominal miscoverage rate $α$ is

C^{*} (x) = [f_{τ_{1}} (x), f_{τ_{2}} (x)],

(11)

where $f_{τ_{1}}$ and $f_{τ_{2}}$ are the conditional quantile functions defined in (1) for quantile levels $τ_{1}$ and $τ_{2},$ respectively. Such a prediction interval with true quantile functions is ideal but cannot be constructed, only a corresponding empirical version can be estimated based on data in practice.

Next, we use the split conformal method vovk2005algorithmic for constructing non-crossing conformal intervals. We split the observations ${(X_{i}, Y_{i})}_{i = 1}^{n}$ into two disjoint sets: a training set $I_{1}$ and a calibration set $I_{2}$ . Non-crossing deep neural estimators of $f_{τ_{1}}$ and $f_{τ_{2}}$ based on the training set $I_{1}$ are given by

({^f}_{1}, {^f}_{2}) \in

arg min f_{1}, f_{2} \in F_{ϕ} R_{I_{1}}^{τ_{1}, τ_{2}} (f_{1}, f_{2}),

(12)

where

Rτ1,τ2I1(f1,f2)\coloneqq∑i∈I1[ρτ1{yi−f1(xi)}+ρτ2{yi−f2(xi)}]+λ∑i=I1V(f1,f2;xi)

and $λ \geq 0$ is a tuning parameter. A key step is to compute the conformity score based on the calibration set $I_{2}$ romano2019conformalized, which is defined as

Ei:=max{^f1(Xi)−Yi,Yi−^f2(Xi)},   i∈I2.

(13)

Let $E = {E_{i}}_{i \in I_{2}}$ . The conformity scores in $E$ quantify the errors incurred by the plug-in prediction interval $[{^f}_{1} (x), {^f}_{2} (x)]$ evaluated on the calibration set $I_{2}$ , and can account for undercoverage and overcoverage.

Finally, for a new input data point $X_{n + 1}$ and $α \in (0, 0.5)$ , the $100 (1 - α) %$ prediction interval for $Y_{n + 1}$ is defined as

(14)

where $Q_{1 - α} (E, I_{2})$ is the $(1 - α) (1 + 1 / | I_{2} |)$ -th empirical quantile of ${E_{i} : i \in I_{2}}$ , namely, the $⌈ (| I_{2} | + 1) (1 + α) ⌉$ -th smallest value in ${E_{i} : i \in I_{2}}$ , and $⌈ a ⌉$ denotes the smallest integer no less than $a$ . Here, $| I |$ is the cardinality of a set $I$ . The empirical quantile $Q_{1 - α} (E, I_{2})$ is a data-driven quantity, which conformalizes the plug-in prediction interval.

In contrast to the conformalized quantile regression procedure in romano2019conformalized, our method produces conformal prediction intervals that avoid the quantile crossing problem. We refer to the proposed non-crossing conformalized quantile regression as NC-CQR, and the corresponding conformal interval (14) as NC-CQR interval.

Remark 3.1

The usual linear quantile regression (QR) model assumes that, for a given $τ \in (0, 1)$ ,

Q_{τ} (Y | X = x) = a_{τ} + b_{τ}^{⊤} x,

(15)

where $a_{τ} \in R$ and $b_{τ} \in R^{d}$ are the intercept and slope parameters. Following romano2019conformalized, a conformal interval based on linear quantile regression can be constructed. Specifically, by splitting the observations ${(X_{i}, Y_{i})}_{i = 1}^{n}$ into two disjoint subsets: a training set $I_{1}$ and a calibration set $I_{2}$ , we can fit model (15) on the training set $I_{1}$ and obtain the estimators for $a_{τ}$ and $b_{τ}$ for a given $τ \in (0, 1)$ , denoted by $({^a}_{τ}, {^b}_{τ})$ . Under model (15), the conformal interval with miscoverage rate $α$ is given by

{^C}_{q r} (x) = [{^f}_{l o} (x) - {^Q}_{q r}, {^f}_{h i} (x) + {^Q}_{q r}],

(16)

where ${^f}_{l o} (x) := x^{⊤} {^b}_{τ_{l o}} + {^a}_{τ_{l o}}$ , ${^f}_{h i} (x) := x^{⊤} {^b}_{τ_{h i}} + {^a}_{τ_{h i}}$ with $τ_{l o} = α / 2$ and $τ_{h i} = 1 - α / 2$ , and ${^Q}_{q r}$ is the $⌈ (| I_{2} | + 1) (1 + α) ⌉$ -th smallest value among the conformity scores $Eqri=max{^flo(x)−yi,yi−^fhi(xi)},i∈I2$ .

We summarize the implementation of NC-CQR interval construction in the following algorithm.

Algorithm Computation of non-crossing conformalized prediction intervals

Input: Observations $S = {(X_{i}, Y_{i})}_{i = 1}^{n}$ and miscoverage level $α \in (0, 0.5) .$

Output: $C (x) = [{^f}_{1} (x) - Q_{1 - α} (E, I_{2}), {^f}_{2} (x) + Q_{1 - α} (E, I_{2})] .$

Steps:

Split $S$ into a training set $S_{1} = {(X_{i}, Y_{i})}_{i \in I_{1}}$ and a calibration set $S_{2} = {(X_{i}, Y_{i})}_{i \in I_{2}}$ .
Fit $S_{1}$ to (12) to obtain $({^f}_{1}, {^f}_{2})$ for $τ_{1} = 1 - α / 2$ and $τ_{2} = α / 2$ .
Compute the conformity score $Ei:=max{^f1(Xi)−Yi,Yi−^f2(Xi)}$ , $i \in I_{2} .$
Find $Q_{1 - α} (E, I_{2})$ , the $⌈ (| I_{2} | + 1) (1 + α) ⌉$ -th smallest value of $E_{i}, i \in I_{2} .$
5. Compute the prediction band $C (x)$ according to (14).

4 Theoretical properties

In this section, we study the theoretical properties of the proposed NC-CQR method. We evaluate NC-CQR using the following two criteria:

Validity: Under proper conditions, a conformal prediction interval $^C$ satisfies that

$P (Y_{n + 1} \in^C (X_{n + 1})) \geq 1 - α .$ (17)
Accuracy: If the validity requirement (17) is satisfied, a conformal prediction interval should be as narrow as possible.

The validity requirement (17) is evaluated based on the finite-sample marginal coverage in (10), which holds in the sense of averaging over all possible test values of $X_{n + 1}$ . The accuracy of a prediction interval is usually measured by the discrepancy defined in (19) between the lengths of the prediction interval and the oracle one.

We assume that the target conditional quantile function $f_{τ}$ defined in (1) is a $β$ -Hölder smooth function with $β > 0$ as stated in condition (C3) below. Let $β = s + r > 0$ , $r \in (0, 1]$ and $s = ⌊ β ⌋ \in N_{0}$ , where $⌊ β ⌋$ denotes the largest integer strictly smaller than $β$ and $N_{0}$ denotes the set of non-negative integers. For a finite constant $B_{0} > 0$ , the Hölder class of functions $f : [0, 1]^{d} \to R$ is defined as

H^{β} ([0, 1]^{d}, B_{0}) = {f : max ∥ α ∥_{1} \leq s ∥ \partial^{α} f ∥_{\infty} \leq B_{0}, max ∥ α ∥_{1} = s sup x \neq y \frac{| \partial^{α} f (y) - \partial^{α} f (x) |}{∥ x - y ∥_{2}^{r}} \leq B_{0}},

(18)

where $\partial^{α} = \partial^{α_{1}} \dots \partial^{α_{d}}$ with $α = (α_{1}, \dots, α_{d})^{⊤} \in N_{0}^{d}$ and $∥ α ∥_{1} = \sum_{i = 1}^{d} α_{i}$ . We assume the following conditions.

(C1) The observations ${(X_{i}, Y_{i})}_{i = 1}^{n}$ are i.i.d. copies of $(X, Y)$ .
(C2) (i) The support of the predictor vector $X$ is a bounded compact set in $R^{d}$ , and without loss of generality we let $X = [0, 1]^{d};$ (ii) Let $ν$ be the probability measure of $X$ . The probability measure $ν$ is absolutely continuous with respect to the Lebesgue measure.
(C3) For any fixed $τ \in (0, 1)$ , the target conditional quantile function $f_{τ}$ defined in (1) is a Hölder smooth function of order $β$ and a finite constant $B_{0}$ .
(C4) There exist constants $γ > 0$ and $κ_{τ} > 0$ such that for any $| δ | \leq γ$ and any $τ \in (0, 1)$ ,

$| F (f_{τ} (x) + δ ∣ x) - F (f_{τ} (x) ∣ x) | \geq κ_{τ} | δ |$

for all $x \in X$ up to a $ν$ -negligible set, where $F (\cdot ∣ x)$
is the conditional cumulative distribution function of
$Y$ given $X = x$ .

Condition (C1) is a basic assumption in conformal inference. The boundedness support assumption in Condition (C2) is made for technical convenience in the proof for deep neural estimation. Condition (C3) is a regular smoothness assumption on the target regression functions so that whose approximation result using deep neural networks can be obtained. Condition (C4) is imposed for self-calibration of the resulting neural estimator.

Theorem 4.1

Under Condition (C1), for a new i.i.d pair $(X_{n + 1}, Y_{n + 1})$ , the proposed NC-CQR interval $^C (X_{n + 1})$ satisfies $P (Y_{n + 1} \in^C (X_{n + 1})) \geq 1 - α .$

Theorem 4.1 shows that the proposed NC-CQR interval enjoys a rigorous coverage guarantee. The proof of this theorem is given in the Appendix.

We now quantify the accuracy of the prediction interval in terms of the difference between the NC-CQR interval $^C (X)$ and the oracle $C^{*} (X)$ in (11) on the support of $X .$ Let $ν$ be the probability measure of $X$ defined in Condition (C2) and define $∥f−g∥Lp(ν)\coloneqq{E|f(X)−g(X)|p}1/p$ for $p \geq 1$ . The difference between the lengths of $^C (X)$ and $C^{*} (X)$ is given by

ΔX(C∗,^C)\coloneqq

E [∥ ({^f}_{2} + Q_{1 - α} (E, I_{2})) - ({^f}_{1} - Q_{1 - α} (E, I_{2})) ∥_{L^{2} (ν)} - ∥ f_{τ_{2}} - f_{τ_{1}} ∥_{L^{2} (ν)}] .

(19)

By the triangle inequality, we have $Δ_{X} (C^{*},^C) \leq I_{1} + I_{2} + I_{3},$ where

I_{1} := E ∥ {^f}_{2} - f_{τ_{2}} ∥_{L^{2} (ν)}, I_{2} := E ∥ {^f}_{1} - f_{τ_{1}} ∥_{L^{2} (ν)} and I_{3} := 2 E ∥ Q_{1 - α} (E, I_{2}) ∥_{L^{2} (ν)} .

To bound $I_{1}$ and $I_{2}$ , we need to bound the error $∥^f - f_{τ} ∥_{L^{2} (ν)} .$ To this end, we first derive bounds for the the excess risk of $({^f}_{1}, {^f}_{2})$ defined as $E{Rτ1,τ2(^f1,^f2)−Rτ1,τ2(fτ1,fτ2)},$ where $R^{τ_{1}, τ_{2}}$ is defined in (5). Without loss of generality, we assume that $| I_{1} | = ⌈ n / 2 ⌉$ and $| I_{2} | = ⌊ n / 2 ⌋$ , where $⌈ a ⌉$ denotes the smallest integer no less than $a$ and $⌊ a ⌋$ denotes the largest integer no greater than $a$ .

Theorem 4.2

Letting $N = 1$ and $L = ⌊ n^{d / 2 (d + β)} ⌋$ , then the width, depth and size of the neural network satisfy

	$W$	$= 114 (⌊ β ⌋ + 1)^{2} d^{⌊ β ⌋ + 1}, D = 21 (⌊ β ⌋ + 1)^{2} ⌈ n^{d / 2 (d + β)} {log}_{2} (8 n^{d / 2 (d + β)}) ⌉,$
	$S$	$= O ((⌊ β ⌋ + 1)^{6} d^{2 ⌊ β ⌋ + 2} ⌈ n^{d / 2 (d + β)} (log n) ⌉) .$

Suppose that Conditions (C1)-(C3) hold. If $λ = log n$ , the non-asymptotic error bound of the excess risk satisfies

E{Rτ1,τ2(^f1,^f2)−Rτ1,τ2(fτ1,fτ2)}≤c1(⌊β⌋+1)8d2(⌊β⌋+1)(logn)5n−βd+β,

where $c_{1} > 0$ is a universal constant independent of $n, d, τ, B, S, W$ and $D$ and $a \lor b := max {a, b}$ .

The convergence rate of the excess risk is $n^{- β / (d + β)}$ up to a logarithmic factor. The next theorem gives an upper bound for the accuracy of the proposed prediction interval defined in (19).

Theorem 4.3

(Non-asymptotic upper bound for prediction accuracy) Suppose that Conditions (C1)-(C4) hold. Let $F_{ϕ} = F_{D, W, U, S, B}$ be a class of ReLU activated feedforward neural networks with width, depth specified as in Theorem 4.2 and let $({^f}_{1}, {^f}_{2}) \in arg {min}_{f_{1}, f_{2} \in F_{ϕ}} R_{n}^{τ_{1}, τ_{2}} (f_{1}, f_{2})$ be the empirical risk minimizer over $F_{ϕ}$ . Then, there exists a constant $N > 0$ , for $n > N$ ,

Δ_{X} (C^{*},^C) \leq c (⌊ β ⌋ + 1)^{8} d^{2 (⌊ β ⌋ + 1)} (log n)^{4} n^{- β / (d + β)},

where $c > 0$ is a constant independent of $n, d, τ, B, S, W$ and $D$ .

Theorem 4.3 gives an upper bound for the difference between the lengths of our proposed prediction interval and the oracle interval. With properly-selected neural network parameters, the oracle band can be consistently estimated by the NC-CQR band.

Finally, we consider the conformal interval based on linear quantile regression defined in (16).

Corollary 4.1

Suppose that at least one of $f_{τ_{1}}$ and $f_{τ_{2}}$ is a non-linear function on a subset of $X$ with non-zero measure. Then for any sample size $n$ , the accuracy of the conformal band ${^C}_{q r}$ defined in (16) is strictly worse than that of the oracle conformal band $C^{*}$ , i.e., there exists an $ϵ > 0$ such that $Δ_{X} (C^{*}, {^C}_{q r}) > ϵ .$

According to Corollary 4.1, a conformal interval based on the linear quantile regression will not reach the oracle accuracy in the presence of nonlinearity.

5 Numerical studies

5.1 Synthetic data

We first simulate a synthetic dataset with a $d$ -dimensional feature $X$ and a continuous response $Y$ from the distributions defined in Section B.1.3. Our method is applied to $5, 000$ independent observations from this distribution, using $2, 000$ of them to train the deep quantile regression estimators and $2, 000$ of them for calibration. The remaining data is for testing. We consider different dimensions $d = 5, 10, 15, 20, 25$ to investigate how the dimensionality of the input affects the overall performance under multivariate input settings. The result is shown in Figure 2. In Figure 1(c), when dimension increases, NC-CQR method performs better than CQR in terms of smaller crossing rate. It shows that the NC-CQR method can mitigate the crossing problem in quantile regression. We also give a 3D visualization of the conformal intervals of our proposed NC-CQR estimation and that of the CQR method when $d = 2$ in Figure 2. One can see from Figure 1(a) that the conformal interval by our proposed NC-CQR method does not have any crossing, while in Figure 1(b), the red region indicates that the lower bound is larger than the upper bound of the interval. More details of the results are given in Section B.1.3.

(a) The crossing rates of two methods as dimension increases.

Our next synthetic example illustrates the adaptive property of the proposed NC-CQR prediction interval. We consider the following model with a discontinuous regression function:

Y = 5 X 1 (X \leq 0.5) + 5 (X - 1) 1 (X > 0.5) + ε .

(20)

A 90% NC-CQR interval is shown in color blue. The green lines represent a 90% conformal prediction intervals based on the standard linear quantile regression. As can be seen from the plots, the NC-CQR intervals automatically adapt to the shape of the regression function and the heteroscedasticity of the model error. Additional examples demonstrating the advantages of NC-CQR prediction intervals are given in the Appendix.

Figure 3: Prediction intervals based on NC-CQR and linear QR for model (20)

5.2 Real data examples

We compute the NC-CQR prediction intervals for the following datasets: bike sharing¹¹1https://archive.ics.uci.edu/ml/datasets/bike+sharing+dataset, house price²²2https://www.kaggle.com/datasets/harlfoxem/housesalesprediction and the air foil³³3https://archive.ics.uci.edu/ml/datasets/airfoil+self-noise data sets. To examine the performance of the ReLU penalty, we compare the NC-CQR with CQR. We compute the conformal prediction intervals based on these two methods using the same deep quantile regression estimator. Their performances are evaluated as in Section 5.1. We subsample $40 %$ of data for training, $40 %$

for calibration, and the remaining data is used for testing. All features are standardized to have 0 mean and unit variance. The nominal coverage rate is set to be

80 %

. Figure 4 shows that the proposed NC-CQR method can mitigate the crossing problem encountered in the CQR estimation. Additional details of the result are given in Section B.2 in the Appendix.

6 Concluding remarks

We have proposed NC-CQR, a penalized deep quantile regression method which avoids the quantile crossing problem via a ReLU penalty. We have derived non-asymptotic upper bounds for the excess risk of the non-crossing empirical risk minimizers. Our error bounds achieve optimal minimax rate of convergence, the with the prefactor of the error bounds depending polynomially on the dimension of the predictor, instead of exponentially. The proposed ReLU penalty for monotonic constraints is applicable to problems with different loss functions. It is applicable to other nonparametric estimation method, including smoothing splines, kernel method and neural network. Therefore, the proposed ReLU penalty method may be of independent interest. Based on the estimated non-crossing conditional quantile curves, we construct conformal prediction intervals that are fully adaptive to heterogeneity. We have also provided conditions under which the proposed NC-CQR conformal prediction bands are valid with the correct coverage probability and achieve optimal accuracy in terms of the width of the prediction band.

A main limitation of the proposed method and the theoretical properties is that they rely on the i.i.d. assumption for the data. There are applications where observations are not independent (e.g., time series data) and there may be distribution shift. It would be interesting to extend the proposed method to deal with such non-i.i.d. settings and establish the corresponding theoretical properties.

Appendix A Appendix: Proofs

In this appendix, we prove the theoretical results stated in Section 4. First, we prove Lemma 2.1. For convenience, we restate this lemma below.

Lemma A.1

Suppose $Z$ is a random variable with c.d.f. $F (\cdot)$ . For a given $τ \in (0, 1)$ , any element in ${z : F (z) = τ}$ is a minimizer of the expected check loss function $E ρ_{τ} (Z - t)$ . Moreover, the pair of the true conditional quantile functions $(f_{τ_{1}}, f_{τ_{2}})$ is the minimizer of the loss function (5), that is,

(f_{τ_{1}}, f_{τ_{2}}) = arg min f_{1}, f_{2} \in F R^{τ_{1}, τ_{2}} (f_{1}, f_{2}),

where $F$ is a class of measurable functions that contains the true conditional quantile functions.

Proof[of Lemma A.1] First, we write

min t E ρ_{τ} (Z - t) = min {\int_{- \infty}^{t} (τ - 1) (z - t) d F (z) + \int_{t}^{\infty} τ (z - t) d F (z)} .

Taking derivative w.r.t $t$ , we obtain

	$0$	$= (1 - τ) \int_{- \infty}^{t} d F (z) - τ \int_{t}^{\infty} d F (z)$
		$= (1 - τ) F (t) - τ (1 - F (t)) = F (t) - τ .$

Since $F (\cdot)$ is a c.d.f., it is monotonic and thus any element of ${z : F (z) = τ}$ minimizes the expected loss $E ρ_{τ} (Z - t)$ . For functions $f_{1}, f_{2}$ , recall that the loss function (5) in the main context is

R^{τ_{1}, τ_{2}} (f_{1}, f_{2}) = E_{Z} [ρ_{τ_{1}} {Y - f_{1} (X)} + ρ_{τ_{2}} {Y - f_{2} (X)}] + λ E_{Z} V (f_{1}, f_{2}; X) .

By the definition of $f_{τ_{k}}$ in (1) in the main context, $f_{τ_{k}}$ satisfies $inf {y : F (y | X = x) \geq τ_{k}}$ for $k = 1, 2$ . Then $f_{τ_{k}}$ minimizes the expected loss $E ρ_{τ_{k}} (Y - f_{k} (x))$ for $k = 1, 2$ . Since the true quantile function $(f_{τ_{1}}, f_{τ_{2}})$ satisfies the monotonicity requirement (2) in the main context, then $E_{Z} V (f_{τ_{1}}, f_{τ_{2}}; X) = 0$ . Therefore, $(f_{τ_{1}}, f_{τ_{2}})$ is the minimizer of the loss function (5) in the main context.

Lemma A.2

Suppose that

(S21)

for some sequence $η_{n}$ and $ξ_{n}$ such that $η_{n} \to 0$ and $ξ_{n} \to 0$ when $n \to \infty$ . Define

An:={x:∣∣^f1(x)−fτ1(x)∣∣>η1/3n}∪{x:∣∣^f2(x)−fτ2(x)∣∣>η1/3n}.

Then, under Conditions (C1)-(C4), for any $X$ independent of $I_{1}$ ,

P [X \in A_{n} | I_{1}] \leq η_{n}^{1 / 3} + ξ_{n} .

Furthermore, if the calibration set $I_{2}$ is partitioned into

I_{2, a} := {i \in {n + 1, \dots, 2 n} : X_{i} \in A_{n}}, I_{2, b} := {i \in {n + 1, \dots, 2 n} : X_{i} \in A_{n}^{c}},

then for any constant $c > 0$ , we have

P [| I_{2, a} | \geq n (η_{n}^{1 / 3} + ξ_{n}) + c \sqrt{n log n}] \leq n^{- 2 c^{2}} .

Proof[of Lemma A.2] By Markov inequality, we have

		$P [X \in A_{n} ∣ I_{1}]$
	$=$
	$\leq$
	$\leq$	$η_{n}^{- 2 / 3} E_{X} [{({^f}_{1} (X) - f_{τ_{1}} (X))}_{1}^{2} ∣ I_{1}] + η_{n}^{- 2 / 3} E_{X} [{({^f}_{2} (X) - f_{τ_{2}} (X))}_{2}^{2} ∣ I_{1}]$
	$\leq$	$η_{n}^{1 / 3} + ξ_{n} .$

When the calibration set $I_{2}$ is partitioned into

I_{2, a} := {i \in I_{2} : X_{i} \in A_{n}}, I_{2, b} := {i \in I_{2} : X_{i} \in A_{n}^{c}},

it follows from Hoeffding’s inequality that, for any $t > 0$ ,

		$P [\| I_{2, a} \| \geq n (η_{n}^{1 / 3} + ξ_{n}) + t]$
	$\leq$	$P [\| I_{2, a} \| \geq n P [X \in A_{n}] + t]$
	$\leq$	$P [\frac{1}{n} 2 n \sum i = n + 1 1 [X_{i} \in A_{n}] \geq P [X_{i} \in A_{n}] + \frac{t}{n}]$
	$\leq$	$P [\frac{1}{n} 2 n \sum i = n + 1 1 [X_{i} \in A_{n}] - P [X_{i} \in A_{n}] \geq \frac{t}{n}]$
	$\leq$	$exp (- \frac{2 t^{2}}{n}) .$

Letting $t = c \sqrt{n log n}$ , the proof of the second inequality in Lemma (A.2) is completed.

Proof[of Theorem 4.1] The proof of Theorem 4.1 is similar to that of Theorem 1 in romano2019conformalized. Let $E_{n + 1}$ be the conformity score

En+1:=max{^f1(Xn+1)−Yn+1,Yn+1−^f2(Xn+1)}

at the test point $(X_{n + 1}, Y_{n + 1}) .$ Recall that the prediction interval

^C (X_{n + 1}) = [{^f}_{1} (X_{n + 1}) - Q_{1 - α} (E, I_{2}), {^f}_{2} (X_{n + 1}) + Q_{1 - α} (E, I_{2})] .

By the construction of the prediction interval, $Y_{n + 1} \in$

Nonparametric Quantile Regression: Non-Crossing Constraints and Conformal Prediction

Estimation of Non-Crossing Quantile Regression Process with Deep ReQU Neural Networks

Sparse Quantile Regression

Quantile index regression

Deep Quantile Regression: Mitigating the Curse of Dimensionality Through Composition

Solution to the Non-Monotonicity and Crossing Problems in Quantile Regression

Prediction intervals with controlled length in the heteroscedastic Gaussian regression

The upper-crossing/solution (US) algorithm for root-finding with strongly stable convergence

1 Introduction

2 Non-crossing nonparametric quantile regression

2.1 Non-crossing quantile regression via ReLU penalty

Lemma 2.1

Remark 2.1

Remark 2.2

Remark 2.3

2.2 ReLU Feedforward neural networks

3 Non-crossing quantile regression for conformal prediction

Remark 3.1

4 Theoretical properties

Theorem 4.1

Theorem 4.2

Theorem 4.3

Corollary 4.1

5 Numerical studies

5.1 Synthetic data

5.2 Real data examples

6 Concluding remarks

Appendix A Appendix: Proofs

Lemma A.1

Lemma A.2

Nonparametric Quantile Regression: Non-Crossing Constraints and Conformal Prediction

Related Research

Estimation of Non-Crossing Quantile Regression Process with Deep ReQU Neural Networks

Sparse Quantile Regression

Quantile index regression

Deep Quantile Regression: Mitigating the Curse of Dimensionality Through Composition

Solution to the Non-Monotonicity and Crossing Problems in Quantile Regression

Prediction intervals with controlled length in the heteroscedastic Gaussian regression

The upper-crossing/solution (US) algorithm for root-finding with strongly stable convergence

1 Introduction

2 Non-crossing nonparametric quantile regression

2.1 Non-crossing quantile regression via ReLU penalty

Lemma 2.1

Remark 2.1

Remark 2.2

Remark 2.3

2.2 ReLU Feedforward neural networks

3 Non-crossing quantile regression for conformal prediction

Remark 3.1

4 Theoretical properties

Theorem 4.1

Theorem 4.2

Theorem 4.3

Corollary 4.1

5 Numerical studies

5.1 Synthetic data

5.2 Real data examples

6 Concluding remarks

Appendix A Appendix: Proofs

Lemma A.1

Lemma A.2