A Selective Overview of Quantile Regression for Large-scale Data

1. Introduction

With the rapid advancement of technology, we are entering the era of large-scale data, characterized by the ability to generate and store vast amounts of information at reduced costs. Large sample sizes and high-dimensional features have become defining characteristics of many modern datasets, and the exponential growth of massive datasets across diverse fields has become increasingly common. For example, in genomics, the development of new technologies enables biologists to generate tens of thousands of datasets [1], and the abundance of genomic sequencing data facilitates the identification of genetic markers for rare diseases [2,3]. In economics and finance, an increasing number of companies are adopting data-driven approaches. By analyzing large datasets, such as stock prices and transaction records, companies can better assess risks and make informed decisions [1]. Social networks like Twitter, Facebook and YouTube also generate large amounts of data, enabling researchers to build predictive models that reveal individual characteristics[1,4,5,6]. Although large-scale data can provide richer information for data analysis in various domains, they also pose significant challenges to traditional statistical models [1]. For massive data with both high dimension and large sample size, here we focus on the main techniques that address the following three major challenges.

The first challenge lies in the larger volume of large-scale data. Datasets can reach exabyte-level sizes, yet the memory capacity of a single computer is limited, making it impossible to store the entire dataset. Moreover, traditional statistical models are predominantly built on the assumption that data can be processed on a single machine, such models could be ineffective for large-scale data [7]. To address this challenge, a fundamental strategy for storing and processing such data is distributed learning, often referred to as the “divide-and-conquer” (dc) approach. The idea of dc-method is to partition a massive dataset into smaller data blocks to mitigate memory constraints. A common implementation of this strategy involves compressing subsets of data on each machine into local estimator and then aggregating them on a central host to obtain the final model estimates [7]. In addition, many researchers have explored dc-strategies for various statistical models, such as linear regression [8,9], logistic regression [10], quantile regression [11,12,13], and references therein.

The second challenge arises from the high-speed generation of large-scale datasets. Learning models from such data in a streaming fashion is becoming increasingly ubiquitous, in which data arrive in streams and chunks (also referred to as streaming datasets) [14]. It is often numerically challenging or sometimes infeasible in such a data environment to update the model estimation with the entire dataset in memory. As a result, traditional offline methods, which require access to the full dataset, become less practical or even obsolete. Instead, online methods are capable of providing a way to analyze the out-of-memory data by constructing and updating the model sequentially at one time without storage requirements from previous raw data. The two key issues of streaming data analysis are (1) the storage challenge that old data are already discarded when new data arrive, and (2) the speed challenge that the analysis must be updated in a very fast way to support real-time information processing. Various online methods have been proposed to date, including aggregated estimate equation [15], cumulatively updating estimate equation [16], stochastic gradient descent (SGD) and its variants [17,18], as well as renewable estimators [19]. In this paper, we assume that at a given time point t, only the current data block $({\mathit{x}}_t,y_t)$ and the summary statistics based on the previous $(t-1)$ data blocks are available (see Figure 1). And we focus on the “renewable estimation” or “online updating", which directly update parameter estimates utilizing summary statistics from both previously processed data and newly arriving batches of data [13].

The third challenge is that large-scale data often show heterogeneity due to heteroskedastic variance or inhomogeneous covariate effects [1]. Traditional statistical models are generally predicated on the assumption that the data are homoscedastic and follow Gaussian distributions. However, these assumptions impose significant limitations when applied to the complex structural characteristics of large-scale data. Moreover, most of the existing methods for large-scale data analysis are focused on mean regression models [19,20,21,22,23]. When random error distributions are highly skewed or exhibit heteroskedasticity, linear mean regression models often yield inefficient estimates and poor predictive performance. In contrast, quantile regression (QR) models [24,25] offer a more robust alternative capable of effectively addressing these challenges. By complementing the exclusive focus of classical least squares regression on the conditional mean, QR offers a systematic strategy to examine how covariates influence the location, scale, and shape of the entire response distribution. It has been actively studied in machine learning, economics, finance, and statistics [26,27,28,29,30,31,32,33]. However, these traditional statistical methods for estimating QR models require processing the entire dataset simultaneously. Due to the infinite size of large-scale datasets and the limited memory capacity of computers, these conventional parameter estimation methods are not applicable [7]. While there has been substantial research on distributed learning algorithms and online updating algorithms for large-scale data [16,19,34,35], the non-differentiability of the objective function in quantile regression makes it challenging to directly adapt these algorithms for estimating quantile regression models.

Therefore, this paper provides a comprehensive methodological framework for QR in the context of large-scale data as shown in Table 1. Specifically, the main contributions of this overview are summarized as follows:

• Overcoming non-differentiability of the QR check function. We introduce common techniques on smooth or surrogate of the QR loss function.
• Parameter Estimation in distributed and streaming settings. We systematically reviewe methods to efficiently compute and aggregate QR estimates when large-scale data is stored across multiple machines or in one machine. Then we summarize recent research progress of renewable estimation for QR regression as the data arrive sequentially.
• Statistical Inference for Large-Scale QR. We briefly introduce statistic theoretic of QR estimator for large-scale data, such as asymptotic distribution and techniques about how to derive empirical variance of estimators.
• Applications. We introduce the applications in meteorological forecasting, demand and price forecasting, geographical distribution data and high-frequency data forecasting.
• Potential research directions. We provide some potential research directions which can promote the development of this field.

Due to space limitations, we primarily focus on the linear QR framework when reviewing existing work on large-scale data. We do not evaluate each method in detail, however, we will offer guidance on the scenarios where specific methods are most effective and discuss the potential for further advancing the existing literature to tackle emerging challenges in large-scale data quantile regression analysis.

The rest of this article is organized as follows. Section 2 briefly reviews the foundations of QR. Section 3 and Section 4 illustrate advances on QR for massive data and streaming data, respectively. Section 5 introduces the applications of QR for large-scale data. Limitations and potential future work are discussed in Section 6.

2. Foundations of Quantile Regression

In this section, we first show the definition and mathematical framework of QR in Section 2.1. Then, we introduce some common estimating techniques in Section 2.2. Finally, we summarize the key differences and advantages compared with mean regression in Section 2.3.

2.1. Definition and Mathematical Framework of QR

A basic tool in different scientific fields for analyzing the impact of a set of regressors X on the distribution of a response Y is QR, and it can be defined with the help of appropriate loss functions. For a fixed $\tau \in (0,1)$, the $\tau$th conditional quantile of Y given X, $q_{\tau }\left(Y\right|X=\mathit{x})$, can be obtained by minimizing asymmetrically weighted mean absolute deviations [24]:

$$\begin{array}{c}\hfill q_{\tau }\left(Y\right|X=\mathit{x})\in \arg \underset{\theta }{min}E\left([{\rho }_{\tau }(Y-\theta )-{\rho }_{\tau }\left(Y\right)]|X=\mathit{x}\right),\end{array}$$

(1)

where ${\rho }_{\tau }\left(u\right)=u\left(\tau -\mathbf{1}(u<0)\right)$ is the check function as shown in Figure 2, and $\mathbf{1}(·)$ is the indicator function. Observably, ${\rho }_{\tau }\left(u\right)$ is the piecewise linear function and ${\rho }_{0.5}\left(u\right)=0.5\left|u\right|$ is the absolute loss function. Note that subtracting ${\rho }_{\tau }\left(Y\right)$ from the expectation as shown in (1) makes the integrand well defined and finite without assuming $E\left(\right|Y\left|\right|X=\mathit{x})<\infty$.

Let $(y_i,{\mathit{x}}_i)$ be i.i.d. samples draw from $(Y,X)$, then the linear QR model follows:

$$\begin{array}{cc}\hfill q_{\tau }\left(y_i\right|{\mathit{x}}_i)& ={\mathit{x}}_i^{\top }\mathit{\beta }\left(\tau \right),i=1,\cdots ,n,\hfill \end{array}$$

(2)

where ${\mathit{x}}_i={(1,x_{i1},\cdots ,x_{ip})}^{\top }$ and $\mathit{\beta }\left(\tau \right)\in {\mathbb{R}}^{p+1}$ is the unknown quantile regression coefficients. The QR estimator $\widehat{\mathit{\beta }}\left(\tau \right)$ can be obtained via minimizing the following objective function:

$$\begin{array}{cc}\hfill \widehat{\mathit{\beta }}\left(\tau \right)& =\underset{\mathit{\beta }\left(\tau \right)}{argmin}\sum _{i=1}^n{\rho }_{\tau }(y_i-{\mathit{x}}_i^{\top }\mathit{\beta }\left(\tau \right)).\hfill \end{array}$$

(3)

Besides estimation, statistical inference for QR is also a critical aspect. Basically, it involves constructing confidence intervals, hypothesis testing and assessing the variability of the estimated parameters, etc. The central work for statistic inference is to derive the asymptotic distribution of the QR estimator as $n\to \infty$ . Assume $F_{ϵ\left(\tau \right)}(·|X)$ and $f_{ϵ\left(\tau \right)}(·|X)$ are the conditional cumulative distribution function and the density function of ${ϵ}_i\left(\tau \right)=y_i-{\mathit{x}}_i^{\top }{\mathit{\beta }}_{\tau }$, which satisfies $F_{ϵ\left(\tau \right)}^{-1}\left(\tau \right|X)=0$. As was shown by Koenker and Bassett [24], under some regularity conditions, $\widehat{\mathit{\beta }}\left(\tau \right)$ is asymptotically normal:

$$\begin{array}{c}\hfill \sqrt{n}(\widehat{\mathit{\beta }}\left(\tau \right)-\mathit{\beta }\left(\tau \right))\stackrel{d}{\to }N\left(0,\tau (1-\tau )D^{-1}{\Omega }_X{D}^{-1}\right),\mathrm{as}n\to \infty \end{array}$$

(4)

where $D=\mathbb{E}\left[f_{ϵ\left(\tau \right)}\left(0\right|X)X{X}^{\top }\right],{\Omega }_X=\mathbb{E}\left(X{X}^{\top }\right)$, and $\stackrel{d}{\to }$ denotes convergence in distribution. Although the asymptotic covariance matrix seems to suggest that estimation of regression quantiles becomes increasingly more precise when very small ($\tau \to 0$) or large ($\tau \to 1$) quantiles are considered (due to the product $\tau (1-\tau )$ in the numerator), in fact this effect is typically dominated by the error $f_{ϵ\left(\tau \right)}\left(0\right|X)$ which will usually be close to zero for such values of $\tau$. Hence, the estimation of interior quantiles close to the median will actually be more precise for most of the standard distributions, while estimation of outer quantiles is associated with larger uncertainty.

Here the sample estimation of the covariance matrix denotes as:

$$\begin{array}{cc}\hfill \widehat{D}=n^{-1}\sum _{i=1}^n{\mathit{x}}_i{\mathit{x}}_i^{\top }{\widehat{f}}_{ϵ\left(\tau \right)}\left(0\right|{\mathit{x}}_i),& {\widehat{\Omega }}_X=n^{-1}\sum _{i=1}^n{\mathit{x}}_i{\mathit{x}}_i^{\top },\hfill \end{array}$$

where ${\widehat{f}}_{ϵ\left(\tau \right)}\left(0\right|{\mathit{x}}_i)$ is the conditional density estimator of $f_{ϵ\left(\tau \right)}(·|{\mathit{x}}_i)$ evaluated at 0. Note that the asymptotic variance of QR estimator involves density function of the error terms, which is unknown generally. Thus, direct implementation via the asymptotic variance is difficult. Generally, an alternative inference would be based on bootstrap method. Please refer to Chapter 3 of Davino et al. [36] for more details. For simplicity, we will omit $\tau$ of $\mathit{\beta }\left(\tau \right)$ or $\widehat{\mathit{\beta }}\left(\tau \right)$ in the following paper whenever there is no ambiguity.

2.2. Common Estimation Techniques

Estimation of QR parameters typically involves solving an optimization problem to minimize the check function in (3). This function is nondifferentiable, which poses challenges for direct computation. Figure 2 illustrates the check function ${\rho }_{\tau }\left(u\right)$ at a specified level $\tau$, highlighting the nondifferentiability of the loss function at $u=0$. The QR estimation can mainly be classified into three types: (1) linear programming (LP) method; (2) smooth QR method, and (3) Bayesian QR method. We summarize them together with the main algorithms in Figure 3, which are briefly discussed below.

Linear programming technique expands the estimation problem by introducing $2n$ auxiliary variables $\mathit{u}={(u_1,\cdots ,u_n)}^{\top }$ and $\mathit{v}={(v_1,\cdots ,v_n)}^{\top }$ with $u_i={(y_i-{\mathit{x}}_i^{\top }\mathit{\beta })}_{+},v_i={({\mathit{x}}_i^{\top }\mathit{\beta }-y_i)}_{+}$, where ${\left(x\right)}_{+}=max(x,0)$ and therefore $\mathit{Y}-\mathit{x}\mathit{\beta }=\mathit{u}-\mathit{v}$ with $\mathit{Y}={(y_1,\cdots ,y_n)}^{\top },\mathit{x}={({\mathit{x}}_1,\cdots ,{\mathit{x}}_n)}^{\top }$. The minimization problem (3) can be rewritten as

$$\underset{\mathit{\beta },\mathit{u},\mathit{v}\in {\mathbb{R}}^{p+1}\times {\mathbb{R}}^{2n}}{min}\left\{\tau {\mathbf{1}}^{\top }\mathit{u}+(1-\tau ){\mathbf{1}}^{\top }\mathit{v}|\mathit{Y}=\mathit{x}\mathit{\beta }+\mathit{u}-\mathit{v}\right\},$$

(5)

where $\mathbf{1}$ is a $n\times 1$ unit vector. This is a constrained minimization problem with polyhedric constraints, i.e., the constraint defines a geometric object with flat faces and straight edges, such that the constrained problem corresponds to the original quantile restriction, but written in terms of the auxiliary variables. After augmenting the auxiliary variables and given the constraints, the minimization problem is now linear in the parameters and can therefore be tackled with linear programming techniques that allow to incorporate the polyhedric constraints. Refer to Lange [37] for a brief introduction to linear programming with a focus on statistics-related applications. By reformulation, (5) can be solved by simplex methods [38] or interior point methods [39].

Algorithms developed for the general LP problems in (5) may not fully deploy the properties of $L_1$ or quantile regression and have their own shortcomings. For example, the interior point algorithm can only give the approximate solutions of the original problem, and rounding is necessary to provide the same accuracy as that of the simplex algorithm. In this case, some heuristic approaches demonstrate the advantages of both speed and accuracy. One such approach is the smoothing method, which aims to approximate the objective function in (3), enabling the use of Newton or quasi-Newton optimization methods. There are three main types of smoothing techniques: 1) Horowitz’s smoothed loss [40,41]; 2) Huber’s smoothed loss [42] and 3) Convolution-type smoothing [43,44,45,46].

More specifically, Horowitz’s smoothed loss approximateds the indicator function $\mathbf{1}(u<0)$ with a smooth function $H(-u/h)$, where $H(·)$ takes values between 0 and 1, and $h\to 0$ is the bandwidth. For example, $H(·)$ can be chosen as the distribution function of standard normal $N(0,1)$. With this approximation, one can replace ${\rho }_{\tau }\left(u\right)$ by

$$l_{h,\tau }^{\mathrm{Horo}}\left(u\right)=u(\tau -H(-u/h)).$$

(6)

Instead, Huber’s smoothed method approximated ${\rho }_{\tau }\left(u\right)$ by the generalized Huber function $l_{\gamma ,\tau }^{\mathrm{Huber}}\left(u\right)$, i.e.,

