{'fraction_non_alphanumeric': 0.04318652702009162, 'fraction_numerical': 0.0144652055731909, 'mean_word_length': 4.262052410482096, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Communication-efficient distributed statistical learning

We present the Communication-efficient Surrogate Likelihood (CSL) framework for solving distributed statistical learning problems. CSL provides a communication-efficient surrogate to the global likelihood that can be used for low-dimensional estimation, high-dimensional regularized estimation and Bayesian inference. For low-dimensional estimation, CSL provably improves upon the averaging schemes and facilitates the construction of confidence intervals. For high-dimensional regularized estimation, CSL leads to a minimax optimal estimator with minimal communication cost. For Bayesian inference, CSL can be used to form a communication-efficient quasi-posterior distribution that converges to the true posterior. This quasi-posterior procedure significantly improves the computational efficiency of MCMC algorithms even in a non-distributed setting. The methods are illustrated through empirical studies.


INTRODUCTION
In modern statistical problems, people may not be able to store large datasets on a single machine, let alone conduct statistical inference via optimization or MCMC. In this big data paradigm, it is natural to consider methods that split the dataset across multiple machines and conduct statistical inference in a distributed manner, in either the frequentist framework of obtaining a pointestimator Duchi et al. [2012], Zhang et al. [2013a], Zhang and Lin [2015], Kannan et al. [2014], Shamir et al. [2014], Lee et al. [2015], or in the Bayesian framework of sampling from the posterior distribution Wang and Dunson [2015], Neiswanger et al. [2015]. Good distributed statistical learning procedures should give rise to statistically optimal estimators with low communication cost.
To achieve low communication cost, one-shot or embarrassingly parallel approaches Zhang et al. [2013a], Lee et al. [2015], Wang and Dunson [2015], Neiswanger et al. [2015] only use one-round of communication by first sending estimators or posterior samples from local machines to a center node, and then combining them to form a global estimator or approximation to the posterior distribution.
In the frequentist framework, most one-shot approaches rely on averaging Zhang et al. [2013a], where the global estimator is the average of the local estimators. Lee et al. [2015] extends this idea from low dimensional regular models Zhang et al. [2013a] to high-dimensional sparse linear regression by combining local debiased Lasso estimates van de Geer et al. [2014]. Recent work by Duchi et al. [2015] show that under certain conditions, these averaging estimators can attain the information-theoretic complexity lower bound-for linear regression, at least O(dk) bits must be communicated in order to attain the minimax rate of parameter estimation, in spite of the existence of low-dimensional structures such as sparsity Braverman et al. [2015], where d is the dimension of the parameter and k is the number of machine. However, these averaging-based one-shot communication approaches tend to suffer from several drawbacks. First, it is not straightforward to conduct statistical inferences, such as creating confidence intervals/regions and doing hypothesis testing, based on the averaging estimator. Second, in order for the averaging estimator to achieve the minimax rate of convergence, each local machine must have access to at least Ω( √ N ) samples, where N is the total sample size. In another word, the number of machines should be much smaller than √ N . This could be restrictive when both N and k are large. Third, when the underlying data generating model is non-linear, our empirical study shows that even for small k of order 10 1 , the averaging estimator only exhibits a slight improvement over the local estimators and tends to have low estimation accuracy.
In the Bayesian framework, embarrassingly parallel approaches run a Markov chain Monte Carlo (MCMC) algorithm across local machines in parallel and transmit local posterior samples to a central node to produce an approximation to the global posterior distribution. Unfortunately, the number of local posterior samples must grow at least exponentially in the dimension d to approximate a local d-dimensional posterior distributions due to the curse of dimensionality. As a consequence, the communication complexity can be prohibitive when d is moderate to high.
When combining local posterior samples in the central node, existing approaches Wang and Dunson [2015], Neiswanger et al. [2015] that approximate the global posterior distribution by a weighted empirical distribution of "averaging draws" tend to suffer from the weight-degeneracy issue (weights collapse to only a few sample) when k is large.
In this paper, we formulate a unified framework for frequentist and Bayesian distributed learning by introducing a Communication-efficient Surrogate likelihood (CSL) function. In the frequentist perspective, CSL can be viewed as a communication-efficient surrogate to the global likelihood function that uses all samples. For example, CSL can replace the role of the likelihood function in forming either the maximum likelihood estimator (MLE) in regular parametric models or the penalized MLE in high-dimensional models. In the Bayesian perspective, CSL can be used to form a quasi-posterior distribution Chernozhukov and Hong [2003] as a surrogate to the full posterior.
CSL can be constructed efficiently by communicating O(dk) number of bits. After its construction, CSL can be efficiently evaluated by using the n samples in a single local machine. Even in a non-distributed Bayesian computation environment, CSL can be used as computationally-efficient surrogate to the likelihood function by pre-dividing the dataset into k subsamples-the computational complexity of one iteration of MCMC is then reduced by a factor of k.
Our CSL-based distributed learning approach overcomes the aforementioned drawbacks suffered by the one-shot and embarrassingly parallel approaches. In the frequentist framework, CSL can be used in the same way as the global likelihood function to either form a point-estimator or conduct statistical inferences. Moreover, a CSL-based estimator can achieve the same rate of convergence as the global likelihood-based estimator and has communication complexity O(dk).
On the computational side, we show that by applying our distributed estimation procedure in an iterative manner, the resulting multi-round algorithm has a geometric convergence rate with contraction factor O(n −1/2 ), where n is the number of samples in each local machine. This O(n −1/2 ) rate of convergence significantly improves the condition number contraction factor 1 attained by the naive distributed gradient method that forms the global gradient in each iteration by combining the local gradients. As an implication, in order to achieve the same accuracy as the global likelihoodbased estimator, n can be independent of the total sample size N as long as O( log N log n ) iterations are applied, which is constant for n > k. In contrast, the aforementioned averaging estimator requires n ≫ √ N . For instance, due to the fast O(n −1/2 ) rate, usually 2 to 3 iterations suffices for our procedure to match the same accuracy of the global likelihood-based estimator even for reasonably large k (See Section 4.1 for more details). Unlike the bootstrap-based approach Zhang et al. [2013a] for boosting the accuracy that has a high computational complexity, the additional complexity of the iterative version of our approach grows only linearly in the number of iterations. Last but not least, our empirical study suggests that a CSL-based estimator may exhibit significant improvement over the averaging estimator in non-linear distributed statistical learning problems, making CSL practically appealing for large distributed learning tasks.
For high-dimensional ℓ 1 -regularized estimator, we provide an algorithm that communicates O(dk) bits that attains the optimal communication/risk trade off [Garg et al., 2014]. This improves over the averaging method of Lee et al. [2015] because it requires p times less computation, and allows for iterative refinement to obtain arbitrarily low optimization error in logarithmic number of rounds. During the preparation of this manuscript, we became aware of the concurrent work of Wang et al. [2016] that also proposes the high-dimensional CSL-based estimator. The theoretical results for ℓ 1 -regularized least squares are equivalent. However, their focus is limited to the highdimensional linear models, and do not consider parametric estimation in low-dimensioanl families, nor do they consider Bayesian posterior estimation.
In the Bayesian framework, our method does not require transmitting local posterior samples and is free from the weight degeneracy issue. This makes the communication complexity of our approach considerably lower than those embarrassingly parallel Bayesian computation approaches.
In the optimization literature, the most related work to ours is the distributed approximate Newton algorithm (DANE) proposed in Shamir et al. [2014]. Although they only rigorously analysed DANE for quadratic objectives, their idea of combining gradient descent with a local Newton method also applies to non-quadratic objectives. Motivated by this optimization idea, we propose the CSL framework that is specifically tailored for distributed statistical learning. CSL is a surrogate likelihood function that can be easily applied to M -estimation, high-dimensional regularized estimation and Bayesian inference. On the theoretical side, the analysis in Shamir et al. [2014] does not imply that DANE can improve over the gradient method for non-quadratic objectives. In contrast, our analysis demonstrates fast convergence rates for a broader class of regular parametric models and high-dimensional models and is not restricted to quadratic objectives. Another related work to DANE and ours is the iterated Hessian sketch (IHS) algorithm Pilanci and Wainwright [2014] that is specifically designed for constrained least-squares optimization. Although DANE applied to quadratic problems can be viewed as a special case of IHS by choosing the sketching matrix as a rescaled sub-sampling matrix, the analysis in Pilanci and Wainwright [2014] only applies to a class of low-incoherence sketching matrices that excludes sub-sampling.
Overall, we build a general framework for statistical inference in the distributed learning setting.
Our iterative distributed statistical learning procedure reveals a trade-off between communication complexity and statistical accuracy. When applied to Bayesian inference, our approach has the same flavour as Pereyra [2015], Amandine et al. [2015] that introduces modern optimization techniques to improve the computational efficiency of sampling-based algorithms.
The remainder of this paper is organized as follows. In Section 2, we informally present the motivation for CSL. Section 3 presents algorithms and theory on three different problems: parameter estimation in low-dimensional regular parametric models (Section 3.1), regularized parameter estimation in the high-dimensional problems (Section 3.2), and Bayesian inference in regular parametric models (Section 3.3). Section 4 shows our experimental results in these three settings. All proofs are provided in the Appendix.

BACKGROUND AND PROBLEM FORMULATION
We begin by setting up our general framework for distributed statistical learning: either in a frequentist or Bayesian framework. Hereafter, we describe the CSL methodology and explain the motivation.
Suppose that the data is stored in a distributed manner where each machine stores a sub-sample of n observations. Let Z j : = {Z ij : i = 1, . . . , n} denote the sub-sample that is stored in the jth machine M j for j = 1, . . . , k. Our goal is to conduct statistical inference on the parameter θ while taking into consideration the communication cost among the machines. For example, we may want to find a point estimator θ and an associated confidence interval (region).
Let L : Θ × Z → R be a twice differentiable loss function such that the true parameter is a minimizer of the population risk L * (θ) : = E θ * [L(θ; Z)], that is Define the local and global loss functions as Here L j (θ) is the loss function evaluated at θ by using the local data stored in machine M j .
Throughout the rest of the paper, the reader may consider the negative log-likelihood function as a typical candidate of loss function L.
For large-scale problems, suppose that we estimate parameter θ directly using the gradient the methodology, and rigorous analysis is provided in Section 3 to follow.
Our motivation is from the Taylor series expansion of L N . Viewing L N (θ) as an analytic function, we can expand it into the infinite series Here θ is any initial estimator of θ, for example, the local empirical loss minimizer arg min θ L 1 (θ) in the first machine M 1 . Because the data is split across machines, evaluating the derivatives ∇ j L N (θ) (j ≥ 1) requires one communication round. However, unlike the d-dim gradient vector ∇L N (θ), the higher-order derivatives require communicating more than O(d 2 ) bits from each machine. This reasoning motivates us to replace the global higher-order derivatives ∇ j L N (θ) (j ≥ 2) with the local ones, leading to the following modification of L N , Therefore, we define the surrogate loss to be (ignoring some additive constants) Comparing expressions (4) and (5), we see that the approximation error is (up to a constant) where we used the fact that ∇ 2 L N (θ) − ∇ 2 L 1 (θ) 2 = O n −1/2 by using matrix concentration, is a good estimator of θ that achieves the root-N rate of convergence under mild conditions. One may also construct confidence regions associated with θ using the sandwiched covariance matrix (for example, see formula (10)). In our distributed learning framework, we propose the following estimator by pretending the surrogate loss function L as the global loss function L N (θ), θ = arg min θ∈Θ L(θ).
In Section 3.1, we show that θ is equivalent to θ up to higher-order terms, and provide two ways to construct confidence regions for θ using local observations stored in machine M 1 . Now we give a heuristic understanding of why θ is a good estimator. For convenience, we may assume that the empirical risk function L N (θ) has a unique minimizer. First look at the global empirical loss minimizer θ. Under our assumption that the loss function is twice-differentiable, θ is the unique solution of equation 2 By solving this equation, we obtain θ − θ * 2 = O p ( ∇L N (θ * ) 2 ) = O p (N −1/2 ), as long as the Hessian matrix ∇ 2 L N (θ * ) is non-singular. Now let is turn to the surrogate loss minimizer θ.
2 There is no Taylor's theorem for vector-valued functions, but we formalize this heuristic in Section 3.1 Similar argument leads to θ − θ * 2 = O p ( ∇ L(θ * ) 2 ) and we only need to show that ∇ L(θ * ) 2 is of order O p (N −1/2 ). In fact, by our construction, which is of order O p (N −1/2 ) as long as θ * − θ 2 = O p k −1/2 with k = N/n the number of machines. For example, this requirement on initial estimator θ is satisfied for the minimizer θ 1 of the sub-sample loss function L 1 (θ) when n > k.
Example (High-dimensional regularized estimator): In the high-dimensional regime where the dimensionality d can be much larger than the sample size N , we need to impose some lowdimensional structural assumption, such as the sparsity assumption, on the unknown true parameter θ * . Under such an assumption, regularized estimators are popular and widely used for estimating θ. For concreteness, we focus on the sparsity assumption that most components of the d-dim vector θ * is zero, and consider the ℓ 1 -regularized estimator θ : = arg min θ∈Θ L N (θ) + λ θ 1 as the benchmark estimator that we want to approximate. Here λ is the regularizing parameter.
In the distributed learning framework, we consider the following estimator using surrogate loss function L, In Section 3.2, we show that θ achieves the same rate of convergence as the benchmark estimator θ under a set of mild conditions. This idea of using the surrogate loss function to approximate the global loss function is general and is applicable to other high-dimensional problems.
Example (Bayesian inference): In Bayesian framework, viewing parameter θ as random, we place a prior distribution π over parameter space Θ. Without confusion, we also use the same notation π(θ) to denote the pdf of the prior distribution at point θ. According to Bayes' rule, the posterior distribution satisfies The loss function L corresponds to the negative log-likelihood function of the statistical model {P θ : θ ∈ Θ} and L N (θ) is the global negative log-likelihood of the observations Z N 1 . The posterior distribution π(θ | Z N 1 ) can be used to conduct statistical inference. For example, we may construct an estimator of θ as the posterior expectation and use the highest posterior region as a credible region for this estimator. Since the additive constant C in expression (7) can be absorbed into the normalizing constant, we may use the surrogate posterior distribution to approximate the global posterior distribution π(θ | Z N 1 ). In Section 3.3, we formalize this argument and show that this surrogate posterior gives a good approximation the global posterior.
From now on, we will refer the methodology of using the surrogate loss function L(·) to approximate the global loss function L(·) for distributed statistical learning as a Communication-efficient Surrogate Likelihood (CSL) method. As a remark, the idea of computing the global likelihoodfunction using sub-samples is useful not only in the distributed learning framework, but also in a common learning framework if the sample size is so large that even a single evaluation of the likelihood function or its gradient can be expensive. Using our surrogate loss function L(θ), we only need one pass over the entire dataset to construct L(θ). After its construction, L(θ) can be efficiently evaluated by using a small subset of the data.

MAIN RESULTS AND THEIR CONSEQUENCES
In this section, we elaborate on the three examples in Section 2.2 of applying the CSL method.
For each of the examples, we provide explicit bound on either the estimation error θ − θ * 2 of the resulting estimator θ or the approximation error π N − π N 1 of the approximated posterior π N (·).

Communication-efficient M -estimators in low-dimensions
In this subsection, we consider low-dimensional parametric family {P θ : θ ∈ Θ}, where the dimensionality d of θ is much smaller than sample size n. Under this setting, the minimizer of the population risk in optimization problem (1) is unique under the set of regularity conditions to follow and θ * is identifiable. As a concrete example, we may consider the negative log-likelihood function ℓ(θ; z) = − log p(z; θ) as the loss function, where p(·; θ) is the pdf of the observation under P θ .
Note that the developments in this subsection can also be extended to misspecified models where the marginal distribution P of the observations is not contained in the model space Under misspecification, we can view the parameter θ * associated with the projection P θ * of the true data generating model P onto the misspecified model space {P θ : θ ∈ Θ} as the "true" parameter.
The results under misspecification are similar to the well-specified case and omitted in this paper.
For low-dimensional parametric models, we impose some regularity conditions on the parameter space, the loss function L and the associated population risk function L * . These conditions are standard in classical statistical analysis of M -estimators. In the rest of the paper, we call a parametric model that satisfies this set of regularity conditions a regular parametric model. Our first assumption describes the relationship of the parameter space Θ and the true parameter θ * .
The second assumption is a local identifiability condition, ensuring θ * to be a local minimum of L * .
Assumption PB (Local convexity): The Hessian matrix I(θ) = ∇ 2 L * (θ) of the population risk function L * (θ) is invertible at θ * : there exists two positive constants (µ − , µ + ), such that Since the loss function is the negative log-likelihood function, the corresponding Hessian matrix is called the information matrix. Our next assumption is a global identifiability condition, which is a standard condition for proving estimation consistency.
Assumption PC (Identifiability): For any δ > 0, there exists ǫ > 0, such that Our last assumption controls moments of higher order derivatives of the loss function, and allows us to obtain high probability bounds on the estimation error using Markov's inequality. Let Assumption PD (Smoothness): There exist constants (G, L) and a function M (z) such that Following the heuristic argument in Section 2.2, we use the surrogate function L defined in (6) as the objective function for constructing an M-estimator in regular parametric models. Our first result shows that under Assumptions PA-PD, given any reasonably good initial estimator θ, any significantly boosts the accuracy in terms of the approximation error θ− θ 2 to the global empirical risk minimizer θ = arg min θ∈Θ L N (θ).
Theorem 3.1. Suppose that Assumptions PA-PD hold and the initial estimator θ lies in the neighbour U (ρ) of θ * . Then any minimizer θ of surrogate loss function L(θ) satisfies with probability at least 1 − C 1 kn −8 . Here constants (C 1 , C 2 ) are independent of (k, n, N ).
Under the conditions of Theorem 3.1, it can be shown that θ − θ * 2 = O p (N −1/2 ) and .1 and inequality (A.7) in Appendix B.1), and therefore as long as θ − θ 2 = O p (n −1/2 ), which is true for θ = θ 1 : = arg min θ L 1 (θ), the empirical risk minimizer in local machine M 1 . To formalize this argument, we have the following corollary that provides an ℓ 2 risk bound for θ.

