A review of distributed statistical inference

Abstract

The rapid emergence of massive datasets in various fields poses a serious challenge to traditional statistical methods. Meanwhile, it provides opportunities for researchers to develop novel algorithms. Inspired by the idea of divide-and-conquer, various distributed frameworks for statistical estimation and inference have been proposed. They were developed to deal with large-scale statistical optimization problems. This paper aims to provide a comprehensive review for related literature. It includes parametric models, nonparametric models, and other frequently used models. Their key ideas and theoretical properties are summarized. The trade-off between communication cost and estimate precision together with other concerns is discussed.

Keywords:

• Distributed computing
• divide-and-conquer
• communication-efficiency
• shrinkage methods
• nonparametric estimation
• principal component analysis
• feature screening
• bootstrap

1. Introduction

With the rapid development of information technology, datasets of massive sizes become increasingly available. E-commerce companies like Amazon have to analyse billions of transaction data for personalized recommendation. Bioinformatics scientists need to locate relevant genes corresponding to some specific phenotype or disease from massive SNPs data. For Internet-related companies, large amounts of text, image, voice, and even video data are in urgent need of effective analysis. Due to the accelerated growth of data sizes, the computing power and memory of one single computer are no longer sufficient. Constraint on network bandwidth and other privacy or security considerations also make it difficult to process the whole data on one central machine. Accordingly, distributed computing systems become increasingly popular.

Similar to parallel computing executed on a single computer, distributed computing is closely related to the idea of divide-and-conquer. Simply speaking, for some statistical problems, we can divide a complicated large task into many small pieces so that they can be tackled simultaneously on multiple CPUs or machines. Their outcomes are then aggregated to obtain the final result. It is conceivable that this procedure can save the computing time substantially if the algorithm can be executed in a parallel way. The main difference between a traditional parallel computing system and a distributed computing system is the way they access memory. For parallel computing, different processors can share the same memory. Consequently, they can exchange information with each other in a super-efficient way. While for distributed computing, distinct machines are physically separated. They are often connected by a network. Accordingly, each machine can only access its own memory directly. Therefore, the inter-machine communication cost in terms of time spending could be significant and thus should be prudently considered.

The rest of this article is organized as follows. Section 2 studies parametric models. Section 3 focuses on nonparametric methods. Section 4 expresses some other related methods. The article is concluded with a short discussion in Section 5.

2. Parametric models

Assume a total of N observations denoted as $Z_i=(X_i^{\mathrm{\top }},Y_i{)}^{\mathrm{\top }}\in {\mathbb{R}}^{p+1}$ with $1\le i\le N$. Here $X_i\in {\mathbb{R}}^p$ is the covariate vector and $Y_i\in \mathbb{R}$ is the corresponding scalar response. Define $\{{\mathbb{P}}_{\mathit{\theta }}:\mathit{\theta }\in \mathrm{\Theta }\}$ to be a family of statistical models parameterized by $\mathit{\theta }\in \mathrm{\Theta }\subset {\mathbb{R}}^p$. We further assume that $Z_i$'s are independent and identically distributed with the distribution ${\mathbb{P}}_{{\mathit{\theta }}^{\ast }}$, where ${\mathit{\theta }}^{\ast }=({\theta }_1^{\ast },\dots ,{\theta }_p^{\ast }{)}^{\mathrm{\top }}$ is the true parameter. Consider a distributed setting, where N sample units are allocated randomly and evenly to K local machines ${\mathcal{M}}_k,1\le k\le K$, such that each machine has n observations. Obviously, we should have $N=nK$. Write $\mathbb{S}=\{1,\dots ,N\}$ as the index set of whole sample. Then, let ${\mathcal{S}}_k$ denote the index set of local sample on ${\mathcal{M}}_k$ with ${\mathcal{S}}_{k_1}\cap {\mathcal{S}}_{k_2}=set$ for any $k_1\ne k_2$. Other than the local machines, there also exists a central machine represented by ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$. A standard architecture should have ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$ to be connected with every ${\mathcal{M}}_k$.