$$l_{\gamma ,\tau }^{\mathrm{Huber}}\left(u\right)=\left\{\begin{array}{cc}u(\tau -1)-{\displaystyle \frac{1}{2}}{(\tau -1)}^2\gamma ,\hfill & u\le (\tau -1)\gamma \hfill \\ {\displaystyle \frac{u^2}{2\gamma }},\hfill & (\tau -1)\gamma \le u\le \tau \gamma \hfill \\ \tau u-{\displaystyle \frac{1}{2}}{\tau }^2\gamma ,\hfill & t\ge \tau \gamma \hfill \end{array}\right.,$$

(7)

where $\gamma >0$ is a threshold. Different from the former, the convolution-type smoothing method approximates ${\rho }_{\tau }\left(u\right)$ with a smooth function $l_{h,\tau }^{\mathrm{conquer}}\left(u\right)$, i.e.,

$$l_{h,\tau }^{\mathrm{conquer}}\left(u\right)=({\rho }_{\tau }\ast K_h)\left(u\right)={\int }_{-\infty }^{\infty }{\rho }_{\tau }\left(\nu \right)K_h(\nu -u)d\nu ,$$

(8)

where ∗ denotes the convolution operator, and throughout the paper, $K_h\left(u\right)=h^{-1}K(u/h)$ with $K(·)$ being a kernel function that integrates to one and $h>0$ being a bandwidth. Following He et al. [44], we adopt the term conquer for convolution-type smoothed quantile regression.

Thus, the corresponding smoothed QR estimator can be obtained via the empirical smoothed loss function in (6)-(8), respectively. In Section 3 and Section 4, we show how these smoothing techniques used in the large-scale setting.

The third option for QR estimation is Bayesian approaches, which utilize the Asymmetric Laplace distribution (ALD) as a likelihood function [47,48]. Let $\mathrm{ALD}(\mu ,\sigma ,\tau )$ be the ALD with location parameter $\mu$, scale parameter $\sigma >0$ and asymmetric (skewness) parameter $\tau$, and the pdf is given by

$$\begin{array}{cc}\hfill f(u;\mu ,\sigma ,\tau )& ={\displaystyle \frac{\tau (1-\tau )}{\sigma }}exp\left\{-{\rho }_{\tau }\left({\displaystyle \frac{u-\mu }{\sigma }}\right)\right\}={\displaystyle \frac{\tau (1-\tau )}{\sigma }}exp\left\{-{\displaystyle \frac{1}{\sigma }}{\rho }_{\tau }\left(u-\mu \right)\right\}.\hfill \end{array}$$

This method imposed the $\mathrm{ALD}(0,\sigma ,\tau )$, on the residual terms, hence the cumulative distribution function of $ϵ$ is $f\left(ϵ\right)=\tau (1-\tau )exp\left\{{\rho }_{\tau }\left(ϵ\right)\right\}$, then $y_i|{\mathit{x}}_i\sim \mathrm{ALD}({\mathit{x}}_i^{\top }\mathit{\beta },\sigma ,\tau )$. Thus, the likelihood is reformatted as:

$$\begin{array}{c}\hfill \mathcal{L}\left(\mathit{Y}\right|\mathit{\beta })={\displaystyle \frac{{\tau }^n{(1-\tau )}^n}{{\sigma }^n}}exp\left\{-{\displaystyle \frac{1}{\sigma }}\sum _{i=1}^n{\rho }_{\tau }(y_i-{\mathit{x}}_{\mathit{i}}^{\top }\mathit{\beta })\right\}.\end{array}$$

(9)

For estimation, given a prior distribution of $\mathit{\beta }$, the Gibbs sampler can be used to derive the Bayesian QR model. In addition, the EM algorithm can also be applied to estimate $\mathit{\beta }$ in (9) due to the fact that a mixture representation was developed for the ALD based on a normal variance–mean mixture with an exponential mixing distribution. More details can be found in Kozumi and Kobayashi [48].

Comparing the aforementioned estimation algorithms, the LP algorithm faces significant computational challenges when applied to large datasets. In contrast, smooth QR and Bayesian QR methods benefit from the continuity of their loss functions, and it is easy to derive iterative algorithm or closed-form estimators. These properties make them more adaptable to modern computational strategies, such as distributed and online algorithms, which are often employed to manage the growing complexity of data in practical applications. Finally, we list some standard QR packages, which can be found in R (with the “quantreg” package), SAS (through the “PROC QUANTREG” procedure), Python (via “statsmodels” library and QuantReg class), Matlab (with the “quantreg” function) and STATA (using the “qreg” command).

2.3. Comparison with Mean Regression: Key Differences and Applications

Quantile regression goes beyond modeling the conditional mean and provides a detailed representation of the entire conditional distribution. This makes QR particularly valuable in a wide range of applications, with its advantages evident in the following key aspects: (1) robustness, (2) interpretability, and (3) flexibility.

More specifically, unlike traditional mean regression, which minimizes the mean squared error and assigns equal weights to all samples, QR is less sensitive to outliers. This robustness arises because QR minimizes the check loss function, which adjusts the influence of data points based on their position relative to the quantile being estimated. Moreover, by analyzing the effects of predictors across different quantiles, QR effectively captures heteroskedasticity, outliers, and skewed distributions in the data. This ability to offer insights into multiple aspects of the data distribution enables QR to deliver a more nuanced and reliable analysis compared to traditional mean-based regression methods. For example, in the medical field, anesthesia research is a significant area of study. However, related variables, such as the length of hospital stay, ventilator dependence, blood loss, and pain scores, often violate the assumptions of normally distributed residuals, homoskedasticity, and independence. Relying solely on linear regression in these cases can result in inaccurate model estimates. To address this issue, Staffa et al. [49] applied multivariable quantile regression, which provides a more comprehensive and accurate understanding of the relationships between risk factors and outcomes.

Besides robustness, interpretability means that QR offers a comprehensive view of the effects of explanatory variables at different quantile levels. Many studies show that QR provides more insights than mean regression. For example, in non-Gaussian data distributions, QR improves model accuracy and is widely used in market return studies [50,51,52]. It also reveals weight change trends among different groups in U.S. girls’ weight and age studies [53]. In medical research, QR uncovers differences in the relationship between ventilator dependence and hospital stay across quantiles [49]. In ecology, QR provides a more complete view of causal relationships, as shown in studies linking trout abundance to stream width-to-depth ratios [54,55]. Empirical results indicate that mean regression might incorrectly suggest no relationship between trout density and stream ratios.

Lastly, in cases where data is missing or observed values are biased from actual values, QR effectively mitigates the impact of missing data on model fitting [56]. For instance, in ecological research, scholars have noted that when ecological limiting factors constrain organisms [55], many other unmeasured ecological factors may also impose limitations on biological responses. This makes mean regression insufficient for capturing the true patterns of biological responses. In contrast, QR can estimate the upper limits of biological response potential without being affected by unmeasured factors. Cade [57] explored in greater detail the multiplicative interactions between measured and unmeasured ecological factors and addressed the aforementioned analytical issues using QR methods.

3. Advances in QR for Massive Data

In this section, we briefly introduce foundational work on QR for massive datasets. In such situations, the size of the dataset can become so large that the storage of it on a single machine is impractical. Additionally, traditional methods often become infeasible due to the significant computational burden. In the following section, we focus on two common approaches for handling massive data: (1) distributed computing and (2) subsampling methods in Section 3.1 and Section 3.2, respectively. We also discuss the scenarios in which each of these methods is most effective. A summary is provided in Table 2.

3.1. Distributed Computing

Distributed computing (dc) is a computing paradigm in which a single system’s tasks are divided across multiple independent computing nodes that work collaboratively to achieve a common goal. These nodes, often geographically dispersed, communicate and coordinate their works through a network. A related concept is the dc framework [22,58,59,60], which breaks down a problem into smaller sub-problems, solves each problem independently, and combines their results to derive the overall solution. And this concept has already been applied in practice by many scholars [23,61,62,63].

For large-scale datasets, we below describe a general dc algorithm for QR estimation. Assume the entire samples ${\{y_i,{\mathit{x}}_i\}}_{i=1}^N$ follows (3), we split the data indices $\{1,2\cdots ,N\}$ into K subsets $\{{\mathcal{D}}_1,\cdots ,{\mathcal{D}}_K\}$ with each size $n_k$, that is ${\sum }_{k=1}^K{n}_k=N$. Correspondingly, the entire data set ${\{y_i,{\mathit{x}}_i\}}_{i=1}^N$ is divided into K batches $\{{\mathcal{M}}_1,{\mathcal{M}}_2,\cdots ,{\mathcal{M}}_K\}$, where ${\mathcal{M}}_k=\{y_i,{\mathit{x}}_i:i\in {\mathcal{D}}_k\}$ for $1\le k\le K$.

During estimation, each data batch is stored on a single machine and a low-dimensional statistic is computed in each subset ${\mathcal{M}}_k$. These statistics are transformed to main machine and then aggregated using an aggregation function, to obtain the final global estimator.

The dc framework is often implemented through two methods, one-shot (os) and communication iterative computation, see Figure 4. Specifically, the os computation refers to a process where tasks are distributed and solved in a single computation, with all nodes working independently before that their results are aggregated [23]. This approach is efficient for problems that do not require intermediate communication or feedback between nodes during processing. For instance, Lee et al. [64] proposed a simple and communication-efficient method for distributed high-dimensional sparse regression based on the os approach.

While the os computation can be efficient for simple or well-structured problems, it lacks efficiency in terms of data utility. If the initial estimation is suboptimal or inaccurate, there is no built-in mechanism for adjustment or improvement. Hence it is highly sensitive to noise, errors, or biases in the input data, which can significantly impact the final result. To improve the statistical efficiency of os estimators, distributed communication iterative methods have been proposed [65,66]. Compared to the os method, communication iterative methods make more comprehensive use of global information, leading to more precise parameter estimation.

In the follow part of this section, we will discuss some seminal research about implementation of QR for distributed computing.

3.1.1. One-Shot Method and Its Variants

os is the simplest distributed method, requiring only that parameters be estimated on local machines and then transmitted back to the master for averaging/aggregate to obtain the final global estimates. Specifically, for each data batch, we first estimate local parameter ${\widehat{\mathit{\beta }}}_{\left(k\right)}$ in (3) using kth batch samples ${\mathcal{M}}_k$. Each local parameter is then sent to the master, which calculates the global estimator by ${\displaystyle \frac{1}{K}}{\sum }_{i=1}^K{\widehat{\mathit{\beta }}}_{\left(k\right)}$, which is denoted as the os estimator ${\widehat{\mathit{\beta }}}_{\mathrm{os}}$.

Extensive research has been proposed in this field for QR. For example, Xu et al. [63] developed a block average approach, employing the os strategy by averaging estimators derived from each local computer.

Chen and Zhou [7] also applied the os strategy by reconstructing the estimator for the entire dataset with an asymptotically negligible approximation error. Moreover, this method not only introduces a resampling technique that simplifies and enhances the efficiency of weight matrix computation, but also yields a global estimator in closed form by combining low-dimensional statistics derived from each data block. This feature further facilitates continuous model updating in streaming data scenarios.

To guarantee statistical efficiency, Volgushev et al. [67] proposed a distributed quantile curve estimation method based on a two-step quantile projection algorithm. First, they estimated specific $\tau$-quantile values (i.e., ${\tau }_m:m=1,\cdots ,M$) in local machine and aggregated them to the master, which outputs a pooled estimator; then, they applied a B-spline method to construct a projection matrix for generating the overall quantile curve. More specifically, let,

$$\begin{array}{cc}\hfill {\widehat{\mathit{\alpha }}}_j& =\underset{\mathit{\alpha }}{argmin}\sum _{m=1}^M{\left({\widehat{\mathit{\beta }}}_{\mathrm{os}}\left({\tau }_m\right)-{\mathit{\alpha }}^{\top }\mathit{B}\right)}^2,\hfill \end{array}$$

where $\mathit{B}={({\mathit{B}}_1,\cdots ,{\mathit{B}}_q)}^{\top }$ is a B-spline basis defined on equidistant knots. Then the quantile curve can be obtained by $\widehat{\mathit{\beta }}\left(\tau \right)={\widehat{\mathsf{\Theta }}}^{\top }\mathit{B}$, where $\widehat{\mathsf{\Theta }}=[{\widehat{\mathit{\alpha }}}_1,\cdots ,{\widehat{\mathit{\alpha }}}_p]$ is the projection matrix. This method can be applied to the inference of conditional distribution functions, enabling researchers to better understand the distributional characteristics of the data. Furthermore, due to its use of subsample estimators for computation, it offers low computational cost, making it more suitable for practical application scenarios.

The drawback of the aforementioned research is that ${\widehat{\mathit{\beta }}}_{\left(k\right)}$ relies solely on local information, which can lead to statistical inefficiency, particularly when the data are not randomly distributed between workers. To address this issue, Pan et al. [68] proposed an updating algorithm. Based on the following asymptotic expression of global estimator,

$$\begin{array}{c}\hfill {\widehat{\mathit{\beta }}}_{\mathrm{os}}={\mathit{\beta }}_0+{\displaystyle \frac{1}{f_{ϵ}\left(0\right)}}{\left(\sum _{i=1}^N{\mathit{x}}_i{\mathit{x}}_i^{\top }\right)}^{-1}\left[\sum _{i=1}^N{\mathit{x}}_i\{\tau -\mathbf{1}({ϵ}_i\le 0)\}\right]+o_p\left({\displaystyle \frac{1}{\sqrt{N}}}\right),\end{array}$$

(10)

where $o_p(·)$ denotes convergence to zero in probability. And the one-step update estimator derived by replace ${\mathit{\beta }}_0,ϵ$ in (10) by a pilot estimator ${\widehat{\mathit{\beta }}}_{\mathrm{pilot}}$ (eg. estimator derived in node ${\mathcal{M}}_1$) and ${\widehat{ϵ}}_i=y_i-{\mathit{x}}_i^{\top }{\widehat{\mathit{\beta }}}_{\mathrm{pilot}}$. That is,

$$\begin{array}{c}\hfill {\widehat{\mathit{\beta }}}_{\mathrm{os}}^{\left(1\right)}={\widehat{\mathit{\beta }}}_{\mathrm{pilot}}+{\displaystyle \frac{1}{{\widehat{f}}_{ϵ}\left(0\right)}}{\left(\sum _{i=1}^N{\mathit{x}}_i{\mathit{x}}_i^{\top }\right)}^{-1}\left[\sum _{i=1}^N{\mathit{x}}_i\{\tau -\mathbf{1}({\widehat{ϵ}}_i\le 0)\}\right],\end{array}$$

(11)

where ${\widehat{f}}_{ϵ}(·)$ is a kernel density estimator.

3.1.2. Communication Iterative Computation

Communication iterative computation, on the other hand, involves multiple iterations of computation and communication. In this approach, nodes repeatedly update their results based on feedback or inputs from other nodes, gradually converging toward the final solution. This strategy utilizes information from the entire dataset, leading to significantly improved accuracy compared to the os strategy. However, communication iterative computation requires data transmission among multiple machines, which can result in high communication costs. Balancing statistical efficiency with reduced communication overhead has therefore become a critical challenge in this area of research. Additionally, in the QR framework, the non-differentiability of its loss function poses another significant challenge. This makes it impossible to directly apply the communication loss function. Consequently, the critical issue in distributed computing for QR is how to approximate or derive the gradient function for a non-differentiable loss function, together with excellent statistical efficiency with well-controlled communication costs. Below we review some related work.

One notable work in this regard is motivated by the communication-efficinet surrogate likelihood method proposed by Jordan et al. [65], and extend it to QR [45,69]. The differences between them are the different techniques for the non-smooth nature of the QR loss function. For example, Hu et al. [69] introduced sub-gradients to construct a surrogate loss function, and established the consistency and asymptotic normality of the proposed methodology. Later, Tan et al. [45] applied the convolution smoothing approach (8) to the local and global objective functions of quantile regression.

Beyond the above work, researchers address the challenges of distributed QR by proposing various loss functions via smoothing techniques introduced in Section 2.2 and iterative algorithms. For example, Chen et al. [41] presented a computationally efficient method based on Horowitz’s smoothed loss in (6), which involves multiple rounds of aggregations. Specifically, for the k-th batch,

$$\begin{array}{cc}\hfill U_k^{\left(t\right)}& =\sum _{i\in {\mathcal{D}}_k}{\mathit{x}}_i\left\{H\left({\displaystyle \frac{y_i-{\mathit{x}}_i^{\top }{\widehat{\mathit{\beta }}}^{(t-1)}}{h}}\right)+\tau -1+{\displaystyle \frac{y_i}{h}}{H}^{\prime }\left({\displaystyle \frac{y_i-{\mathit{x}}_i^{\top }{\widehat{\mathit{\beta }}}^{(t-1)}}{h}}\right)\right\},\hfill \\ \hfill V_k^{\left(t\right)}& =\sum _{i\in {\mathcal{D}}_k}{\mathit{x}}_i{\mathit{x}}_i^{\top }{\displaystyle \frac{1}{h}}{H}^{\prime }\left({\displaystyle \frac{y_i-{\mathit{x}}_i^{\top }{\widehat{\mathit{\beta }}}^{(t-1)}}{h}}\right).\hfill \end{array}$$

(12)

Here $H^{\prime }(·)$ is the derivative of $H(·)$, and ${\widehat{\mathit{\beta }}}^{(t-1)}$ is the estimator from the $(t-1)$-th iteration, and update the estimator at the $\left(t\right)$-th iteration as:

$${\widehat{\mathit{\beta }}}^{\left(t\right)}={\left(\sum _{k=1}^K{V}_k^{(t-1)}\right)}^{-1}\left(\sum _{k=1}^K{U}_k^{(t-1)}\right).$$

(13)

After limited q iterations, the authors show that the statistical efficiency of the final estimator becomes the same as the one computed on the whole data. They also designed a new one-pass algorithm for streaming data. For more details, please refer to Chen et al. [41].

Note that $l_{h,\tau }^{\mathrm{horo}}\left(u\right)$ gains smoothness at the cost of convexity, He et al. [44] and Tan et al. [46] employed the convolution-type smoothing loss in (8), which transforms the loss function into a convex and twice-differentiable form. Their work focuses on statistical inference for QR in the context of the “increasing dimension". For massive data, Jiang and Yu [70] proposed smoothing QR for a distributed system based on the conquer approach and its Taylor expression. Specifically, via (8), the empirical smoothed loss function becomes,

$${\mathcal{S}}_h\left(\mathit{\beta }\right)={\displaystyle \frac{1}{N}}\sum _{i=1}^N{l}_{h,\tau }^{\mathrm{conquer}}(y_i-{\mathit{x}}_i^{\top }\mathit{\beta })={\displaystyle \frac{1}{N}}\sum _{i=1}^N{\int }_{-\infty }^{\infty }{\rho }_{\tau }\left(\nu \right)K_h(\nu -(y_i-{\mathit{x}}_i^{\top }\mathit{\beta }))d\nu ,$$

(14)

and the smoothing estimator is given by ${\widehat{\mathit{\beta }}}_h^{\mathrm{conquer}}={\widehat{\mathit{\beta }}}_h^{\mathrm{conquer}}\left(\tau \right)={argmin}_{\mathit{\beta }}{\mathcal{S}}_h\left(\mathit{\beta }\right)$.

Note that the entire data are located in K local machines, each machine cannot communicate ${\mathcal{S}}_h\left(\mathit{\beta }\right)$ to the master machine, then Jiang and Yu [70] adopted a Taylor expansion of ${\mathcal{S}}_h\left(\mathit{\beta }\right)$ around an initial estimator ${\widehat{\mathit{\beta }}}^{\left(0\right)}$, yielding:

$$\begin{array}{cc}\hfill {\mathcal{S}}_h\left(\mathit{\beta }\right)\approx {\tilde{\mathcal{S}}}_h\left(\mathit{\beta }\right)=& {\mathcal{S}}_h\left({\mathit{\beta }}_0\right)+{\left(\mathit{\beta }-{\widehat{\mathit{\beta }}}^{\left(0\right)}\right)}^{\top }{\mathcal{S}}_h^{\left(1\right)}\left({\widehat{\mathit{\beta }}}_0\right)+\hfill \\ \hfill & {\displaystyle \frac{1}{2}}{\left(\mathit{\beta }-{\widehat{\mathit{\beta }}}^{\left(0\right)}\right)}^{\top }{\mathcal{S}}_h^{\left(2\right)}\left({\widehat{\mathit{\beta }}}_0\right)\left(\mathit{\beta }-{\widehat{\mathit{\beta }}}^{\left(0\right)}\right),\hfill \end{array}$$

(15)

where ${\mathcal{S}}_h^{\left(1\right)}\left(\mathit{\beta }\right)$ and ${\mathcal{S}}_h^{\left(2\right)}\left(\mathit{\beta }\right)$ are the gradient and Hessian matrix:

$$\begin{array}{c}\hfill {\mathcal{S}}_h^{\left(1\right)}\left(\mathit{\beta }\right)=N^{-1}\sum _{i=1}^N{\mathit{x}}_i\left\{\tilde{K}(-{ϵ}_i/h)-\tau \right\},{\mathcal{S}}_h^{\left(2\right)}\left(\mathit{\beta }\right)=N^{-1}\sum _{i=1}^N{\mathit{x}}_i{\mathit{x}}_i^{\top }{K}_h(-{ϵ}_i/h),\end{array}$$

where $\tilde{K}\left(t\right)={\int }_{-\infty }^{t}K\left(u\right)du$. Note that (15) is a convex quadratic loss function, thus under a distributed system, ${\mathcal{S}}_h^{\left(1\right)}\left({\widehat{\mathit{\beta }}}_0\right)$ and ${\mathcal{S}}_h^{\left(2\right)}\left({\widehat{\mathit{\beta }}}_0\right)$ can be rewritten as:

$$\begin{array}{cc}\hfill {\mathcal{S}}_h^{\left(1\right)}\left({\widehat{\mathit{\beta }}}_0\right)& ={\displaystyle \frac{1}{N}}\sum _{k=1}^K{n}_k{\mathcal{S}}_{k,h}^{\left(1\right)}\left({\widehat{\mathit{\beta }}}_0\right)={\displaystyle \frac{1}{N}}\sum _{k=1}^K\sum _{i\in {\mathit{D}}_i}{\mathit{x}}_i\left\{\tilde{K}(-{ϵ}_i^{\left(0\right)}/h)-\tau \right\},\hfill \\ \hfill {\mathcal{S}}_h^{\left(2\right)}\left({\widehat{\mathit{\beta }}}_0\right)& ={\displaystyle \frac{1}{N}}\sum _{k=1}^K{n}_k{\mathcal{S}}_{k,h}^{\left(2\right)}\left({\widehat{\mathit{\beta }}}_0\right)={\displaystyle \frac{1}{N}}\sum _{k=1}^K\sum _{i\in {\mathit{D}}_i}{\mathit{x}}_i{\mathit{x}}_i^{\top }{K}_h(-{ϵ}_i^{\left(0\right)}/h),\hfill \end{array}$$

where ${ϵ}_i^{\left(0\right)}=y_i-{\mathit{x}}_i^{\top }{\widehat{\mathit{\beta }}}^{\left(0\right)}$. To avoid communication cost, similar to Jordan et al. [65], one can replace ${\mathcal{S}}_h^{\left(2\right)}\left({\widehat{\mathit{\beta }}}_0\right)$ with ${\mathcal{S}}_{1,h_1}^{\left(2\right)}\left({\widehat{\mathit{\beta }}}_0\right)$, and the convolution-type smoothing QR can be derived by Newton-Raphson algorithm, i.e., ${\widehat{\mathit{\beta }}}^{(t+1)}={\widehat{\mathit{\beta }}}^{\left(t\right)}-{\left\{{\mathcal{S}}_{h_1,1}^{\left(2\right)}\left({\widehat{\mathit{\beta }}}^{\left(t\right)}\right)\right\}}^{-1}{\mathcal{S}}_h^{\left(1\right)}\left({\widehat{\mathit{\beta }}}^{\left(t\right)}\right)$, where ${\mathcal{S}}_{h_1,1}^{\left(2\right)}(·)$ is the Hessian matrix with bandwidth $h_1$ which selected base on the first machine. Moreover, Jiang and Yu [70] also extended the above distributed smoothing QR estimator in high dimensions by incorporating $l_1$ penalty. For more details, please refer to Jiang and Yu [70].

3.1.3. Distributed High-Dimensional QR

Modern data acquisitions have facilitated the collection of massive and high-dimensional data, here we briefly summarize some work related QR under massive and high-dimensional data [67,70,71,72,73,74]. For example, Volgushev et al. [67] studied a os QR method with fixed or growing dimension; Zhao et al. [71] proposed a bias-correct method [75] to obtain a de-biased estimator before taking the averaging. However, the de-biasing step is computationally time-consuming and suffers the drawbacks of the os method. Wang and Lian [72] extended the communication-efficient surrogate likelihood method proposed by Jordan et al. [65] for distributed QR, while Chen et al. [73] transformed QR loss to a least squares loss and provided a distributed estimator that is both computational and communicatively efficient. Recently, motivated by Chen et al. [73], Wang and Shen [74] extended Tan et al. [45] to high-dimensional setting, proposed a novel method for distributed high-dimensional QR. Overall, different from distributed estimation in Section 3.1.1 and Section 3.1.2, support recovery or variable selection is also a task. One notable practice is to incorporate the penalty function, such as $l_1$ proposed by Tibshirani [76] and smooth clipped absolute deviation (SCAD) proposed by Fan and Li [77].

3.1.4. Distributed Statistical Inference for QR

In addition to estimation, other statistical inference tools, such as hypothesis testing and confidence intervals, also play a crucial role in scientific research and data analysis. These inference tools allow researchers to quantify the uncertainty of the estimators obtained from the data, and then help practitioners interpret the results appropriately. In the above discussion of this section, various distributed estimation methods have been introduced. Thus, we briefly summarize some attempts on the statistical inference for distributed QR.

In fact, there have been potential challenges for the statistical inference of the standard QR. Classical approaches to inference are typically based on estimating its asymptotic variance from the data directly, or conducting bootstrap to approximate the asymptotic distribution. First, estimating the limiting variance requires the choice of a bandwidth parameter, and existing research indicates that classical rules for bandwidth selection need to be adjusted in a dc setting [78]. For example, Pan et al. [68] proved that ${\widehat{\mathit{\beta }}}_{\mathrm{os}}^{\left(1\right)}$ in (11) shared the same asymptotic distribution with the global estimator, i.e., (4). Thus, they performed valid statistical inference with non-parametric estimator of $f_{ϵ}\left(0\right|X)$ and a standard moment estimator of ${\Omega }_x$, which can be estimated easily on a distributed system. Second, the bootstrap approach requires repeatedly computing QR estimates up to thousands of times, and therefore is unduly expensive for massive data inevitably. To solve this problem, Volgushev et al. [67] proposed several simple inference procedures which directly utilize the fact that in a dc setting, estimators from subsamples are available. Their procedures are very easy to implement and require only a very small amount of computation on the central computer without additional communication costs. Please refer to Volgushev et al. [67] for more details.

There are also some work related to large-scale statistical inference for QR without in the distributed framework. For example, He et al. [44] propose a bootstrapped procedure to perform hypothesis testing and confidence estimation with larger dimension and larger sample size. Later, Lee et al. [79] developed a fast inference procedure for QR with ultra-large dataset based on stochastic subgradient descent (S-subGD) updates. The inference procedure handles cross-sectional data sequentially: (i) updating the parameter estimate with each incoming “new observation", (ii) aggregating it as a Polyak–Ruppert average, and (iii) computing a pivotal statistic for inference using only a solution path.

3.2. Subsampling-Based Method for QR

Distributed algorithms provide an effective solution for statistical modeling with massive datasets. However, they face several challenges in practical applications: 1) high cost of distributed computing systems and 2) high learning cost. As an alternative approach, subsampling repeatedly generates small subsamples for parameter estimation and statistical inference, which has been widely used to reduce the computational burden when handling massive data. It performs analysis on a small subsample drawn from the full data and provides a practical solution to extracting information from massive data with limited computing power. Consequently, various subsampling methods are proposed, and has been applied for linear regression [80,81,82], logistic regression [83,84] and generalized linear regression [85,86]. Existing optimal subsampling techniques require a differentiable objective function or loss function, hence traditional subsampling strategies can not be applied directly to QR. Below we describe some work on subsampling-based methods for QR [21,86,87,88].