Corollary 3.2. Under conditions of Theorem 3.1, we have
and C is some constant independent of (n, k, N ).
Note that the Hájek-Le Cam minimax theorem guarantees that for any estimator θ N based on Therefore, the estimator θ is (first-order) minimax-optimal and achieves the Cramér-Rao lower bound when the loss function L is the negative log-likelihood function.
One-step approximation: The computational complexity of exactly minimizing the surrogate loss L(θ) in (8) can be further reduced by using a local quadratic approximation to L. In fact, we have by Taylor's theorem that Similar to before, we replace the global gradient ∇L N (θ) with the local one ∇L 1 (θ), which leads to the quadratic surrogate loss, Because surrogate loss functions L H and L agree up to the second-order Taylor expansion, they behave similarly when used as objective functions for constructing M -estimators. This motivates the estimator, which can be computed in closed-form using matrix inversion. The next theorem shows that θ H satisfies the same estimation bound as θ. Unlike the classical one-step MLE that requires the initial estimator to be within an O(N −1/2 ) neighbourhood of the truth θ * , we only require θ − θ * 2 to be O(n −1/2 ).
Theorem 3.3. Suppose that Assumptions PA-PD hold and the initial estimator θ satisfies θ − The analogue of Corollary 3.2 can also be stated for θ H .
Iterative local estimation algorithm: It is worth mentioning that Theorem 3.1 (Theorem 3.3) suggests that an Iterative Local Estimation Algorithm (abbreviated as ILEA, see Algorithm 1) reduces the approximation error θ − θ 2 by a factor of n −1/2 in each iteration as long as the initial in each iteration of ILEA, we set θ as the current iterate θ (t) , construct the surrogate loss function L (t) (θ), and then solve the next iterate θ (t+1) by either exactly minimizing the surrogate loss as or by the local one-step quadratic approximation as Theorem 3.1 (or Theorem 3.3) guarantees, with high probability, the error bound where C 3 is positive constant independent of (n, k, N ). If the desired accuracy is the statistical accuracy θ−θ * 2 of the MLE and our initial estimator is O(n −1/2 ) consistent, then we need proceed at most ⌈ log k log n ⌉ iterations. ILEA interpolates between the gradient method and Newton's algorithm. When n is large relative to k, then ILEA behaves like Newton's algorithm, and we achieve the optimal statistical accuracy in one iteration. If n is a fixed constant size, then ILEA degenerates to a preconditioned gradient method. By appropriately choosing the sub-sample size n, ILEA achieves a trade-off among storage, communication, and computational complexities, depending on specific constraints of computing resources.
Do one of the following in Machine M 1 :

