Quantile Regression Under Memory Constraint

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.

89,018 characters · 16 sections · 53 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

frontmatter\runtitle{Quantile Regression Under Memory Constraint} \begin{aug} , , \thankstext{t1}{Research supported by NSFC, Grant No.11431006 and 11690013, the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning, Youth Talent Support Program, 973 Program (2015CB856004) and a grant from Australian Research Council.} \runauthor{X. Chen, W. Liu and Y. Zhang} \address{X. Chen\\ Y. Zhang\\ Information, Operations and Management Sciences\\ Stern School of Business\\ New York University\\ New York, New York 10012\\ USA\\ \printead{e1}\\ \phantom{E-mail:\ }\printead*{e3}} \address{W. Liu\\ Department of Mathematics\\ Institute of Natural Sciences and Moe-LSC\\ Shanghai Jiao Tong University\\ Shanghai\\ China\\ \printead{e2}} \end{aug} \begin{abstract} This paper studies the inference problem in quantile regression (QR) for a large sample size $n$ but under a limited memory constraint, where the memory can only store a small batch of data of size $m$. A natural method is the na\"ive divide-and-conquer approach, which splits data into batches of size $m$, computes the local QR estimator for each batch, and then aggregates the estimators via averaging. However, this method only works when $n=o(m^2)$ and is computationally expensive. This paper proposes a computationally efficient method, which only requires an initial QR estimator on a small batch of data and then successively refines the estimator via multiple rounds of aggregations. Theoretically, as long as $n$ grows polynomially in $m$, we establish the asymptotic normality for the obtained estimator and show that our estimator with only a few rounds of aggregations achieves the same efficiency as the QR estimator computed on all the data. Moreover, our result allows the case that the dimensionality $p$ goes to infinity. The proposed method can also be applied to address the QR problem under distributed computing environment (e.g., in a large-scale sensor network) or for real-time streaming data. \end{abstract} \begin{keyword}[class=MSC] \kwd[Primary ]{62F12} \kwd[; secondary ]{62J02} \end{keyword} \begin{keyword} \kwd{Quantile regression} \kwd{sample quantile} \kwd{divide-and-conquer} \kwd{distributed inference} \kwd{streaming data} \end{keyword}

Introduction

The development of modern technology has enabled data collection of unprecedented size, which leads to large-scale datasets that cannot be fit into memory or are distributed in many machines over limited memory. For example, the memory of a personal computer only has a storage size in GBs while the dataset on the hard disk could have a much larger size. In addition, in a sensor network, each sensor is designed to collect and store a limited amount of data, and computations are performed via communications and aggregations among sensors (see, e.g., Wang17sensor). Other examples include high-speed data streams that are transient and arrive at the processor at a high speed. In online streaming computation, the memory is usually limited as compared to the length of the data stream Gama2013,zhangqi2007. Under memory constraints in all these scenarios, classical statistical methods, which are developed under the assumption that the memory can fit all the data, are no longer applicable; thus, many estimation and inference methods need to be re-investigated. For example, suppose that there are $n$ samples for some very large $n$, a fundamental question in data analysis is as follows:

{5pt}

boxedminipage{.9\textwidth} How to calculate the sample quantiles of $n$ samples when the memory can only store $m$ samples with $n\gg m$?

As one of the most popular interview questions from high-tech companies, this problem has attracted much attention from computer scientists over the last decade; see manku1998, greenwald2004power, zhangqi2007, guha2009 and the references therein. However, this is mainly a computation problem with a fixed dataset, which does not involve any statistical modeling.

Motivated by this sample quantile calculation problem, we study a more general problem of quantile regression (QR) under memory constraints. Quantile regression, which models the conditional quantile of the response variable given covariates, finds a wide range of applications to survival analysis (e.g., Wang:14:Sur,Xu:16:Sur), health care (e.g., Sherwood:13,Luo:13:gap), and economics (e.g., Belloni:17). In the classical QR model, assume that there are $n$ i.i.d. samples $\{(\boldsymbol{X}_i, Y_i)\}$ from the following model

equation[equation omitted — 135 chars of source]

where $\boldsymbol{X}_i'=(1, X_{i1}, \ldots, X_{ip})$ is the random covariate vector with the dimension $p+1$ drawn from a common population $\boldsymbol{X}$. The error $\epsilon_i$ is an unobserved random variable satisfying $\mathbb{P}(\epsilon_i\leq 0|\boldsymbol{X}_i)=\tau$ for some specified $0<\tau<1$ (known as the quantile level). In other words, $\boldsymbol{X}_i'\boldsymbol{\beta}(\tau)$ is the $\tau$-th quantile of $Y_i$ given $\boldsymbol{X}_i$. When all the $n$ samples can be fit into memory, one can estimate $\boldsymbol{\beta}(\tau)$ via the classical QR estimator koenker2005quantile,

equation[equation omitted — 200 chars of source]

where $\rho_{\tau}(x)=x(\tau-I\{x\leq 0\})$ is the asymmetric absolute deviation function (a.k.a. check function) and $I(\cdot)$ is the indicator function. However, when samples are distributed across many machines or the sample size $n$ is extremely large and thus the samples cannot be fit into memory, it is natural to ask the following question:

{5pt}

boxedminipage{.9\textwidth} How to estimate and conduct inference about $\boldsymbol{\beta}(\tau)$ when the memory can only store $m$ samples with $n\gg m$?

The divide-and-conquer (DC), as one of the most important algorithms in computer science, has been commonly adopted to deal with this kind of big data challenge. We below describe a general DC algorithm for statistical estimation. Specifically, we split the data indices $\{1,2,...,n\}$ into $N$ subsets $\mathcal{H}_{1},...,\mathcal{H}_{N}$ with equal size $m$ and $N=n/m$. Correspondingly, the entire dataset $\{Y_{i},\boldsymbol{X}_{i},1\leq i\leq n\}$ is divided into $N$ batches $\mathcal{D}_{1},\ldots,\mathcal{D}_{N}$, where $\mathcal{D}_k=\{Y_i,\boldsymbol{X}_i,i\in \mathcal{H}_k\}$ for $1 \leq k \leq N$. By swapping each batch of data $\mathcal{D}_k$ into the memory, one constructs a low dimension statistic $T_{k}=g_k(\mathcal{D}_{k})$ for $\mathcal{D}_k$ with some function $g_k(\cdot)$. Then the estimator $\widehat{\boldsymbol{\beta}}$ is obtained by the aggregation of $\{T_{k}\}_{k=1}^N$ (i.e., $\widehat{\boldsymbol{\beta}}=G(T_{1},\ldots,T_{N})$ for some aggregation function $G(\cdot)$). In recent years, this DC framework has been widely adopted in distributed statistical inference (see, e.g., li2013statistical,chen2014split,battey2015distributed,zhao2016partially,shi2016massive,banerjee2016divide,volgushev2017distributed and Section (ref) for detailed descriptions).

In addition to memory-constrained estimation on a single machine (where the size of the dataset is much larger than memory size), another natural situation for using the DC framework comes from the application of large-scale wireless sensor networks (see, e.g., shrivastava2004medians,greenwald2004power, wain2006, huang2011sampling,Wang17sensor). In a sensor network with $N$ sensors, the data are collected and stored in different sensors. Moreover, due to limited energy carried by sensors, communication cost is one of the main concerns in data aggregation. The samples are not transferred to the base station or neighboring sensors directly. Instead, each sensor first summarizes the samples into a low dimensional statistic $T_{k}$, which can be transferred with a low communication cost. Figure (ref) visualizes a typical sensor network with data flows as a routing tree with the base station as the root. An internal sensor node in the $i$-th layer receives statistics $T_{(\cdot)}$ from its children nodes in $(i+1)$-th layer, and then combines received statistics with its own $T_{(\cdot)}$ and sends the resulting statistic to its parent in the $(i-1)$-th layer. The final estimator in the base station (or central node) can be computed by $\widehat{\boldsymbol{\beta}}=G(T_{1},\ldots,T_{N})$.

figure[figure omitted — 478 chars of source]

A critical problem in statistical DC framework is how to construct local statistics $g_k(\cdot)$ and aggregation function $G(\cdot)$. In many existing studies, a typical choice of $g_k(\cdot)$ is to use the same estimator as the one designed for the estimation from the entire data. For example, in QR, one may choose

eqnarray[eqnarray omitted — 236 chars of source]