For example, Wang and Ma [87] studied optimal subsampling procedures and derived the asymptotic distribution of a general subsampling estimator and two versions of optimal subsampling probabilities. Specically, taking a random subsample using sampling with replacement from the full data ${\{y_i,{\mathit{x}}_i\}}_{i=1}^N$ satisfying (3) according to the probabilities ${\pi }_i(i=1,\cdots ,N)$, such that ${\sum }_{i=1}^N{\pi }_i=1$. Denote the subsample as ${\{y_i^{★},{\mathit{x}}_i^{★}\}}_{i=1}^n$, with associated subsampling probabilities ${\pi }_i(i=1,\cdots ,n)$. The loss function of subsampling QR denote as:

$$\begin{array}{c}\hfill Q_n^{\ast }\left(\mathit{\beta }\right)={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{{\rho }_{\tau }(y_i^{\ast }-{\mathit{x}}_i^{\ast \top }\mathit{\beta })}{N{\pi }_i^{\ast }}}.\end{array}$$

(16)

The subsample estimator, denoted as $\tilde{\mathit{\beta }}$, is the minimizer of $Q_n^{\ast }\left(\mathit{\beta }\right)$ The asymptotic normality of $\tilde{\mathit{\beta }}$ has been established, which can be used to identify the optimal probabilities by minimizing the trace of the asymptotic variance-covariance matrix, either for a linearly transformed parameter estimator, or the original parameter estimator. See Theorem 1-3 in Wang and Ma [87]. Later on, subsequent studies have also validated the effectiveness of this approach [89], which proposed a weighted least squares resampling approach for estimating the variance-covariance matrix of QR subsample-based estimator.

