Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
87,134 characters · 16 sections · 79 citation commands
Sparse Quantile Regression
\onehalfspacing
Quantile regression has been increasingly popular since the seminal work of Koenker1978. See Koenker2005 for a classic and comprehensive text on quantile regression and Koenker17 for a review of recent developments. This paper is concerned with estimating a sparse high-dimensional quantile regression model:
where $Y\in \mathbb{R}$ is the outcome of interest, $X\in \mathbb{R}^{p}$ is a $p$-dimensional vector of covariates, $\theta_\ast$ is the vector of unknown parameters, and $U$ is a regression error. Let $Q_{\tau }(U|X)$ denote the $\tau $-th quantile of $U$ conditional on $X$. Assume that $ Q_{\tau }(U|X)=0$ almost surely for a given $\tau \in (0,1)$ and that the data consist of a random sample of $n$ observations $\left( Y_{i},X_{i}\right) _{i=1}^{n}$. As usual, $p$ can be much larger than $n$; however, sparsity $s$, the number of nonzero elements of $\theta_\ast$, is less than $n$.
To date, an $\ell _{1}$-penalized approach to estimating (ref) has been predominant in the literature mainly thanks to its computational advantages. See e.g., belloni2011, Wang13, belloni2014uniform, Belloni2019, zheng2015, Lee2018oracle, Lv2018, Wang18wild and wang2019wp among many others. The $ \ell _{1}$-penalized quantile regression ($\ell _{1}$-PQR hereafter) is akin to the well known approach of Lasso tibshirani1996. Smooth yet non-convex penalized estimation approaches have also been proposed as alternatives to $\ell _{1}$-PQR. These include methods of adaptive Lasso (adaptive $\ell _{1}$-) and non-convex penalized quantile regressions Wu:Liu:09,Wang12,Fan14adaptive,Fan2014,Peng:Wang:15. See also wang2022 for the state-of-the-art theoretical analysis of $\ell _{1}$-based and non-convex penalized quantile regressions.
Recently, there is emerging interest in adopting an $\ell _{0}$-based approach since the latter is regarded as a more direct solution to estimation problem under sparsity. For instance, bertsimas2016 took an $\ell _{0}$-constrained approach in order to solve the best subset selection problem in linear regression models. Huang:2018 proposed a scalable computational algorithm for $\ell _{0} $-penalized least squares solutions. chen2018, chen2018arXiv studied the $\ell _{0}$ -constrained and $\ell _{0}$-penalized empirical risk minimization approaches to high dimensional binary classification problems. Bertsimas:Parys:2020 and Hazimeh:Mazumder:2020 made further advances in $\ell _{0}$-based methods for mean regression models. DHM:jmlr:2021 considered an $\ell _{0}$-regularized approach for a hinge-loss based classification problem.
In this paper, we pursue an $\ell _{0}$-based approach to estimating sparse quantile regression. We are inspired by bertsimas2016, who provided a piece of numerical evidence---without theoretical analysis---that the $\ell _{0}$-constrained least absolute deviation (LAD) estimator outperforms $\ell _{1}$-penalized LAD estimator in terms of both sparsity and predictive accuracy. That is, bertsimas2016 made a convincing case for adopting an $\ell _{0}$ -based approach in median regression. In convex optimization, a constrained approach is equivalent to a penalized method boyd2004convex. For non-convex problems, both are distinct and it is unclear which method is better. Therefore, in the paper, we consider both $\ell _{0}$-constrained and $\ell_{0}$-penalized quantile regression ($ \ell_0 $-CQR and $\ell_0$-PQR hereafter).
The main contributions of this paper are twofold. First, we derive an exponential inequality on the tail probability of the excess quantile predictive risk and apply it to obtain non-asymptotic upper bounds on a triplet of population quantities for the $\ell _{0}$-PQR estimator: the mean excess predictive risk, the mean-square regression function estimation error, and the mean-square parameter estimation error. The resulting rates of convergence for the triplets are at the order of $s \ln p /n$, which is the same as those of $\ell_1$-PQR and non-convex penalized quantile regression belloni2011,wang2019wp,wang2022 and nearly matches the minimax lower bound $s \ln (p/s) /n$ obtained in Theorem 4.1 of wang2019wp. However, the optimal tuning parameter $\lambda$ in $\ell_1$-PQR is of order $ \sqrt{\ln p/n}$, whereas it is of order ${\ln p/n}$ in $\ell_0$-PQR. We also characterize expected Hamming loss for the $\ell_0$-penalized estimator. In a nutshell, $\ell_0$-PQR produces a sparser estimator than $\ell_1$-PQR, while maintaining the same level of prediction and estimation errors. In addition, we establish analogous results for the $\ell_0$-CQR estimator under the assumption that the imposed sparsity is at least as large as true sparsity. Our non-asymptotic results build on bousquet2002 and massart2006 and are applicable for $\ell_0$-based, general $M$ -estimation with a Lipschitz objective function that includes sparse logistic regression as a special case. Therefore, our theoretical results may be of independent interest beyond quantile regression.
The second contribution is computational. Both $\ell_0$-CQR and $\ell_0$-PQR estimation problems can be equivalently reformulated as mixed integer linear programming (MILP) problems. This reformulation enables us to employ efficient mixed integer optimization (MIO) solvers to compute exact solutions to the $\ell _{0}$-based quantile regression problems. However, the method of MIO is concerned with optimization over integers, which could be computationally challenging for large scale problems. To scale up $\ell _{0}$-based methods, bertsimas2016 developed fast first-order approximation methods for both $\ell _{0}$-constrained least squares and absolute deviation estimators. Huang:2018 also proposed a scalable computational algorithm for approximating the $\ell _{0}$ -penalized least squares solutions. Building on these papers, we propose a new first-order computational approach, which can deliver high-quality approximate $\ell _{0}$-PQR solutions and thus can be used as a warm-start strategy for boosting the computational performance of the MILP based implementation approach. As a standalone algorithm, our first-order approach renders the $\ell _{0}$-PQR computationally as scalable as commonly used $\ell_1$-PQR.
As an illustrative application, we consider conformal prediction of birth weights and have a horse race among $\ell _{0}$-CQR, $\ell _{0}$-PQR, $\ell _{1}$-PQR, adaptive Lasso and non-convex penalized quantile regressions with $n\approx 1000$ and $p$ ranging from $p\approx 20$ to $p\approx 1600$. Recently, romano2019conformalized combined conformal prediction with quantile regression and proposed conformalized quantile regression that rigorously ensures a non-asymptotic, distribution-free coverage guarantee, independent of the underlying regression algorithm. When we implement conformal prediction using competing estimation methods, we find that both $\ell _{0}$ -CQR and $\ell _{0}$-PQR are capable of delivering much sparser solutions than $\ell _{1}$-PQR, while maintaining tighter width yet comparable coverage of prediction confidence intervals in the high dimensional settings. Adaptive Lasso and non-convex penalized estimation approaches do improve the performance of $\ell _{1}$-PQR, but they also tend to select more covariates than $\ell _{0}$-CQR and $\ell _{0}$-PQR and perform poorly especially when the covariates are highly dependent on each other. Furthermore, we obtain similar results in Monte Carlo experiments. Therefore, $\ell _{0}$-CQR and $\ell _{0}$-PQR are worthy competitors to $ \ell _{1}$-PQR, adaptive Lasso and non-convex penalized quantile regression approaches ---superior if a researcher prefers sparsity---as supported by non-asymptotic theory, a real-data application and Monte Carlo experiments.
The rest of this paper is organized as follows. In Section (ref), we set up the sparse quantile regression model and present the $\ell _{0}$-based approaches. In Section (ref), we establish non-asymptotic statistical properties of the proposed $\ell _{0}$ -PQR and $\ell _{0}$-CQR estimators. In Section (ref), we provide both MILP- and first-order (FO)-based computational approaches for solving the $\ell _{0}$-PQR problems. In Section (ref), we perform a simulation study on the finite-sample performance of our proposed estimators. In Section (ref), we illustrate our method in a real data application concerning conformal prediction of birth weights. We then conclude the paper in Section (ref). Appendix (ref) collates proofs of all theoretical results of the paper and an online appendix contains further details on the variable selection results of our empirical study.
Let $\Vert a\Vert _{0}$ be the $\ell _{0}$ norm of a vector $a$, which is the number of nonzero components of $a$. The usual $\ell _{1}$ and $\ell _{2} $ norms are denoted by $\Vert \cdot \Vert _{1}$ and $\Vert \cdot \Vert _{2}$, respectively. For any $t$ and $u$, let
Let $\Theta $ denote a parameter space, which is assumed to be a compact subspace of $\mathbb{R}^{p}$. Define
We first define $\ell_0$-CQR. For any given sparsity $q\geq 0$, let $ \widetilde{\theta }$ denote an $\ell _{0}$-constrained quantile regression ($ \ell _{0}$-CQR) estimator, which is defined as a solution to the following minimization problem:
In practice, choosing $q$ is important: $\ell_0$-CQR will result in selecting many more (or far fewer) covariates if the imposed sparsity is too large (or too small).
To mitigate the issue of unknown true sparsity $s$, we now focus on $\ell _{0}$-PQR. Let $\widehat{\theta }$ denote an $\ell _{0}$-PQR estimator, which is defined as a solution to the following minimization problem:
where $\lambda $ is a nonnegative tuning parameter and $k_{0}$ is a fixed upper bound for the true sparsity $s$. In $\ell _{0}$-PQR, the main tuning parameter is $\lambda $. To adapt to an unknown $s$, we rely on $\ell _{0}$ -penalization that is controlled by $\lambda $. The sparsity bound $k_{0}$ is different from $q$ in $\ell _{0}$-CQR. The latter acts as a tuning parameter, which will be calibrated to maximize the predictive performance, whereas the former is predetermined and imposed throughout the implementation of $\ell_{0}$-PQR. We will set $k_{0}$ with a large value in numerical exercises.
To make our proposed estimators operational, we follow the standard machine learning approach. That is, we first randomly split the dataset into three samples: training, validation and test samples. For each candidate value of the tuning parameter $q$ or $\lambda$, we estimate the model using the training sample. Then, the tuning parameter is selected based on the quantile prediction risk using the validation sample. Finally, out-of-sample performance is evaluated using the test sample.
We provide general regularity conditions that include quantile regression as a special case. Define $S(\theta )\equiv \mathbb{E}\left[ \rho (Y,X^{\top }\theta )\right] $.
Note that for quantile regression,
where $F_{U|X}(\cdot |x)$ is the cumulative distribution function of $U$ conditional on $X=x$ and $F_{X}$ is the cumulative distribution function of $ X$. Thus, Assumption (ref) is satisfied.
Assumption (ref) is satisfied for quantile regression with $L=1$. For any two real numbers $x$ and $y$, let $x\vee y\equiv \max \{x,y\}$ and $ x\wedge y\equiv \min \{x,y\}$.
Assumption (ref) requires that each component of $X$ and that of $\theta $ be bounded by a universal constant. This condition could be restrictive yet is commonly adopted in the literature. For example, zheng2018 assumed the uniform boundedness of each of the covariates, citing the literature that points out that \textquotedblleft a global linear quantile regression model is most sensible when the covariates are confined to a compact set\textquotedblright\ to avoid the problem of quantile crossing.
Assumption (ref) is very mild. A similar condition is assumed in massart2006 to avoid measurability issues and to use the concentration inequality by bousquet2002. By Assumption (ref) and taking $\Theta ^{\prime }=\Theta \cap \mathbb{Q}^{p}$, Assumption (ref) (i) holds by the denseness of the set of rational numbers and the continuity of the function $\rho $. Suppose that, for some non-negative random variable $Z$ with $\mathbb{E}\left( Z\right) <\infty $, and $\rho (Y,X^{\top }\theta )\leq Z$ holds with probability 1 for every $\theta \in \Theta $. Then using this condition together with Assumption (ref) (i), we can also deduce from the dominated convergence theorem that $ S(\theta _{\ast })=\inf_{\theta \in \Theta ^{\prime }}S(\theta )$ and thus Assumption (ref) (ii) also holds. In the quantile regression case, we can take the dominating variable $Z$ to be $\left\vert Y\right\vert +pB^{2}$, which has finite mean provided that the mean of $|Y|$ is also finite. The requirement that $\mathbb{E}|Y|<\infty $ is not strictly necessary because we can redefine the quantile regression objective function by $\rho (Y,X^{\top }\theta )-\rho (Y,X^{\top }\theta _{\ast })$, whose magnitude is uniformly bounded above by $2pB^{2}$.
For each $\theta $, define $R(\theta )\equiv \mathbb{E}[|X^{\top }(\theta -\theta _{\ast })|^{2}],$ which is the the expected squared difference of the true quantile regression function $X^{\top }\theta _{\ast }$ and a linear fit evaluated at a given parameter vector $\theta $.
Assumption (ref) relates $R(\theta )$ to the difference of their corresponding quantile predictive risks. Given Assumption (ref) , if, for some $k\geq k_{0}$, the distribution $F_{U|X}\left( z|x\right) $ admits a Lebesgue density $f_{U|X}\left( z|x\right) $ that is bounded below by a positive constant $c_{u}$ for all $z$ in an open interval containing $ \left[ -B^{2}\left( k+s\right) ,B^{2}\left( k+s\right) \right] $ and for all $x$ in the support of $X$, then Assumption (ref) holds with $ \kappa _{0}=\sqrt{c_{u}/2}.$
For any subset $J\subset \{1,...,p\}$, let $X_{J}$ denote the $\left\vert J\right\vert $-dimensional subvector of $X\equiv (X^{(1)},\ldots ,X^{(p)})^{\top }$ formed by keeping only those elements $X^{(j)}$ with $ j\in J$. Suppose that, for some $k\geq k_{0}$ and for any subset $J\subset \{1,...,p\}$ such that $\left\vert J\right\vert \leq k+s$, the smallest eigenvalue of $\mathbb{E}\left( X_{J}X_{J}^{\top }\right) $ is bounded below by a positive constant $\omega $. Since $R(\theta )=(\theta -\theta _{\ast })^{\top }\mathbb{E}\left[ XX^{\top }\right] (\theta -\theta _{\ast })$ and $ \Vert \theta -\theta _{\ast }\Vert _{0}\leq k+s$ for $\theta \in \mathbb{B} (k),$ it then follows that Assumption (ref) holds with $\kappa _{1}=\sqrt{\omega }$. This assumption is related to the sparse eigenvalue condition used in the high dimensional regression literature (see, e.g. Raskutti2011). For example, if $X$ is a random vector with mean zero and the covariance matrix $\Sigma $ whose $(i,j)$ component is $\Sigma _{i,j}=r^{|i-j|}$ for some constant $r>0$, then the smallest eigenvalue of $ \Sigma $ is bounded away from zero where the lower bound is independent of the dimension $p$ van2009 and thus $\mathbb{E} \left( X_{J}X_{J}^{\top }\right) $ is bounded below by a universal positive constant for every $J\subset \{1,...,p\}$.
The following theorem is the key step to the main results of this section for $\ell_0$-PQR.
Results ((ref)) and ((ref)) of Theorem (ref) are non-asymptotic and establish exponential inequalities on the tail probabilities of the excess quantile predictive risk $S(\widehat{\theta })-S\left( \theta _{\ast }\right) $ as well as the $ \ell _{0}$-distance between the $\ell _{0}$-PQR estimator and the true parameter value. Applying inequality ((ref)), we can obtain non-asymptotic upper bounds on a triplet of population quantities: (i) the mean excess predictive risk $\mathbb{E}[S(\widehat{\theta } )-S(\theta _{\ast })]$; (ii) the mean-square regression function estimation error $\mathbb{E}[R(\widehat{\theta })]$; (iii) the mean-square parameter estimation error $\mathbb{E}[\Vert \widehat{\theta }-\theta _{\ast }\Vert _{2}^{2}]$. The results concerning these bounds are given in the next theorem.
If $k_0/s$ is bounded by a fixed constant, we can deduce from Theorem (ref) that
which suggests that the optimal $\lambda $ be of the following form:
where $C_{\lambda }$ is a positive constant that needs to be chosen by a researcher. Under (ref) and the side condition that $k_0/s \leq C_k$ for some fixed constant $C_k$, we have that
We now specialize Theorem (ref) to quantile regression. The following corollary provides the main results for $\ell_0$-PQR.
Corollary (ref) provides non-asymptotic bounds on the mean-square regression function and parameter estimation errors as well as the excess quantile prediction risk. The resulting rates of convergence are of order $ s\ln p/n$, which is the same as those of $\ell _{1}$-PQR and non-convex penalized quantile regression (see, e.g., belloni2011 and Wang13 for earlier results and wang2019wp and wang2022 for the latest results). These are nearly minimax optimal rates of convergence because it is shown in wang2019wp that the minimax lower bound for $\mathbb{E}[\Vert \widehat{\theta }-\theta _{\ast }\Vert _{2}^{2}]$ is of order $s\ln (p/s)/n$ . The optimal tuning parameter $\lambda $ in $\ell _{1}$-PQR is of order $ \sqrt{\ln p/n}$, whereas it is of order ${\ln p/n}$ in $\ell _{0}$-PQR.
Using the method for proving Corollary (ref), we can obtain the following result for $\ell_0$-CQR.
Corollary (ref) shows that the $\ell_0$-CQR estimator is also nearly minimax optimal, provided that the imposed sparsity $q$ is at least as large as the true sparsity $s$ and that $q/s$ is bounded by a fixed constant. Therefore, our theory predicts that $\ell_0$-PQR and $\ell_0$-CQR would perform similarly in applications.
Applying ((ref)) of Theorem (ref), we now derive a theoretical result regarding the $\ell _{0}$-PQR in terms of expected Hamming loss. Specifically, the following theorem presents an upper bound on the expectation of the Hamming distance between $\widehat{\theta }$ and $\theta _{\ast }$.
Note that $D_{H}(\widehat{\theta },\theta _{\ast })\equiv s^{-1}\Vert \widehat{\theta }-\theta _{\ast }\Vert _{0}$ is the Hamming distance---normalized by diving it by $s$---between $\widehat{\theta }$ and $ \theta _{\ast }$, that is, $s^{-1}$ times the number of elements of the $ \ell _{0}$-PQR estimator that are different from the corresponding elements of the true parameter vector. Theorem (ref) shows that $\mathbb{ E}[D_{H}(\widehat{\theta },\theta _{\ast })]$ can be bounded by a constant that is slightly larger than 4, provided that the tuning parameter $\lambda $ is suitably chosen. Note that
Thus, we do not expect that $\mathbb{E}[\Vert \widehat{\theta }-\theta _{\ast }\Vert _{0}]$ can be small since it is impossible to make $\mathbb{P}( \widehat{\theta }\neq \theta _{\ast })$ small. Instead, what we obtain in Theorem (ref) is that $\mathbb{E}[\Vert \widehat{\theta } -\theta _{\ast }\Vert _{0}]$ is bounded by $(4+\nu )s$, independent of $p$.
In view of Theorem (ref), Theorem (ref) suggests that the estimated sparsity and the selected set of covariates of $\ell _{0}$-PQR cannot be too distinct from $s$ and the true set of nonzero elements of $\theta _{\ast }$. By Theorem (ref) , the resulting sparsity of $\ell _{0}$-PQR is likely to be substantially smaller than $k_{0}$ with a suitable choice of $\lambda $ and $ k_{0}$; therefore, we expect that the constraint $\theta \in \mathbb{B} (k_{0})$ in (ref) will not be binding in practice. Moreover, since the choice of $C_{\lambda }$ in Theorem (ref) is independent of $(s,n,p)$, the minimax optimal rates are still intact.
The MIO approach is useful for solving variable selection problems with $ \ell _{0}$-norm constraints or penalties bertsimas2016,chen2018, chen2018arXiv. Assume that the parameter space $\Theta $ takes the form $\Theta =\prod\nolimits_{j=1}^{p}\left[ \underline{\theta }_{j},\overline{\theta }_{j}\right] ,$ where $\underline{ \theta }_{j}$ and $\overline{\theta }_{j}$ are lower and upper parameter bounds such that $-\infty <\underline{\theta }_{j}\leq \theta _{j}\leq \overline{\theta }_{j}<\infty $ for $j\in \{1,...,p\}$. We now present an implementation of $\ell _{0}$-PQR, which builds on the method of mixed integer linear programming (MILP). Specifically, the $\ell _{0}$-penalized minimization problem ((ref)) can be equivalently reformulated as the following MILP problem:
We now explain the equivalence between ((ref)) and ((ref)). If we remove from the problem ((ref)) the second term of the objective function as well as all the $(d_{1},...,d_{p})$ control variables together with their constraints ((ref)) and ((ref)), the resulting minimization problem reduces to the linear programming reformulation of the standard linear quantile regression problem Koenker2005. In the presence of the penalty term and the $(d_{1},...,d_{p})$ controls, the inequality and dichotomization constraints ((ref)) and ((ref)) ensure that, whenever $d_{j}=0$, the value $\theta _{j}$ must also be zero and the sum $\sum\nolimits_{j=1}^{p}d_{j}$ thus captures the number of non-zero components of the vector $\theta $. The last constraint (ref) imposes that the estimated sparsity is at most $k_0$. As a result, both minimization problems ((ref)) and ((ref)) are equivalent. This equivalence enables us to employ modern MIO solvers to solve $\ell _{0}$-PQR problems.
The MIO formulation ((ref)) is concerned with optimization over integers, which could be computationally challenging for large scale problems. bertsimas2016 have developed discrete first-order algorithms enabling fast computation of near optimal solutions to $\ell _{0}$-constrained least squares and least absolute deviation estimation problems. Huang:2018 have also proposed fast and scalable algorithms for computing approximate solutions to $\ell _{0}$-penalized least squares estimation problems. These algorithms build on the necessary conditions for optimality in the $\ell _{0}$-constrained or penalized optimization problems. Motivated from these papers, in this subsection, we present a first-order approximation algorithm that can be used as either a standalone solution algorithm or a warm-start strategy for enhancing the computational performance of our MIO approach to the $\ell _{0} $-PQR problem.
For $\tau \in (0,1)$, the quantile regression objective function ((ref)) can be equivalently expressed as
The function $S_{n}(\theta )$ is nonsmooth. Following nesterov2005, we can construct a smooth approximation of $S_{n}(\theta )$ by
where $w$ denote the vector of controls $\left( w_{1},...,w_{n}\right) $ in the maximization problem ((ref)). Note that nesterov2005's smoothing method is different from a convolution-based smoothing method for quantile regression by FGH:2021 and HE2023.
Assume that the parameter space $\Theta $ is of an equilateral cube form $ \Theta =\left[ -B,B\right] ^{p}$ for some $B>0$. Let $t$ be any given vector in $\mathbb{R}^{p}$. Let $\widehat{\beta }$ be a solution to the following $ \ell _{0}$-penalized minimization problem:
where $\lambda $ is a non-negative penalty tuning parameter. It is straightforward to see that the solution $\widehat{\beta }$ can be computed as follows. Let $\widetilde{\beta }$ be a $p$ dimensional vector given by
for $j\in \{1,...,p\}$. Then the solution $\widehat{\beta }=$ $\widetilde{ \beta }$ if $\Vert \widetilde{\beta }\Vert _{0}\leq k_{0}$. Otherwise, letting $S(t)$ denote the set of $k_{0}$ indices that keep track of the largest $k_{0}$ components of $t$ in absolute value, we set $\widehat{\beta } _{j}=\widetilde{\beta }_{j}$ for $j\in S(t)$ and $\widehat{\beta }_{j}=0$ for $j\notin S(t)$. Therefore, the problem ((ref)) admits a simple closed-form solution. We will exploit this fact and develop a first-order approximation algorithm.
Define
For any vector $t\in \mathbb{R}^{p}$, suppose we can construct a quadratic envelope of $S_{n}(\theta ;\delta )$ with respect to the vector $t$ in the sense that
for some non-negative real scalar $l$, which does not depend on the parameter vector $\theta $. Note that ((ref)) holds whenever the gradient function $\bigtriangledown _{\theta }S_{n}(\cdot ;\delta )$ is Lipschitz continuous such that
for some Lipschitz constant $h$, which does not depend on $t$ and $t^{\prime }$. By the envelope theorem,
where $(\widehat{w}_{1,\delta },...,\widehat{w}_{n,\delta })$ is the solution to the minimization problem ((ref)). Using nesterov2005, we can deduce that ((ref)) holds with
and hence ((ref)) holds for every $l\geq h$.
Define
Note that $\widetilde{Q}_{n}(\theta ;t,\delta ,l)$ is an upper envelope of $ Q_{n}(\theta ;\delta )$ around the vector $t$ with the property that $ \widetilde{Q}_{n}(t;t,\delta ,l)=Q_{n}(t;\delta )$.
For $t\in \mathbb{R}^{p}$, define the mapping
Arranging the terms, we can easily deduce
We say that a point $t\in \mathbb{R}^{p}$ is a stationary point of the mapping $H_{\delta ,l}$ if $t\in H_{\delta ,l}(t)$. For each given value of $ \delta $, let $\widehat{\theta }_{\delta }$ denote a solution to the problem of minimizing $Q_{n}(\theta ;\delta )$ over $\theta \in \mathbb{B}(k_{0})$. We propose to approximate $\widehat{\theta }_{\delta }$ by solving for the stationary point of the mapping $H_{\delta ,l}$. This can be justified by the following proposition, which is a straightforward extension of Theorem 3.1 of bertsimas2016.
Proposition (ref) implies that any solution to the minimization of $Q_{n}(\theta ;\delta )$ over $\theta \in \mathbb{B}(k_{0})$ is also a stationary point of the mapping $H_{\delta ,l}$. Moreover we can solve for a stationarity point by iterating until convergence. Result ((ref)) indicates that the convergence rate is $O(N^{-1})$, where $N$ is the number of performed iterations. Note that we can use ((ref)) to obtain a closed-form solution to the $\ell _{0}$ -penalized minimization problem ((ref)) for every $t\in \mathbb{R}^{p}$ and therefore solving for a stationary point of $H_{\delta ,l}$ would incur relatively little computational cost.
We now turn to the $\ell _{0}$-PQR problem ((ref)). The next proposition builds on the results of nesterov2005 concerning the uniform approximation bound of the smooth function to the non-smooth quantile loss function. Let $c_{\tau }\equiv \tau ^{2}\vee \left( 1-\tau \right) ^{2}$.
Given a tolerance level $\epsilon $, Proposition (ref) implies that, for any given $\delta \leq 2\epsilon c_{\tau }^{-1}$, if we solve for the minimization of $Q_{n}(\theta ;\delta )$ over $\theta \in \mathbb{B} (k_{0})$, the resulting solution $\widehat{\theta }_{\delta }$ is an $ \epsilon $-level approximate $\ell _{0}$-PQR estimator in the sense that
This thus yields the following algorithm for computing a near optimal solution to the $\ell _{0}$-PQR problem ((ref)).
We end this section with a remark that unlike MIO, the first-order approximation method only delivers a feasible solution to the minimization of $Q_{n}(\theta ;\delta )$ over $\theta \in \mathbb{B}(k_{0})$ and this solution does not necessarily coincide with a global optimal solution to the $\ell _{0}$-PQR problem ((ref)).
In this section, we perform Monte Carlo simulation experiments to evaluate the performance of our $\ell _{0}$-based quantile regression approaches. We consider the following data generating setup. Let $Z=(Z_{1},...,Z_{p-1})$ be a $p-1$ dimensional multivariate normal random vector with mean zero and covariance matrix $\Sigma $ with its element $\Sigma _{i,j}=\left( 0.5\right) ^{\left\vert i-j\right\vert }$. Let $X=(X_{1},...,X_{p})$ be a $p$ dimensional covariate vector with its components $X_{1}=1$ and $ X_{j}=Z_{j-1}1\{\left\vert Z_{j-1}\right\vert \leq 6\}$ for $j\in \{2,..,p\}$ . The outcome $Y$ is generated according to the model:
where $\varepsilon $ is a random disturbance which is independent of $X$ and follows the univariate normal distribution with mean zero and standard deviation $0.25$. We considered two configurations of the true parameter vector $\theta _{\ast }$. For configuration (i), we set the sparsity $s=5$ and the true parameter value $\theta _{\ast ,j}=1$ for $s$ equispaced values. For configuration (ii), we employed a more challenging case where $ s=20$ and all the $s$ nonzero components of $\theta _{\ast }$ occurred at equispaced indices between 1 and $p$ with its first 5 components equal to 1 and remaining $\left( s-5\right) $ nonzero components being set to be $ \left( 2^{-1},2^{-2},...,2^{-(s-5)}\right) $, which decreased exponentially to zero.
We compared the finite-sample performance among the $\ell _{0}$-PQR and $ \ell _{0}$-CQR of the present paper, the $\ell _{1}$-PQR of belloni2011, the adaptive Lasso penalized quantile regression of Fan14adaptive and the nonconvex penalized quantile regression of Wang12. In each simulation repetition, we generated a training sample of $n=100$ observations for estimating the parameter vector $\theta $ and another independent validation sample of $100$ observations for calibrating the tuning parameters of these estimation approaches. Moreover, we also generated a test sample of $5000$ observations for evaluating the out-of-sample predictive performance.
We focused our simulation study on median regression ($\tau =0.5$). To implement the $\ell _{1}$-PQR approach, we used the $\ell _{1}$-penalized quantile regression estimator of belloni2011 with the penalty level given by
where $\Lambda \left( 1-\alpha |X\right) $ is the $\left( 1-\alpha \right) $ level quantile of the random variate $\Lambda $, which is defined in belloni2011, conditional on the covariate vector $ X $. Following belloni2011, we set $\alpha =0.1$. Moreover, we calibrated the optimal tuning value $c_{BC}$ from a set of candidate values $ \mathcal{S}$ using the aforementioned validation sample in the setup with $ p\geq 100$. For the low dimensional setup with $p<100$, we performed this calibration over an expanded set $\mathcal{S\cup \{}0\}$, thereby allowing for an estimating model that did not penalize any parameter. For simulations under parameter configuration (i), we set $\mathcal{S}=\{0.1,0.2,...,1.9,2\}$ . Under parameter configuration (ii), which is a more difficult case for estimation, we further enlarged the tuning value search space by taking $ \mathcal{S}$ to be $\{0.01,0.02,...1.99,2\}$.
To implement the $\ell _{0}$-CQR method, we solved over the training sample the $\ell _{0}$-constrained estimation problem ((ref)) for sparsity level $q$ ranging from $1$ up to $p\wedge 25$. To solve ((ref)) with $\tau =0.5$ for a given value of $q$, following bertsimas2016, we used the MIO-based, $\ell _{0}$ -constrained LAD approach with a warm-start strategy by supplying the MIO solver an initial guess computed via the discrete first-order approximation algorithms. We then calibrated the optimal sparsity level among this set of $ q$ values using the calibration sample. The resulting $\ell _{0}$-CQR estimator was then constructed based on the model associated with the calibrated optimal sparsity level.
For the $\ell _{0}$-PQR method, noting that the scale of the quantile regression objective function $S_{n}(\theta )$ varies whenever that of $Y$ changes, to relate the penalty term to the scale of $Y$ and to the derived rate (ref), we adopted the following simple rule:
which is proportional to the sample average of the absolute value of $Y$.
For a given value of $c$ in ((ref)), we solved the problem ((ref)) with $k_{0}=100$ using our MIO computational approach of Section (ref), where we warm-started the MIO solver by supplying as an initial guess the approximate solution obtained through the first-order method of Section (ref). As in the $\ell _{1}$-PQR case, we calibrated the optimal tuning scalar $c$ over the set $\mathcal{S}$ using the calibration sample in the setup with $p\geq k_{0}$ and over the expanded set $\mathcal{S\cup \{} 0\} $ in the setup with $p<k_{0}$.
We provided further details here on the implementation of our first-order approximation procedure in Algorithm (ref). We set the tolerance level $\epsilon $ to be $2\cdot 10^{-4}$ and parameter $l$ of the quadratic envelope in ((ref)) to be $2h$, where $h$ is the Lipschitz constant given by ((ref)). Note that Algorithm (ref) also requires an initial guess. We therefore ran it for $T=50$ times, each of which was performed with a different initial guess and used the output that delivered the best penalized objective function value in ((ref)) as the resulting first-order approximate solution. We chose these $T$ initial guesses sequentially where the first one was the $\ell _{1} $-PQR solution of belloni2011 implemented with its tuning value $c_{BC}$ set to be identical to the given value $c$ in ((ref)) whereas, for $t\in \{2,...,T\}$, the $t$-th initial guess was subsequently constructed as the solution to the standard quantile regression of the outcome $Y$ on those covariates selected in the output of Algorithm (ref) which was initiated with the $\left( t-1\right) $-th initial guess. We found this implementation procedure worked very well in both our simulation study here and the empirical application of Section (ref).
We specified the parameter space $\Theta $ to be $[-10,10]^{p}$ for the MIO computation of both the $\ell _{0}$-PQR and $\ell _{0}$-CQR estimators. Throughout this paper, we used the MATLAB implementation of the Gurobi Optimizer (version 8.1.1) to solve all the MIO problems. Moreover, all numerical computations were done on a desktop PC (Windows 7) equipped with 128 GB RAM and a CPU processor (Intel i9-7980XE) of 2.6 GHz. To reduce computation cost in all MIO computations associated with the covariate configuration of $p=500$, we set the MIO solver time limit to be 10 minutes beyond which we forced the solver to stop early and used the best discovered feasible solution to construct the resulting $\ell _{0}$-PQR and $\ell _{0}$ -CQR estimators.
We also compared our $\ell _{0}$-based estimators with the adaptive Lasso quantile regression approach of Fan14adaptive. Specifically, the latter approach seeks to minimize the following weighted $\ell _{1}$ -penalized quantile regression objective function
where $\widehat{\theta }^{ini}=(\widehat{\theta }_{1}^{ini},...,\widehat{ \theta }_{p}^{ini})$ is an initial high dimensional quantile regression estimator and $g_{\mu }$ is a penalty weight function. For implementation, we set $\widehat{\theta }^{ini}$ to be the $\ell _{1}$-PQR estimator and considered the following two choices for the penalty weight $g_{\mu }$. The first choice, which is based on the derivative of the smoothly clipped absolute deviation (SCAD) penalty function fan2001, is given by
for some parameters $a>2$ and $\mu \geq 0$. For ease of reference, we use AL-SCAD as shorthand for the adaptive Lasso quantile regression approach implemented with the SCAD based penalty weight ((ref)). For the second choice, we specify $g_{\mu }$ to be the derivative of the minimax concave penalty (MCP) function zhang2010, which is given by
for some parameters $a>1$ and $\mu \geq 0$. We refer to AL-MCP as shorthand for the estimation approach based on the minimization of ((ref)) with the MCP based penalty weight ((ref)).
Finally, we considered the approach of nonconvex penalized quantile regression of Wang12 implemented with either the SCAD or MCP penalty. Specifically, this approach is based on minimizing the penalized objective function
where $G_{\mu }$ is either the SCAD or the MCP penalty function, both of which are nonconvex. Following zou2008 and Wang12, we adopted the local linear approximation algorithm to solve this nonconvex minimization problem. This algorithm proceeds as follows. Let $\widehat{ \theta }^{(0)}=(\widehat{\theta }_{1}^{(0)},...,\widehat{\theta }_{p}^{(0)})$ be an initial estimator, which we take as the $\ell _{1}$-PQR estimator. Given an estimator $\widehat{\theta }^{(m)}$ at the $m$th iteration stage, we solve for $\widehat{\theta }^{(m+1)}$ by minimizing ((ref) ) with $\widehat{\theta }^{ini}$ being replaced by $\widehat{\theta }^{(m)}$ and $g_{\mu }$, which is the derivative of $G_{\mu }$, taking the form ((ref)) if $G_{\mu }$ is the SCAD penalty or ((ref)) if $ G_{\mu }$ is the MCP penalty function. We then iterate this process until the vector of weight differences $\left[ g_{\mu }\left( \left\vert \widehat{ \theta }_{j}^{(m+1)}\right\vert \right) -g_{\mu }\left( \left\vert \widehat{ \theta }_{j}^{(m)}\right\vert \right) \right] $ converges in $\ell _{2}$ -norm within a numerical tolerance of $10^{-4}$. We refer to QR-SCAD and QR-MCP as shorthand for the nonconvex penalized quantile approaches implemented respectively with the SCAD and MCP penalties. Throughout the implementation of all the methods that use ((ref)) and ((ref)), we set $a=3.7$ and focused on the calibration of the tuning parameter $\mu $, which was performed over the set $\mathcal{S}$ using the calibration sample in the setup with $p\geq 100$ and over the expanded set $ \mathcal{S\cup \{}0\}$ in the setup with $p<100$.
We reported performance results based on 100 simulation repetitions. We considered the following performance measures. Abusing the notation a bit, let $\widehat{\theta }$ denote the estimated parameters under a given quantile regression approach. To assess the predictive performance, we reported the relative risk, which is the ratio of the median predictive risk evaluated at the estimate $\widehat{\theta }$ over that evaluated at the true value $\theta _{\ast }$. We approximated the out-of-sample predictive risk using the generated 5000-observation test sample. Let $in\_RR$ and $ out\_RR$ respectively denote the average of in-sample and that of out-of-sample relative risks over the simulation repetitions.
We also reported the estimation performance in terms of both the average parameter estimation error defined as $\mathbb{E}[ \Vert \widehat{\theta } -\theta _{\ast }\Vert _{2}] $ and the average regression function estimation error defined as $\mathbb{E}[|X^{\top }(\widehat{\theta }-\theta _{\ast })|^{2}]$. Finally, we examined the variable selection performance. We say that a covariate $X_{j}$ is effectively selected if and only if the magnitude of $\widehat{\theta }_{j}$ is larger than a small tolerance level (e.g., $10^{-5}$ as used in our numerical study) which is distinct from zero in numerical computation. Let $Avg\_sparsity$ denote the average number of effectively selected covariates. Let $Corr\_sel$ be the proportion of the truly relevant covariates being effectively selected. Let $Orac\_sel$ be the proportion of obtaining an oracle variable selection outcome where the set of effectively selected covariates coincides exactly with that of the truly relevant covariates. Finally, let $Num\_irrel$ denote the average number of effectively selected covariates whose true regression coefficients are zero.
For parameter configuration (i), we performed simulations with $p\in \{10,500\}$ to assess the performance in both the low and high dimensional settings. The results for these two settings are presented respectively in Tables (ref) and (ref). For $\ell _{0}$-PQR, we report performance measures for both the implementation based on the first-order (FO) approximation and that based on the MIO, which was warm-started by using the FO solutions as initial guesses. We find that, regarding the predictive performance, all the competing approaches performed comparably well for both in-sample and out-of-sample relative risks under the low dimensional covariate design. By contrast, for the high dimensional design, $ \ell _{1}$-PQR was considerably dominated by all the other approaches in terms of out-of-sample predictive performance.
Turning to the variable selection results, we note that all the eight estimation approaches had perfect $Corr\_sel$ rates and hence were effective for selecting the relevant covariates. However, superb $Corr\_sel$ performance might just be a consequence of overfitting, which may result in excessive selection of irrelevant covariates and adversely impact on the out-of-sample predictive performance. From the results on the variable selection performance measures, we note that the number of irrelevant variables selected under $\ell _{1}$-PQR was quite large relatively to those under the other seven approaches in the high dimensional setup even though all of the considered estimation approaches exhibited the effect of reducing the covariate space dimension. This echoes with the finding in the setup of $ p=500$ that $\ell _{1}$-PQR had far better in-sample fit in terms of $in\_RR$ yet worse out-of-sample fit in terms of $out\_RR$ relatively to all the other quantile regression approaches. Besides, while we could observe nonzero and high values of $Orac\_Sel$ for the $\ell _{0}$-PQR, $\ell _{0}$ -CQR, adaptive Lasso based and other nonconvex penalized estimation approaches in both the low and high dimensional setups, the $\ell _{1}$-PQR approach could rarely induce oracle variable selection outcome in these simulations. We also find, except for $\ell _{1}$-PQR, which performed relatively poorly, all the other approaches performed comparably well in the parameter and regression function estimation performances. Finally, for the $ \ell _{0}$-PQR approach, the FO-based algorithm as a standalone solution algorithm also performed very well. In fact, we find that in the high dimensional setup, all the MIO-based $\ell _{0}$-PQR computations could not converge within the 10-minute computational time limit and the discovered solutions upon early stopping coincided with the FO-based solutions which were used to warm-start the MIO solver. This indicates that the first-order algorithm already located high-quality $\ell _{0}$-PQR solutions upon which further improvements through the global optimization solver of MIO could not be obtained within the given computational time constraint. These simulation results thus shed lights on the usefulness of our first-order approximation approach for solving $\ell _{0}$-PQR problems in the presence of computational resource constraints.
In this section we report results of simulations conducted with $p=500$ under parameter configuration (ii) where $s=20$ and the true nonzero regression coefficients can decay exponentially to zero. Under this simulation design, it turns out that none of the estimation approaches could yield nonzero $Corr\_sel$ and $Orac\_Sel$ values. This result is not surprising as the true coefficients of some relevant covariates are of very small magnitudes so that none of the estimation methods in this simulation study could effectively select these variables. We now summarize in Table 3 the results of remaining performance measures. From this table, it is evident that $\ell _{1}$-PQR tended to select far more irrelevant covariates and thus performed quite poorly in terms of out-of-sample prediction as well as parameter and regression function estimation performances. The apparent overfitting of $\ell _{1}$-PQR could be mitigated through employing instead the adaptive Lasso based approaches of AL-SCAD and AL-MCP or the nonconvex penalized estimation approaches of QR-SCAD and QR-MCP. Compared to these four alternative penalized quantile regression approaches, it is worth noting that both the $\ell _{0}$-PQR and $\ell _{0}$-CQR approaches were capable of further reducing substantially the incidence of selecting irrelevant covariates while retaining comparable performance in prediction, parameter and regression function estimation, and the selection of relevant covariates. Finally, we also note in this simulation that the FO-based implementation of $\ell _{0}$-PQR delivered similar pattern of estimated sparsity to that of the MIO-based $\ell _{0}$-PQR and its predictive and estimation performances did not fall much behind those of the MIO-based $ \ell _{0}$-PQR and $\ell _{0}$-CQR approaches. These results also suggest that our FO-based $\ell _{0}$-PQR method can be a useful standalone approach to sparse estimation of high dimensional quantile regression models.
In this section, we compare $\ell _{0}$-PQR and $\ell _{0}$-CQR with $\ell _{1}$-PQR and the other four alternative penalized quantile regression approaches (AL-SCAD, AL-MCP, QR-SCAD, QR-MCP) of Section (ref) via a real data application to conformal prediction of birth weights. In particular, we employ conformalized quantile regression romano2019conformalized to construct prediction intervals for birth weights.
We now describe the split conformalized quantile regression procedure romano2019conformalized. First, we split the data into a proper training set, indexed by $\mathcal{I}_{1}$, and a calibration set, indexed by $\mathcal{I}_{2}$. For each quantile regression algorithm, we use the proper training set $\mathcal{I}_{1}$ to obtain the estimates of two conditional quantile functions $\widehat{Q}_{\alpha /2}(Y|X=x)$ and $ \widehat{Q}_{1-\alpha /2}(Y|X=x)$ for a given level $\alpha \in (0,0.5)$. Then, the following scores are evaluated on the calibration set $\mathcal{I} _{2}$ as
for each $i\in \mathcal{I}_{2}$. Finally, given new covariates $X_{n+1}$, construct the prediction interval for $Y_{n+1}$ as\
where $Q_{1-\alpha }(E,\mathcal{I}_{2})$ is the $(1-\alpha )(1+1/|\mathcal{I} _{2}|)$-th empirical quantile of $\{E_{i}:i\in \mathcal{I}_{2}\}$. Remarkably, Theorem 1 of romano2019conformalized guarantees that the prediction interval (ref) satisfies the marginal, distribution-free, finite-sample coverage in the sense that
provided that the data $\{(Y_{i},X_{i}):i=1,\ldots ,n+1\}$ are exchangeable.
We look at the dataset on birth weights originally analyzed by almond2005costs. We use the excerpt from cattaneo2010efficient available at \url{http://www.stata-press.com/data/r13/cattaneo2.dta}. The sample size is 4642 and the outcome of interest is infant birth weight measured in kilograms. The basic covariates include 20 variables concerning mother's age, mother's years of education, father's age, father's years of education, number of prenatal care visits, trimester of first prenatal care visit, birth order of an infant, months since last birth, an indicator variable whether a newborn died in previous births, mother's smoking behavior during pregnancy, mother's alcohol consumption during pregnancy, mother's marital status, mother's and father's hispanic status and race (being white or not), an indicator variable whether a mother was born abroad, and three dummy variables indicating seasons of the birth.
We consider the following five different dictionary specifications. The first one is concerned with a covariate vector of $p=21$ that includes a regression intercept together with the aforementioned 20 basic explanatory variables. The second specification modifies the first by discretizing both the maternal and paternal years of education into four categories indicating whether the schooling year is less than 12, exactly 12, between 12 and 16, or at least 16. In addition, we replace number of prenatal care visits, months since last birth and both parents' ages by cubic B-spline terms using 4 interior B-spline knots. These allow us to approximate smooth functions of these variables in the quantile regression analysis. We exclude the B-spline intercept terms so that the resulting covariate vector for the second specification has dimension $p=49$. The third covariate specification consists of all variables in the second specification together with those obtained by interacting the B-spline expansion terms with the other explanatory variables. This then renders $p=609$ in the third covariate specification scenario. Both the fourth and fifth specifications are constructed using the same procedure as for the third case except that we enlarge the covariate dimensions by using respectively 12 and 16 interior B-spline knots for these two high dimensional scenarios, each of which comprises 1281 and 1617 covariates respectively. Finally, for each covariate specification, to put the variables on a similar scale, all stochastic covariates thus constructed are further standardized to have mean zero and variance unity.
We conduct conformal prediction of infant birth weights with nominal level $ \alpha =0.1$. We split the sample randomly into four subsets of about equal size: $\mathcal{I}_{1}$, $\mathcal{I}_{2}$, $\mathcal{I}_{3}$ and $\mathcal{I }_{4}$. As described above, the set $\mathcal{I}_{1}$ is the training sample for the estimation of conditional quantile functions. We perform this estimation respectively using the $\ell _{0}$-PQR, $\ell _{0}$-CQR, $\ell _{1}$-PQR, AL-SCAD, AL-MCP, QR-SCAD and QR-MCP approaches. We calibrate the tuning parameters in these competing quantile regression approaches using the validation sample $\mathcal{I}_{2}$. Let $\mathcal{S} =\{0.1,0.2,...,1.9,2\}\cup \{0.1(0.7)^{s}:s=1,...,8\}$. For penalized estimation approaches, we calibrate the tuning constants $c_{BC}$ of ((ref)), $c$ of ((ref)), and $ \mu $ of ((ref)) and ((ref)) over the set $\mathcal{S}\cup \{0\}$ for the estimation with the first two covariate specifications where $p\in \{21,49\}$ and over the set $\mathcal{S} $ for the high dimensional estimation case where $p\in \{609,1281,1617\}$. Except for these modifications, implementation of all penalized quantile regression approaches and estimation and calibration for the $\ell _{0}$-CQR are performed in the same fashion\ as described in the simulation study.
With the calibrated optimal tuning parameter value, we use the set $\mathcal{ I}_{3}$ to estimate the out-of-sample quantile prediction risk and conformalize quantile regression estimates by constructing $\{E_{i}:i\in \mathcal{I}_{3}\}$. We then evaluate the coverage performance of the prediction interval (ref) over the test sample $\mathcal{I}_{4}$. We carry out 10 replications of such random splitting exercises and report the averages of estimated sparsity, out-of-sample prediction risk as well as length and coverage of the prediction interval across these replications. To mitigate the computational cost, every MIO computation in this empirical study is conducted under a 5-minute computational time constraint.
We summarize in Figures (ref) and (ref) statistical performances under the aforementioned competing quantile regression approaches. For the basic covariate specification with $p=21$, we also juxtapose and compare the performance results with those obtained through standard quantile regression (Std-QR) of Koenker1978. As the Std-QR approach does not incur any tuning parameter, we conduct the Std-QR based conformal prediction by splitting the dataset evenly into three subsamples of which we use the first for parameter estimation, the second for estimating the out-of-sample quantile prediction risk and conformalizing the quantile regression estimates, and the third for estimating the coverage of the conformalized prediction interval. We also perform 10 replications of this sample-splitting procedure and report the average performance results under the the Std-QR approach. Moreover, for each regularized estimation method, we report in Figure (ref) its computational performance, which is measured by the employed CPU seconds that are averaged over the range of tuning parameter values and across the random splitting replications. See also online appendix of the paper for further details on the variable selection results of this empirical study.
From Figure (ref), we find that $\ell _{1}$-PQR tended to induce a far denser estimating model than all the other competing regularized estimation approaches across all covariate specifications. At the 5% quantile level, the average number of selected covariates under $\ell _{1}$-PQR could go over 30 when $p>10^{3}$, whereas that quantity did not exceed 12 for all the other high dimensional quantile regression approaches. At the 95% quantile level, the average estimated sparsity for $\ell _{1}$-PQR was around 12 when $p=21$. Yet, this figure rose quickly as the covariate dimension increased. It reached around 30, then moved up to 72, and eventually escalated toward 155 as $p$ increased from 609 to 1617. While this excessive sparsity pattern could be curbed under both our $\ell _{0}$-based approaches and the approaches of AL-SCAD, AL-MCP, QR-SCAD and QR-MCP, we note that, for the high dimension scenario with $p=1617$, the average estimated sparsities under AL-SCAD and AL-MCP still went over 60 and those under QR-SCAD and QR-MCP were smaller yet remained above 30. By contrast, none of the $\ell _{0}$-based approaches selected more than 9 variables across all the covariates specifications.
Concerning the predictive performance, at the lower quantile level, all the estimation approaches exhibited commensurate out-of-sample quantile prediction risks in the cases where $p\in\{21,49\}$. For $p=609$, both $\ell _{0}$-CQR and $\ell _{1}$-PQR performed slightly better than the other approaches. However, in the high dimensional scenarios where $ p\in\{1281,1617\}$, both $\ell _{0}$-PQR and $\ell _{0}$-CQR had similar prediction performances, which clearly dominated those of $\ell _{1}$-PQR and the other four alternative penalized quantile regression approaches. At the higher quantile level, all the estimation approaches also performed comparably well in the low dimensional cases. However, except for the $\ell _{0}$-based approaches whose out-of-sample prediction risks were all capped below around 6.5%, those risks under all the other high dimensional estimation approaches went over 60% at $p=1281$ and surged toward 90% when $p$ reached 1617.
For the conformalized prediction intervals, Figure (ref) indicates that coverages of these intervals were quite similar across all the estimation approaches and on average dovetailed well with the nominal size. However, lengths of the prediction intervals varied across both the estimation methods and the covariate specifications. For $p\in\{21,49\}$, all methods performed comparably well though QR-SCAD and QR-MCP appeared to produce slightly wider prediction intervals at $p=49$. For the case with $ p=609$, QR-MCP delivered the tightest prediction interval with length 1.67; those of the other approaches had lengths ranging from 1.75 to 1.84. For $ p=1281$, the prediction interval under AL-SCAD was the tightest whereas that under AL-MCP was the widest. The maximal difference in these interval lengths was about 0.18 in this covariate specification. Yet, this maximal difference reached 0.33 and hence nearly doubled at $p=1617$ where the prediction interval lengths under the $\ell _{0}$-based approaches were all around 1.79 yet those under the other regularized estimation approaches were at least around 2 and this could exceed 2.11 for the case of $\ell _{1}$ -PQR. On the whole, the statistical performance results in Figures (ref) and (ref) reveal that the $\ell _{0}$-based approaches were capable of delivering sparser solutions than all the other competing estimation approaches whilst maintaining quite favorable performances in the prediction accuracy.
We now turn attention to computational performance of the considered estimation approaches in this study. From Figure (ref), it is evident that $\ell _{1}$-PQR enjoyed the best computational performance with its average computation time being capped below around 31 CPU seconds across all the estimation scenarios. FO-based $\ell _{0}$-PQR also performed very well. It was initially slower than both AL-SCAD and AL-MCP in the low dimensional estimation cases where its computational time remained below 24 CPU seconds. Yet, its computational performance became quite competitive relatively to those of the adaptive Lasso based approaches in the high dimensional cases with $p\in\{1281,1617\}$ where the average computation time for both QR-SCAD and QR-MCP could go well over 110 seconds.
Relative to the other approaches, the two methods, $\ell _{0}$-CQR and MIO-based $\ell _{0}$-PQR, were observed to be far more computationally intensive. Across all estimation scenarios, Figure (ref) indicates a substantial computational performance difference between the MIO and non-MIO-based approaches. This could be anticipated because of the high computational complexity in the estimation problems for both $\ell _{0}$-CQR and MIO-based $\ell _{0}$-PQR. Based on these empirical results, we find that the FO-based implementation of $\ell _{0}$-PQR could strike a good balance between statistical and computational performances and thus be a valuable standalone estimation approach. On the whole, we note that the $ \ell _{0}$-based approaches could be competitive alternatives to $\ell _{1}$ -PQR, adaptive Lasso based and the other nonconvex penalized estimation approaches employed in this numerical study.
In this paper, we study estimation of a sparse high dimensional quantile regression model. The main contributions of this paper are twofold. First, we derive non-asymptotic expectation bounds on the excess quantile prediction risk as well as the mean-square parameter and regression function estimation errors of both the $\ell _{0}$-PQR and $\ell _{0}$-CQR estimators. These theoretical results imply the near minimax optimal rates of convergence. Moreover, we characterize expected Hamming loss for the $ \ell _{0}$-penalized estimator. The second contribution is computational. We provide an exact computation approach for $\ell _{0}$-PQR through the method of mixed integer optimization. We also develop a first-order approximation algorithm for solving large scale $\ell _{0}$-PQR problems. Through Monte Carlo simulations and a real-data application, we find that both $\ell _{0}$ -PQR and $\ell _{0}$-CQR perform fairly well and produce much sparser solutions than $\ell _{1}$-PQR does and also outperform the adaptive Lasso and non-convex penalized quantile regression approaches. Our theoretical and numerical results suggest that the $\ell _{0}$-based approaches are worthy competitors to the $\ell _{1}$-based and non-convex penalized estimation methods in sparse quantile regression. Recently, Hazimeh:Mazumder:2020 developed fast computational methods for $\ell _{0}$-penalized least squares with an additional $\ell _{1}$- or $\ell _{2}$-penalty term. It is an interesting future research topic to extend their approach to quantile regression and investigate its statistical properties.
We are indebted to the editor, Elie Tamer, an associate editor and two anonymous referees for constructive comments and suggestions. We would like to thank Roger Koenker, Rahul Mazumder, Guillaume Pouliot and participants at 2019 Optimization-Conscious Econometrics Conference in Chicago, 2020 Econometric Society World Congress, and 2021 ASSA Annual Meeting for helpful comments. We are also grateful to Rahul Mazumder for providing us his code. This work was supported in part by the Ministry of Science and Technology, Taiwan (MOST109-2410-H-001-027-MY2), Academia Sinica (AS-CDA-106-H01), the European Research Council (ERC-2014-CoG-646917-ROMIA), and the UK Economic and Social Research Council (ESRC) through research grant (ES/P008909/1) to the CeMMAP.