11
(1). Update θ (t+1) ∈ arg min θ∈Θ L t (θ) ; // Exactly minimizing surrogate function For example, when L corresponds to the negative log-likelihood function, Σ = I(θ * ) −1 will be the inverse of the information matrix. It is easy to see that the plug-in estimator is a consistent estimator of the asymptotic covariance matrix Σ, that is, Σ → Σ in probability as N → ∞. Based on the limiting distribution of √ N ( θ − θ * ) and the plug-in estimator Σ, we can conduct statistical inference, for example, constructing confidence intervals for θ.
The following corollary shows that for any reasonably good initial estimator θ, the asymptotic distribution of either the minimizer θ of the surrogate function L(θ) or the local one-step quadratic approximated estimator θ matches that of the global empirical risk minimizer θ. Moreover, we also have a consistent estimator Σ of Σ only using the local information in Machine M 1 . Therefore, we can conduct statistical inference locally without access to the entire data while achieving the same power as global statistical inference procedures.
Moreover, the following plug-in estimator is a consistent estimator for Σ as n → ∞. If we also have k → ∞, then the plug-in estimator is also a consistent estimator for Σ as (n, k) → ∞. Similar results hold for the local one-step quadratic approximated estimator θ H under the assumptions of Theorem 3.3.
Corollary 3.4 illustrates that we may pretend L(θ) as the global loss function and use it for statistical inference-Σ is precisely the plug-in estimator of the sandwiched covariance matrix using surrogate loss function L(θ) (cf. equation (11)). In the special case when L(θ) is the negative loglikelihood function, we may instead use ∇ 2 L( θ) −1 as our plug-in estimator for Σ . Σ ′ tends to be a better estimator than Σ when k ≫ n, since variance O(k −1 ) of the middle term in equation (13) is much smaller than variance O(n −1 ) of the middle term in equation (12). See Section 4.1 for an empirical comparison of using Σ and Σ ′ for constructing confidence intervals.
3.2 Communication-efficient regularized estimators with ℓ 1 -regularizer In this subsection, we consider high-dimensional estimation problems where the dimensionality d of parameter θ can be much larger than the sample size n. Although the development here can apply to a broader class of problems, we focus on ℓ 1 -regularized procedures. ℓ 1 -regularized estimators work well under the sparsity assumption that most components of the true parameter θ * is zero. Let S = supp(θ * ) be a subset of {1, . . . , d} that encodes the sparsity pattern of θ * and s = |S| = d j=1 I(θ * j = 0). Applying the same idea of using the surrogate loss function L(θ) as a proxy to the global likelihood function in ℓ 1 -regularized estimation procedures, we propose the following communication-efficient regularized estimator θ ∈ arg min θ∈Θ L(θ) + λ θ 1 .
We study the statistical precision of this estimator in the high-dimensional regime.
We first present a theorem on the statistical error bound θ − θ * 2 of estimator θ for general loss function L. Then we will verify the conditions in the theorem for high-dimensional linear models and generalized linear models. To begin with, we state our assumptions.
Assumption HA (Restricted strongly convexity): Local loss function L 1 (θ) at machine L 1 is restricted strongly convex over S: for all δ ∈ C(S) : = {v : v S 1 ≤ 3 v S c 1 }, where δ is some positive constant independent of n.
As the name suggested, restricted strongly convexity requires the global loss function L n (θ) to be a strongly convex function when restricted to the cone C(S).
Assumption HB (Restricted Lipschitz Hessian): Both local and global loss function L 1 (θ) and L N (θ) have restricted Lipschitz Hessian at radius R: for all δ ∈ C(S) ∩ B R (θ * ), where M is some positive constant independent of N .
The restricted Lipschitz Hessian condition is always satisfied for linear models where the Hessian Theorem 3.5. Suppose that Assumption HA and Assumption HB at radius The lower bound condition on the regularization parameter λ for θ is slightly stronger than that for the estimator θ based on the global loss function, which is λ ≥ 2 ∇L N (θ * ) ∞ . Since the estimation error upper bound provided by Theorem 3.5 is proportional to the regularization parameter, it is reasonable to expect that θ will have a slight larger error than θ, depending on how good the initial estimator θ is. For example, in generalized linear models, if smallest regularization parameters λ are chosen for θ and θ, then the estimation error of θ will be greater than that of θ by an amount of As long as θ − θ * 1 and θ − θ * 2 are sufficiently small, this difference will be negligible with respect to the estimation error bound of θ, which is s log d N . For example, we may choose θ to be the local ℓ 1 regularized estimator θ 1 : = arg min θ L 1 (θ) + λ 1 θ with estimation error s log d n , so that We may also consider an iterative estimation procedure analogous to Algorithm 1 for improving the higher-order estimation accuracy of the communication-efficient regularized estimator θ. The convergence rate can be analyzed by inducting on Theorem 3.5. Now we apply Theorem 3.5 to two examples.