The subsampling methods mentioned above primarily focus on the weak convergence of the resulting estimator, which often fails to effectively capture the convergence rate. Consequently, the theoretical framework does not offer clear guidance on determining the appropriate subsample size. To address this, Ai et al. [86] proposed a Poisson sampling QR method for large-scale data, where subsamples are selected through independent Bernoulli trials to approximate the full dataset. This method eliminates the need to calculate sampling probabilities in RAM all at once, thereby significantly reducing computational costs when dealing with extremely large datasets. Beside, when the dataset cannot be stored on a single machine, a single optimal subsampling procedure becomes inadequate for addressing such distributed scenarios. To overcome this issue, Chao et al. [88] proposed a distributed QR subsampling procedure, and developed the optimal subsampling probabilities and subset sizes selection criteria simultaneously via a two-step algorithm. The theoretical results, such as consistency and asymptotic normality of resultant estimators are established.

4. Advances in Quantile Regression for Streaming Data

In contrast to Section 3, here we focus on studies for QR for streaming data. As mentioned in the Section 1, streaming data consist of data blocks arriving sequentially over time, and it requires frequent updates to the analysis result as new data is collected. However, due to the sheer volume of data, it also becomes challenging to store all the information from each batch. As a result, renewable estimation or online updating methods have been proposed, as shown in Figure 1. Here we first introduce some notations used in this section. Assume $D_{\left(b\right)}=\{D_t:t=1,2,\cdots ,b\}$ are aggregated streaming data up to batch b, where $D_t={\left\{(y_{ti},{\mathit{x}}_{ti})\right\}}_{i=1}^{n_i}$ are the streaming data of t-batch with sample size $n_t$ satisfying (3), i.e.,