Let $\mathcal{L}:\mathrm{\Theta }\times {\mathbb{R}}^{p+1}\mapsto \mathbb{R}$ be the loss function. Assume that the true parameter ${\mathit{\theta }}^{\ast }$ minimizes the population risk ${\mathcal{L}}^{\ast }\left(\mathit{\theta }\right)=\mathbb{E}\left[\mathcal{L}\right(\mathit{\theta };Z\left)\right]$, where $\mathbb{E}$ stands for expectation with respect to ${\mathbb{P}}_{{\mathit{\theta }}^{\ast }}$. Define the local loss on the kth machine as ${\mathcal{L}}_k\left(\mathit{\theta }\right)=n^{-1}\sum _{i\in {\mathcal{S}}_k}\mathcal{L}(\mathit{\theta };Z_i)$. Correspondingly, define the global loss function based on the whole sample as $\mathcal{L}\left(\mathit{\theta }\right)=N^{-1}\sum _{i\in \mathbb{S}}\mathcal{L}(\mathit{\theta };Z_i)=K^{-1}\sum _{k=1}^K{\mathcal{L}}_k\left(\mathit{\theta }\right)$, whose minimizer is $\hat{\mathit{\theta }}={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{\mathit{\theta }\in \mathrm{\Theta }}\mathcal{L}\left(\mathit{\theta }\right)$. In most cases, the whole sample estimator $\hat{\mathit{\theta }}$ should be $\sqrt{N}$-consistent and asymptotically normal (Lehmann & Casella, 2006). If N is small enough so that the whole sample $\mathbb{S}$ can be read into the memory of one single computer, then $\hat{\mathit{\theta }}$ can be easily computed. The entire computation can be executed in the memory of this computer. On the other hand, if N is too large so that the whole sample $\mathbb{S}$ cannot be placed on any single computer, then a distributed system must be used. In this case, the whole sample estimator $\hat{\mathit{\theta }}$ is no longer computable (or at least very difficult to compute) in practice. Then, how to develop novel statistical methods for distributed systems becomes a problem of great importance.

2.1. One-shot approach

To solve the problems, various methods have been proposed. They can be roughly divided into two classes. The first class contains so-called one-shot methods. They are to be reviewed in this subsection. The other class contains various iterative methods. They are to be reviewed in the next subsection.

The basic idea of the one-shot approach is to calculate some relevant statistics on each local machine. Subsequently, these statistics are sent to a central machine, where they are assembled into the final estimator. The most popular and direct way of aggregation is simple average. Specifically, for each $1\le k\le K$, machine ${\mathcal{M}}_k$ uses local sample ${\mathcal{S}}_k$ to compute the local empirical minimizer ${\hat{\mathit{\theta }}}_k={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{\mathit{\theta }\in \mathrm{\Theta }}{\mathcal{L}}_k\left(\mathit{\theta }\right)$. These local estimates (i.e., ${\hat{\mathit{\theta }}}_k$'s) are then transferred to the centre machine ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$, where they are averaged as $\overline{\mathit{\theta }}=K^{-1}\sum _{k=1}^K{\hat{\mathit{\theta }}}_k$. This leads to the final simple averaging estimator $\overline{\mathit{\theta }}$ (see Figure 1(a)).

Figure 1. Illustrations of the two different approaches. (a) one-shot approach and (b) iterative approach.

Obviously, the one-shot style of distributed framework is highly communication-efficient. Because it requires only one single round of communication between each ${\mathcal{M}}_k$ and ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$. Hence, the communication cost is of the order $O\left(Kp\right)$, where p is the dimension of each estimate $\hat{\mathit{\theta }}$. Theoretical properties of simple averaging estimator were also studied in the literature. For example, it was shown in Zhang et al. (2013, Corollary  2) that, under appropriate regularity conditions, (1) $\mathbb{E}\Vert \overline{\mathit{\theta }}-{\mathit{\theta }}^{\ast }{\Vert }_2^2\le \frac{C_1}{N}+\frac{C_2}{n^2}+O(\frac{1}{Nn}+\frac{1}{n^3}),$(1) where $C_1,C_2$ are some positive constants. If n is sufficiently large such that $n^{-2}=o\left(N^{-1}\right)$, then the dominant term in (1(1) $\mathbb{E}\Vert \overline{\mathit{\theta }}-{\mathit{\theta }}^{\ast }{\Vert }_2^2\le \frac{C_1}{N}+\frac{C_2}{n^2}+O(\frac{1}{Nn}+\frac{1}{n^3}),$(1) ) becomes $C_1/N$, and is of the order $O\left(N^{-1}\right)$. It is the same as that of the whole sample estimator. This also implies that, in order to obtain the global convergence rate, we should not divide the whole sample into too many parts. A further improved theoretical result was obtained by Rosenblatt and Nadler (2016). They showed that the one-shot estimator is the first order equivalent to the whole sample estimator. However, the second-order error terms of $\overline{\mathit{\theta }}$ can be non-negligible for nonlinear models. Similar observation was also obtained by Huang and Huo (2015). The work of Duchi et al. (2014) revealed that the minimal communication budget to attain the global estimation error for linear regression is $\mathcal{O}\left(Kp\right)$ bits up to a logarithmic factor under some conditions. This result matches the simple averaging procedure and confirms the sharpness of the bound in (1(1) $\mathbb{E}\Vert \overline{\mathit{\theta }}-{\mathit{\theta }}^{\ast }{\Vert }_2^2\le \frac{C_1}{N}+\frac{C_2}{n^2}+O(\frac{1}{Nn}+\frac{1}{n^3}),$(1) ). To further reduce the bias, a novel subsampling method was developed by Zhang et al. (2013). By this technique, the error bound is improved to be $O(N^{-1}+n^{-3})$, which relaxes the restriction on the number of machines.

Instead of the linear combination of local maximum likelihood estimates (MLEs) as simple average, Liu and Ihler (2014) proposed a KL-divergence based combination method. The final estimator is computed by ${\hat{\mathit{\theta }}}_{\mathrm{K}\mathrm{L}}={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{\mathit{\theta }\in \mathrm{\Theta }}\sum _{k=1}^K\mathrm{KL}(p\left(x|{\hat{\mathit{\theta }}}_k\right)\Vert p\left(x|\mathit{\theta }\right)),$where $p\left(x|\mathit{\theta }\right)$ is the probability density of ${\mathbb{P}}_{\mathit{\theta }}$ with respect to some proper measure μ, and KL-divergence is defined by $\mathrm{KL}\left(p\right(x)\Vert q(x\left)\right)={\int }_{\mathcal{X}}p\left(x\right)\log \left\{p\left(x\right)/q\left(x\right)\right\}\mathrm{d}\mu \left(x\right)$. It was shown that ${\hat{\mathit{\theta }}}_{\mathrm{K}\mathrm{L}}$ is exactly the global MLE $\hat{\mathit{\theta }}$ if $\{{\mathbb{P}}_{\mathit{\theta }}:\mathit{\theta }\in \mathrm{\Theta }\}$ is a full exponential family (defined in their paper). This sheds light on the inference about generalized linear models (GLMs) based on exponential likelihood.

In many cases, some local machines might suffer from data of poor quality. This could lead to abnormal local estimates, which further degrade the statistical efficiency of the final estimator. To fix the problem, Minsker (2019) devised a robust assembling method. It leads to an estimator as ${\hat{\mathit{\theta }}}_{\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{u}\mathrm{s}\mathrm{t}}={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{\mathit{\theta }\in \mathrm{\Theta }}\sum _{k=1}^K\rho (|\mathit{\theta }-{\hat{\mathit{\theta }}}_k|)$, where $\rho (\cdot )$ is a robust loss function satisfying some conditions. For example, when $\rho \left(u\right)=u$ and $p=1$ (univariate case), ${\hat{\mathit{\theta }}}_{\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{u}\mathrm{s}\mathrm{t}}$ is the median of ${\hat{\mathit{\theta }}}_k$'s. It should be more robust against outliers compared with the simple average. Under some regularity conditions, they showed that ${\hat{\mathit{\theta }}}_{\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{u}\mathrm{s}\mathrm{t}}$ achieves the same convergence rate as the whole sample estimator provided $K\le O\left(\sqrt{N}\right)$.

2.2. Iterative approach

Although one-shot approach involves little communication cost, it suffers from several disadvantages. First, the local machines need to have sufficient amount of data (e.g., $n\gg \sqrt{N}$). Otherwise the aggregated estimator cannot reach the convergence rate as the global estimator. This prevents us from utilizing many machines to speed up the computation process (Jordan et al., 2019; Wang et al., 2017). Second, the simple averaging estimator is often poor in performance for nonlinear models (Huang & Huo, 2015; Jordan et al., 2019; Rosenblatt & Nadler, 2016). Last, when p is diverging with N, the situation could be even worse (Lee et al., 2017; Rosenblatt & Nadler, 2016). This suggests that carefully designed algorithms allowing a reasonable number of iterations should be useful for a distributed system.

Inspired by the one-step method in the M-estimator theory, Huang and Huo (2015) proposed an one-step refinement of the simple averaging estimator. Let us recall that $\overline{\mathit{\theta }}$ is the one-shot averaging estimator. To further improve its statistical efficiency, it should be broadcast to each local machine. Next, local gradient $\mathrm{\nabla }{\mathcal{L}}_k\left(\overline{\mathit{\theta }}\right)$ and local Hessian ${\mathrm{\nabla }}^2{\mathcal{L}}_k\left(\overline{\mathit{\theta }}\right)$ can be computed on each ${\mathcal{M}}_k$. Then, they are reported to ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$ to form the whole sample gradient $\mathrm{\nabla }\mathcal{L}\left(\overline{\mathit{\theta }}\right)=K^{-1}\sum _{k=1}^K\mathrm{\nabla }{\mathcal{L}}_k\left(\overline{\mathit{\theta }}\right)$ and Hessian ${\mathrm{\nabla }}^2\mathcal{L}\left(\overline{\mathit{\theta }}\right)=K^{-1}\sum _{k=1}^K{\mathrm{\nabla }}^2{\mathcal{L}}_k\left(\overline{\mathit{\theta }}\right)$. Thus an one-step updated estimator can be constructed on ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$ as (2) ${\hat{\mathit{\theta }}}^{\left(1\right)}=\overline{\mathit{\theta }}-\left[{\mathrm{\nabla }}^2\mathcal{L}\right(\overline{\mathit{\theta }}\left){]}^{-1}\mathrm{\nabla }\mathcal{L}\right(\overline{\mathit{\theta }}).$(2) Compared with one-shot estimator, ${\hat{\mathit{\theta }}}^{\left(1\right)}$ involves one more round of communication cost. Nevertheless, the statistical efficiency of the resulting estimator could be well improved. In fact, Huang and Huo (2015) showed that $\mathbb{E}\Vert {\hat{\mathit{\theta }}}^{\left(1\right)}-{\mathit{\theta }}^{\ast }{\Vert }_2^2\le \frac{C_1}{N}+O(\frac{1}{n^4}+\frac{1}{N^2}),$where $C_1>0$ is some constant. Obviously, this is a lower upper bound of mean squared error than that in (1(1) $\mathbb{E}\Vert \overline{\mathit{\theta }}-{\mathit{\theta }}^{\ast }{\Vert }_2^2\le \frac{C_1}{N}+\frac{C_2}{n^2}+O(\frac{1}{Nn}+\frac{1}{n^3}),$(1) ). To attain the global convergence rate, the local sample size needs to satisfy $n^{-4}=o\left(N^{-1}\right)$, which is a much milder condition. Furthermore, they showed that ${\hat{\mathit{\theta }}}^{\left(1\right)}$ also has the same asymptotic efficiency as the whole sample estimator $\hat{\mathit{\theta }}$ under some regularity conditions.

A natural idea to further extend the one-step estimator is to allow the iteration (2(2) ${\hat{\mathit{\theta }}}^{\left(1\right)}=\overline{\mathit{\theta }}-\left[{\mathrm{\nabla }}^2\mathcal{L}\right(\overline{\mathit{\theta }}\left){]}^{-1}\mathrm{\nabla }\mathcal{L}\right(\overline{\mathit{\theta }}).$(2) ) to be executed many times. Specifically, let ${\hat{\mathit{\theta }}}^{\left(t\right)}$ be the estimator of the tth iteration. Then, we can use (2(2) ${\hat{\mathit{\theta }}}^{\left(1\right)}=\overline{\mathit{\theta }}-\left[{\mathrm{\nabla }}^2\mathcal{L}\right(\overline{\mathit{\theta }}\left){]}^{-1}\mathrm{\nabla }\mathcal{L}\right(\overline{\mathit{\theta }}).$(2) ) by replacing $\overline{\mathit{\theta }}$ with ${\hat{\mathit{\theta }}}^{\left(t\right)}$ to generate the next step estimator ${\hat{\mathit{\theta }}}^{(t+1)}$ (see Figure 1(b)). However, this requires a large number of Hessian matrices to be computed and transferred. If the parameter dimension p is relatively high, this will lead to a significant communication cost of the order $O\left(K{p}^2\right)$. To fix the problem, Shamir et al. (2014) proposed an approximate Newton method, which conducts Newton-type iteration distributedly without transferring the Hessian matrices. Following this strategy, Jordan et al. (2019) developed an approximate likelihood approach. Their key idea is to update Hessian matrix on one single machine (e.g., ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$) only. Then, (2(2) ${\hat{\mathit{\theta }}}^{\left(1\right)}=\overline{\mathit{\theta }}-\left[{\mathrm{\nabla }}^2\mathcal{L}\right(\overline{\mathit{\theta }}\left){]}^{-1}\mathrm{\nabla }\mathcal{L}\right(\overline{\mathit{\theta }}).$(2) ) can be revised to be ${\hat{\mathit{\theta }}}^{(t+1)}={\hat{\mathit{\theta }}}^{\left(t\right)}-[{\mathrm{\nabla }}^2{\mathcal{L}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}\left({\hat{\mathit{\theta }}}^{\left(t\right)}\right){]}^{-1}\mathrm{\nabla }\mathcal{L}\left({\hat{\mathit{\theta }}}^{\left(t\right)}\right),$where ${\mathrm{\nabla }}^2{\mathcal{L}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$ is the Hessian matrix computed on the central machine. By doing so, the communication cost due to transmission of Hessian matrices can be saved. Under some conditions, they showed that (3) $\Vert {\hat{\mathit{\theta }}}^{(t+1)}-\hat{\mathit{\theta }}{\Vert }_2\le \frac{C_1}{\sqrt{n}}\Vert {\hat{\mathit{\theta }}}^{\left(t\right)}-\hat{\mathit{\theta }}{\Vert }_2,\mathrm{f}\mathrm{o}\mathrm{r}t\ge 0,$(3) holds with high probability, where $C_1>0$ is some constant. By the linear convergence formula (3(3) $\Vert {\hat{\mathit{\theta }}}^{(t+1)}-\hat{\mathit{\theta }}{\Vert }_2\le \frac{C_1}{\sqrt{n}}\Vert {\hat{\mathit{\theta }}}^{\left(t\right)}-\hat{\mathit{\theta }}{\Vert }_2,\mathrm{f}\mathrm{o}\mathrm{r}t\ge 0,$(3) ), we can see that it requires $[\log K/\log n]$ iterations to achieve the $\sqrt{N}$-consistency as the whole sample estimator $\hat{\mathit{\theta }}$, provided ${\hat{\mathit{\theta }}}^{\left(0\right)}$ is $\sqrt{n}$-consistent. Note that if $n=K=\sqrt{N}$, one iteration suffices to attain the optimal convergence rate. However, the satisfactory performance of this method relies on a good choice of the machine, on which the Hessian needs to be updated (Fan, Guo et al., 2019a). To fix the problems, Fan, Guo et al. (2019a) added an extra regularized term to the approximate likelihood used in Jordan et al. (2019). With this modification, the performance of the resulting estimator can be well improved. Theoretically, they showed a similar linear convergence rate under some more general conditions, which require no strict homogeneity of the local loss functions.

2.3. Shrinkage methods

We study the shrinkage methods for sparse estimation in this subsection. For a high-dimensional problem, especially when the dimension of ${\mathit{\theta }}^{\ast }$ is larger than the sample size N, it is difficult to estimate ${\mathit{\theta }}^{\ast }$ without any additional assumptions (Hastie et al., 2015). A popular constraint for tackling these problems is sparsity, which assumes only a subset of the entries in ${\theta }^{\ast }$ is non-zero. The index of non-zero entries is called the support of ${\mathit{\theta }}^{\ast }$, that is $\mathrm{supp}\left({\mathit{\theta }}^{\ast }\right)={\mathcal{A}}^{\ast }=\{1\le j\le p:{\theta }_j^{\ast }\ne 0\}.$To induce a sparse solution, an additional regularization term of $\mathit{\theta }$ is often introduced in the loss function. Specifically, we need to solve the shrinkage regression problem as $\underset{\mathit{\theta }\in \mathrm{\Theta }}{min}\left\{\mathcal{L}\right(\mathit{\theta })+\sum _{j=1}^p{\rho }_{\lambda }(|{\theta }_j|\left)\right\}$, where ${\rho }_{\lambda }(\cdot )$ is a penalty function with a regularization parameter $\lambda >0$. Popular choices are LASSO (Tibshirani, 1996), SCAD (Fan & Li, 2001) and others discussed in Zhang and Zhang (2012). For simplicity, we consider the LASSO estimator in the framework of the linear regression problem. Specifically, the whole sample estimator is computed as ${\hat{\mathit{\theta }}}_{\lambda }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{\mathit{\theta }\in \mathrm{\Theta }}\{\frac{1}{N}\Vert \mathit{Y}-\mathit{X}\mathit{\theta }{\Vert }_2^2+\lambda \Vert \mathit{\theta }{\Vert }_1\},$where $\mathit{Y}=(Y_1,\dots ,Y_N{)}^{\mathrm{\top }}\in {\mathbb{R}}^N$ is the vector of response, $\mathit{X}=(X_1,\dots ,X_N{)}^{\mathrm{\top }}\in {\mathbb{R}}^{N\times p}$ is the design matrix, and $\Vert \mathit{\theta }{\Vert }_1=\sum _{j=1}^p|{\theta }_j|$ denotes the $l_1$-norm of $\mathit{\theta }$. It is known that the LASSO procedure would produce biased estimators for the large coefficients. This is undesirable for the simple average procedures, since average cannot eliminate the systematic bias. To reduce bias, Javanmard and Montanari (2014) proposed a debiasing technique for the lasso estimator, that is (4) ${\hat{\mathit{\theta }}}_{\lambda }^{\left(\mathrm{d}\right)}={\hat{\mathit{\theta }}}_{\lambda }+\frac{1}{N}\mathit{M}{\mathit{X}}^{\mathrm{\top }}(\mathit{Y}-\mathit{X}{\hat{\mathit{\theta }}}_{\lambda }),$(4) where $\mathit{M}\in {\mathbb{R}}^{p\times p}$ is an approximation to the inverse of $\hat{\mathrm{\Sigma }}={\mathit{X}}^{\mathrm{\top }}\mathit{X}/N$. It appears that when $\hat{\mathrm{\Sigma }}$ is invertible (e.g., when $N\gg p$), setting $\mathit{M}={\hat{\mathrm{\Sigma }}}^{-1}$ gives ${\hat{\mathit{\theta }}}_{\lambda }^{\left(\mathrm{d}\right)}=({\mathit{X}}^{\mathrm{\top }}\mathit{X}{)}^{-1}{\mathit{X}}^{\mathrm{\top }}\mathit{Y}$, which is the ordinary least squares estimator and obviously unbiased. Hence, procedure (4(4) ${\hat{\mathit{\theta }}}_{\lambda }^{\left(\mathrm{d}\right)}={\hat{\mathit{\theta }}}_{\lambda }+\frac{1}{N}\mathit{M}{\mathit{X}}^{\mathrm{\top }}(\mathit{Y}-\mathit{X}{\hat{\mathit{\theta }}}_{\lambda }),$(4) ) compensates for the bias incurred by ${\ell }_1$ regularization in some sense.

By this debiasing technique, Lee et al. (2017) developed an one-shot type estimator for the LASSO problem. Specifically, let ${\hat{\mathit{\theta }}}_{k,\lambda }^{\left(\mathrm{d}\right)}$ be the debiased LASSO estimator computed on ${\mathcal{M}}_k$. Then an averaging estimator can be constructed on ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$ as ${\overline{\mathit{\theta }}}_{\lambda }=K^{-1}\sum _{k=1}^K{\hat{\mathit{\theta }}}_{k,\lambda }^{\left(\mathrm{d}\right)}$. Unfortunately, the sparsity level can be seriously degraded by averaging. For this reason, a hard threshold step often comes as a remedy. It was noticed that the debiasing step is computationally expensive. Hence an improved algorithm was also proposed to alleviate the computational cost of this step. Under certain conditions, they showed that the resulting estimator has the same convergence rate as the whole sample estimator. Battey et al. (2018) investigated the same problem with additional study on hypothesis testing. Furthermore, a refitted estimation procedure was used to preserve the global oracle property of the distributed estimator. An extension to high dimensional GLMs can also be found in Lee et al. (2017) and Battey et al. (2018). For this model, Chen and Xie (2014) implemented a majority voting method to aggregate the regularized local estimates. For the low dimensional sparse problem with smooth loss function (e.g., GLMs, Cox model), Zhu et al. (2019) developed a local quadratic approximation method with an adaptive-LASSO type penalty. They showed rigorously that the resulting estimator can be as good as the global oracle estimator.

Intuitively, above one-shot methods may need a stringent condition on the local sample size to meet the global convergence rate due to the limited communication. In fact, the simple averaging estimator requires $n\ge O(K{s}^2\log p)$ to match the oracle rate in the context of sparse linear model (Lee et al., 2017), where $s=|{\mathcal{A}}^{\ast }|$ is the number of non-zero entries of ${\mathit{\theta }}^{\ast }$. For this problem, Wang et al. (2017) and Jordan et al. (2019) independently proposed a communication-efficient iterative algorithm, which constructs a regularized likelihood by using local Hessian matrix. As demonstrated by Wang et al. (2017), an one-step estimator ${\hat{\mathit{\theta }}}^{\left(1\right)}$ suffices to achieve the global convergence rate if $n\ge O(K{s}^2\log p)$ (the condition used in  Lee et al., 2017). Furthermore, if multi-round communication is allowed, ${\hat{\mathit{\theta }}}^{(t+1)}$ (i.e., estimator of the $(t+1)$th iteration) can match the estimator based on the whole sample as long as $n\ge O(s^2\log p)$ and $t\ge O(\log K)$, under some certain conditions.

2.4. Non-smooth loss based models

The methods we described above typically require the loss function $\mathcal{L}$ to be sufficiently smooth, although a non-smooth regularization term is permitted (see e.g., Jordan et al., 2019; Wang et al., 2017; Zhang et al., 2013; Zhu et al., 2019). However, there are also some useful methods involving non-smooth loss functions, such as quantile regression and support vector machine. It is then of great interest to develop distributed methods for these methods.

We first focus on the quantile regression (QR) model. The QR model has a widespread use in statistics and econometrics, and performs more robustly against the outliers than the ordinary quadratic loss (Koenker, 2005). Specifically, a QR model assumes $Y_i=X_i^{\mathrm{\top }}{\mathit{\theta }}_{\tau }^{\ast }+{ϵ}_i,i\in \mathbb{S}$, where $X_i\in {\mathbb{R}}^p$ is the covariate vector, $Y_i$ is the corresponding response, ${\mathit{\theta }}_{\tau }^{\ast }\in {\mathbb{R}}^p$ is the true regression coefficient, and ${ϵ}_i$ is the random noise satisfying $\mathbb{P}({ϵ}_i\le 0|X_i)=\tau$, where $\tau \in (0,1)$ is a known quantile level. It is known that ${\mathit{\theta }}_{\tau }^{\ast }$ is the minimizer of $\mathbb{E}\left[{\rho }_{\tau }\right(Y_i-X_i^{\mathrm{\top }}\mathit{\theta }\left)\right]$. Here ${\rho }_{\tau }\left(u\right)=u(\tau -\mathbf{1}\{u\le 0\left\}\right)=u\left(\mathbf{1}\right\{u>0\}+\tau -1)$ is the non-differentiable check-loss function, where $\mathbf{1}\{\cdot \}$ is the indicator function. When data size N is moderate, we can estimate ${\mathit{\theta }}_{\tau }^{\ast }$ by ${\hat{\mathit{\theta }}}_{\tau }=\underset{\mathit{\theta }\in \mathrm{\Theta }}{min}{N}^{-1}\sum _{i\in \mathbb{S}}{\rho }_{\tau }(Y_i-X_i^{\mathrm{\top }}\mathit{\theta })$ on one single machine. However, when N is very large, a distributed system has to be used. Accordingly, distributed estimators have to be developed.

In this regard, Volgushev et al. (2019) studied the one-shot averaging type estimator. Specifically, a local estimator ${\hat{\mathit{\theta }}}_{k,\tau }$ is first computed on each local machine ${\mathcal{M}}_k$. Then, the averaging estimator is assembled as ${\overline{\mathit{\theta }}}_{\tau }=K^{-1}\sum _{k=1}^K{\hat{\mathit{\theta }}}_{k,\tau }$ on the central machine ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$. They further investigated the theoretical properties of the averaging estimator in detail. It was shown that the if the number of machines satisfies $K=o(\sqrt{N}/\log N)$, then ${\overline{\mathit{\theta }}}_{\tau }$ should be as efficient as the whole sample estimator ${\hat{\mathit{\theta }}}_{\tau }$ under some regularity conditions. Chen and Zhou (2020) proposed an estimating equation based one-shot approach for the QR problem. The asymptotic equivalence between the resulting estimator and the whole sample estimator was also established under $K=o\left(N^{1/4}\right)$ and some other conditions. It can be seen that the performance of one-shot approaches relies more on the local sample size. In fact, Volgushev et al. (2019) showed that $K=o\left(\sqrt{N}\right)$ is a necessary condition for the global efficiency of the simple averaging estimator ${\overline{\mathit{\theta }}}_{\tau }$. To remove the constraint $K=o\left(\sqrt{N}\right)$ on the number of machines, Chen et al. (2019) proposed an iterative approach. Their key idea is to approximate the check-loss function by a smooth alternative. More specifically, they approximated $\mathbf{1}\{u>0\}$ by a smooth function $H\left(u/h\right)$, where $H(\cdot )$ is a smooth cumulative distribution function and $h>0$ is the tuning parameter controlling the approximation goodness (see Figure 2(a)). With this modification, the algorithm can update the estimates by (5) ${\hat{\mathit{\theta }}}_{\tau }^{(t+1)}=[\mathit{V}\left({\hat{\mathit{\theta }}}_{\tau }^{\left(t\right)}\right){]}^{-1}\mathit{U}\left({\hat{\mathit{\theta }}}_{\tau }^{\left(t\right)}\right),$(5) where $\mathit{U}\left(\mathit{\theta }\right)=\sum _{k=1}^K{\mathit{U}}_k\left(\mathit{\theta }\right)$, $\mathit{V}\left(\mathit{\theta }\right)=\sum _{k=1}^K{\mathit{V}}_k\left(\mathit{\theta }\right)$, and ${\mathit{U}}_k\in {\mathbb{R}}^p,{\mathit{V}}_k\in {\mathbb{R}}^{p\times p}$ depend only on the bandwidth h and local sample ${\mathcal{S}}_k$. It was shown that a constant number of rounds of iteration suffices to match the convergence rate of the whole sample estimator. Thus, the communication cost is roughly of the order $O\left(K{p}^2\right)$. But this is not applicable when p is very large. Thus, for the high dimensional QR problem, Zhao et al. (2014) and Zhao et al. (2019) adopted an one-shot averaging method based on the debiased local estimates as that in (4(4) ${\hat{\mathit{\theta }}}_{\lambda }^{\left(\mathrm{d}\right)}={\hat{\mathit{\theta }}}_{\lambda }+\frac{1}{N}\mathit{M}{\mathit{X}}^{\mathrm{\top }}(\mathit{Y}-\mathit{X}{\hat{\mathit{\theta }}}_{\lambda }),$(4) ). Accordingly, Chen et al. (2020) proposed a communication-efficient multi-round algorithm inspired by the approximate Newton method (Shamir et al., 2014). This iterative approach removes the restriction on the number of machines. A revised divide-and-conquer stochastic gradient descent method for QR and other models with diverging dimension can be found in Chen, Liu et al. (2021b).

Figure 2. Approximation of two non-smooth loss functions. (a) QR loss with $\tau =0.6$ and (b) hinge loss.

We next consider the support vector machine (SVM), which is one of the most successful statistical learning methods (Vapnik, 2013). The classical SVM is aimed at the binary classification problem, i.e., the response variable $Y_i\in \{-1,1\}$. Formally, a standard linear SVM solves the problem $\underset{\mathit{\theta }\in \mathrm{\Theta }}{min}{N}^{-1}\sum _{i\in \mathbb{S}}(1-Y_i{X}_i^{\mathrm{\top }}\mathit{\theta }{)}_{+}+\lambda \Vert \mathit{\theta }{\Vert }_2^2$, where $(u{)}_{+}=u\mathbf{1}(u>0)$ is the hinge loss, and $\lambda >0$ is the regularization parameter. By the same smooth technique used in Chen et al. (2019), i.e., replacing the hinge loss with a smooth alternative (see Figure 2(b)), Wang, Yang et al. (2019) proposed an iterative algorithm like (5(5) ${\hat{\mathit{\theta }}}_{\tau }^{(t+1)}=[\mathit{V}\left({\hat{\mathit{\theta }}}_{\tau }^{\left(t\right)}\right){]}^{-1}\mathit{U}\left({\hat{\mathit{\theta }}}_{\tau }^{\left(t\right)}\right),$(5) ). To reduce the communication cost incurred by transferring matrices, they further employed the approximate Newton method (Shamir et al., 2014). Theoretically, they showed the asymptotic normality of the estimator, which can be used to construct confidence interval. For the ultra-high dimensional SVM problem, Lian and Fan (2018) studied the one-shot averaging method with debiasing procedure similar to (4(4) ${\hat{\mathit{\theta }}}_{\lambda }^{\left(\mathrm{d}\right)}={\hat{\mathit{\theta }}}_{\lambda }+\frac{1}{N}\mathit{M}{\mathit{X}}^{\mathrm{\top }}(\mathit{Y}-\mathit{X}{\hat{\mathit{\theta }}}_{\lambda }),$(4) ).

3. Nonparametric models

Different from parametric models, a nonparametric model typically involves infinite-dimensional parameters. In this section, we focus mainly on the nonparametric regression problems. Specifically, consider here a general regression model as $Y_i=f^{\ast }\left(X_i\right)+{ϵ}_i,i\in \mathbb{S}$, where $f^{\ast }(\cdot )$ is an unknown but sufficiently smooth function and ${ϵ}_i$ is the random noise with zero mean. The aim of nonparametric regression is to estimate function $f^{\ast }\in \mathcal{F}$, where $\mathcal{F}$ is a given nonparametric class of functions.

3.1. Local smoothing

One way to estimate $f^{\ast }(\cdot )$ is to fit a locally constant model by kernel smoothing (Fan & Gijbels, 1996). More concretely, the whole sample estimator is given by ${\hat{f}}_h\left(x\right)=\sum _{i\in \mathbb{S}}{W}_{h,X_i}\left(x\right)Y_i,$where the $W_{h,X_i}\left(x\right)\ge 0$ is the local weight at $X=x$ satisfying $\sum _{i\in \mathbb{S}}{W}_{h,X_i}\left(x\right)=1$. Specifically, for a Nadaraya–Watson kernel estimator, we should have $W_{h,X_i}\left(x\right)=K\left(\right(X_i-x\left)/h\right)/\sum _{i^{\prime }\in \mathbb{S}}K\left(\right(X_{i^{\prime }}-x\left)/h\right)$, where $K(\cdot )$ is a kernel function and $h>0$ is the bandwidth. In the univariate case ($p=1$), classical theory stated that the mean squared error of ${\hat{f}}_h\left(x\right)$ is of the order $O(h^4+(Nh{)}^{-1})$ (Fan & Gijbels, 1996). Thus, the optimal rate $O\left(N^{-4/5}\right)$ can be achieved by choosing bandwidth $h=O\left(N^{-1/5}\right)$.

For a distributed kernel smoothing, an one-shot estimator can also be constructed. Let ${\hat{f}}_{k,h}\left(x\right)$ be the local estimator computed on ${\mathcal{M}}_k$. Then an averaging estimator can be obtained as ${\overline{f}}_h\left(x\right)=K^{-1}\sum _{k=1}^K{\hat{f}}_{k,h}\left(x\right)$. Chang, Lin, Wang (2017) studied the theoretical properties of ${\overline{f}}_h\left(x\right)$ in a specific function space $\mathcal{F}$. They established the same minimax convergence rate of ${\overline{f}}_h\left(x\right)$ as that of the whole sample estimator. However, they found that a strict restriction on the number of machines K is needed to achieve this optimal rate. To fix the problem, two solutions were provided. They are, respectively, a date-dependent bandwidth selection algorithm and an algorithm with a qualication step.

Nearest neighbours method can be regarded as another local smoothing method. Qiao et al. (2019) studied the Nearest neighbours classification in a distributed setting, where the optimal number of neighbours to achieve the optimal rate of convergence was derived. Li et al. (2013) discussed the problem of density estimation for scattered datasets. Kaplan (2019) focused on the choice of bandwidth for nonparametric smoothing techniques. All the works above in this subsection indicate that the bandwidth (or local smoothing parameter) used in the distributed setting should be adjusted according to the whole sample size N, other than the local sample size n.

3.2. RKHS methods

We next discuss another popular nonparametric regression method. This is reproducing kernel Hilbert space (RKHS) method. An RKHS $\mathcal{H}$ can be induced by a continuous, symmetric and positive semi-definite kernel function $\mathcal{K}(\cdot ,\cdot ):{\mathbb{R}}^p\times {\mathbb{R}}^p\mapsto \mathbb{R}$. Two typical examples are: the polynomial kernel $\mathcal{K}(x_1,x_2)=(x_1^{\mathrm{\top }}{x}_2+1{)}^d$ with an integer $d\ge 1$, and the radical kernel $\mathcal{K}(x_1,x_2)=\exp (-\gamma \Vert x_1-x_2{\Vert }_2^2)$ with $\gamma >0$. Refer to, for example, Berlinet and Thomas-Agnan (2011); Wahba (1990) for more details about RKHS. Then, our target is to find an $\hat{f}\in \mathcal{H}$ so that the following penalized empirical loss can be minimized. That leads to the whole sample estimator as (6) ${\hat{f}}_{\lambda }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{f\in \mathcal{H}}\left\{\frac{1}{N}\sum _{i\in \mathbb{S}}\right(Y_i-f\left(X_i\right){)}^2+\lambda \Vert f{\Vert }_{\mathcal{H}}^2\},$(6) where $\Vert \cdot {\Vert }_{\mathcal{H}}$ is the norm associated with the RKHS $\mathcal{H}$ and $\lambda >0$ is the regularization parameter. This problem is also known as kernel ridge regression (KRR). By the representer theorem for the RKHS (Wahba, 1990), any solution to the problem (6(6) ${\hat{f}}_{\lambda }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{f\in \mathcal{H}}\left\{\frac{1}{N}\sum _{i\in \mathbb{S}}\right(Y_i-f\left(X_i\right){)}^2+\lambda \Vert f{\Vert }_{\mathcal{H}}^2\},$(6) ) must have the linear form as ${\hat{f}}_{\lambda }\left(x\right)=\sum _{i\in \mathbb{S}}{\alpha }_i\mathcal{K}(X_i,x)$, where ${\alpha }_i\in \mathbb{R}$ for each $i\in \mathbb{S}$. By this property, we can treat the KRR as a parametric problem with unknown parameter $\alpha =({\alpha }_1,\dots ,{\alpha }_N{)}^{\mathrm{\top }}\in {\mathbb{R}}^N$. The error bounds of the whole sample estimator ${\hat{f}}_{\lambda }$ has been well established in the existing literature (e.g., Steinwart et al., 2009; Zhang, 2005). However, a standard implementation of the KRR involves inverting a kernel matrix in ${\mathbb{R}}^{N\times N}$ (Saunders et al., 1998). Therefore, when N is extremely large, it is time consuming or even computationally infeasible to process the whole sample on a single machine. Thus, we should consider a distributed system.

In this regard, Zhang et al. (2015) studied the distributed KRR by taking the one-shot averaging approach. Specifically, each machine ${\mathcal{M}}_k$ computes local KRR estimate ${\hat{f}}_{k,\lambda }$ by (6(6) ${\hat{f}}_{\lambda }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{f\in \mathcal{H}}\left\{\frac{1}{N}\sum _{i\in \mathbb{S}}\right(Y_i-f\left(X_i\right){)}^2+\lambda \Vert f{\Vert }_{\mathcal{H}}^2\},$(6) ) based on local sample ${\mathcal{S}}_k$. Then the central machine ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$ averages them to obtain final estimator ${\overline{f}}_{\lambda }=K^{-1}\sum _{k=1}^K{\hat{f}}_{k,\lambda }$. Theoretically, they established the optimal convergence rate of mean squared error for ${\overline{f}}_{\lambda }$ with different types of kernel functions, under some regularity conditions. Lin et al. (2017) derived a similar optimal error bound under some relaxed conditions. Xu et al. (2016) extended the loss function in (6(6) ${\hat{f}}_{\lambda }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{f\in \mathcal{H}}\left\{\frac{1}{N}\sum _{i\in \mathbb{S}}\right(Y_i-f\left(X_i\right){)}^2+\lambda \Vert f{\Vert }_{\mathcal{H}}^2\},$(6) ) to a further general form. Some related works on the distributed KRR problem by one-shot averaging approach can be found in Shang and Cheng (2017), Lin and Zhou (2018), Guo et al. (2019); Mücke and Blanchard (2018) and Wang (2019) and many others. It was noted that these one-shot approaches require the number of machines diverges in a relative slow speed to meet the global convergence rate. To fix the problem, Chang, Lin, Zhou (2017) proposed a semi-supervised learning framework by utilizing the additional unlabelled data (i.e., observations without response $Y_i$). Latest work of Lin et al. (2020) allowed communication between machines to improve the performance. In order to choose an optimal tuning parameter λ in (6(6) ${\hat{f}}_{\lambda }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{f\in \mathcal{H}}\left\{\frac{1}{N}\sum _{i\in \mathbb{S}}\right(Y_i-f\left(X_i\right){)}^2+\lambda \Vert f{\Vert }_{\mathcal{H}}^2\},$(6) ), Xu et al. (2018) proposed a distributed generalized cross-validation method.

For semiparametric models, Zhao et al. (2016) considered a partially linear model with heterogeneous data in a distributed setting. Specifically, they assumed the following model (7) $Y_i=X_i^{\mathrm{\top }}{\mathit{\theta }}_{\left(k\right)}^{\ast }+f^{\ast }\left(W_i\right)+{ϵ}_i,i\in {\mathcal{S}}_k,$(7) where $W_i\in \mathbb{R}$ is an additional covariate, $f^{\ast }(\cdot )$ is the unknown function, and ${\mathit{\theta }}_{\left(k\right)}^{\ast }\in {\mathbb{R}}^p$ is the true linear coefficient associated with the data on ${\mathcal{M}}_k$ for $1\le k\le K$. In other words, the local data on different machines are assumed to share the same nonparametric part, but are allowed to have different linear coefficients. To estimate the unknown function and coefficients, they extended the classical RKHS theory to cope with the partially linear function space. Under some regularity conditions, the resulting estimator of the nonparametric part is shown to be as efficient as the whole sample estimator, provided the number of machines does not grow too fast. The case with high dimensional linear part was also investigated. For example, under the homogeneity assumption (i.e., the linear coefficients ${\mathit{\theta }}_{\left(k\right)}^{\ast }$'s are assumed to be identical to ${\mathit{\theta }}^{\ast }$ across different machines), Lv and Lian (2017) adopted the one-shot averaging approach with debiasing technique analogous to (4(4) ${\hat{\mathit{\theta }}}_{\lambda }^{\left(\mathrm{d}\right)}={\hat{\mathit{\theta }}}_{\lambda }+\frac{1}{N}\mathit{M}{\mathit{X}}^{\mathrm{\top }}(\mathit{Y}-\mathit{X}{\hat{\mathit{\theta }}}_{\lambda }),$(4) ) to estimate the linear coefficient. Lian et al. (2019) considered the same heterogeneous model as in (7(7) $Y_i=X_i^{\mathrm{\top }}{\mathit{\theta }}_{\left(k\right)}^{\ast }+f^{\ast }\left(W_i\right)+{ϵ}_i,i\in {\mathcal{S}}_k,$(7) ), but the linear part is assumed in a high dimensional setting (i.e., $p>N$). For this model, they proposed a novel projection approach to estimate the common nonparametric part (not in an RKHS framework). Theoretically, the asymptotic normality of the one-shot averaging estimator for the nonparametric function was established under some certain conditions.

4. Other related works

4.1. Principal component analysis

Principal component analysis (PCA) is a common procedure to reduce the dimension of the data. It is widely used in the practical data analysis. Unlike the regression problems, PCA is an unsupervised method, which does not require a response variable Y. To conduct a PCA, a covariance matrix $\hat{\mathrm{\Sigma }}$ needs to be constructed as $\hat{\mathrm{\Sigma }}=N^{-1}\sum _{i\in \mathbb{S}}{X}_i{X}_i^{\mathrm{\top }}$, where $X_i$'s are assumed to be centralized already. Next, a standard singular value decomposition (SVD) is applied to $\hat{\mathrm{\Sigma }}$. That leads to $\hat{\mathrm{\Sigma }}=\hat{V}\hat{D}{\hat{V}}^{\mathrm{\top }}$, where $\hat{D}$ is a diagonal matrix of eigenvalues and $\hat{V}$ is an orthogonal matrix of eigenvectors. Then, the columns of $\hat{V}$ are the principal component directions that we need.

In a distributed setting, simple average of the eigenvectors estimated locally cannot give a valid result. To solve the problem, Fan, Wang et al. (2019b) developed a divide-and-conquer algorithm for estimating eigenspaces. It involves only one single round of communication. This algorithm is quite easy to implement as well. We state it as follows (Fan, Wang et al., 2019b, Algorithm 1):

1. For each $k=1,\dots ,K$, machine ${\mathcal{M}}_k$ computes d leading eigenvectors of the local sample covariance matrix ${\hat{\mathrm{\Sigma }}}_k=n^{-1}\sum _{i\in {\mathcal{S}}_k}{X}_i{X}_i^{\mathrm{\top }}$, denoted by ${\hat{v}}_{1,k},\dots ,{\hat{v}}_{d,k}\in {\mathbb{R}}^p$. Next, they are arranged by columns in ${\hat{V}}_k=({\hat{v}}_{1,k},\dots ,{\hat{v}}_{d,k})\in {\mathbb{R}}^{p\times d}$, which is then sent to the central machine ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$.

2. The central machine ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$ averages K local projection matrices to obtain $\tilde{\mathrm{\Sigma }}=K^{-1}\sum _{k=1}^K{\hat{V}}_k{\hat{V}}_k^{\mathrm{\top }}$. Then it computes d leading eigenvectors of $\tilde{\mathrm{\Sigma }}$, denoted by ${\tilde{v}}_1,\dots ,{\tilde{v}}_d\in {\mathbb{R}}^p$. ${\tilde{v}}_1,\dots ,{\tilde{v}}_d$ are the estimators of the first d principal component directions that we need.

It is noticeable that the communication cost of above one-shot algorithm is of the order $O\left(Kdp\right)$. This can be considered to be communication-efficient since d is usually much smaller than p in practice. Fan, Wang et al. (2019b) showed that, under some appropriate conditions, the distributed estimator achieves the same convergence rate as the global estimator. The cases with heterogeneous local data were also investigated in their work. To further remove the restriction on the number of machines, Chen, Lee et al. (2021) proposed a communication-efficient multi-round algorithm based on the approximate Newton method (Shamir et al., 2014).

4.2. Feature screening

Massive datasets often involve ultrahigh dimensional data, for which feature screening is critically important (Fan & Lv, 2008). To fix the idea, consider a standard linear regression model as $Y_i=X_i^{\mathrm{\top }}{\mathit{\theta }}^{\ast }+{ϵ}_i,i\in \mathbb{S}$, where $X_i\in {\mathbb{R}}^p$ is the covariate vector, $Y_i$ is the corresponding response, ${\mathit{\theta }}^{\ast }\in {\mathbb{R}}^p$ is the true parameter, and ${ϵ}_i$ is the random noise. To screen for the most promising features, the seminal method of sure independence screening (SIS) has been proposed by Fan and Lv (2008). Specifically, let ${\mathcal{A}}^{\ast }=\{1\le j\le p:{\theta }_j^{\ast }\ne 0\}$ be the true sparse model. Let ${\omega }_j$ be the Pearson correlation between jth feature and response Y. Then, SIS screens features by a hard threshold procedure as ${\hat{\mathcal{A}}}_{\gamma }=\{1\le j\le p:|{\hat{\omega }}_j|>\gamma \}$, where γ is a prespecified threshold and ${\hat{\omega }}_j$ is the whole sample estimator of ${\omega }_j$. Under some specific conditions, Fan and Lv (2008) showed the sure screening property for SIS, that is, $\mathbb{P}\left({\mathcal{A}}^{\ast }\subset {\hat{\mathcal{A}}}_{\gamma }\right)\to 1\mathrm{a}\mathrm{s}N\to \mathrm{\infty }.$However, the estimator ${\hat{\omega }}_j$ is usually biased for many correlation measures. This indicates that a direct one-shot averaging approach is unlikely to be the best practice for the distributed system. To fix the problem, Li et al. (2020) proposed a novel debiasing technique. They found that many correlation measures can be expressed as ${\omega }_j=g({\nu }_1,\dots ,{\nu }_s)$, including Pearson correlation used above, Kendall τ rank correlation, SIRS correlation (Zhu et al., 2011), etc. Therefore, they used U-statistics to estimate the components ${\nu }_q,1\le q\le s$ on local machines. Then, these unbiased estimators of ${\nu }_q$'s given by local machines are averaged on the central machine ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$. Consequently, ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$ can construct distributed estimator ${\tilde{\omega }}_j$ by the averaging estimators of the components in the known function g. Finally, they showed the sure screening property of ${\hat{\mathcal{A}}}_{\gamma }$ based on the distributed estimators under some regularity conditions. When the feature dimension is much larger than the sample size (i.e., $p\gg N$), another distributed computing strategy is to partition the whole data by features, other than by samples. Refer to, for example, Song and Liang (2015); Yang et al. (2016) for more details.

4.3. Bootstrap

Bootstrap and related resampling techniques provide a general and easily implemented procedure for automatically statistical inference. However, these methods are usually computationally expensive. Especially when sample size N is very large, it would be even practically infeasible to conduct. To mitigate this computing issue, various alternative methods have been proposed, such as subsamping approach (Politis et al., 1999) and ‘m-out-of-n’ bootstrap (Bickel et al., 2012). Their key idea is to reduce the resample size. However, due to the difference between the size of whole sample and resample, an additional correction step is generally required to rescale the result. This makes these methods less automatic.

To solve this problem, Kleiner et al. (2014) proposed the bag of little bootstraps (BLB) method. It integrates the idea of subsampling and can be computed distributedly without a correction step. Suppose that N sample units have been randomly and evenly partitioned to K machines. Consider that we want to assess the accuracy of the point estimator for some parameter $\mathit{\theta }$. Then we summarize their algorithm as follows.

1. For each $1\le k\le K$, machine ${\mathcal{M}}_k$ draws r samples of size N (instead of n) from ${\mathcal{S}}_k$ with replacement. Then it computes r estimates of $\mathit{\theta }$ based on the r resamples drawn above, respectively. After that, each ${\mathcal{M}}_k$ computes some accuracy measure (e.g., variance, confidence region) by the r estimates above, denoted by ${\hat{\xi }}_k$. Finally, all of the local machines send ${\hat{\xi }}_k$'s to the central machine ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$.

2. The central machine ${\mathcal{M}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$ aggregates these ${\hat{\xi }}_k$'s by $\overline{\xi }=K^{-1}\sum _{k=1}^K{\hat{\xi }}_k$. And $\overline{\xi }$ is the final accuracy measure that we need.

It is remarkable that one does not need to process datasets of size N on local machines actually, although the nominal size of resample is N. This is because each machine contains at most n sample units. In fact, randomly generating some certain weight vectors of length n suffices to approximate the resampling process.

5. Future study

To conclude the article, we would like to discuss here a number of interesting topics for future study. First, for datasets with massive sizes, a distributed system is definitely needed. Obviously, there could be no place to store the data. On the other hand, for datasets with sufficiently small sizes, traditional memory based statistical methods can be immediately used. Then, there leaves a big gap between the big and small datasets. Those middle-sized data are often of sizes much larger than the computer memory but smaller than the hard drive. Consequently, they can be comfortably placed on a personal computer, but can hardly be processed by memory as a whole. For those datasets, their sizes are not large enough to justify an expensive distributed system. They are also not small enough to be handled by traditional statistical methods. How to analyse datasets of this size seems to be a topic worth studying. Second, when the whole data are allocated to local machines randomly and evenly, the data on different machines are independent and identically distributed and balanced. Then, all of the methods discussed above can be safely implemented. However, when the data on different machines are collected from (for example) different regions, the homogeneity of the local data would normally be hard to satisfy. The situation could be even worse if the sample sizes allocated to different local machines are very different. How to cope with these heterogeneous and unbalanced local data is a problem of great importance (Wang et al., 2020). The idea of meta analysis may be applicable to these situations (Liu et al., 2015; Xu & Shao, 2020; Zhou & Song, 2017). Finally, in the era of big data, personal privacy is under unprecedented threat. How to protect users' private information during the learning process deserves urgent attention. In this regard, differential privacy (DP) provides a theoretical approach for privacy-preserving data analysis (Dwork, 2008). Some related works associated with distributed learning are Agarwal et al. (2018), Truex et al. (2019) and Wang, Ishii et al. (2019) and many others. Although it is a hot research area recently, there are still many open challenges. Thus, it is of great interest to study the privacy-preserving distributed statistical learning problems practically and theoretically.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work is supported by National Natural Science Foundation of China (No. 11971171), the 111 Project (B14019) and Project of National Social Science Fund of China (15BTJ027). Weidong Liu's research is supported by National Program on Key Basic Research Project (973 Program, 2018AAA0100704), National Natural Science Foundation of China (No. 11825104, 11690013), Youth Talent Support Program, and a grant from Australian Research Council. Hansheng Wang's research is partially supported by National Natural Science Foundation of China (No. 11831008, 11525101, 71532001). It is also supported in part by China's National Key Research Special Program (No. 2016YFC0207704).

Notes on contributors

Yuan Gao

Mr. Yuan Gao is a Ph.D. candidate in school of statistics at East China Normal University.

Weidong Liu

Dr. Weidong Liu is the Distinguished Professor in school of mathematical sciences at Shanghai Jiao Tong University.

Hansheng Wang

Dr. Hansheng Wang is a professor in Guanghua School of Management at Peking University.

Xiaozhou Wang

Dr. Xiaozhou Wang is an assistant professor in school of statistics at East China Normal University.

Yibo Yan

Mr. Yibo Yan is a Ph.D. candidate in school of statistics at East China Normal University.

Riquan Zhang

Dr. Riquan Zhang is a professor in school of statistics at East China Normal University.