and a simple averaging function $G(T_{1},\ldots,T_{N})= \frac{1}{N} \sum_{k=1}^N T_k$, where $T_k=g_k(\mathcal{D}_{k})$ is the local statistic. We call this kind of DC methods the na\"ive-DC algorithm where the estimator is denoted by $\widehat{\boldsymbol{\beta}}_{ndc}=\sum_{k=1}^{N}\widehat{\boldsymbol{\beta}}_{QR,k}/N$. Despite its popularity, the na\"ive-DC algorithm might fail when $n/m$ is large. For example, in a special case of quantile estimation (i.e., $p=0$), it is straightforward to show that $\sqrt{n}|\widehat{\beta}_{ndc}-\beta(\tau)|\rightarrow \infty$ in probability when $n/m^{2}\rightarrow \infty$ (see Theorem (ref) in Appendix). Similar phenomenon occurs for general $p$, see volgushev2017distributed. In fact, the local QR estimators $\widehat{\boldsymbol{\beta}}_{QR,k}$ are biased estimators with the bias $O(1/m)$. Although the averaging aggregation is useful for variance reduction, it is unable to reduce the bias, which makes the na\"ive-DC fail when $n$ is large as compared to $m$. In the DC framework, bias reduction in $T_{k}$ is more critical than the variance reduction. This is a fundamental difference from many classical inference problems that require to balance the variance and bias (cf. nonparametric estimation). Furthermore, in the na\"ive-DC algorithm, we need to solve $N=n/m$ optimization problems, which could be computationally expensive.

The deficiency of the na\"ive-DC approach calls for a new DC scheme to achieve the following two important goals in distributed inference.

enumerate• The obtained estimator $\widehat{\boldsymbol{\beta}}$ should achieve the same statistical efficiency as merging all the data together under a weak condition on the sample size $n$ as a function of $m$. More precisely, it is desirable to remove the constraint $n=o(m^2)$ in na\"ive-DC so that the procedure can be applied to situations such as large-scale sensor networks, where the number of sensors $N=n/m$ is excessively large. • The second goal is on the computational efficiency. For example, the na\"ive-DC requires solving a non-smooth optimization for computing the local statistic $T_{k}=g_k(\mathcal{D}_k)$. Since there is no explicit formula for $g_k(\cdot)$, the computation is quite heavy (especially considering each sensor has a limited computational power).

This paper develops new constructions of $g_k(\cdot)$ and $G(\cdot)$ with a multi-round aggregation scheme in Algorithm (ref), which simultaneously achieves the two goals referred above. Our method is applicable to both scenarios of small memory on a single machine and large-scale sensor networks. Instead of using the local QR estimation as in (ref), we adopt a smoothing technique in the literature (see, e.g., horowitz1998bootstrap, whang2006smoothed, wang2012corrected, pang2012variance, wu2015smoothed), and propose a new estimator for QR called linear estimator of QR (LEQR), which serves as the cornerstone for our DC approach. Our LEQR has an explicit formula in the form of direct sums of the transformation of $\{Y_{i}, \boldsymbol{X}_i\}$, which is quite different from the optimization-based QR estimator in (ref). It is also worthwhile noting that the linearity is the most desirable property in the DC framework for both theoretical development and computation efficiency.

The high-level description of the proposed multi-round DC approach is provided as follows. Our method only needs to compute an initial QR estimator $\widehat{\boldsymbol{\beta}}_0$ based on a small part of samples (e.g., $\mathcal{D}_{1}$). Based on $\widehat{\boldsymbol{\beta}}_0$, for each batch of data $\mathcal{D}_k$, we compute local statistics $\{T_k\}$ using the proposed LEQR. The local statistics are in a simple form of weighted sums of $\boldsymbol{X}_i$ and $\boldsymbol{X}_i \boldsymbol{X}_i'$. The aggregation function is constructed by adding up the local statistics and then solving a linear system, which gives the first-round estimator $\widehat{\boldsymbol{\beta}}^{(1)}$. Now, we can repeat this DC algorithm using $\widehat{\boldsymbol{\beta}}^{(1)}$ as the initial estimator. After $q$ iterations, we denote our final estimator by $\widehat{\boldsymbol{\beta}}^{(q)}$. Theoretically, under some conditions on the growth rate of $p \rightarrow \infty$ as a function $m$ and $n$, we first establish the Bahadur representation of $\widehat{\boldsymbol{\beta}}^{(q)}$ and show that the Bahadur remainder term achieves a nearly optimal rate (up to a log-factor) when $q$ satisfies some mild conditions (see (ref) and Theorem (ref)). Furthermore, as long as $n=o(m^{A})$ for some constant $A$, the final estimator $\widehat{\boldsymbol{\beta}}^{(q)}$ achieves the same asymptotic efficiency as the QR estimator (ref) computed on the entire data (see Theorem (ref)).

The new DC approach is particularly suitable for QR in sensor networks in which communication cost is one of the major concerns. The proposed procedure only requires $O(p^2)$ bits communication between any two sensors. We also highlight two other important applications of our method:

enumerate• Our method can be adapted to make inference for online streaming data, which arrives at the processor at a high speed. Our method provides a sequence of successively refined estimators of $\boldsymbol{\beta}(\tau)$ for streaming data and can deal with an arbitrary length of data stream. The online quantile estimation problem (which is a special case of QR with $p=0$) for streaming data has been extensively studied in computer science literature (see, e.g., munro1980, zhangqi2007, guha2009, wang2003exp). However, these works mainly focus on developing approximations to the sample quantile, which are insufficient to obtain limiting distribution results for the purpose of inference. We extend the quantile estimation to the more general QR problem and provide the asymptotical normality result for the proposed online estimator (see Theorem (ref)). • Our method also serves as an efficient optimization solver for classical QR on a single machine. As compared to the standard interior-point method for solving the QR estimator in (ref) that requires the computational complexity of $O(n^{1.25} p^3 \log n)$ portnoy1997, our approach requires $O(m^{1.25} p^3\log m + np^2+p^{3})$ since it only solves an optimization on a small batch of data for construction initial estimator. Therefore, our method is computationally more efficient.

We will illustrate them in Section (ref) after we provide the detailed description of the method.

Organization and notations

The rest of the paper is organized as follows. In Section (ref), we review the related literature on recent works on distributed estimation and inference. Section (ref) describes the proposed inference procedure for QR under memory constraints. Section (ref) presents the theoretical results. In Section (ref), we demonstrate the performance of the proposed inference procedure by simulated experiments, followed by conclusions in Section (ref). The proofs and additional experimental results are provided in Appendix.

In the QR model in (ref), let $F(\cdot|\boldsymbol{x})$ and $f(\cdot|\boldsymbol{x})$ denote the CDF and PDF of $\epsilon$ conditioning on $\boldsymbol{X}=\boldsymbol{x}$, respectively, throughout the paper. Then, for any $\boldsymbol{x}$, we have $F(0|\boldsymbol{x})=\tau$. For two sequences of real numbers $f(n)$ and $g(n)$, let $f(n)=\Omega (g(n))$ denote that $f$ is bounded below by $g$ (up to constant factor) asymptotically. For a set of random variables $X_n$ and a corresponding set of constants $a_n$, $X_n=O_p(a_n)$ means that $X_n/a_n$ is stochastically bounded and $X_n=o_p(a_n)$ means that $X_n/a_n$ converges to zero in probability as $n$ goes to infinity. For a real number $c$, we will use $\lfloor c \rfloor$ to denote largest integer less than or equal to $c$. Finally, denote the Euclidean norm for a vector $\boldsymbol{x}\in\mathbb{R}^p$ by $\|\boldsymbol{x}\|_2$, and denote the spectral norm for a matrix $\boldsymbol{X}$ by $\|\boldsymbol{X}\|$.

Related works

The explosive growth of data presents new challenges for many classical statistical problems. In recent years, a large body of literature has emerged on studying estimation and inference problems under memory constraints or in distributed environments (please see the references described below as well as other works such as kleiner2014scalable,wang2013parallelizing,Wang:15:MCMC). Examples include but are not limited to density parameter estimation li2013statistical, generalized linear regression with non-convex penalties chen2014split, kernel ridge regression zhang2015divide, high-dimensional sparse linear regression lee2017communication, high-dimensional generalized linear models battey2015distributed, semi-parametric partial linear models zhao2016partially, QR processes volgushev2017distributed, $M$-estimators with cubic rate shi2016massive, and some non-standard problems where rates of convergence are slower than $n^{1/2}$ and limit distributions are non-Gaussian banerjee2016divide. All these results rely on averaging, where the global estimator is the average of the local estimators computed on each batch of data. For the averaging estimators to achieve the same asymptotic distribution for inference as pooling the data together, it usually requires the number of batches (i.e., the number of machines) to be $o(m)$ (i.e., $n=o(m^2)$). However, in some applications such as sensor networks and streaming data, the number of batches can be large.