$$\begin{array}{c}\hfill q_{\tau }\left(y_{ti}\right|{\mathit{x}}_{ti})={\mathit{x}}_{ti}^{\top }\mathit{\beta }\left(\tau \right),i=1,\cdots ,n_t.\end{array}$$

(17)

And the total sample size is $N_b={\sum }_{t=1}^b{n}_t$. For simplicity, in this section, we omit $\tau$ whenever there is no ambiguity, and use “$˜$” to denote the online estimator (i.e., ${\tilde{\mathit{\beta }}}_b$ is the renewable estimator obtained via data up to batch b, $D_{\left(b\right)}$) and “$^$” to denote the estimates derived from data at stationary time points (i.e, ${\widehat{\mathit{\beta }}}_t$ obtained via data from $D_t$), respectively.

The core challenge of streaming QR lies in constructing a renewable estimation framework that allows parameters to be updated incrementally as new data arrive, without the need to store or reprocess the entire dataset. Existing methods for renewable estimation for QR can be categories into three types: (1) the one-shot divide-and-conquer algorithms (os-dc) as well as its variants when viewing stream data as data blocks partitioned over the time domain; (2) sequential QR (sqr) algorithm partly motivated by Bayesian QR; and (3) the cumulative estimating equation together with a Taylor expansion method (cee-taylor). Table 3 summarizes these three types of renewable algorithm, and we will introduce the above three strategies in detail in the following sections.

4.1. os-dc Methods and Its Variants

Note that most of these techniques in Section 3.1.1 and Section 3.1.2 are only applicable to static big data but not stream data except the os-dc algorithms when viewing stream data as data blocks partitioned over the time domain [23,41]. For example, Chen et al. [41] applied Horowitz’s smoothing QR for streaming data by using the dc strategy via a one-pass algorithm (see Algorithm 2 in Chen et al. [41]) to provide a sequence of successively refined estimators. While the os-dc algorithms can be naturally extended to stream data by implementing the aggregation in a block-by-block manner, it often requires the block size to grow to infinity to establish their asymptotic theory. However, the block size in stream data is generally determined by hardware capacity and hence fixed. In addition, updates based on direct applications of os-dc algorithms may not be efficient enough to support the frequent update in the stream data context.

As a straightforward extension of the os-dc technique, Wang et al. [11] developed an updating estimator based on the normal approximation of the QR estimator from a single batch of data. Based on its asymptotic distribution, the local QR estimator $\{{\widehat{\mathit{\beta }}}_1,\cdots ,{\widehat{\mathit{\beta }}}_b\}$ can be treated as sample generated from the following multivariate distribution:

$$\begin{array}{c}\hfill \left\{N({\mathit{\beta }}_0,{\displaystyle \frac{\Delta }{n_1}}{\Sigma }_1^{-1}),\cdots ,N({\mathit{\beta }}_{\mathbf{0}},{\displaystyle \frac{\Delta }{n_b}}{\Sigma }_b^{-1})\right\},\end{array}$$

where $\Delta ={\displaystyle \frac{\tau (1-\tau )}{f_{ϵ}^2\left(0\right)}}$ with $f_{ϵ}(·)$ being the pdf of error term, and ${\displaystyle \frac{1}{n_t}}{\sum }_{i=1}^{n_t}{\mathit{x}}_{\mathit{ti}}{\mathit{x}}_{\mathit{ti}}^{\top }\stackrel{p}{\to }{\Sigma }_t$ for $t=1,\cdots ,b$. Then, based on the multivariate normal distribution, the renewable estimator ${\tilde{\mathit{\beta }}}_{\left(b\right)}$ can be formulated as:

$$\begin{array}{c}\hfill {\tilde{\mathit{\beta }}}_{\left(b\right)}={\left(\sum _{t=1}^b{\mathit{X}}_t^{\top }{\mathit{X}}_t\right)}^{-1}\left(\sum _{t=1}^b{\mathit{X}}_t^{\top }{\mathit{X}}_t{\widehat{\mathit{\beta }}}_t\right),\end{array}$$

(18)

where ${\mathit{x}}_t={({\mathit{x}}_{t1},\cdots ,{\mathit{x}}_{t{n}_t})}^{\top }$. Obviously, in order to run the renewable algorithm (18), one only needs to store $\{{\mathit{x}}_t^{\top }{\mathit{x}}_t,{\widehat{\mathit{\beta }}}_t\}$ for each batch of data, while the historical data are not completely involved. The idea of Wang et al. [11] has been extended to other complex regression problems, such as expectile regression [96], Huber robust regression [97,98], unconditional quantile regression [99] and composite quantile regression [95].

Unlike Wang et al. [11], Chen and Zhou [7] also proposed a simple and efficient dc-based strategy for the estimation of the QR model with massive datasets, which only keeps compact statistics of each data block and uses them to reconstruct the estimator of the entire data, with an asymptotically negligible approximation error. This property makes the proposed method particularly appealing whenever the data blocks are retained in multiple servers or received in the form of data stream. Specifically, for the data up to batch b, the proposed renewable estimator ${\tilde{\mathit{\beta }}}_{\left(b\right)}$ can be formulated in the closed-form:

$$\begin{array}{c}\hfill {\tilde{\mathit{\beta }}}_{\left(b\right)}={\left(\sum _{t=1}^b{\displaystyle \frac{n_t}{N_b}}{\widehat{A}}_t\right)}^{-1}\left(\sum _{t=1}^b{\displaystyle \frac{n_t}{N_b}}{\widehat{A}}_t{\widehat{\mathit{\beta }}}_t\right),\end{array}$$

(19)

where ${\widehat{A}}_t$ is the estimator of $A_t=E\left[X{X}^{\top }{f}_{Y|X}\left(X^{\top }{\widehat{\mathit{\beta }}}_t\right|X)\right]$, which is estimated via a resampling method inspired by Zeng and Lin [100]. Here $f_{Y|X}\left(y\right|\mathit{x})$ is the conditional density function of Y given X.

4.2. sqr Algorithm

Note that local processing in Zhang et al. [23] and Wang et al. [11] requires solving a QR problem, which results in updates that are too slow for real-time processing of stream data. Recall that Bayesian QR uses a linear model with an ALD error term to achieve equivalence between the original QR estimator and the Bayesian maximum-a-posteriori (MAP) estimator [47]. Additionally, the scale-mixture representation of the ALD holds [48]. Motivated by these foundations, Fan and Lin [12] proposed a sqr algorithm, where sequential updates occur naturally within the Bayesian setup. Specifically, the posterior distribution obtained from the past analysis is used as the prior information for later analysis.

The data block arriving at time point $t(t=1,2,\cdots )$ is denoted as $D_t$. They first converted the check loss optimization based on each data block to a least squares problem to expedite local processing to accommodate the frequent update in stream data, resulting in ${\widehat{\mathit{\beta }}}_t$; the sqr estimator is then updated up to batch b, denoted as ${\tilde{\mathit{\beta }}}_{\left(b\right)}$. This estimator is calculated by solving a nonzero centered ridge regression problem [101], which has a closed-form solution given as follows:

$$\begin{array}{c}\hfill {\tilde{\mathit{\beta }}}_{\left(b\right)}={\displaystyle \frac{n_b}{{\sum }_{t=1}^b{n}_t}}{\widehat{\mathit{\beta }}}_b+{\displaystyle \frac{{\sum }_{t=1}^{b-1}{n}_t}{{\sum }_{t=1}^b{n}_t}}{\tilde{\mathit{\beta }}}_{(b-1)}.\end{array}$$

(20)

From a theoretical perspective, the proposed estimator is asymptotically normal and achieves a linear convergence rate under mild conditions. Additionally, the computational complexity is significantly reduced from $O\left(p^2\right)$ to $O\left(p\right)$, offering a notable improvement over distributed algorithms. The general idea of sqr has also been extended to composite quantile regression [102]. Please see Algorithm 1 in Fan and Lin [12] for more details.

Similarly, Tian et al. [90] also proposed a Bayesian QR approach for streaming data, in which the posterior distribution was updated by utilizing the aggregated statistics of current and historical data. Leveraging the posterior distribution of historical data as the prior for current data stream, a posterior distribution is established without any loss of information, which ensures that the proposed Bayesian streaming QR estimator is equivalent to the global oracle estimator calculated based on the entire data.

4.3. cee-taylor-Based Renewable Estimation

For massive data, the cumulative estimating equations estimator was proposed first by Lin and Xi [15] and subsequently adapted to the sequential estimation setting by Schifano et al. [16], yet their estimation consistency may not be guaranteed in the situation where streaming datasets arrive perpetually with $b\to \infty$. To overcome this unnatural limitation, Luo and Song [19] presented an incremental updating algorithm to analyze streaming datasets using generalized linear models. The main procedure can be summarized as follows.