Example (Sparse linear regression):
In sparse linear regression, observations {z ij = (x ij , y ij ) : where x ij is a d-dimensional covariate vector, y ij is the response and β ∈ R d is the unknown regression coefficient to be estimated. Recall the sparsity assumption that s = d j=1 I(θ * j = 0) = o(n). For linear regression, the global loss function takes the form as We consider random design where x ij is i.i.d. A-sub-Gaussian, that is, for all α ∈ R d , be the covariance matrix of the design. For this class of design, it is known that Assumption HA is satisfied with high probability as long as Σ is strictly positive definite and n ≥ C 0 s log d for some constant C 0 > 0 depending on the minimal eigenvalue of Σ Raskutti et al.
[2010]. For linear model, the Lipschitz constant M in Assumption HB is zero and therefore HB is also satisfied.
Theorem 3.6. If x ij is A-sub-Gaussian, Σ is strictly positive definite and n ≥ C 0 s log d, then with probability at least 1 − c 1 exp{−c 2 n}, it holds that If the initial estimator satisfies θ − θ * 1 ≤ C 2 s log d n , then with the same probability, it holds that The constants (c 1 , c 2 , C 0 , C 1 , C 2 , C 3 ) are independent of (n, k, d, s).
For sparse linear regression under the sparsity condition, the minimax rate of estimating θ is s log d s N . Therefore, Theorem 3.6 shows that our approximated estimator θ is nearly minimaxoptimal if n ≥ Cs √ N log d for some constant C > 0. When this lower bound on the local sample size n fails, we may still apply the iterative estimation procedure (Algorithm 1) to boost the estimation accuracy and obtain a minimax-optimal estimator as we remarked after Theorem 3.5.
Example (Generalized linear models): In this section, we apply Theorem 3.5 to generalized linear models with ℓ 1 -regularizer. We begin with some background on generalized linear models.
Recall that the data z ij = (x ij , y ij ), where y ij is the response and x ij is the d-dim covariate vector.
A generalized linear model assumes the conditional distribution of y ij given x ij to be where σ is a scalar parameter, θ is the unknown d-dim parameter to be estimated and φ is a link function, for example, e.g. φ(x) = log(1 + e x ) in Logistic regression, and φ(x) = e x in Poisson regression. We still assume sparsity that s = d j=1 I(θ * j = 0) = o(n). Now the global loss function and its gradient is given by Under a random design assumption, we verify Assumptions HA and HB, and obtain the following result.

Communication-efficient Bayesian inference
In this subsection, we consider a Bayesian framework of distributed statistical learning for regular parametric models. In a standard Bayesian framework, one specifies a prior distribution π over the parameter space Θ and then conducts statistical inference by sampling from the global posterior where L : Θ × Z → R is the negative log-likelihood function and D is the normalizing constant. In the rest of this subsection, we tacitly assume that the loss function L introduced in Section 2.1 is the negative log-likelihood function. Extensions to the Gibbs posterior Bissiri Pier and Stephen [2013] where L is replaced with a generic loss function L in posterior (14)  [2015] in distributed Bayesian learning framework utilizes the decomposition so that the global posterior π(θ | Z N 1 ) can be written as the product of sub-sample posteriors where the prior is raised to power k −1 so that it is appropriately weighted in product (15) and D j is the normalizing constant. In a typical MapReduce framework, separate Markov chains are run in machines {M j } k j=1 based on the local data in an embarrassingly parallel manner so that the draws can be collected from each sub-sample posterior. After running individual Markov chains, all local posterior draws are transmitted to a central node. In the central node, an approximation π N (θ) to the global posterior π N (θ) : = π(θ | Z N 1 ) is formed. A main drawback of these approaches is that the communication cost can be extremely high-for example, exponentially large in the dimension d-since the number of draws from each local posterior must be large enough to be representative of the local posterior distribution.
Following the heuristic in Section 2.2, we propose a different sampling scheme for distributed Bayesian inference using the surrogate loss function L(θ). Our sampling scheme is communication efficient and requires running one single Markov chain in a local machine. Here is an outline of the algorithm: 1. Compute a good initial estimate θ, e.g. the one-step estimate θ H in Section 3.1.
3. Transmit all local gradients to Machine M 1 and form the global gradient ∇L N (θ) = 1 k k j=1 ∇L j (θ).
5. Machine M 1 runs a Markov chain to sample from the surrogate posterior π N (θ) ∝ exp − N L(θ) π(θ), and uses the draws to conduct statistical inference.
The following result shows that the surrogate posterior π N (·) is close to the global posterior π(· | Z N 1 ) as long as the initial estimator θ is reasonably close to θ * .
If we use the local one-step estimator θ H as the initial estimator θ, then the approximation error becomes This illustrates that we may choose k = N/n up to o N 1/2 (log N ) −1 why still maintaining π N − π N 1 = o p (1). The overall communication of this procedure requires two passes over the entire dataset (one for computing θ H and one for constructing L(θ)). To allow larger k, we may apply the iterative algorithm in Section 3.1 to improve the accuracy of the initial estimator θ. As our theory only covers regular parameter models in low dimensions, it is still an open problem of how to design communication-efficient Bayesian procedures for high-dimensional problems with provably theoretical guarantees.

SIMULATIONS
In this section, we conduct simulations for the three examples of using CSL methodology in Section 2.2.

Distributed M -estimation in logistic regression
In logistic regression, i.i.d. observations Z N 1 = {Z ij = (X ij , Y ij ) : i = 1, . . . , n. j = 1, . . . , k} are generated from the model In our simulation, the true regression coefficient θ * is a d-dim vector with d ∈ {2, 10, 50} and the d-dim covariate vector X ij is independently generated from N (0, I d ). For each replicate of the simulation, we uniformly sample the parameter θ * from the d-dim unit cube [0, 1] d .
We implement the one-step Newton-Raphson estimator θ (1) with the averaging estimator θ A Zhang et al.
[2013b] as our initial estimator θ. We also implement the iterative local estimation algorithm to produce 2-step and 3-step estimators θ (2) and θ (3) by iteratively applying the one-step estimation procedure. We compare the our communication estimators with the (optimal) global M -estimator θ global and the sub-sample estimator θ sub that only uses the local data in Machine M 1 . Two different regimes are considered: 1. the total sample size N is fixed at N = 2 19 ≈ 10 6 , and the local sample size n varies from 10 2 to 10 4 ; 2. the local sample size n is fixed at 64 (d = 2), 256 (d = 10) or 2048 (d = 50), and the number of machines k varies from 10 2 to 10 4 .        have the best performance. When n is sufficiently large, the 1-step, 2-step and 3-step estimators have almost the same performance as θ global . However, as n becomes small, further application of the iterative local estimation procedure in Algorithm 1 does not improve the statistical accuracy. This is in fact consistent with Theorem 3.3-the contraction coefficient is dominated by the sum of two terms: the initial estimation error θ (t) − θ global 2 and the local Hessian approximation error |||∇L 1 (θ * ) − ∇L N (θ * )||| 2 . Even though the initial estimation error can be reduced to a small level, the local Hessian approximation error still persists for small n and prevents further improvement from applying the iterative procedure. We remark that the condition that the local size n should exceed a d-dependent threshold is a mild requirement in practice.
Indeed, the local machine storage limit in reality is often large enough to ensure n ≫ d. Even under the scenario (small n) where our theory fails to predict, the 1-step, 2-step and 3-step estimators still have better performance than θ A and θ sub . In plots (b), (d) and (e), we fix the local sample size n under different d such that n exceeds the d-dependent threshold, and gradually increase the number of machines k. In our regime, k is comparable or even much larger than n, and therefore the averaging estimator θ A does not improve as more data is available. This is consistent with the theoretical results in Zhang et al. [2013b] that requires k ≫ n for θ A to have comparable performance as θ global . By using our approach, even a single step application of Algorithm 1 significantly improves the accuracy of θ A . Moreover, θ (2) and θ (3) achieve almost the same accuracy as θ global .
Consistent with our theory, for a fixed number of steps t, the t-step estimate θ (t) tends to have larger estimation error than θ global as k grows. In plot (d), even for k as large as 10 4 (much larger than the local sample size n ∼ 10 2 ), the 2-step estimate θ (2) already achieves the same level of estimation accuracy as the global estimator θ global . These simulation results illustrate the potential of our communication-efficient estimation procedures for large-scale distributed statistical learning. Now we assess the performance of the inference procedures based on the plug-in estimators Σ and Σ ′ under the logistic model (16). We use Σ or Σ ′ and the 3-step estimator θ (3) to construct a 95% confidence interval (CI) for the first component θ 1 of θ as The coverage of the CI based on 100 trials is calculated. Figure 2 shows the result. In plot (a), coverage based on both plug-in estimators are low at n = 2 7 because the sample size is so small  In plots (b), we change the number of machines k while fixing the local sample size n (total sample size N = nk).
that the center θ (3) of the CI has a large bias (see Figure 1 (c)). In plot (b), the CI based on Σ ′ has low coverage when the number k of machines is small, which is consistent with our theory. In all other regimes of (n, k), both CI have coverage that is close to the nominal level 95%. Moreover, the CI based on Σ ′ is slightly better than the one based on Σ for large k, which empirically justifies our intuition in the discussion right after Corollary 3.4.

Distributed sparse linear regression
We evaluate the CSL estimator on the sparse linear regression problem. The data is generated as N (0, 1), the noise ǫ ij is i.i.d. N (0, 1), and θ * is s-sparse with signal-to-noise ratio |θ i | σ = 5. In the first experiment, we keep the total data size N fixed, and increase the number of machines k. This corresponds to each machine having a smaller local sample size n as k increases. We observe that the one-step CSL estimator has nearly constant error, even though each machine has less local data. In fact at k = 30, the local data size is n = 720, which is much smaller than d, yet the CSL estimator achieves the same mean-square error as lasso on all N points. The error of the averaging estimator increases dramatically as n decreases, since the mean-squared error is s log d n , showing that the averaging algorithm is not suitable in the distributed setting.
In the second experiment, we keep n fixed and increase k and N . As predicted by our theory, the one-step CSL estimator has error that is linear on the log-log scale because the mean-squared error scales as s log d nk . The averaging estimator only has error that slowly decreases with the increased sample size, due to the bias induced by regularization. The averaging estimator does not attain mean-square error of s log d nk .

Distributed Bayesian inference
Our synthetic dataset is generated from the logistic model (16) under dimension d ∈ {2, 10, 50}. We use the 3-step estimator θ (3) in Section 3.1 as the initial estimator θ and implement the Bayesian procedures based on the (approximated) posterior distribution π n (θ) and π N (θ) by sampling using Markov Chain Monte Carlo algorithm. We use the Metropolis algorithm, where at each iteration the proposal distribution for θ is a d-dim Gaussian distribution centered at the current iterate θ (t) .
In each case, we run the Markov chain for 20000 iterations and treat the first half as burn-in.     θ under different (d, n, k) combinations (n is chosen so that θ (3) is a good approximation to the global estimator θ, see Figure 1). Consistent with our theoretical prediction, π N (θ) provides a good approximation to π N (θ) as long as the initial estimator θ is sufficiently close to θ, even when k is much larger than n (see plot (b)). Since the computation of the approximate poster distribution π N (θ) only uses the local data in Machine M 1 , the computation of the acceptance ratio using π N (θ) is k times as fast as that using the full data posterior π N (θ) in each iteration of the Metropolis algorithm.

REAL DATA APPLICATION
We apply distributed logistic regression to the skin dataset Rajen and Abhinav [2012]. The goal is to predict whether a given color sample described by its B, G, R values (each ranges from 0 − 255) corresponds to a skin sample or non-skin sample. The dataset is generated using skin textures from face images of diversity of age, gender, and race people. The total sample size is 245, 057, out of which 50, 859 is the skin samples and 194, 198 is non-skin samples. The skin dataset contains three features-B, G, R values of the color, and a 0-1 response variable indicates whether the sample is non-skin (0) or skin (1). We randomly split the dataset into a training set of size N = 200, 000 and a testing set N 0 = 45057, and use B-spline transforms (df= 15) for each feature as predictors to allow a non-linear dependence between the response and features. Therefore, the dimension of the covariate X is d = 45.
We randomly split the entire training set into k 0 = 100 subsets, each of size n = 2000.
We apply our distributed M -estimation method for logistic regression to a training set with k ∈ {20, 40, 60, 80, 100} subsets, and test the fitted model to the testing set. Now, we exactly minimize the surrogate function to form the 1-step estimator θ (1) with the averaging estimator θ A Zhang et al.
[2013b] as our initial estimator. We also implement the iterative local estimation algorithm to produce 2-step and 3-step estimators θ (2) and θ (3) by iteratively applying the onestep estimation procedure. Figure 6 plots the misclassification rate versus the number of subsets used. As we can see, the 1-step estimator gains significant improvement on the prediction performance over the initial averaging estimator, and both the 2-step and 3-step estimators have similar prediction performance as the 1-step estimator. This suggests that for the skin dataset and our split setting, the one-step approximation of the likelihood function already has reasonably good performance.

CONCLUDING REMARKS
In this paper, we presented the Communication-efficient Surrogate Likelihood (CSL) framework for solving distributed statistical learning problems. We applied this methodology to three examples: low-dimensional M -estimation, high-dimensional regularized estimation and low-dimensional Bayesian inference. We believe that the general idea of constructing a surrogate function to the negative log-likelihood function (or general loss function) that can be efficient evaluated using only a sub-sample of the observations is also useful in big-data problems where the sample size is so large that even the calculation of the likelihood function is expensive. One future direction is to investigate how this framework can be applied to high-dimensional Bayesian inference with strong theoretical supports. Another direction is to find a sharp theoretical lower bound on the local sample size n for the final estimator to remain optimal, for example, in the minimax sense.

APPENDIX A. PROOFS OF MAIN RESULTS
A.1 Proof of Theorem 3.1 For j = 1, . . . , k, let M j = 1 n n i=1 M (z ij ) and δ ρ = min{ρ, ρµ − /4M }. Consider the following "good events": , and Before proving the claimed error bound for θ, we state two auxiliary results that are used in the proof. The first result provides a control on the probability of bad event k j=0 E c j , which is proved Here c j (j = 1, 2, 3) are constants independent of (n, k, N, d, G, L).
The second result characterizes the error bound θ− θ 2 in terms of the gradient norm ∇ L( θ) 2 at θ, which formalizes the heuristic argument in Section 2.2. Its proof is provided in Appendix B.2.
Lemma A.2. Suppose that Assumptions PA-PD hold. Then under event E 0 ∩ E 1 we have Therefore, it remains to prove a high probability upper bound of the gradient norm ∇ L( θ) 2 .
Simple calculation yields By the optimality of the global empirical risk minimizer θ, we have ∇L N ( θ) = 0.
By adding and subtracting ∇ N ( θ) in equation (A.1), we obtain By the integral form of Taylor's expansion, we have that for any j ∈ {1, . . . , k}, under event E j . Combining the three preceding displays, we obtain that under event k j=0 E k , Combining Lemma A.1 and the above display yields the claimed error bound on θ − θ 2 .
A.2 Proof of Corollary 3.2 Recall definitions of events {E j } k j=0 in Section A.1. In the remaining of this proof, we use C to denote some constant independent of (n, k, N ), whose magnitude may change from line to line. We need the following auxiliary result, whose proof is provided in Appendix B.3.
Combining this Lemma, inequality (A.7) in Appendix B.1 and Theorem 3.1, we obtain that Now by applying Hölder's inequality and Lemma B.1 in Appendix B.1, we obtain Combining this with bound (A.8) in Appendix B.1 on P( k j=0 E c j ), we obtain that under Assumption PA, A.3 Proof of Theorem 3.3 Before analysing the one-step Newton-Raphson estimator θ H , we prove some auxiliary results.
Use λ min (A) to denote the minimal eigenvalue of a symmetric matrix A.
The error can be decomposed as We analyze the two terms respectively. The first term can be expressed as which yields the bound The second term can be analysed by using Theorem 5.3 in Bubeck [2014], which guarantees that under the assumption θ − θ 2 ≤ µ N 2M N , it holds that where µ N : = λ min [∇ 2 L N ( θ)] and M N is the Lipschitz constant of the Hessian ∇ 2 L N (θ), that is ∇ 2 L N (θ 1 ) − ∇ 2 L N (θ 2 ) 2 ≤ M N θ 1 − θ 2 2 for all θ 1 , θ 2 ∈ U (ρ). Putting pieces together, we Now the claimed bounds on θ H − θ 2 is a direct consequence of the preceding display and Lemma A.4.
A.4 Proof of Corollary 3.4 The proof of the second part on the consistency of plug-in estimators for Σ is standard by using the consistency of θ implied by the first part, the central limit theorem and Slutsky's theorem.
Therefore we only prove the first part on the asymptotic expansion of θ. Based on Theorem 3.1, we only need to establish the asymptotic expansion (10) of the global empirical risk minimizer θ.
By the integral form of Taylor's expansion, we have Then simple linear algebra yields . Then, the claimed expansion is an easy consequence of inequality (A.7) and Assumption D.

A.5 Proofs for regularized M-estimators
Proof of Theorem 3.5. This theorem follows from applying Corollary 1 of Negahban et al. [2012] to the objective F (θ). We check that L(θ) satisfies the restricted strong convexity condition.
The restricted strong convexity of L is implied by the same property of L 1 , since We can upper bound ∇ L(θ * ) ∞ as follows: Using Assumption HB, Proof of Theorem 3.6. To apply Theorem 3.5, we have to compute By applying the sub-exponential concentration inequality , Pr( 1 where c Σ is a constant that depends on Σ. By union bound over all (j, k) pairs, .
Proof of Theorem 3.7. To apply Theorem 3.5, we need to compute verify Assumptions HA and HB.
The restricted strong convexity of L is verified in Proposition 1 of Negahban et al. [2012]. Next we verify Assumption HB. Thus, Thus M = L φ BL, where L φ is a local upper bound on φ ′′′ , L is the upper restricted eigenvalue of X, and B = max x ∞ .
We also need to compute an upper bound on where we used the same argument as in the proof of Theorem 3.6.
By Lemma 6 of Negahban et al.
[2012], we know ∇L N (θ * ) ∞ ≤ C log d N . Thus by Theorem 3.5, we have shown A.6 Proof of Theorem 3.8 Recall the definition of the "good events" E j for j = 1, . . . , k in Section A.1 as

Moreover, we define events
Then under the assumptions of the theorem and our previous developments in Section A.1, we have To prove the claimed result, we need three auxiliary lemmas. The first lemma provides the local expansions of global loss function L N (θ) and surrogate function L(θ) around the global empirical loss minimizer θ. The proof is provided in Appendix B.5.
Lemma A.5. Under event k j=1 E j , we have that for all θ ∈ U (ρ), Our second lemma shows that the global identifiability assumption PC for L 1 (θ) implies the identifiability for the surrogate loss L(θ). The proof is provided in Appendix B.6.
Our last lemma shows that if the results in the previous two lemma holds, then we a Bernsteinvon mises result for the approximated posterior π N . The proof is provided in Appendix B.7.
where N d (µ, Σ) (·) is the pdf of a d-dim Gaussian distribution with mean vector µ and covariance matrix Σ, and the remainder term Here C is a constant independent of (n, k, N ).
Now we apply Markov's inequality, Jensen's inequality and the union bound to obtain that there exists constants c 1 , c 2 , c 3 independent of (n, k, N, d, G, L) such that P k j=0 E c j ≤ c 1 + c 2 (log 2d) 16 L 16 + c 3 G 16 k n 8 . We will apply Lemma 6 in Zhang et al. [2013b] with θ * = θ and F 1 = L in the notation therein.
The first inequality is true since under event E 0 ∩ E 1 and Assumption D, we have To prove the second inequality, we apply the integral form of Taylor's expansion to obtain that where matrix H 1 = 1 0 ∇ 2 L 1 θ * + t( θ − θ * ) dt satisfies under event E 1 . Therefore, triangle's inequality yields that under event E 0 ∩ E 1 , This proves the second inequality and therefore the claimed result.

B.3 Proof of Lemma B.3
The claimed inequality is a immediate consequence of equation ( Under Assumption D and event k j=0 E ′ j , we can bound ∇ 2 L N ( θ) as This proves the first claimed inequality.
Putting pieces together, we can prove (A.13) and therefore the claimed Bernstein-von Mises result for π N .