To address the challenge, instead of using one-shot aggregation via averaging, recent works by jordan2016communication and wang2016efficient proposed iterative methods with multiple rounds of aggregations, which relaxes the condition $n=o(m^2)$. These methods have been applied to $M$-estimator and Bayesian inference. Their framework is based on an approximate Newton algorithm Shamir:14 and thus requires the twice-differentiability of the loss function. However, the QR loss is non-differentiable and thus their approach cannot be utilized. More detailed discussions on these two works will be provided in Remark (ref). wain2006 proposed a multi-round decentralized quantile estimation algorithm under a restrictive communication-constrained setup. However, their method cannot be applied to solve QR problems.

There is a large body of literature on estimation and inference for QR and its variant (e.g., censored QR). We will not be able to provide a detailed survey here and we refer the readers to koenker2005quantile,quantilehandbook for more background knowledge and recent development of QR. However, it is worth noting that the smoothing idea in QR literature has been adopted for developing our linear estimator for QR (LEQR in Section (ref)), which serves as the cornerstone of our method. The idea of smoothing the non-smooth QR objective goes back to horowitz1998bootstrap, where he studied the bootstrap refinement in inference in quantile models. Since the smoothing idea overcomes the difficulty in higher-order expansion of the scores associated with the QR objective, it plays an important role in solving various QR problems. For example, wang2012corrected and wu2015smoothed proposed different smoothed objectives to determine the corrected scores under the presence of covariate measurement errors for QR and censored QR. Galvao16smooth proposed a fix-effects estimator for the smoothed QR in linear panel data models and derived the corresponding limiting distribution. pang2012variance proposed an induced-smoothing idea for estimating the variance of inverse-censoring-probability weighted estimator in Bang02median. whang2006smoothed considered the problem of inference using the empirical likelihood mehod for QR and demonstrated that the smoothed empirical likelihood can help achieve higher-order refinements (i.e., $O(n^{-1})$ of the coverage error). The smoothing idea can also facilitate the computation, especially for the first-order optimization methods (see, e.g., zheng2011gradient). We also adopt the smoothing idea for constructing our LEQR estimator (see Section (ref)), which heavily relies on the first-order optimality condition of the objective. Instead of using the smoothing technique for computing a one-stage estimator as in existing literature, our use of the smoothing technique enables successive refinement of the LEQR estimator (see Propositions (ref) and (ref)).

Methodology

In this section, we introduce the proposed method. We start with a new linear type estimator for quantile regression, which serves as an important building block for our inference approach.

A linear type estimator of quantile regression

We first propose a linear type estimator for quantile regression, which is named as LEQR (Linear Estimator for Quantile Regression). Recall the classical quantile regression estimator from (ref). Note that for the quantile regression, the loss function $\rho_{\tau}(x)=x(\tau-I\{x\leq 0\})=x(I\{x> 0\}+\tau-1)$ is non-differentiable. Using the smoothing idea (see the literature surveyed in Section (ref)), we approximate the indicator factor $I\left\{x>0\right\}$ with a smooth function $H(x/h)$, where $h \rightarrow 0$ is the bandwidth. With this approximation, we replace $\rho_{\tau}(x)$ in quantile regression by $K_h(x)=x(H(x/h)+\tau-1)$ and define

eqnarray[eqnarray omitted — 187 chars of source]

Now the right hand side in ((ref)) is differentiable and we note that $\frac{\mathrm{d} K_h(x)}{\mathrm{d} x}= H(x/h)+\tau-1+(x/h)H'(x/h)$. Here the function $H(u)$ is a smooth approximation of the indicator factor $I\{x>0\}$ satisfying $H(u)=1$ when $u\geq 1$ and $H(u)=0$ when $u\leq -1$ (see more details on the condition of $H(\cdot)$ in Condition (C3) on page (ref)). For example, one may choose $H(u)$ for $u \in[-1,1]$ to be the integral of a smooth kernel function with support on $[-1,1]$, e.g., a biweight (or quartic) kernel (see (ref)). We further note that (ref) is a non-convex optimization and thus it is difficult to compute the minimizer $\widehat{\boldsymbol{\beta}}_h$. However, it is not a concern since (ref) is only introduced for the motivation purpose. The proposed LEQR estimator, which is explicitly defined in (ref) below, does not require to solve (ref).

Since $H(\cdot)$ is a smooth function, by the first-order optimality condition, the solution $\widehat{\boldsymbol{\beta}}_{h}$ in (ref) satisfies (see Theorem 2.6 in Beck:14)

equation[equation omitted — 305 chars of source]

From (ref) we can express $\widehat{\boldsymbol{\beta}}_h$ by {

eqnarray[eqnarray omitted — 480 chars of source]

}

However, there is no closed-form expression of $\widehat{\boldsymbol{\beta}}_{h}$ from this fixed point equation. Instead of using $\widehat{\boldsymbol{\beta}}_{h}$ on the right hand side of ((ref)), we replace $\widehat{\boldsymbol{\beta}}_{h}$ by a consistent initial estimator $\widehat{\boldsymbol{\beta}}_{0}$, which leads to the proposed LEQR:

eqnarray[eqnarray omitted — 474 chars of source]

Allowing the choice of initial estimator $\widehat{\boldsymbol{\beta}}_{0}$ is crucial for our inference procedure under the memory constraint. While it is difficult to solve the fixed point equation in (ref) when all data cannot be loaded into memory, we are still able to compute an initial estimator using only a batch of samples, e.g., $\mathcal{D}_{1}$. Given the initial estimator $\widehat{\boldsymbol{\beta}}_{0}$, the LEQR in (ref) only depends on sums of $\boldsymbol{X}_i$ and $\boldsymbol{X}_i\boldsymbol{X}_i'$, which can be easily implemented via a divide-and-conquer scheme (see Section (ref) for details). Comparing to the na\"ive-DC method that needs to solve $N$ optimization problems on each batch of data, LEQR only needs to solve one optimization, and thus is computationally more efficient. Moreover, our estimator in (ref) is essentially solving a linear equation system, and we do not need to explicitly compute a matrix inversion. There are a number of efficient methods for solving a linear system numerically, such as conjugate gradient method hestenes1952methods and stochastic variance reduced gradient method johnson2013accelerating. In our simulation studies, we use the conjugate gradient method, and due to space limitations, more detailed explanations of this method are provided in Appendix (ref).

Divide-and-conquer LEQR

Based on LEQR, we now introduce a divide-and-conquer LEQR for estimating $\boldsymbol{\beta}(\tau)$. For each batch of data $\mathcal{D}_k$ for $1\leq k\leq N$, let us define the following quantities. {

eqnarray[eqnarray omitted — 466 chars of source]

\normalfont } The inference procedure is presented in Algorithm (ref). Using our theory in Theorem (ref) and (ref), for the $g$-th iteration, we choose the bandwidth $h=h_g= \max(\sqrt{p/n},(p/m)^{2^{g-2}})$ for $1\leq g\leq q$.

algorithm[algorithm omitted — 1,915 chars of source]

Algorithm (ref) can not only be used under the memory constraint, but also be applied to distributed setting, to reduce computational cost for classical quantile regression, and to deal with streaming data. We illustrate these important applications as follows.