Denote ${\mathcal{L}}_{N_b}\left(\mathit{\beta }\right)$ and ${\mathcal{L}}_{n_b}\left(\mathit{\beta }\right)$ as the general aggregated loss function based on data up to batch b (i.e., ${\mathcal{D}}_{\left(b\right)}$) and data at batch b (i.e., ${\mathcal{D}}_b$), respectively. First, decompose ${\mathcal{L}}_{N_b}\left(\mathit{\beta }\right)$ into a combination of the historical data and the arrival of new data, i.e.,

$$\begin{array}{c}\hfill {\mathcal{L}}_{N_b}\left(\mathit{\beta }\right)={\displaystyle \frac{1}{N_b}}\left(N_{b-1}{\mathcal{L}}_{N_{b-1}}\left(\mathit{\beta }\right)+n_b{\mathcal{L}}_{n_b}\left(\mathit{\beta }\right)\right).\end{array}$$

(21)

Then, by applying Taylor expansion to the first part of (21) around the online estimates up to $b-1$ batches, ${\tilde{\mathit{\beta }}}_{(b-1)}$, ${\mathcal{L}}_{N_b}\left(\mathit{\beta }\right)$ can be asymptotically equivalent to (omitting terms that do not contain the unknown parameter):

$$\begin{array}{c}\hfill {\tilde{\mathcal{L}}}_{N_b}\left(\mathit{\beta }\right)={\displaystyle \frac{1}{N_b}}\left\{{\displaystyle \frac{1}{2}}{(\mathit{\beta }-{\tilde{\mathit{\beta }}}_{b-1})}^{\top }\left[\sum _{t=1}^{b-1}{n}_t{\nabla }^2{\mathcal{L}}_{n_t}\left({\tilde{\mathit{\beta }}}_t\right)\right](\mathit{\beta }-{\tilde{\mathit{\beta }}}_{b-1})+n_b{\mathcal{L}}_{n_b}\left(\mathit{\beta }\right)\right\},\end{array}$$

(22)

where ${\nabla }^2{\mathcal{L}}_{n_t}(·)$ is the second order function of ${\mathcal{L}}_{n_t}(·)$ of t-th batches. Then the renewable estimation based on data up to batch b is defined as a minimizer of (22). We adopt the term cee-taylor-based renewable estimator for the above procedure.

Extensive researches have extended the idea of cee-taylor to QR for streaming data by incorporating the smoothing techniques to ensure that the objective function remains both differentiable and convex [91,92,93,94]. For example, both Jiang and Yu [91] and Xie et al. [93] proposed a renewable method for QR via the convolution-type smoothing objective function:

$${\mathcal{S}}_{N_b,h}\left(\mathit{\beta }\right)={\displaystyle \frac{1}{N_b}}\sum _{t=1}^b{n}_t{\mathcal{S}}_{n_t,h}\left(\mathit{\beta }\right)={\displaystyle \frac{1}{N_b}}\sum _{t=1}^b\sum _{i=1}^{n_t}{l}_{h,\tau }^{\mathrm{conquer}}(y_i-{\mathit{x}}_i^{\top }\mathit{\beta }),$$

where ${\mathcal{S}}_{n_t,h}\left(\mathit{\beta }\right)={\displaystyle \frac{1}{n_t}}{\sum }_{i=1}^{n_t}{l}_{h,\tau }^{\mathrm{conquer}}(y_i-{\mathit{x}}_i^{\top }\mathit{\beta })$ with $l_{h,\tau }^{\mathrm{conquer}}(·)$ defined in (8). Then the proposed renewable estimator at bacth b can be obtained according to (21) and (22) with ${\mathcal{L}}_{N_b}\left(\mathit{\beta }\right)$ replaced by ${\mathcal{S}}_{N_b,h}\left(\mathit{\beta }\right)$. That is, the the renewable estimator up to batch b, ${\tilde{\mathit{\beta }}}_b$, is given by:

$$\begin{array}{c}\hfill {\tilde{\mathit{\beta }}}_b=\underset{\mathit{\beta }}{argmin}{\displaystyle \frac{1}{N_b}}\left\{{\displaystyle \frac{1}{2}}{(\mathit{\beta }-{\tilde{\mathit{\beta }}}_{b-1})}^{\top }\left[\sum _{t=1}^{b-1}{n}_t{\nabla }^2{\mathcal{S}}_{n_t,h}\left({\tilde{\mathit{\beta }}}_t\right)\right](\mathit{\beta }-{\tilde{\mathit{\beta }}}_{b-1})+n_b{\mathcal{S}}_{n_b,h}\left(\mathit{\beta }\right)\right\}.\end{array}$$

(23)

Moreover, Xie et al. [93] proposed a online high-dimensional smoothed-type QR algorithm. For variable selection, this paper replace ${\tilde{\mathit{\beta }}}_t$ in (23) with the following lasso estimator:

$$\begin{array}{c}\hfill {\tilde{\mathit{\beta }}}_b\left({\lambda }_b\right)=\underset{\mathit{\beta }}{argmin}\left[{\displaystyle \frac{1}{N_b}}\left\{{\displaystyle \frac{1}{2}}{(\mathit{\beta }-{\tilde{\mathit{\beta }}}_{b-1})}^{\top }\left[\sum _{t=1}^{b-1}{n}_t{\nabla }^2{\mathcal{S}}_{n_t,h}\left({\tilde{\mathit{\beta }}}_t\right)\right](\mathit{\beta }-{\tilde{\mathit{\beta }}}_{b-1})+n_b{\mathcal{S}}_{n_b,h}\left(\mathit{\beta }\right)\right\}+{\lambda }_b\parallel \mathit{\beta }\parallel \right],\end{array}$$

where ${\lambda }_b$ is the regularization parameter adaptively selected after the arrival of the b-th batch data. Additionally, an online debiased lasso procedure is developed for inference, where the debiasing matrix is computed using only the current data and summary statistics from previous batches. For theoretical respect, the proposed online estimators enable efficient estimation and valid inference. For further details, refer to Xie et al. [93].

In addition to the convolution-type smoothing objective function, Horowitz’s smoothing method could also be adopted for QR renewable estimation. Note that $l_{h,\tau }^{\mathrm{horo}}\left(u\right)$ achieves smoothness at the cost of convexity. To address this issus, Sun et al. [92] introduced a convex smooth quantile loss, which is infinitely differentiable and converges uniformly to the quantile loss. Specifically, they modified $l_{h,\tau }^{\mathrm{horo}}\left(u\right)$ by adding $hg(-u/h)$and proposed the following loss

$${\tilde{l}}_{h,\tau }^{\mathrm{horo}}\left(u\right)=l_{h,\tau }^{\mathrm{horo}}\left(u\right)+hg(-u/h),$$

(24)

where $g(·)$ is nonnegative bounded function satisfying ${\displaystyle \frac{u}{h}}{H}^{\prime }(-u/h)=g^{\prime }(-u/h)$. Then following the same procedure in CEE-Taylor-based renewable estimation described in (21) and (22) with ${\mathcal{L}}_{N_b}\left(\mathit{\beta }\right)$ replaced by the smoothing loss function in (24). Later, Peng and Wang [94] also adopted the Horowitz’s smoothing method to propose a two-stage online debiased lasso estimation and statistical inference method for high-dimensional QR models in the presence of streaming data. For more details, please refer to Sun et al. [92] and Peng and Wang [94].

4.4. Renewable Estimation for Streaming Quantile Eegression Under Heterogeneity Scenario

Most existing online learning methods are based on homogeneity assumptions, which require the samples in a sequence to be independent and identically distributed (i.i.d.). However, as data are generated over time, maintaining homogeneity between different batches becomes increasingly challenging, especially as the growth of number of batches and intervals between data collections. This suggests that different batches may not necessarily originate from the same underlying model. In fact, batch correlation and dynamically evolving batch-specific effects are key characteristics of real-world streaming data, such as electronic health records and mobile health data. As a result, research into online updating methods that address heterogeneity has become an emerging and significant area of study.

For example, Luo and Song [103] studied multivariate online regression analysis with heterogeneous streaming data for linear regression, while Wei et al. [104] proposed a new online estimation procedure to analyze the constantly incoming streaming datasets in the framework of generalized linear models with dynamic regression coefficients to address potential heterogeneity. Recently, Ding et al. [105] studied renewable risk assessment of heterogeneous streaming time-to-event cohorts, where risk estimates for each cohort of interest are continuously updated using the latest observations and historical summary statistics. They introduced parameters to quantify the potential discrepancy between batch-specific effects from adjacent cohorts and applied penalized estimation techniques to identify nonzero discrepancy parameters, adaptively adjusting risk estimates based on current data and historical trends.

Similarly, for QR with heterogeneous streaming data, Chen and Yuan [13] developed an online renewable QR method with streaming datasets, respectively under the homogeneity and heterogeneity assumptions. The renewable strategy is based on the cee-taylor approach, and a homogeneity detection method was proposed to identify potential homogeneous parameters that are the same for all data batches through penalized estimation.

5. Application

In the previous sections, we introduce some seminal work for QR for large-scale data, including the model setup, parameter estimation algorithm, and statistical inference. In this section, we aim to introduce the application fields for large-scale QR in detail and introduce some benchmark datasets.

5.1. Meteorological Forecasting

One key application field of large-scale QR is predicting air quality or pollutant emissions. Meteorological data typically exhibit the following characteristics: First, such data are often collected from various monitoring stations, resulting in decentralized storage. Second, because sensors are generally placed in diverse geographical locations, the collected data tend to exhibit heterogeneity due to natural geographical differences. Lastly, meteorological data are typically collected at high frequencies, such as hourly or even by the minute, leading to a massive volume of data. Consequently, retaining all historical data may sometimes be infeasible, creating a streaming data scenario.

Chen and Zhou [7] applied distributed computing methods to study the factors influencing PM2.5 emissions in Beijing. Their dataset contained hundreds of thousands of observations. The empirical results indicated that distributed algorithms could achieve the same accuracy as oracle (pooled) estimation while significantly reducing data transmission and computational costs. Furthermore, QR enabled the model to better capture the heterogeneous effects of independent variables on the response variable. In contrast, Wang et al. [11] employed a renewable approach to build a QR model for forecasting greenhouse gas (GHG) emissions. Below we list a dataset, sourced from the University of California, Irvine Machine Learning Repository (UCI/GHG), contained over a million data points. Due to large volume of dataset, storing all historical data for subsequent analysis was impractical. Empirical results demonstrated that the online updating method performed comparably to the oracle QR.

5.2. Demand and Price Forecasting

The second empirical research scenario for large-scale QR is in the consumer market forecasting. Compared to other fields, the consumer market has the following characteristics: First, demand in consumer markets is often highly asymmetric, significantly influenced by extreme events such as promotions and holidays, which can lead to substantial fluctuations. Second, consumer market data typically comes from multiple sources, including online transactions, offline sales, social media feedback, and market surveys, resulting in enormous data volumes. Lastly, consumer markets are susceptible to outliers, such as seasonal changes or sudden events, which can significantly affect forecast results. As mentioned earlier, QR can capture different quantiles of distributions, providing a more precise description of demand fluctuations. Hence, it is generally more robust compared to conditional mean regression.

Pan et al. [68] applied subsampling method to predict vehicle prices based on QR. The dataset is from the Kaggle US Used Cars dataset (Kaggle). In the empirical analysis, 28 factors were selected to assess car prices and evaluate computational costs. The results showed that even with millions of data points transmitted during the analysis, the communication cost remained acceptable. Fan and Lin [12] applied the online updating algorithm to forecast prices using the same dataset. The results were consistent with reality experiences, cars with higher horsepower and lower mileage are more likely to be sold at higher prices. Additionally, by comparing the computational speed between online updating method and distributed computing, the online updating method reduced computation time by $90\%$. Xie et al. [93] used the Seoul Bike Sharing Demand dataset for demand forecasting with renewable estimation algorithm (UCI/bike). The study found that temperature had a negative effect on bike rentals during the winter (lower quantiles) but a positive effect in other seasons. This is expected because people are less likely to rent bikes when temperatures are relatively low in winter.

5.3. Geographical Distribution Data

Another common application of large-scale QR is geographical data, where data storage distributed by geographical divisions. Such data typically comes from multiple sources, and with distributed algorithms research can build model in parallel, improving computational efficiency. Additionally, because of the natural geographical segmentation, data collection often faces uneven spatial distribution and heterogeneity. QR can significantly enhance the robustness of model fitting.

Huang et al. [106] used an iterative communication distributed algorithm to conduct large-scale QR on wine quality data (UCI/wine) and geological data from Xuchang. The study found that, due to the heterogeneity of data from different geographical locations, the iterative distributed quantile regression performed significantly better than the one-shot method and achieved prediction accuracy comparable to the oracle method.

5.4. High-Frequency Data Forecasting

The fourth application of large-scale QR is in high-frequency data, characterized by large volumes, rapid changes, and high real-time demands. This type of data is commonly observed in aviation, financial markets, and electrical load monitoring. For instance, aviation data includes real-time flight positions, speeds, weather conditions, and delays, all of which require timely analysis for effective decision-making in the aviation industry. In financial markets, stock prices, trading volumes, foreign exchange rates, and futures prices can fluctuate on a millisecond basis. Financial institutions rely on analyzing such data for high-frequency trading, risk management, and market monitoring. These data types often exhibit significant uncertainty and volatility, which traditional mean regression methods fail to adequately capture. In contrast, QR provides a more comprehensive view of the data distribution, allowing it to account for the effects of outliers and extreme events more effectively.

Sun et al. [92] built a prediction model on household electricity usage prediction, with data sourced from UCI/electric. The empirical results showed that the prediction results with renewable algorithm achieved same accuracy as oracle estimation and even with very small batch sizes, the final model prediction accuracy remained largely unaffected. Xie et al. [93] applied the proposed method to assess the stock returns of the CSI 300 Index. The prediction results showed that the renewable estimator produced stable results and led to reliable conclusions. The study also suggested that, in practice, financial managers may track stocks over a longer period and cautiously construct portfolios to ensure stability, and that the online updating method could be an effective tool for analyzing stock datasets. Chen and Yuan [13] studied the influence factors of flight delays, the data available at www.transtats.bts.gov. The empirical study demonstrated that the online updating algorithm effectively reduced the computational costs of parameter estimation. Additionally, the proposed heterogeneity detection method was able to capture dynamic parameter changes between different batches, better uncovering the heterogeneous information in the data.

6. Concluding Remarks and Future Research

This paper provides a selective review of QR methods for large-scale data, addressing both static datasets with massive size (massive data) and streaming data. First, we revisited some foundational work on QR, focusing on the common smoothing techniques to overcome the non-differentiability of the QR loss function. Next, we reviewed the development of QR for massive data, which includes one-shot method and its variant, communicative iterative methods and subsampling-based methods. In addition to estimation techniques, we also summarized developments in distributed high-dimensional QR and distributed statistical inference. Subsequently, we discussed three types of online updating techniques of QR for streaming data, covering the os-dc methods and its variants, the sqr algorithm and the cee-taylor-based renewable estimations. Finally, we introduced application areas of QR and highlighted some commonly used benchmark datasets. It is important to note that this review primarily focuses on linear QR models. Our primary concern is the implementation of parameter estimation for large-scale data, rather than the complexity of model structures.

To conclude this paper, we discuss some limitations about the existing work and potential future research directions.

Firstly, while distributed computing can alleviate the computational burden of central estimation, it faces significant scalability challenges as the growth of data size and the increasing number of local machines. A central server may require a high-quality communication network to manage limited communication resources, resulting in higher communication and infrastructure costs between the server and local workers. Additionally, maintaining direct connections between the central server and each machine becomes impractical as the number of local machines escalate. A promising solution to this challenge is to optimize the distributed network. For instance, Elgabli et al. [107] proposed a decentralized framework for distributed computing. The key advantage of this algorithm is that the competition for communication resources is limited to only half of the workers. Instead of communicating directly with the central server, each worker exchanges its locally trained model with only two neighboring workers, thereby significantly reducing communication costs.

Secondly, most existing studies focus on scenarios with a fixed p. However, in the era of big data, datasets with extremely high dimensions have been becoming increasingly common, which exceeds the storage capacity of a single machine. In such situations, traditional penalized regression methods may fail to perform effectively. Wang et al. [108] proposed a distributed algorithm designed for feature-partitioned distributed systems. This approach first identifies a decorrelated structure in the data using a dc framework and subsectionly fits a penalized regression model to the decorrelated data. This method can be extended to QR, providing a more robust solution for managing extremely high-dimensional, large-scale datasets.

The third concern involves data privacy. In a distributed system relying on a central server, information stored on local nodes may be vulnerable to attackers targeting the central server. Therefore, the developement of reliable distributed systems for data protection is essential. One potential solution involves encoding or encrypting the information or statistics shared by local nodes, which plays a critical role in enhancing data security. Another promising approach is the application of transfer learning [109] or federated learning [106]. These methods can either identify informative patterns similar to the target or develop models independently on local nodes to enhance the target model. By avoiding direct data transmission, these techniques could significantly reduce privacy risks.

In the fourth place, the growing complexity of data structures presents another significant challenge. The studies reviewed in this paper predominantly address linear QR. However, modern datasets often exhibit intricate and non-linear structures that exceeds the capability of traditional linear approaches, underscoring the need for more advanced statistical methodologies to effectively capture and analyze such complexities. Wang et al. [110] proposed a nonparametric QR method based on random features. Specifically, distributed data is directly mapped into a randomized lower-dimensional feature space using a kernel function that captures the complex structures of the model. This transformation enables the use of standard linear QR with a ridge penalty for subsequent analysis. Moreover, the inherent complex structure of the data, such as the censored observations, presents additional challenges for linear QR in large-scale settings. Sit and Xing [111] proposed a censored QR method designed for distributed settings. This approach combines the linear QR estimator [41] with inverse probability weighting (IPW) to handle distributed censored data. However, a limitation of this method is that the IPW weights are calculated using only local data, which could potentially reduce the efficiency as it fails to incorporate global information.

Finally, while QR effectively captures heterogeneity in data structures and renewable estimation facilitates the analysis of dynamic data sets, challenges remain to address the heterogeneity of model structures across different data batches, as discussed in Section 4.4. For instance, Chen and Yuan [13], Wei et al. [104] utilized adaptive lasso penalization to identify nonzero components between two sequential streaming estimators, thereby highlighting the presence of dynamic coefficients. Similarly, Luo and Song [103] developed a real-time regression analysis method based on the Kalman filter, which simultaneously updates point estimates and detects dynamic hidden effects across data batches. This strategy can be further extended to dynamic regression methods for large-scale QR problems.