enumerate• {\em Quantile regression in large-scale sensor networks.} Algorithm (ref) is directly applicable to distributed sensor network with $N$ sensors, where each sensor collects a batch of data $\mathcal{D}_k$ (for $k=1,\ldots,N$). The base station first broadcasts the initial estimator $\widehat{\boldsymbol{\beta}}_{0}$ computed on $\mathcal{D}_1$ in (ref) to all sensors (see Figure (ref) for an illustration). Then each sensor computes $(\boldsymbol{U}_k, \boldsymbol{V}_k)$ locally, which will be transferred from bottom to the base station. In particular, each sensor $k$ only keeps the summation of $(\boldsymbol{U}_{\cdot}, \boldsymbol{V}_{\cdot})$ from all its children nodes and its own $(\boldsymbol{U}_k, \boldsymbol{V}_k)$ and then transfers the summed statistics to its parent. After receiving $\left(\sum_{k=1}^{N}\boldsymbol{U}_{k}, \sum_{k=1}^{N}\boldsymbol{V}_{k}\right)$, the base station will compute $\widehat{\boldsymbol{\beta}}^{(g)}$ for $g=1$ in (ref). Then this distributed procedure can be repeated for $g=2,\ldots, q$. • {\em Computational reduction of quantile regression.} Algorithm (ref) can also be utilized as an efficient solver for classical quantile regression on a single machine. For the ease of illustration, let us assume the quantile regression estimator in (ref) (or (ref)) is solved by the standard interior-point method. Using the standard interior-point method, the initial estimator requires the computational complexity $O(m^{1.25} p^3\log m)$ portnoy1997. The computation of $\sum_{k=1}^N \boldsymbol{U}_k$ and $\sum_{k=1}^N \boldsymbol{V}_k$ require $O(np)$ and $O(np^{2})$, respectively. Therefore, the computational complexity for $\widehat{\boldsymbol{\beta}}^{(q)}$ is at most $O(m^{1.25} p^3\log m + np^2 +p^{3})$, where $O(p^3)$ comes from the inversion of $\sum_{k=1}^N \boldsymbol{V}_k$. This greatly saves the computational cost as compared to the interior-point method for computing quantile regression estimator on the entire data, which requires a complexity of $O(n^{1.25} p^3\log n)$. \begin{algorithm}[!t] \caption{{ Online LEQR}} {Initialization: For the first $m$ samples $\{Y_{i}, \boldsymbol{X}_{i}, -m+1\leq i\leq 0\}$, compute the initial standard quantile regression estimator $\widehat{\boldsymbol{\beta}}[{0}]$. Then constructs the initial $(\boldsymbol{U}(0), \boldsymbol{V}(0))$ based on $\widehat{\boldsymbol{\beta}}[{0}]$ according to (ref) with $h=(p/m)^{1/2}$. Parameter Setup: Define $a_{2k-1}=2^{k-1}+1/2$ and $a_{2k}=2^{k-1}+3/4$ for $k\geq 1$ and $a_{0}=-\infty$. Let $s_l=\lfloor m^{a_{l-1}}\rfloor+1$ and $r_l=\lfloor m^{a_{l}}\rfloor$ for any $l\geq 1$ and let $r_0=0$. Define a sequence of bandwidths $h_{1}=(p/m)^{1/2}$ and $h_{l}=(p/m^{a_{l-1}})^{1/2}$ for $l\geq 2$. } \begin{algorithmic}[1] \FOR{each interval $l=1,2\ldots, $} \FOR{indices in the interval $l$, $j=s_l, s_l+1, \ldots, r_l$} \STATE Receive an online sample $(Y_j, \boldsymbol{X}_j)$. \STATE Compute \begin{align*} \widetilde{\boldsymbol{U}}(j) & =\boldsymbol{X}_j\Bigl\{H\Big(\dfrac{Y_j-\boldsymbol{X}'_j\widehat{\boldsymbol{\beta}}[r_{l-1}])}{h_{l}}\Big)+\tau-1+\dfrac{Y_j}{h_{l}}H'\Big(\dfrac{Y_j-\boldsymbol{X}'_j\widehat{\boldsymbol{\beta}}[r_{l-1}]}{h_{l}}\Big)\Bigr\} \\ \widetilde{\boldsymbol{V}}(j) & =\boldsymbol{X}_j\boldsymbol{X}_j'\frac{1}{h_{l}}H'\Big(\dfrac{Y_j-\boldsymbol{X}'_j\widehat{\boldsymbol{\beta}}[r_{l-1}]}{h_{l}}\Big), \end{align*} where $\widehat{\boldsymbol{\beta}}[r_{l-1}]$ is the estimator computed upto the end of $(l-1)$-th interval. \STATE Update $(\boldsymbol{U}(j),\boldsymbol{V}(j))$ by \begin{align*} \boldsymbol{U}(j) & \triangleq \sum_{i=s_l}^j \widetilde{\boldsymbol{U}}(j) = \begin{cases} \widetilde{\boldsymbol{U}}(j) & \quad if\; j=s_l\\ \boldsymbol{U}(j-1)+ \widetilde{\boldsymbol{U}}(j) & \quad if\; j> s_l \end{cases}\\ \boldsymbol{V}(j) & \triangleq \sum_{i=s_l}^j \widetilde{\boldsymbol{V}}(j) = \begin{cases} \widetilde{\boldsymbol{V}}(j) & \quad if\; j=s_l\\ \boldsymbol{V}(j-1)+ \widetilde{\boldsymbol{V}}(j) & \quad \text{if}\; j> s_l \end{cases} \end{align*} \STATE Calculate $(\boldsymbol{U}(r_{l-1})+\boldsymbol{U}(j), \boldsymbol{V}(r_{l-1})+\boldsymbol{V}(j))$ and compute the online quantile regression estimator at time $j$: \begin{eqnarray} \widehat{\boldsymbol{\beta}}[j]=(\boldsymbol{V}(r_{l-1})+\boldsymbol{V}(j))^{-1}(\boldsymbol{U}(r_{l-1})+\boldsymbol{U}(j)). \end{eqnarray} \STATE Remove $(\widetilde{\boldsymbol{U}}(j),\widetilde{\boldsymbol{V}}(j)), (\boldsymbol{U}(j-1),\boldsymbol{V}(j-1))$ from the memory. \ENDFOR \STATE Remove $(\boldsymbol{U}(r_{l-1}),\boldsymbol{V}(r_{l-1}))$ from the memory and only keep $\widehat{\boldsymbol{\beta}}[r_{l}]$, $(\boldsymbol{U}(r_{l}), \boldsymbol{V}(r_{l}))$ in the memory. \ENDFOR \end{algorithmic} \end{algorithm} • {\em Online quantile regression for streaming data.} To deal with online streaming data, it is critical to design a one-pass algorithm since streaming data are transient. To this end, based on Algorithm (ref), we develop a new one-pass algorithm (see Algorithm (ref)) that provides a sequence of successively refined estimators. For streaming data, we divide the data into intervals $\{(s_l,r_l)\}_{l=1}^\infty$. The starting and ending positions of the $l$-th interval are chosen as \[ s_l=\lfloor m^{a_{l-1}}\rfloor+1,\text{ and }r_l=\lfloor m^{a_{l}}\rfloor \] for $l\geq 1$ where $a_{2k-1}=2^{k-1}+1/2$ and $a_{2k}=2^{k-1}+3/4$ for $k\geq 1$ and $a_{0}=-\infty$. The intervals are chosen to ensure that the sample size of the $l$-th interval $n_l$ is approximately $n_{l-2}^2$. As we will show in the proof of Theorem (ref), if an initial estimator is computed from $m$ samples, there will be no improvement of the online LEQR estimator after $m^2$ fresh samples and this is the ending point of an interval where we compute a new initial estimator. In Algorithm (ref), for each interval $l$ and each $j$ such that $s_l \leq j\leq r_l$, the memory only maintains $\widehat{\boldsymbol{\beta}}[r_{l-1}]$, $(\boldsymbol{U}(r_{l-1}), \boldsymbol{V}(r_{l-1}))$, $(\boldsymbol{U}(j), \boldsymbol{V}(j))$, and $\widehat{\boldsymbol{\beta}}[j]$. We note that $(\boldsymbol{U}(r_{l-1}), \boldsymbol{V}(r_{l-1}))$ are the weighted sums of $\boldsymbol{X}_i$ and $\boldsymbol{X}_i\boldsymbol{X}_i'$ for $s_{l-1}\leq i \leq r_{l-1}$ and $(\boldsymbol{U}(j), \boldsymbol{V}(j))$ can be easily updated from $(\boldsymbol{U}(j-1), \boldsymbol{V}(j-1))$ in an online fashion. Therefore, except for an $O(m)$ space for deriving the initial estimator $\widehat\boldsymbol{\beta}[0]$, the online LEQR only requires $O(p^2)$ memory, which is independent on $n$. In Theorem (ref), we will show that the online LEQR algorithm achieves the same statistical efficiency for any $l$ and $j$ as the standard quantile regression estimator when merging all the streaming data together. Also, the asymptotic normality of $\sqrt{m+j}(\widehat{\boldsymbol{\beta}}[j]-\boldsymbol{\beta}(\tau))$ holds uniformly in $1\leq j\leq m^{A}$ for any constant $A>0$ (note that $m+j$ is the sample size until time $j$).

Theoretical results

In this section, we provide a Bahadur representation of $\widehat{\boldsymbol{\beta}}^{(q)}$, based on which we derive the asymptotic normality result for DC LEQR $\widehat{\boldsymbol{\beta}}^{(q)}$ and the online LEQR $\widehat{\boldsymbol{\beta}}[j]$. We also discuss adaptive choices of bandwidth and extensions to heterogeneous settings.

Asymptotics for DC LEQR

We note that (ref) can be equivalently written as:

equation[equation omitted — 130 chars of source]

where

eqnarray[eqnarray omitted — 469 chars of source]

We first state some regularity conditions for our theoretical development and then give Propositions (ref)-(ref) on the expansions of $\boldsymbol{A}_{n,h}$ and $\boldsymbol{D}_{n,h}$. We assume the model (ref) with $n$ i.i.d. samples $\{(\boldsymbol{X}_i, Y_i)\}$ (the non-i.i.d. case will be discussed in the next subsection). Let $f(\cdot|\boldsymbol{X})$ be the conditional density function of the noise $\epsilon$ given $\boldsymbol{X}$. Further, we define $ \boldsymbol{D}=\mathbb{E}(\boldsymbol{X}\boldsymbol{X}'f(0|\boldsymbol{X})). $ {

enumerate• The conditional density function $f(\cdot|\boldsymbol{X})$ is Lipschitz continuous (i.e., $|f(x|\boldsymbol{X})-f(y|\boldsymbol{X})|\leq L|x-y|$ for any $x,y\in\mathbb{R}$ and some constant $L>0$). Also, assume that $0<c_{1}\leq \lambda_{\min}(\boldsymbol{D})\leq\lambda_{\max}(\boldsymbol{D})\leq c_{2}<\infty$ for some constants $c_{1},c_{2}$. • Let the smoothing function $H(\cdot)$ satisfy $H(u)=1$ if $u\geq 1$ and $H(u)=0$ if $u\leq -1$. Further, suppose $H(\cdot)$ is twice differentiable and its second derivative $H^{(2)}(\cdot)$ is bounded. Moreover, we assume the bandwidth $h=o(1)$. • Assume that $p=o\left(nh/(\log n)\right)$ and $\sup_{\|\boldsymbol{\theta}\|_{2}=1}\mathbb{E} e^{\eta(\boldsymbol{\theta}'\boldsymbol{X})^{2}}<\infty$ for some $\eta>0$. • Assume for some $\kappa>0$, $p=o\left((n^{1-4\kappa}h/\log n)^{1/3}\right)$. Suppose that $\sup_{j}\mathbb{E} |X_{1,j}|^{a}<\infty$ for some $a\geq 2/\kappa$ and $\sup_{\|\boldsymbol{\theta}\|_2=1}\mathbb{E} (\boldsymbol{\theta}{'}\boldsymbol{X})^{4}<\infty$.

Condition (C1) contains a standard eigenvalue condition related to covariates $\boldsymbol{X}$ and the smoothness of the conditional density function $f$. Condition (C2) is a mild condition on $H$ for smooth approximation.

Condition (C3) and (C3$^*$) illustrate the relationship between the dimension $p$ and sample size $n$ and the moment condition on covariates $\boldsymbol{X}$. Either one of them leads to our theoretical results in Propositions (ref)-(ref). As compared to Condition (C3$^*$), Condition (C3) is weaker in terms of the relationship of $p$ and $n$, but requires a stronger moment condition on $\boldsymbol{X}$.

Under these conditions, we have the following Proposition (ref) and (ref) for the asymptotic behavior of $\boldsymbol{A}_{n,h}$ and $\boldsymbol{D}_{n,h}$.

propositionSuppose we have an initial estimator $\hat{\boldsymbol{\beta}}_{0}$ with $\|\hat{\boldsymbol{\beta}}_{0}-\boldsymbol{\beta}(\tau)\|_2=O_{\mathbb{P}}(a_{n})$ with $a_{n}=O(h)$. Assume that (C1), (C2), and (C3) (or (C3$^{*}$)) hold. We have \begin{eqnarray*} \Bigl\|\boldsymbol{A}_{n,h}-\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{X}_{i}(I\{\epsilon_{i}\geq 0\}+\tau-1)\Bigr\|_2=O_{\mathbb{P}}\Big{(}\sqrt{\frac{ph\log n}{n}}+a^{2}_{n}+h^{2}\Big{)}. \end{eqnarray*}

propositionSuppose the conditions in Proposition (ref) hold. We have \begin{eqnarray*} \|\boldsymbol{D}_{n,h}-\boldsymbol{D}\|=O_{\mathbb{P}}\Big{(}\sqrt{\frac{p\log n}{nh}}+a_{n}+h\Big{)}. \end{eqnarray*}

Combining Propositions (ref) and (ref) with (ref) and with some algebraic manipulations, we have

eqnarray[eqnarray omitted — 193 chars of source]

with

eqnarray[eqnarray omitted — 159 chars of source]

By choosing the bandwidth $h$ shrinking at an appropriate rate, the dominating term of (ref) is $a_n^2$, which means that one round of aggregation enables a refinement of the estimator with its bias reducing from $a_n$ to $a_n^2$ (note that $\|\hat{\boldsymbol{\beta}}_{0}-\boldsymbol{\beta}(\tau)\|_2=O_{\mathbb{P}}(a_{n})$). Therefore, an iterative refinement of the initial estimator will successively improve the estimation accuracy. The effect of bias reduction is mainly due to the term $\frac{Y_i}{h}H'\Big(\frac{Y_i-\boldsymbol{X}'_i\widehat{\boldsymbol{\beta}}_0}{h}\Big)$ in (ref), which is induced by the smoothing technique (please see more details in the proof of Proposition (ref)).

The previous discussions only involve one round of aggregation. } Now we are ready to present the theoretical results for our DC LEQR $\widehat\boldsymbol{\beta}^{(q)}$ in Algorithm (ref) with multiple rounds of aggregations. By a recursive argument based on (ref), we establish the following Bahadur representation.

theoremAssume the initial estimator $\hat{\boldsymbol{\beta}}_{0}$ in (ref) satisfies $\|\hat{\boldsymbol{\beta}}_{0}-\boldsymbol{\beta}(\tau)\|_2=O_{\mathbb{P}}(\sqrt{p/m})$. Let $h_{g}= \max(\sqrt{p/n},(p/m)^{2^{g-2}})$ for $1\leq g\leq q$. Assuming that (C1), (C2), and (C3) (or (C3$^{*}$)) hold with $h=h_q$ and $p$ also satisfies $p=O(m/(\log n)^2)$, then we have \begin{eqnarray} \hat{\boldsymbol{\beta}}^{(q)}-\boldsymbol{\beta}(\tau)=\frac{\boldsymbol{D}^{-1}}{n}\sum_{i=1}^{n}\boldsymbol{X}_{i}(I\{\epsilon_{i}\geq 0\}+\tau-1)+\boldsymbol{r}_{n} \end{eqnarray} with \begin{eqnarray} \|\boldsymbol{r}_{n}\|_2=O_{\mathbb{P}}\Bigl(\sqrt{\frac{ph_{q}\log n}{n}}+\Big{(}\frac{p}{m}\Big{)}^{2^{q-1}}\Bigr). \end{eqnarray}

The classical initial estimator based on QR in (ref) will satisfy $\|\widehat{\boldsymbol{\beta}}_0-\boldsymbol{\beta}(\tau)\|_2=O_{\mathbb{P}}(\sqrt{p/m})$ under some regularity conditions; see HE2000120. According to our choice of bandwidth, we have $h_1=\sqrt{p/m}$ and thus the convergence rate of the initial estimator is $O(h_1)$, which satisfies the condition in Proposition (ref). Furthermore, we note that any initial estimator $\hat{\boldsymbol{\beta}}_{0}$ with $\|\hat{\boldsymbol{\beta}}_{0}-\boldsymbol{\beta}(\tau)\|_2=O_{\mathbb{P}}(\sqrt{p/m})$ can be used in the first iteration and the same Bahadur representation in Theorem (ref) holds.

The condition $p=O(m/(\log n)^2)$ in Theorem (ref) ensures that $\sqrt{p/m}=o(1)$, which implies the consistency of the initial estimator (and the $1/(\log n)^2$ factor is used for balancing the terms in $\|\boldsymbol{r}_{n}\|_2$ in (ref)). This condition on $p$ cannot be implied by either (C3) or (C3$^{*}$); and we choose to present (C3) (or (C3$^{*}$)) at the beginning since they are the minimum requirements to obtain Propositions (ref)--(ref). On the other hand, the condition $p=o(nh_q/\log n)$ in (C3) will be satisfied if $p=o(n/(\log n)^2)$ since $h_q \geq \sqrt{p/n}$. Therefore, we can simply impose a stronger condition $p=o(m/(\log n)^{2})$ that unifies the condition $p=O(m/(\log n)^2)$ and that in (C3).

rem[Nearly optimal rate of the Bahadur remainder term] From Theorem (ref), as long as \begin{equation} q\geq 2 +\log\{\log (\sqrt{p/n})/\log (p/m)\}/\log 2, \end{equation} the bandwidth for the $q$-th iteration is $h_q= \sqrt{p/n}$. Then the first term in the right hand side of (ref) is $(p/n)^{3/4} \sqrt{\log n}$ and the second term is bounded by $p/n$, which is dominated by the first term. Therefore, the Bahadur remainder term $\boldsymbol{r}_n$ of our method achieves a nearly optimal rate \begin{equation} \|\boldsymbol{r}_n\|_2=O_{\mathbb{P}}\left((p/n)^{3/4}(\log n)^{1/2}\right). \end{equation} In fact, for classical QR estimator $\widehat{\boldsymbol{\beta}}_{QR}$ and fixed $p$, it is known that the rate $n^{-3/4}$ can not be improved except for a $\log n$ term koenker2005quantile. Note that in a common scenario when $n=O(m^A)$ and $p=O(m^{\delta})$ for some constants $A \geq 1$ and $0<\delta<1$, the right hand side of ((ref)) is bounded by a constant. Therefore, a constant number of rounds of aggregations is sufficient to obtain a nearly optimal rate in Bahadur representation.

Applying the central limit theorem to Theorem (ref), we obtain the asymptotic distribution of $\widehat{\boldsymbol{\beta}}^{(q)}-\boldsymbol{\beta}(\tau)$ as follows.

theoremSuppose that all the conditions in Theorem (ref) hold. Further, assume that $n=O(m^{A})$ for some constant $A \geq 1$, $p=o(\min\{n^{1/3}/(\log n)^{2/3},m^{\delta}\})$ for some $0<\delta<1$ and the number of iterations $q$ satisfies (ref). By choosing the bandwidth sequence $h_{g}= \max(\sqrt{p/n},(p/m)^{2^{g-2}})$ for $1 \leq g \leq q$, for any $\boldsymbol{v}\in \mathbb{R}^{p}$ with $\boldsymbol{v}\neq 0$, \begin{equation} \frac{n^{1/2}\boldsymbol{v}'\left(\hat{\boldsymbol{\beta}}^{(q)}-\boldsymbol{\beta}(\tau)\right)}{\sqrt{\boldsymbol{v}'D^{-1}\mathbb{E} [\boldsymbol{X}\boldsymbol{X}']\boldsymbol{D}^{-1}\boldsymbol{v}}}\Rightarrow N(0,\tau(1-\tau)) \end{equation} as $m, n\rightarrow \infty$.

Note that to establish central limit theorem, we need $p=o(n^{1/3}/(\log n)^{2/3})$ in Theorem (ref), which ensures that $\|\boldsymbol{r}_n\|_2 = o(1/\sqrt{n})$ (see (ref)). Therefore, as the first term in (ref) is an average of $n$ i.i.d. zero-mean random vectors, the reminder term $\boldsymbol{r}_n$ is dominated by each coordinate of the first term in (ref). For the classical QR estimator $\hat{\boldsymbol{\beta}}_{QR}$ in (ref) (assuming all data is pooled together), the corresponding condition should be $p=o(n^{1/3}/(\log n)^{2/3})$, see HE2000120. This is the same with our condition except for the term $m^{\delta}$ which is required for the consistency of the initial estimator in our method. We also note that the conditions $n=O(m^{A})$ and $p=o(m^{\delta})$ ensure that the number of required iterations $q$ from (ref) is a constant (see Remark (ref)). Therefore, we only need to perform a constant number of aggregations as $m,n \rightarrow \infty$.

Theorem (ref) shows that $\widehat{\boldsymbol{\beta}}^{(q)}$ achieves the same asymptotic efficiency as $\widehat{\boldsymbol{\beta}}_{QR}$ in ((ref)) computed directly on all the samples. When $p$ is fixed, as compared to the na\"ive-DC that also achieves (ref) but under the condition $n=o(m^2)$, our approach removes the restriction on the relationship of $m$ and $n$ by applying multiple rounds of aggregations. It is also important to note that the required number of rounds $q$ in (ref) is usually quite small even with a large dimension $p$.

Given (ref), we only need consistent estimators of $\boldsymbol{D}$ and $\mathbb{E} [\boldsymbol{X}\boldsymbol{X}']$ to construct confidence interval of $\boldsymbol{v}'\boldsymbol{\beta}(\tau)$ for any given $\boldsymbol{v}$. It is natural to use $\boldsymbol{D}_{n,h}$ and $\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{X}_{i}\boldsymbol{X}'_{i}$ to estimate $\boldsymbol{D}$ and $\mathbb{E} [\boldsymbol{X}\boldsymbol{X}']$, respectively. These estimators can be easily implemented under memory constraint by averaging the local sample estimators on each batch of data. The proofs of Propositions (ref), (ref) and Theorems (ref), (ref) are provided in Appendix (ref).

rem[Data-adaptive choices of bandwidth] In practice, one could use the bandwidth $h_{g}=c_g\max(\sqrt{p/n},(p/m)^{2^{g-2}})$ with a scaling constant $c_g$ to further improve the empirical performance. An intuitive data-adaptive way of choosing $c_g$ is provided as follows. Given a set of candidate choices for $c_g$ (e.g., $\{c_1,\ldots, c_L\}$), we choose the best $c_g$ by minimizing \begin{equation} S(c):=\Big{\|}\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{X}_{i}(I\{Y_i-\boldsymbol{X}_{i}'\hat{\boldsymbol{\beta}}^{(g)}_{c}\geq 0\}+\tau-1)\Big{\|}_2, \end{equation} where $\hat{\boldsymbol{\beta}}^{(g)}_{c}$ denotes the estimator in the $g$-th round of aggregation with the constant $c$ in the bandwidth. That is, $c_g= \mathop{\rm arg\min}_{c\in \{c_1,\ldots, c_L\}} S(c)$. In a distributed setting, the method only requires a small amount of communication. More specifically, given $\hat{\boldsymbol{\beta}}^{(g)}_{c}$, each machine $k$ returns $\sum_{i \in \mathcal{H}_k}\boldsymbol{X}_{i}(I\{Y_i-\boldsymbol{X}_{i}'\hat{\boldsymbol{\beta}}^{(g)}_{c}\geq 0\}+\tau-1)$ (i.e., an $O(p)$ vector) to the center for computing $S(c)$. We also evaluate the performance of our algorithm with the use of data-adaptive bandwidth in Section (ref). It is worthwhile noting that the bandwidth tuning is not a critical issue for our algorithm (in contrast to many other smoothed QR estimators) since our estimator is constructed via multiple rounds of aggregations. Even using an inaccurate constant in bandwidth (as long as the bandwidths for different rounds shrink at the right rate of $(p/m)^{2^{g-2}}$), our method can achieve good performance by simply performing more rounds of aggregations. In Section (ref), we will provide simulation studies to show that our algorithm is insensitive to the scaling constant in the bandwidth.

rem[Discussions with related literature] Note that our estimator $\hat{\boldsymbol{\beta}}$ can be written as { \begin{eqnarray*} \hat{\boldsymbol{\beta}}=\hat{\boldsymbol{\beta}}_{0}+\boldsymbol{D}^{-1}_{n,h} \left[\frac{1}{n}\sum\limits_{i=1}^n \boldsymbol{X}_i\Big\{H\Big(\dfrac{Y_i-\boldsymbol{X}'_i\widehat{\boldsymbol{\beta}}_0}{h}\Big)+\tau-1+\dfrac{Y_i-\boldsymbol{X}'_i\hat{\boldsymbol{\beta}}_{0}}{h}H'\Big(\dfrac{Y_i-\boldsymbol{X}'_i\widehat{\boldsymbol{\beta}}_0}{h}\Big)\Big\}\right]. \end{eqnarray*} \normalfont } This formula is closely related to the estimator for quantile regression considered in pang2012variance, where they introduced the estimator (non-censored version) \begin{eqnarray*} \hat{\boldsymbol{\beta}}_{PLW}=\hat{\boldsymbol{\beta}}_{0}+\boldsymbol{A}^{-1}_{n} \left[\frac{1}{n}\sum\limits_{i=1}^n \boldsymbol{X}_i\Big\{H\Big(\dfrac{Y_i-\boldsymbol{X}'_i\widehat{\boldsymbol{\beta}}_0}{\sqrt{\boldsymbol{X}'_{i}W\boldsymbol{X}_{i}}}\Big)+\tau-1\Big\}\right], \end{eqnarray*} where $ \boldsymbol{A}_{n}=\frac{1}{n}\sum\limits_{i=1}^n \frac{\boldsymbol{X}_i\boldsymbol{X}_i'}{\sqrt{\boldsymbol{X}'_{i}W\boldsymbol{X}_{i}}}H'\Big(\frac{Y_i-\boldsymbol{X}'_i\widehat{\boldsymbol{\beta}}_0}{\sqrt{\boldsymbol{X}'_{i}W\boldsymbol{X}_{i}}}\Big) $ and $W$ is a weight matrix with order $O(1/n)$. pang2012variance applied the smoothing to the original score function $\sum_{i=1}^{n}\boldsymbol{X}_{i}(I\{Y_{i}-\boldsymbol{X}'_{i}\widehat{\boldsymbol{\beta}}_{0}\geq 0\}+\tau-1)$, while our estimator comes from the smoothing to the loss function $\rho_{\tau}(\cdot)$ and hence we have an additional term $\frac{Y_i-\boldsymbol{X}'_i\hat{\boldsymbol{\beta}}_{0}}{h}H'\Big(\frac{Y_i-\boldsymbol{X}'_i\widehat{\boldsymbol{\beta}}_0}{h}\Big)$. Note that this term plays a key role in reducing the bias induced by the initial estimator $\hat{\boldsymbol{\beta}}_{0}$ to $\|\hat{\boldsymbol{\beta}}_{0}-\boldsymbol{\beta}(\tau)\|^{2}_{2}$. As pointed above, this allows our estimator to have a successive improvement on the estimation accuracy by iteratively updating the initial estimator. Moreover, the recent work by jordan2016communication and wang2016efficient also proposed iterative approaches in distributed setting for successive refinement of an estimator. However, there are a few key differences between our DC LEQR and their approaches. First, their results require the loss function to have Lipschitz continuous second-order derivative, which is not satisfied by the original quantile loss function. Even if we replace the indicator function $I\{x\geq 0\}$ in the quantile loss by the smoothed version $H(x/h)$, the second derivative of the loss function will not satisfy their conditions (e.g., Assumption PD in jordan2016communication requires that the “Lipschitz constant” of the second derivative has a uniform upper bound). Furthermore, our results allow $p\rightarrow \infty$ in the inference problem without $\ell_1$-regularization.

rem[General heterogenous case] We further consider a more general heterogenous case where $\{\boldsymbol{X}_i,\epsilon_i\}$'s are independent, but not identically distributed from the model (ref). Due to space limitations, this heterogenous case is relegated to Appendix (ref).

Asymptotics for Online LEQR

The next theorem gives the limiting behavior of the online LEQR in Algorithm (ref).

theoremSuppose that (C1)-(C3) hold and $p=o(m/(\log m)^{2})$. We have for any $A>0$ and uniformly in $1\leq j\leq m^{A}$, \begin{eqnarray} \hat{\boldsymbol{\beta}}[j]-\boldsymbol{\beta}(\tau)&=&\frac{\boldsymbol{D}^{-1}}{m+j}\sum_{i=-m+1}^{j}\boldsymbol{X}_{i}\Big{(}I\{\epsilon_{i}\geq 0\}+\tau-1\Big{)}+\boldsymbol{r}_{m,j}, \end{eqnarray} where \begin{eqnarray} \|\boldsymbol{r}_{m,j}\|_2=O_{\mathbb{P}}\Big(\sqrt{\frac{p}{m+j}}\Big{\\Big{(}\frac{p}{m}\Big{)}^{1/4}\sqrt{\log m}+\frac{\sqrt{p}}{m^{1/4}}\Big{\}}\Big). \end{eqnarray} Furthermore, when $p=o(m^{1/4})$, we have for any $\boldsymbol{v}\in \mathbb{R}^{p}$ with $\boldsymbol{v}\neq 0$, \begin{equation} \frac{(m+j)^{1/2}\boldsymbol{v}'\left(\hat{\boldsymbol{\beta}}[j]-\boldsymbol{\beta}(\tau)\right)}{\sqrt{\boldsymbol{v}'D^{-1}\mathbb{E} [\boldsymbol{X}\boldsymbol{X}']\boldsymbol{D}^{-1}\boldsymbol{v}}}\Rightarrow N(0,\tau(1-\tau)) \end{equation} as $m\rightarrow \infty$.

We note that $m+j$ is the total number of used samples (including the samples for initialization) up to time $j$. From ((ref)), we have $\|\widehat{\boldsymbol{\beta}}[j]-\boldsymbol{\beta}(\tau)\|_2=O_{\mathbb{P}}(\sqrt{p/(m+j)})$ for any $1\leq j\leq m^{A}$ when $\boldsymbol{r}_{m,j}$ is dominated by the first term. By Theorem 1 in siegmund1969, we can further obtain $\|\widehat{\boldsymbol{\beta}}[j]-\boldsymbol{\beta}(\tau)\|_2=O_{\mathbb{P}}\Big{(}\sqrt{\frac{p\log\log m}{m+j}}\Big{)}$ uniformly for $1\leq j\leq m^{A}$. To establish the asymptotic distribution in (ref), we need $p=o(m^{1/4})$, which ensures that $\|\boldsymbol{r}_{m,j}\|_2 = o(1/\sqrt{m+j})$ (see (ref)).

Simulations

In this section, we provide simulation studies to illustrate the performance of DC LEQR for constructing confidence intervals for QR problems. We generate data from a linear regression model

equation[equation omitted — 108 chars of source]

where $\boldsymbol{X}_i=(1, X_{i1}, \ldots, X_{ip})'\in \mathbb{R}^{p+1}$ is a random covariate vector. Here, $(X_{i1}, \ldots, X_{ip})'$ follows a multivariate uniform distribution $\mathrm{Unif}([0,1]^p)$ with $\text{Corr}(X_{ij},X_{ik})=0.5^{|j-k|}$ for $1\leq j\neq k\leq p$. Please refer to falk1999simple for the construction of such a multivariate uniform distribution. The regression coefficient vector $\boldsymbol{\beta}=\mathbf{1}_{p+1}$. The errors $\epsilon_i$'s are generated independently from the following three distributions:

enumerate[(1)] • homoscedastic normal, $\epsilon_i\sim N(0,1)$; • heteroscedastic normal, $\epsilon_i\sim N(0,(1+0.3X_{i1})^2)$; • exponential, $\epsilon_i\sim\mathrm{Exp}(1)$.

For each quantile level $\tau$, we compute the corresponding true QR coefficient vector $\boldsymbol{\beta}(\tau)$ in the QR model (ref) by shifting $\epsilon_i$ such that $\mathbb{P}(\epsilon_i\leq 0|\boldsymbol{X}_i)=\tau$:

enumerate[(1)] • homoscedastic normal, $\boldsymbol{\beta}(\tau)=\boldsymbol{\beta}+\Phi^{-1}(\tau)\boldsymbol{e}_1$; • heteroscedastic normal, $\boldsymbol{\beta}(\tau)=\boldsymbol{\beta}+\Phi^{-1}(\tau)\boldsymbol{e}_1+0.3 \Phi^{-1}(\tau)\boldsymbol{e}_2$; • exponential, $\boldsymbol{\beta}(\tau)=\boldsymbol{\beta}+F_{\mathrm{exp}}^{-1}(\tau)\boldsymbol{e}_1$;

Here $\Phi$ and $F_{\mathrm{exp}}$ are the cumulative distribution function of standard normal distribution and exponential distribution with parameter $1$. The vector $\boldsymbol{e}_i$ (for $i=1,\ldots, p+1$) is the $(p+1)$-dimensional canonical vector with the $i$-th element being one and all the other elements being zero.

We use the integral of a biweight (or quartic) kernel as the smoothing function $H$:

equation[equation omitted — 208 chars of source]

It is easy to see that it satisfies Condition (C2).

Coverage rates

table[table omitted — 2,756 chars of source]

table[table omitted — 2,758 chars of source]

table[table omitted — 2,727 chars of source]

We compute the DC LEQR $\widehat{\boldsymbol{\beta}}^{(q)}$ in (ref) and measure the performance of $\widehat{\boldsymbol{\beta}}^{(q)}$ in terms of the statistical inference. In particular, we report average coverage rates of the confidence interval of $\boldsymbol{v}_0'\boldsymbol{\beta}(\tau)$, where $\boldsymbol{v}_0=(p+1)^{-1/2}\mathbf{1}_{p+1}$. We set the nominal coverage probability $1- \alpha_0$ to $95\%$. From Theorem (ref), an oracle $(1- \alpha_0)$-th confidence interval for $\boldsymbol{v}_0'\boldsymbol{\beta}(\tau)$ is given by

equation[equation omitted — 253 chars of source]

where $z_{\alpha_0/2}$ is the $(1-\alpha_0/2)$-quantile of the standard normal distribution. To construct the confidence interval, we estimate $\boldsymbol{D}$ and $\mathbb{E}\left[\boldsymbol{X}\boldsymbol{X}'\right]$ by $\boldsymbol{D}_{n,h}$ and $\frac{1}{n} \sum_{i=1}^n \boldsymbol{X}_i \boldsymbol{X}_i'$, respectively. There are two major advantages of this approach. First, since $\boldsymbol{D}_{n,h}$ has already been obtained in computing DC LEQR, we estimate $\boldsymbol{D}$ without any extra computation. Second, both $\boldsymbol{D}_{n,h}$ and $\frac{1}{n} \sum_{i=1}^n \boldsymbol{X}_i \boldsymbol{X}_i'$ are in the form of summation over $n$ terms, which can be easily computed in a distributed setting with little communication cost. As we will show in Table (ref), the proposed estimator is very close to the truth. The scaling constant in bandwidth is simply set to one (as in our theorems) and more detailed experiments on the sensitivity analysis of the scaling constant is provided in Section (ref). We report empirical coverage rates as an average of 1000 independent runs of the simulations.

In Tables (ref)--(ref), we present the empirical coverage rates of our DC LEQR estimator, the na\"ive-DC estimator, and the oracle QR estimator in (ref) computed on all data points (denoted as QR All) for three different noise models. More precisely, we generate the error from one of three distributions (i.e., homoscedastic normal for Table (ref), heteroscedastic normal for Table (ref), and exponential for Table (ref)) and consider three different quantile levels $\tau=0.1,0.5,0.9$. In our experiment, we set $m=100$, $p=15$ and vary $n$ form $m^{1.6}$ to $m^3$ (i.e., $\log_m(n)$ from 1.6 to 3). From (ref), it is easy to see that we need number of aggregations $q \geq 4$. Thus, we report the performance of DC LEQR $\widehat{\boldsymbol{\beta}}^{(q)}$ for $q=4$ and $q=5$. We also report the case of $p=3$ in Appendix (ref).

As one can see from Tables (ref)--(ref), for most of the settings, the coverage rates of our DC LEQR are close to the nominal level of $95\%$ after $4$ rounds of aggregations ($q=4$). The coverage performance becomes quite stable for $q=5$ iterations. On the other hand, for the na\"ive-DC estimator, the coverage rates are quite low in most settings, especially when $n$ is larger than $m^2$.

Note that for na\"ive-DC and QR All, we use the same estimator of the limiting variance in (ref) as in our DC LEQR when $q=4$. More precisely, we use $\boldsymbol{D}_{n,h}$ computed in the $4$-th iteration to estimate $\boldsymbol{D}$ in (ref) when constructing the confidence intervals of na\"ive-DC and QR All estimators. We will show in Table (ref) below that the proposed estimator of the limiting variance performs well. In fact, we also use the true limiting variance to construct the confidence interval and the coverage rates for all the methods are almost the same.

In addition, we also report the simulation study with large dimension $(p= 1000)$. The results and analysis are relegated to Appendix (ref). From the results, we can infer that the coverage rates get better as the iterative refinement proceeds. In particular, the coverage rates are close to the nominal level 95% after $4$ iterations when the dimension $p=1000$. In summary, when the dimension $p$ is large, the proposed DC LEQR algorithm still achieves desirable performance with a small number of iterations.

Bias and variance analysis

To see the improvement of DC LEQR over na\"ive-DC when $n$ is excessively larger than the subset size $m$, we also report the mean bias and variance of our proposed DC LEQR, na\"ive-DC, and QR All in Tables (ref)--(ref). The mean bias and variance of $\boldsymbol{v}_0'\widehat{\boldsymbol{\beta}}$ are based on 1000 independent runs of simulations.

From Tables (ref)--(ref), the bias of our method is quite small while the na\"ive-DC approach has a much larger bias regardless of the sample size $n$. For the variance, it decays with the rate $1/n$ as $n$ goes large for all methods. For most cases of using na\"ive-DC, as $\log_m(n)$ exceeding $2$, the squared bias becomes comparable or larger than the variance, which explains the reason of the failure of na\"ive-DC when $n$ is large as compared to $m$. On the other hand, the bias of our proposed DC LEQR is similar to that of QR All and much smaller than that of na\"ive-DC.

Sensitivity analysis and data-adaptive choice of the bandwidth

In this section, we show the empirical performance of the data-adaptive choice of bandwidth in Remark (ref) and the sensitivity of the scaling constant in bandwidth. Due to space limitations, we report $\tau=0.1$, homoscedastic normal noise case as an example. More noise cases (e.g., heteroscedastic normal and exponential cases) are relegated to Appendix (ref), and observations are similar to the homoscedastic normal case.

table[table omitted — 1,794 chars of source]

Table (ref) shows coverage rates of the DC LEQR with $q=1,2,3,4,5$ iterations. Similar to the setting in Tables (ref)--(ref), we choose $m=100$, $p=15$, and $n$ varies from $m^{1.6}$ to $m^3$. We report the performance of DC LEQR using different fixed constants $c=1,3,5,10$ in bandwidth $h_g$ for $1\leq g \leq q$ (using the same constant in all iterations) as well as our data-adaptive choice of bandwidth. For the data-adaptive bandwidth, we choose the best scaling constant from a list of 1000 equally spaced constants from a very small number (0.1) to a large one (100) according to Remark (ref). Note that different scaling constants will be chosen for different iterations. As one can see, for $q=1$ and $q=2$, the adaptive method indeed achieves better coverage than other scaling constants. On the other hand, when $q \geq 3$, all different choices of scaling constants lead to coverage rates close to the nominal level of 95%. This experiment suggests that for our proposed iterative aggregation approach, even when a sub-optimal scaling constant is used, one can still achieve good performance by performing more iterations.

table[table omitted — 1,739 chars of source]

Moreover, to investigate the sensitivity of the scaling constant in terms of variance estimation, we also present the square root of the ratio of the estimated variance of our approach versus the true limiting variance, i.e.,

equation[equation omitted — 333 chars of source]

where $\boldsymbol{D}_{n,h}$ is computed for iterations $q=1,2,3,4,5$ and for each fixed $q$, the bandwidth $h_g=c \max(\sqrt{p/n},(p/m)^{2^{g-2}})$ for $1 \leq g \leq q$. In Table (ref), we report the performance of the variance estimation with different choices of the scaling constant $c$ of the bandwidth $h$ in constructing $\boldsymbol{D}_{n,h}$ .

From Table (ref), the ratio is very close to $1$ when $q=4$ or $5$. Therefore, when $q$ is large, the proposed variance estimator is a reliable one in the distributed setting. Moreover, we notice that the ratio is very stable for different choices of the scaling constant $c$, which illustrates the robustness of the estimator.

table[table omitted — 2,105 chars of source]

Computation efficiency

We further conduct experiments to illustrate the computation efficiency of our algorithm for different $m$ and $n$ with $\tau=0.1$, $p=15$ and $\epsilon_i \sim N(0,1)$. We compare the computation time of DC LEQR versus that of na\"ive-DC as well as QR All in Table (ref). We report the bias and variance and the coverage rates for reference of the performance of the estimators.

First of all, the computation time of DC LEQR is about as twice faster than that of the QR All, especially when $n$ is large. It is also faster than the na\"ive-DC and with a much better coverage when $n$ is much larger than $m$. Moreover, the time of our algorithm grows almost linearly in both $n$ and $q$, which is consistent with the computation time analysis in Section (ref). In contrast, we observe that the time of QR grows faster than a linear function in the sample size $n$. We also observe that, for each fixed $n$, the value $m$ has little effect on the computation time of DC LEQR. For na\"ive-DC, the squared bias in these cases dominate the variance, so the coverage rates are far below the nominal level. In the meantime, DC LEQR has around $95\%$ coverage in all four cases after $2$ iterations, and shows a similar behavior of bias and variance as in Table (ref).

Recently, researchers have developed new optimization techniques based on alternating direction method of multiplier (ADMM) for solving QR problems (see, e.g., yu2017parallel,gu2017admm). We further conduct comparisons to the ADMM approach and the details are provided in Appendix (ref).

Finally, we also conduct simulation studies for online LEQR and the results are presented in Appendix (ref).

Conclusions and future works

In this paper, we propose a novel inference approach for quantile regression under the memory constraint. The proposed method achieves the same asymptotic efficiency as the quantile regression estimator using all the data. Furthermore, it allows a weak condition on the sample size $n$ as a function of memory size $m$ and is computationally attractive. One key insight from this work is that na\"ively splitting data and averaging local estimators could be sub-optimal. Instead, the iterative refinement idea can lead to much improved performance for some inference problems in distributed environments.

In some applications, one would expect a weaker assumption on the distribution of data where the data could be correlated. It would be an interesting future direction to study the problem of inference for correlated data in distributed settings. Moreover, for the online problem, it is also interesting to consider the case where the model (e.g., $\boldsymbol{\beta}(\tau)$) is evolving over time. In this case, some exponential decaying techniques to down weight historical data might be useful.

In the future, we would also like to further explore this idea to other QR problems under memory constraints or in a distributed setup, e.g., $\ell_1$-penalized high-dimensional quantile regression (see, e.g., belloni2011l1,wang2012quantile,fan2016multitask) and censored quantile regression (see, e.g., wang2009locally,kong2013global,volgushev2014censored,Leng:14,He:17:censored).

Acknowledgements

The authors are very grateful to three anonymous referees and the associate editor for their detailed and constructive comments that considerably improved the quality of this paper.