# Weighted Distributed Estimation under Heterogeneity

Jia Gu
Center for Statistical Science, Peking University, and
Song Xi Chen
School of Mathematical Science and Guanghua School of Management
Peking University corresponding author

###### Abstract

This paper considers distributed M-estimation under heterogeneous
distributions among distributed data blocks. A weighted distributed estimator
is proposed to improve the efficiency of the standard ”Split-And-Conquer”
(SaC) estimator for the common parameter shared by all the data blocks. The
weighted distributed estimator is shown to be at least as efficient as the
would-be full sample and the generalized method of moment estimators with the
latter two estimators requiring full data access. A bias reduction is
formulated to the WD estimator to accommodate much larger numbers of data
blocks than the existing methods without sacrificing the estimation
efficiency, and a similar debiased operation is made to the SaC estimator. The
mean squared error (MSE) bounds and the asymptotic distributions of the WD and
the two debiased estimators are derived, which shows advantageous performance
of the debiased estimators when the number of data blocks is large.

Keyword: Bias Correction; Distributed Inference; Federated Learning;
Heterogeneity; Split and Conquer Method; Weighted Estimation.

##  1 Introduction

Modern big data have brought new challenges to statistical inference. One such
challenge is that despite the shear volume of the data, a full communication
among the data points may not be possible due to the expensive cost of data
communication or the privacy concern. The distributed or the ”split-and-
conquer” (SaC) method has been proposed to divide the full data sample to
smaller size data blocks as data communication is too costly to attain an
estimation task in a timely fashion. The SaC estimator is also suited to the
situations where the data are naturally divided to data blocks and data
communication among the data blocks are prohibited due to privacy concern. The
SaC estimation had been considered in [17] for the U-statistics, [30] for the
M-estimation, [5] for the generalized linear models, [27] and [4] for the
quantile regression, while [2] studied the high dimensional testing and
estimation with sparse penalties. [3] studied the estimation efficiency and
asymptotic distributions for the general asymptotic symmetric statistics [14]
and found differences in the efficiency and the asymptotic distributions
between the non-degenerate and degenerate cases.

Bootstrap resampling-based methods had been introduced to facilitate
statistical inference. [13] proposed the bag-of-little bootstrap (BLB) method
for the plug-in estimators by making up economically the full sample for the
distributed inference. [23] suggested a sub-sampled double bootstrap method
designed to improve the computational efficiency of the BLB. [3] proposed the
distributed and the pseudo-distributed bootstrap methods with the former
conducted the resampling within each data block while the latter directly
resampled the distributed statistics.

Privacy has been a major concern in big data applications where people are
naturally reluctant to share the raw data to form a pool of big data as
practised in the traditional full sample estimation. However, the data holders
may like to contribute summary statistics without having to give away the full
data information. Federated Learning or the distributed inference with a
central host has been proposed to accommodate such reality [19, 28, 16, 12],
where summary statistics of the data blocks or the gradients of the objective
functions associated with the private data blocks are submitted to a central
host for forming aggregated estimation or computation.

Homogeneous distribution among the data blocks are assumed in majority of the
distributed inference studies with only a few exceptions [32, 6]. Federated
Learning, on the other hand, was introduced to mitigate many challenges
arising from classical distributed optimization. In particular, heterogeneous
or Non-IID distributed data across different data blocks is one of the
defining characteristics and challenges in the Federated Learning [16, 12].
Indeed, it is natural to expect the existence of heterogeneity, especially for
data stored in different locations or generated by different stochastic
mechanism, for instance mobile phones of different users. However, there has
been little published works on the statistical properties of estimators
considered in the Federated Learning.

This paper considers distributed estimation under heterogeneous distributions
among the data blocks, which is closely related to the Federated Learning and
especially the multi-task learning (MTL) [31]. We consider distributed
M-estimation where there is a common parameter shared by the distributions of
the data blocks and data-block specific heterogeneous parameters. Our
treatment of the heterogeneity is made by explicit parameterization, which is
different from the MTL where the heterogeneity is regularized by penalty
terms. It is noted that [6] considered a heterogeneous setting, but under a
fully parametric likelihood framework. Our study reveals that in the presence
of the heterogeneity the full sample M-estimator of the common parameter
obtained by requiring full data communication, can be less efficient than the
SaC estimator. However, this phenomenon disappears if the objective function
of the M-estimation satisfies a generalized second-order Bartlett’s identity,
which are satisfied by the parametric and quasi likelihoods, and the least
square estimation in the parametric regression.

We propose a weighted distributed (WD) estimator, which is asymptotically at
least as efficient as the full sample and the SaC estimator when the number of
data blocks $K=o(N^{1/2})$ where $N$ is the full sample size. The mean-squared
error (MSE) bound and the asymptotic distribution of the proposed WD estimator
are derived, as well as the asymptotic equivalence between the WD and the
generalized method of moment (GMM) estimator. We propose a debiased weighted
distributed (dWD) estimator with a data splitting mechanism on each data block
to remove the correlation between the empirical bias correction and the
weights used to tackle the heterogeneity. The dWD is asymptotically as
efficient as the WD estimator, but with a more relaxed constraint of
$K=o(N^{2/3})$. The bias-correction is also applied to the SaC formulation
leading to a more communication-efficient dSaC estimator, which is shown to be
more accurate than the subsampled average mixture estimator (SAVGM) [30] in
the homogeneous case.

The paper is organized as follows. The estimation framework and necessary
notations for the study are outlined in Section 2. The relative efficiency
between the full sample and the SaC estimators under the heterogeneity is
discussed in Section 3 to motivate the construction of the weighted
distributed (WD) estimator. The WD estimator is introduced in Section 4 along
with its efficiency, asymptotic distribution and MSE bound. Statistical
properties of two debiased estimators dSaC and dWD are revealed in Section 5.
Section 6 provides numerical verification to the theoretical results. Section
7 concludes with a discussion. Technical details are reported in the
supplementary materials (SM).

##  2 Preliminaries

Suppose that there is a large data sample of size $N$, which is divided into
$K$ data blocks of sizes $\\{n_{k}\\}_{k=1}^{K}$ such that
$N=\sum_{k=1}^{K}n_{k}$ and let $n=NK^{-1}$ be the average sample size of the
data blocks. For the relative sample size among data blocks, we assume the
following assumption.

######  Assumption 1.

There exist constants $0<c<1<C$ such that
$c\leq\frac{n_{k_{1}}}{n_{k_{2}}}\leq C$ for all pairs of $(k_{1},k_{2})$, and
if $K$ is a fixed constant we further assume that
$\frac{n_{k}}{N}\rightarrow\gamma_{k}\in(0,1)$ for a set of constants
$\\{\gamma_{k}\\}_{k=1}^{K}$.

The $k$-th data block consists of a sub-sample $\\{X_{k,i}\\}_{i=1}^{n_{k}}$
which are independent and identically distributed (IID) random vectors from a
probability space $(\Omega,\mathcal{F},P)$ to
$(\mathbb{R}^{d},\mathcal{R}^{d})$ with $F_{k}$ as the distribution. The K
distributions $\\{F_{k}\\}$ share a common parameter
$\phi\in\mathbb{R}^{p_{1}}$, while each $F_{k}$ has another parameter
$\lambda_{k}\in\mathbb{R}^{p_{2}}$ specific to $F_{k}$ of the $k$-th data
block. There are maybe other hidden parameters which define $F_{k}$, which are
however not directly involved in the semi-parametric M-estimation, and thus
are not of interest in the study.

The parameters of interests in the $k$-th block are
$\theta_{k}=(\phi^{T},\lambda_{k}^{T})^{T}$, and the overall parameters of
interests are
$\theta=(\phi^{T},\lambda_{1}^{T},\lambda_{2}^{T},...,\lambda_{K}^{T})^{T}\in\mathbb{R}^{p_{1}+Kp_{2}}$.
Suppose there is a common objective function $M(X;\phi,\lambda_{k})$ that is
convex with respect to the parameter $(\phi,\lambda_{k})$ and facilitates the
M-estimation of the parameters in each data block. In general, the criteria
function can be made block specific, say $M_{k}$ function. Indeed, the
presence of the heterogeneous local parameters $\\{\lambda_{k}\\}_{k=1}^{K}$
leads to different $M_{k}(x,\phi)=M(x,\phi,\lambda_{k})$ for the inference on
$\phi$, which connects to the multi-task learning (MTL).

In the $k$-th data block the true parameter
$\theta_{k}^{*}=(\phi^{*T},\lambda_{k}^{*T})^{T}$ is defined as the unique
minimum of the expected objective function, namely

$\theta_{k}^{*}=(\phi^{*T},\lambda_{k}^{*T})^{T}=\underset{\theta_{k}\in\Theta_{k}}{argmin}\quad\mathbb{E}_{F_{k}}M(X_{k,1};\phi,\lambda_{k}).$
(1)

The true common parameter $\phi^{\ast}$ appears in all $\theta_{k}^{*}$, and
the block-specific $\\{\lambda_{k}^{*T}\\}_{k=1}^{K}$ may differ from each
other. The entire set of true parameters
$\theta^{*}=(\phi^{*T},\lambda_{1}^{*T},\cdots,\lambda_{K}^{*T})^{T}$, can be
also identified as

$\theta^{*}=\underset{\theta\in\Theta}{argmin}\sum_{k=1}^{K}\gamma_{k}\mathbb{E}_{\theta_{k}^{*}}M(X_{k,1};\phi,\lambda_{k}).$
(2)

If the data could be shared across the data blocks, we would attain the
conventional full sample M-estimator

$\hat{\theta}_{full}=\underset{\theta\in\Theta}{argmin}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}M(X_{k,i};\phi,\lambda_{k}),$
(3)

which serves as a benchmark for the distributed estimators. Let
$\psi_{\phi}(X_{k,i};\phi,\lambda_{k})=\frac{\partial
M(X_{k,i};\phi,\lambda_{k})}{\partial\phi}$ and
$\psi_{\lambda}(X_{k,i};\phi,\lambda_{k})=\frac{\partial
M(X_{k,i};\phi,\lambda_{k})}{\partial\lambda_{k}}$ be the score functions. The
estimating equations for the full sample M-estimators are

$\begin{cases}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}\psi_{\phi}(X_{k,i};\phi,\lambda_{k})=0,\\\
\sum_{i=1}^{n_{k}}\psi_{\lambda}(X_{k,i};\phi,\lambda_{k})=0\quad
k=1,...,K.\end{cases}$ (4)

The above full sample estimation is not attainable for the distributed
situations due to privacy or the costs associated with the data
communications. The distributed estimation first conducts local estimation on
each data block, namely the local M-estimator

$\hat{\theta}_{k}=(\hat{\phi}_{k},\hat{\lambda}_{k})=\underset{\theta_{k}\in\Theta_{k}}{argmin}\sum_{i=1}^{n_{k}}M(X_{k,i};\theta_{k})$

with the corresponding estimating equations

$\begin{cases}\sum_{i=1}^{n_{k}}\psi_{\phi}(X_{k,i};\phi_{k},\lambda_{k})=0,\\\
\sum_{i=1}^{n_{k}}\psi_{\lambda}(X_{k,i};\phi_{k},\lambda_{k})=0.\end{cases}$
(5)

Then, the ”split-and-conquer” (SaC) estimator for the common parameter $\phi$
is

$\hat{\phi}^{SaC}=\frac{1}{N}\sum_{k=1}^{K}n_{k}\hat{\phi}_{k}.$ (6)

The heterogeneity among the distributions and the inference models among the
data blocks bring new dimensions to the discussion of the relative efficiency
and the estimation errors, which are the focus of this paper. We are to show
that the conventionally weighted SaC estimator (6) may not be the best
formulation for the estimation of $\phi$. Throughout this paper, unless
otherwise stated, $\|\cdot\|_{2}$ and
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}$
represent the $L_{2}$ norm of a vector and a matrix, respectively. Besides, we
will use $C$ and $C_{i}$ to denote absolute positive constants independent of
$(n_{k},K,N)$.

An important question is the efficiency and the estimation errors of the SaC
estimator $\hat{\phi}^{SaC}$ relative to the full sample estimator
$\hat{\phi}_{full}$. For the homogeneous case, Chen and Peng (2021) [3] found
that for the asymptotic symmetric statistics, the SaC estimator (6) attains
the same efficiency of the full sample estimator in the non-degenerate case,
but encounters an efficiency loss in the degenerate case due to a lack of
communications among different data blocks. Zhang et al. (2013) [30] derived
the mean square error (MSE) bound for the SaC estimator in the homogeneous
case and showed that whenever $K\leq\sqrt{N}$, the SaC estimator achieves the
best possible rate of convergence when all $N$ samples are accessible.

Consider the simultaneous estimating equations of the full sample M-estimation

$\Psi_{N}(\mathbf{X};\theta)=\begin{pmatrix}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}\psi_{\phi}(X_{k,i};\phi,\lambda_{k})\\\
\sum_{i=1}^{n_{1}}\psi_{\lambda}(X_{1,i};\phi,\lambda_{1})\\\ \vdots\\\
\sum_{i=1}^{n_{K}}\psi_{\lambda}(X_{K,i};\phi,\lambda_{K})\end{pmatrix}.$ (7)

Define

$\displaystyle\Psi_{\theta}(\theta_{k})$ $\displaystyle=$
$\displaystyle(\Psi_{\phi}(\theta_{k})^{T},\Psi_{\lambda}(\theta_{k})^{T})^{T}=\mathbb{E}\nabla_{\theta_{k}}M(X_{k,1};\theta_{k}),$
$\displaystyle\Psi_{\theta}^{\theta}(\theta_{k})$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\Psi_{\phi}^{\phi}(\theta_{k})&\Psi_{\phi}^{\lambda}(\theta_{k})\\\
\Psi_{\lambda}^{\phi}(\theta_{k})&\Psi_{\lambda}^{\lambda}(\theta_{k})\end{pmatrix}=\mathbb{E}\nabla_{\theta_{k}}^{2}M(X_{k,1};\theta_{k}),$
$\displaystyle J_{\phi|\lambda}(\theta_{k})$ $\displaystyle=$
$\displaystyle\Psi_{\phi}^{\phi}(\theta_{k})-\Psi_{\phi}^{\lambda}(\theta_{k})\Psi_{\lambda}^{\lambda}(\theta_{k})^{-1}\Psi_{\lambda}^{\phi}(\theta_{k})\quad\hbox{and}$
$\displaystyle S_{\phi}(X_{k,i};\theta_{k})$ $\displaystyle=$
$\displaystyle\psi_{\phi}(X_{k,i};\theta_{k})-\Psi_{\phi}^{\lambda}(\theta_{k})\Psi_{\lambda}^{\lambda}(\theta_{k})^{-1}\psi_{\lambda}(X_{k,i};\theta_{k}).$

Then we can apply Taylor’s expansion and obtain (see Section 1.1 in SM for
details)

$\hat{\phi}_{full}-\phi^{*}=-\\{\sum_{k=1}^{K}\frac{n_{k}}{N}J_{\phi|\lambda}(\theta_{k}^{*})\\}^{-1}\frac{1}{N}\big{\\{}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}S_{\phi}(X_{k,i};\theta_{k}^{*})\big{\\}}+o_{p}(N^{-1/2}),$
(8)

For the local estimator $(\hat{\phi}_{k},\hat{\lambda}_{k})$ based on the
$k$-th data block that solves (5), by replicating the same derivation leading
to (8), we have

$\begin{cases}\hat{\phi}_{k}-\phi^{*}&=-n_{k}^{-1}J_{\phi|\lambda}(\theta_{k}^{*})^{-1}\sum_{i=1}^{n_{k}}S_{\phi}(X_{k,i};\theta_{k}^{*})+o_{p}(n_{k}^{-1/2}),\\\
\hat{\lambda}_{k}-\lambda_{k}^{*}&=-n_{k}^{-1}J_{\lambda|\phi}(\theta_{k}^{*})^{-1}\sum_{i=1}^{n_{k}}S_{\lambda}(X_{k,i};\theta_{k}^{*})+o_{p}(n_{k}^{-1/2}),\end{cases}$
(9)

where

$\displaystyle J_{\lambda|\phi}(\theta_{k})$ $\displaystyle=$
$\displaystyle\Psi_{\lambda}^{\lambda}(\theta_{k})-\Psi_{\lambda}^{\phi}(\theta_{k})\Psi_{\phi}^{\phi}(\theta_{k})^{-1}\Psi_{\phi}^{\lambda}(\theta_{k})\quad\hbox{
and }$ $\displaystyle S_{\lambda}(X_{k,i};\theta_{k})$ $\displaystyle=$
$\displaystyle\psi_{\lambda}(X_{k,i};\theta_{k})-\Psi_{\lambda}^{\phi}(\theta_{k})\Psi_{\phi}^{\phi}(\theta_{k})^{-1}\psi_{\phi}(X_{k,i};\theta_{k}).$
(10)

The distributed inference setting is closely related to the Multi-Task
Learning (MTL) which fits separate local parameters
$\phi_{k}\in\mathbb{R}^{p}$ to the data of different data blocks (tasks)
through convex loss functions $\\{\ell_{k}\\}$. In particular, the MTL is
formulated as [24]:

$\underset{\Phi,\Omega}{min}\quad\bigg{\\{}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}\ell_{k}(\phi_{k}^{T}X_{k,i},Y_{k,i})+\mathcal{R}(\Phi,\Omega)\bigg{\\}},$
(11)

where $\\{(X_{k,i},Y_{k,i}),i=1,2,\cdots,n_{k}\\}$ are data in the $k$-th
block, $\Phi$ is the matrix with $\\{\phi_{k}\\}_{k=1}^{K}$ as column vectors,
$\Omega\in\mathbb{R}^{K\times K}$ and $\mathcal{R}(\cdot,\cdot)$ measures the
extent of the heterogeneity among different data blocks. Choices of
$\mathcal{R}(\cdot,\cdot)$ include
$\mathcal{R}(\Phi,\Omega)=\delta_{1}tr(\Phi\Omega\Phi^{T})+\delta_{2}\|\Phi\|_{F}^{2}$
for $\delta_{1},\delta_{2}>0$ and $\Omega=I_{K\times
K}-\frac{1}{K}1_{K}1_{K}^{T}$ such that
$tr(\Phi\Omega\Phi^{T})=\sum_{k=1}^{K}\|\phi_{k}-\bar{\phi}_{K}\|_{2}^{2}$
where $\bar{\phi}_{K}=\frac{1}{K}\sum_{k=1}^{K}\phi_{k}$, which leads to the
mean-regularized MTL [7]. The second term of $\mathcal{R}$ performs
regularization on each local model, trying to control the magnitude of the
estimates of $\phi_{k}$.

The distributed framework is well connected to the MTL in two key aspects. One
is that despite we use the same objective (loss) function $M$ over the data
blocks, the heterogeneity induced by local parameters
$\\{\lambda_{k}\\}_{k=1}^{K}$ and the distributions effectively define
$M_{k}(\phi,x)=M(x,\phi,\lambda_{k})$, which is equivalent to the block
specific loss functions $\ell_{k}$ used in MTL. Another aspect is that
although the MTL assumes different parameters $\\{\phi_{k}\\}$ over the data
blocks, it regularizes them toward a common one. In contrast, we assume there
is a common parameter $\phi$ shared by the heterogeneous distributions.

##  3 Full Sample versus SaC Estimation

It is naturally expected that the full sample estimator $\hat{\phi}_{full}$
should be at least as efficient as the distributed SaC estimator
$\hat{\phi}^{SaC}$ since the former utilizes the full sample information
including the communications among different data blocks. However, we are to
show that this is not necessarily true in the presence of heterogeneity. To
appreciate this point, we first list more regularity conditions needed in the
analysis.

######  Assumption 2.

(Identifiability) The parameters $\theta_{k}^{*}=(\phi^{*},\lambda_{k}^{*})$
is the unique minimizer of $M_{k}(\theta_{k})=\mathbb{E}M(X_{k,1};\theta_{k})$
for $\theta_{k}\in\Theta_{k}$.

######  Assumption 3.

(Compactness) The parameter space $\Theta_{k}$ is a compact and convex set in
$\mathbb{R}^{p}$ and the true parameter $\theta_{k}^{*}$ is an interior point
of $\Theta_{k}$ and
$\underset{\theta_{k}\in\Theta_{k}}{sup}\|\theta_{k}-\theta_{k}^{*}\|_{2}\leq
r$ for all $k\geq 1$ and some $r>0$. The true common parameter $\phi^{*}$ is
an interior point of a compact and convex set $\Phi\subset\Theta_{k}$.

######  Assumption 4.

(Local strong convexity) The population objective function on the k-th data
block $M_{k}(\theta_{k})=\mathbb{E}M(X_{k,1};\theta_{k})$ is twice
differentiable, and there exists a constant $\rho_{-}>0$ such that
$\nabla^{2}_{\theta_{k}}M_{k}(\theta_{k}^{*})\succeq\rho_{-}I_{p\times p}$.
Here $A\succeq B$ means $A-B$ is a positive semi-definite matrix.

These three assumptions are standard ones on the parameter space and
population objective functions as those in Zhang et al. (2013) [30] and Jordan
et al. (2019) [11] for the homogeneous case. In the heterogeneous setting,
Duan et al. (2021) [6] only requires the parameter space for the common
parameter to be bounded, i.e. $\|\phi-\phi^{*}\|\leq r$ under a fully
parametric setting, while in our assumption, we need the overall parameter
space to be bounded. This stronger assumption is needed since we do not fully
specify the distributions $\\{F_{k}\\}_{k=1}^{K}$ of the random variables and
will be useful when we derive the MSE bound for the weighted distributed
estimator which will be proposed in Section 4.

######  Assumption 5.

(Smoothness) There are finite positive constants $R,L,v$ and $v_{1}$ such that
for all $k\geq 1$,
$\mathbb{E}\|\nabla_{\theta_{k}}M(X_{k,1};\theta_{k}^{*})\|_{2}^{2v_{1}}\leq
R^{2v_{1}}\text{ and
}\mathbb{E}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\nabla_{\theta_{k}}^{2}M(X_{k,1};\theta_{k}^{*})-\nabla_{\theta_{k}}^{2}M_{k}(\theta_{k}^{*})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}^{2v}\leq
L^{2v}.$ In addition, for any $x\in\mathbb{R}^{d}$,
$\nabla_{\theta_{k}}^{2}M(x;\theta_{k})$ and
$\nabla_{\theta_{k}}M(x;\theta_{k})\nabla_{\theta_{k}}M(x;\theta_{k})^{T}$ are
$G(x)-$ and $B(x)-$Lipschitz continuous, respectively, in the sense that

${\left|\kern-0.96873pt\left|\kern-0.96873pt\left|\nabla_{\theta_{k}}^{2}M(x;\theta_{k})-\nabla_{\theta_{k}}^{2}M(x;\theta_{k}^{{}^{\prime}})\right|\kern-0.96873pt\right|\kern-0.96873pt\right|}_{2}\leq
G(x)\|\theta_{k}-\theta_{k}^{{}^{\prime}}\|_{2},$
${\left|\kern-0.96873pt\left|\kern-0.96873pt\left|\nabla_{\theta_{k}}M(x;\theta_{k})\nabla_{\theta_{k}}M(x;\theta_{k})^{T}-\nabla_{\theta_{k}}M(x;\theta_{k}^{{}^{\prime}})\nabla_{\theta_{k}}M(x;\theta_{k}^{{}^{\prime}})^{T}\right|\kern-0.96873pt\right|\kern-0.96873pt\right|}_{2}\leq
B(x)\|\theta_{k}-\theta_{k}^{{}^{\prime}}\|_{2},$

for all $\theta_{k},\theta_{k}^{{}^{\prime}}\in
U_{k}:=\\{\theta_{k}|\|\theta_{k}-\theta_{k}^{*}\|_{2}\leq\rho\\}$ for some
$\rho>0$, and $\mathbb{E}G(X_{k,1})^{2v}\leq
G^{2v},\mathbb{E}B(X_{k,1})^{2v}\leq B^{2v}$ for some positive constants $G$
and $B$.

The Lipschitz continuity of the outer product of the first-order derivative is
required to control the estimation error when we estimate the asymptotic
covariance matrix of the local estimator $\hat{\theta}_{k}$, and it can be
directly verified under the logistic regression case; see Section 1.2 in the
SM for details.

######  Proposition 1.

Under Assumptions 1 \- 4 and Assumption 5 with $v,v_{1}\geq 1$, and if $K$ is
fixed, then $\hat{\theta}_{k}\overset{P}{\rightarrow}\theta_{k}^{*}$ and
$\hat{\theta}_{full}\overset{P}{\rightarrow}\theta^{*}$;
$\hat{\phi}^{SaC}=\frac{1}{N}\sum_{k=1}^{K}n_{k}\hat{\phi}_{k}$ and
$\hat{\phi}_{full}$ are consistent to $\phi^{\ast}$.

######  Theorem 1.

Under Assumptions 1 \- 4 and Assumption 5 with $v,v_{1}\geq 2$, if $K$ is a
fixed constant, the SaC estimator $\hat{\phi}^{SaC}$ and the full sample
estimator $\hat{\phi}_{full}$ satisfy

$\displaystyle\sqrt{N}(\hat{\phi}^{SaC}-\phi^{*})$
$\displaystyle\overset{d}{\rightarrow}\mathcal{N}(\mathbf{0},\sum_{k=1}^{K}\gamma_{k}J_{\phi|\lambda}(\theta_{k}^{*})^{-1}\Sigma_{k}(\theta_{k}^{*})J_{\phi|\lambda}(\theta_{k}^{*})^{-1}),$
(12a) $\displaystyle\sqrt{N}(\hat{\phi}_{full}-\phi^{*})$
$\displaystyle\overset{d}{\rightarrow}\mathcal{N}(\mathbf{0},(\sum_{k=1}^{K}\gamma_{k}J_{\phi|\lambda}(\theta_{k}^{*}))^{-1}(\sum_{k=1}^{K}\gamma_{k}\Sigma_{k}(\theta_{k}^{*}))(\sum_{k=1}^{K}\gamma_{k}J_{\phi|\lambda}(\theta_{k}^{*}))^{-1}),$
(12b)

where
$J_{\lambda|\phi}(\theta_{k}^{*})=\Psi_{\lambda}^{\lambda}(\theta_{k}^{*})-\Psi_{\lambda}^{\phi}(\theta_{k}^{*})\Psi_{\phi}^{\phi}(\theta_{k}^{*})^{-1}\Psi_{\phi}^{\lambda}(\theta_{k}^{*})$
and $\Sigma_{k}=Var\\{S_{\phi}(X_{k,1};\theta_{k}^{*})\\}$.

Define $V(\Sigma,A)=(A^{T})^{-1}\Sigma A^{-1}$ as a mapping from
$\mathbb{S}_{++}^{p_{1}\times p_{1}}\times GL(\mathbb{R}^{p_{1}})$ to
$\mathbb{S}_{++}^{p_{1}\times p_{1}}$, where $\mathbb{S}_{++}^{p_{1}\times
p_{1}}$ and $GL(\mathbb{R}^{p_{1}})$ denote the symmetric positive definite
matrices and invertible real matrices of order $p_{1}$, respectively. Since
$\Sigma_{k=1}^{K}\gamma_{k}=1$ and $\gamma_{k}>0$, the asymptotic variance of
$\hat{\phi}^{SaC}$ can be interpreted as a convex combination of function
values
$\\{V(\Sigma_{k}(\theta_{k}^{*}),J_{\phi|\lambda}(\theta_{k}^{*}))\\}_{k=1}^{K}$
and that of $\hat{\phi}_{full}$ can be seen as
$V(\sum_{k=1}^{K}\gamma_{k}\Sigma_{k}(\theta_{k}^{*}),\sum_{k=1}^{K}\gamma_{k}J_{\phi|\lambda}(\theta_{k}^{*}))$.
However, $V(\cdot,\cdot)$ is not convex with respect to its arguments
$(\Sigma,A)$, which means that the inequality

$\\{\sum_{k=1}^{K}\gamma_{k}J_{\phi|\lambda}(\theta_{k}^{*})\\}^{-1}\\{\sum_{k=1}^{K}\gamma_{k}\Sigma_{k}(\theta_{k}^{*})\\}\\{\sum_{k=1}^{K}\gamma_{k}J_{\phi|\lambda}(\theta_{k}^{*})\\}^{-1}\preceq\sum_{k=1}^{K}\gamma_{k}J_{\phi|\lambda}(\theta_{k}^{*})^{-1}\Sigma_{k}(\theta_{k}^{*})J_{\phi|\lambda}(\theta_{k}^{*})^{-1}$

does not always hold. In other words, $\hat{\phi}_{full}$ is not necessarily
more efficient than $\hat{\phi}^{SaC}$.

To gain understanding of Theorem 1 and to motivate the weighted distributed
estimator, we consider the errors-in-variables model. Suppose that one
observes $K$ blocks of independent data samples
$\\{(X_{k,i},Y_{k,i})\\}_{i=1}^{n}$ for $k=1,2...,K$ and $N=nK$, where
$(X_{k,i},Y_{k,i})$ are IID and generated from the following model:

$\begin{cases}X_{k}=Z_{k}+e_{k},\\\
Y_{k}=\phi^{*}+\lambda_{k}^{*}Z_{k}+f_{k},\\\ \end{cases}$ (13)

where $\\{Z_{k}\\}_{k=1}^{K}$ are random variables whose measurements
$\\{(X_{k},Y_{k})\\}_{k=1}^{K}$ are subject to errors
$\\{(e_{k},f_{k})\\}_{k=1}^{K}$, and $(e,f)$ is bivariate normally distributed
with zero mean and covariance matrix $\sigma^{2}I_{2}$ and is independent of
$Z_{k}$. Here, $\phi^{*}$ is the common parameter across all data blocks while
$\lambda_{k}^{*}(\lambda_{k}^{*}>0)$ represents the block specific parameter.
We assume that $Var(e)=Var(f)$ to avoid any identification issue arisen when
$Z$ is also normally distributed [20]. There is a considerable literature on
the regression problem with measurement errors, as summarised in [8, 22].

We consider the approach displayed in Example 5.26 of [26] which constructs a
kind of marginal likelihood followed by centering to make a bona fide score
equation, as detailed in Section 1.3 of the SM. The M-function is

$M(X_{k},\theta_{k})=\frac{1}{2\sigma^{2}(1+\lambda_{k}^{2})}(\lambda_{k}X_{k}-(Y_{k}-\phi))^{2},$
(14)

with the score equation satisfying $\mathbb{E}\nabla
M(X_{k,1},Y_{k,1}|Z_{k,1},\theta_{k}^{*})=\mathbf{0}_{2\times 1}$.

For simplicity we assume $K=2$, then from Theorem 1 we have

$\begin{cases}Var(\hat{\phi}_{full})\approx\\{\frac{\sigma^{2}\mathbb{E}Z^{2}}{var(Z)}\frac{2}{\frac{1}{1+\lambda_{1}^{*2}}+\frac{1}{1+\lambda_{2}^{*2}}}+\frac{\sigma^{4}(\mathbb{E}Z)^{2}}{var^{2}(Z)}\frac{\frac{2}{(1+\lambda_{1}^{*2})^{2}}+\frac{2}{(1+\lambda_{2}^{*2})^{2}}}{(\frac{1}{1+\lambda_{1}^{*2}}+\frac{1}{1+\lambda_{2}^{*2}})^{2}}\\}\frac{1}{N},\\\
Var(\hat{\phi}^{SaC})\approx\\{\frac{\sigma^{2}\mathbb{E}Z^{2}}{var(Z)}\frac{(1+\lambda_{1}^{*2})+(1+\lambda_{2}^{*2})}{2}+\frac{\sigma^{4}(\mathbb{E}Z)^{2}}{var^{2}(Z)}\\}\frac{1}{N}.\end{cases}$
(15)

Note that the coefficients to $\frac{\sigma^{2}\mathbb{E}Z^{2}}{var(Z)}$ in
the first terms of the variances are harmonic and arithmetic means of
$\\{1+\lambda_{1}^{*2},1+\lambda_{2}^{*2}\\}$, respectively. By the mean
inequality the coefficient in the first term of $Var(\hat{\phi}^{SaC})$ is
larger than that in $Var(\hat{\phi}_{full})$. The second term of the variances
involves $(\mathbb{E}Z)^{2}$ as a multiplicative factor. Thus, if the
unobserved $Z$ has zero mean, the full-sample estimator would be at least as
good as the SaC estimator in terms of variation when the sample size goes to
infinity. However, the story may change when $\mathbb{E}Z\not=0$, because the
second term of $Var(\hat{\phi}_{full})$ has a factor which is the square of a
ratio between the quadratic mean and the arithmetic mean of
$(\frac{1}{1+\lambda_{1}^{*2}},\frac{1}{1+\lambda_{2}^{*2}})$. The factor is
larger than or equal to 1 if and only if $\lambda_{1}^{*}=\lambda_{2}^{*}$
namely the homogeneous case. In the heterogeneous case, by adjusting
$\frac{\sigma^{4}(\mathbb{E}Z)^{2}}{var^{2}(Z)}/\frac{\sigma^{2}\mathbb{E}Z^{2}}{var(Z)}$,
we can find cases such that $\lambda_{1}^{*}\neq\lambda_{2}^{*}$ such that the
full sample estimator has a larger variance than the SaC estimator. Simulation
experiments presented in Section 6 display such cases.

##  4 Weighted Distributed Estimator

The previous section shows that the full sample estimator $\hat{\phi}_{full}$
under heterogeneity may be less efficient than the simple averaged
$\hat{\phi}^{SaC}$. This phenomenon suggests that the conventional wisdom in
the homogeneous context case may not be applicable to the heterogeneous case.
One may also wonder if the simple SaC estimator can be improved under the
heterogeneity. Specifically, how to better aggregate the local estimator
$\hat{\phi}_{k}$ for more efficiency estimation to the common parameter $\phi$
is the focus of this section.

###  4.1 Formulation and Results

Consider a class of estimators formed by linear combinations of the local
estimators $\\{\hat{\phi}_{k}\\}$:

$\\{\hat{\phi}_{w}^{SaC}|\hat{\phi}_{w}^{SaC}=\sum_{k=1}^{K}W_{k}\hat{\phi}_{k},W_{k}\in\mathbb{R}^{p_{1}\times
p_{1}},\sum_{k=1}^{K}W_{k}=I_{p_{1}}\\}.$

We want to minimize the asymptotic variance of $\hat{\phi}_{w}^{SaC}$ with
respect to $\\{W_{k}\\}_{k=1}^{K}$. According to a generalization of Theorem 1

$AsyVar(\hat{\phi}_{w}^{SaC})=\sum_{k=1}^{K}n_{k}^{-1}W_{k}A_{k}^{-1}\Sigma_{k}(A_{k}^{T})^{-1}W_{k}^{T},$
(16)

where $A_{k}=J_{\phi|\lambda}(\theta_{k}^{*})$ and
$\Sigma_{k}=Var\\{S_{\phi}(X_{k,i};\theta_{k}^{*})\\}$. It is noted that the
asymptotic variance is defined via the asymptotic normality of the
M-estimation. For the time being, $A_{k}$ and $\Sigma_{k}$ are assumed known
and denote $H_{k}=A_{k}^{-1}\Sigma_{k}(A_{k}^{T})^{-1}$. We choose the trace
operator as a measure on the size of the asymptotic covariance matrix and this
leads to the minimization problem

$\underset{W_{k}}{Minimize}\quad
tr\bigg{(}\sum_{k=1}^{K}n_{k}^{-1}W_{k}H_{k}W_{k}^{T}\bigg{)}\quad
s.t.\quad\sum_{k=1}^{K}W_{k}=I_{p_{1}},$ (17)

which is a convex optimization problem. It can be solved via the Lagrangian
multiplier method which gives
$W_{k}^{*}=(\sum_{s=1}^{K}n_{s}H_{s}^{-1})^{-1}n_{k}H_{k}^{-1}$. If we replace
the trace with the Frobenius norm in the objective function (17), the same
solution is attained as shown in Section 1.4 of the SM. The SaC estimator
under the optimal weights $W_{k}^{*}$ is called the weighted distributed (WD)
estimator and denoted as $\hat{\phi}^{WD}$. By construction, the WD estimator
is at least as efficient as the SaC estimator (6). To compare the relative
efficiency between $\hat{\phi}_{full}$ and $\hat{\phi}^{WD}$, we note that

$\displaystyle AsyVar(\hat{\phi}_{full})$ $\displaystyle=$
$\displaystyle\big{\\{}(\sum_{k=1}^{K}n_{k}A_{k})^{T}(\sum_{k=1}^{K}n_{k}\Sigma_{k})^{-1}(\sum_{k=1}^{K}n_{k}A_{k})\big{\\}}^{-1}\quad\hbox{and}$
$\displaystyle AsyVar(\hat{\phi}^{WD})$ $\displaystyle=$
$\displaystyle\left(\sum_{k=1}^{K}n_{k}A_{k}^{T}\Sigma_{k}^{-1}A_{k}\right)^{-1}.$
(18)

Define $F(\Sigma,A)=A^{T}\Sigma^{-1}A$. If we can show the convexity of $F$,
an application of Jensen’s inequality will establish the relative efficiency
of the two estimators. In fact, we have the following lemma.

######  Lemma 1.

Suppose $H$ and $K$ are positive definite matrices of order $p$, and $X$ and
$Y$ are arbitrary $p\times m$ matrices. Then,

$Q=X^{T}H^{-1}X+Y^{T}K^{-1}Y-(X+Y)^{T}(H+K)^{-1}(X+Y)\succeq\mathbf{0}.$

The lemma implies that

$(\sum_{k=1}^{K}n_{k}A_{k})^{T}(\sum_{k=1}^{K}n_{k}\Sigma_{k})^{-1}(\sum_{k=1}^{K}n_{k}A_{k})\preceq\sum_{k=1}^{K}n_{k}A_{k}^{T}\Sigma_{k}^{-1}A_{k},$

which means that the WD estimator is at least as efficient as the full sample
estimator, and can be more efficient than $\hat{\phi}_{full}$. That is to say,
the simultaneous estimating equations (7), which are obtained from the first-
order derivative of the the simple summation of local objectives
$\sum_{i=1}^{n_{k}}M(X_{k,i};\theta_{k})$, are not the best formulation of the
M-estimation problem, since the formulation itself does not utilize the
heterogeneity existed in the data blocks. In contrast, the WD estimator
exploits the potential efficiency gain from the heterogeneity by re-weighting
of the local estimators, and this is why the full sample estimator may not be
as efficient as the WD estimator.

###  4.2 Likelihood and Quasi-likelihood

The above results lead us to wonder whether we can attain more efficient
distributed estimators than the full sample estimator under the heterogeneity
if we restrict to a fully parametric setting. When the distribution of
$X_{k,i}$ is fully parametric with density function
$f(\cdot;\phi,\lambda_{k})$, the Fisher information matrix in the $k$-th data
block is

$I(\theta_{k})=I(\phi,\lambda_{k})=\begin{pmatrix}I_{\phi\phi}&I_{\phi\lambda_{k}}\\\
I_{\lambda_{k}\phi}&I_{\lambda_{k}\lambda_{k}}\end{pmatrix}=-\mathbb{E}\begin{pmatrix}\frac{\partial^{2}}{\partial\phi^{2}}logf(X_{k,1};\theta_{k})&\frac{\partial^{2}}{\partial\phi\partial\lambda^{T}}logf(X_{k,1};\theta_{k})\\\
\frac{\partial^{2}}{\partial\lambda\partial\phi^{T}}logf(X_{k,1};\theta_{k})&\frac{\partial^{2}}{\partial\lambda^{2}}logf(X_{k,1};\theta_{k})\end{pmatrix},$

and the partial information matrix as
$I_{\phi|\lambda_{k}}=I_{\phi\phi}-I_{\phi\lambda_{k}}I_{\lambda_{k}\lambda_{k}}^{-1}I_{\lambda_{k}\phi}$.
Now, the objective function for the M-estimation (also the maximum likelihood
estimation (MLE)) is $M(X_{k,i};\phi,\lambda_{k})=-\log
f(X_{k,i};\phi,\lambda_{k})$. Routine derivations show that
$\Sigma_{k}=Var\\{S_{\phi}(X_{k,1};\theta_{k}^{*})\\}=I_{\phi|\lambda_{k}}$
and $A_{k}=J_{\phi|\lambda}(\theta_{k}^{*})=I_{\phi|\lambda_{k}}$. Thus,

$\displaystyle AsyVar(\hat{\phi}_{full})$ $\displaystyle=$ $\displaystyle
AsyVar(\hat{\phi}^{WD})=\left(\sum_{k=1}^{K}n_{k}I_{\phi|\lambda_{k}}\right)^{-1}\quad\hbox{
and }$ $\displaystyle AsyVar(\hat{\phi}^{SaC})$ $\displaystyle=$
$\displaystyle\frac{1}{N^{2}}\sum_{k=1}^{K}n_{k}I_{\phi|\lambda_{k}}^{-1}.$

A direct application of Lemma 1 shows that

$AsyVar(\hat{\phi}_{full})=AsyVar(\hat{\phi}^{WD})\preceq
AsyVar(\hat{\phi}^{SaC}).$ (19)

Thus, the full sample MLE can automatically adjust for the heterogeneity and
is at least as efficient as SaC estimator $\hat{\phi}^{SaC}$. Besides, the
weighted distributed estimators $\hat{\phi}^{WD}$ can fully recover the
efficiency gap of the SaC estimator.

The same relationship among $\hat{\phi}_{full},\hat{\phi}^{SaC}$ and
$\hat{\phi}^{WD}$ also holds for the maximum quasi-likelihood estimator (MQLE)
with independent observations (see Section 1.5 in the SM for details). If one
looks into the asymptotic variances of the MLE and MQLE, it can be found that
the underlying reason for (19) is that the two special M-estimation functions
satisfy the second order Bartlett’s identity [1, 18]:

$\mathbb{E}\nabla M(X_{k},\theta_{k}^{*})\nabla
M(X_{k},\theta_{k}^{*})^{T}=\mathbb{E}\nabla^{2}M(X_{k},\theta_{k}^{*}).$

By the variance formula of the asymptotic distribution of the M-estimator and
Lemma 1, we readily have that the Bartlett’s identity can be relaxed by
inserting a factor $\gamma\not=0$ such that

$\mathbb{E}\nabla
M(X_{k},\theta_{k}^{*})\nabla^{T}M(X_{k},\theta_{k}^{*})=\gamma\mathbb{E}\nabla^{2}M(X_{k},\theta_{k}^{*}).$
(20)

An important example for such a case is the least square estimation for the
parametric regression with homoscedastic and non-autocorrelated disturbances (
see Section 1.6 in the SM for details). Otherwise the full sample least square
estimator may not be efficient and there is an opportunity for the weighted
distributed least square estimation. In summary, as long as the objective
function $M(x_{k},\theta_{k})$ satisfies (20), then $\hat{\phi}_{full}$
attains the same asymptotic efficiency as $\hat{\phi}^{WD}$, and
$\hat{\phi}^{SaC}$ is at most as efficient as the former two estimators.

###  4.3 Relative to Generalized Method of Moment Estimation

To provide a benchmark on the efficiency of the WD estimation, we consider the
generalized method of moment (GMM) estimator [9]. The GMM estimator possess
certain optimal property for semiparametric inference that the weighted
distributed estimation can compare with, despite the GMM requires more data
sharing than the distributed inference would require.

The score functions of the M-estimation on each data block can be aggregated
and combined to form the moment equations

$\begin{cases}\sum_{i=1}^{n_{k}}\psi_{\phi}(X_{k,i};\phi,\lambda_{k})=0,\\\
\sum_{i=1}^{n_{k}}\psi_{\lambda}(X_{k,i};\phi,\lambda_{k})=0,\quad
k=1,...,K.\end{cases}$ (21)

There are $pK$ estimating equations, where the dimension of $\theta^{*}$ is
$pK-(K-1)p_{1}$. Thus, the parameter is over-identified which offers potential
in efficiency gain for the GMM [9]. The GMM estimation based on the moment
restrictions (21) is asymptotically equivalent to solving the following
problem:

$\hat{\theta}_{GMM}=\underset{\theta_{k}=(\phi,\lambda_{k})\in\Theta_{k},1\leq
k\leq
K}{argmin}\quad\tilde{\psi}_{N}^{T}(\theta)W_{0}\tilde{\psi}_{N}(\theta),$
(22)

where $W_{0}=Var(\tilde{\psi}_{N}(\theta^{*}))^{-1}$ is the optimal weighting
matrix [9, 29] and

$\tilde{\psi}_{N}(\theta)=(\sum_{i=1}^{n_{1}}\psi_{\phi}(X_{1,i};\theta_{1})^{T},\sum_{i=1}^{n_{1}}\psi_{\lambda}(X_{1,i};\theta_{1})^{T},\cdots,\sum_{i=1}^{n_{K}}\psi_{\phi}(X_{K,i};\theta_{K})^{T},\sum_{i=1}^{n_{K}}\psi_{\lambda}(X_{K,i};\theta_{K})^{T})^{T}$

The asymptotic variance of the GMM estimator [9] is
$AsyVar(\hat{\theta}_{GMM})=(G_{0}^{T}W_{0}G_{0})^{-1}$, where
$G_{0}^{T}=\mathbb{E}\\{\frac{\partial\tilde{\psi}_{N}^{T}(\theta^{*})}{\partial\theta}\\}$.
A derivation given in Section 1.7 of the SM shows that

$AsyVar(\hat{\phi}_{GMM})=\\{\sum_{k=1}^{K}n_{k}J_{\phi|\lambda}\Sigma_{k}^{-1}J_{\phi|\lambda}\\}^{-1}.$
(23)

Thus, the weighted distributed estimator’s efficiency is the same as that of
the GMM estimator. This is very encouraging to the proposed WD estimator as it
attains the same efficiency as the GMM without requiring much data sharing
among the blocks, which avoids the expenses of the data transmission and
preserves the privacy of the data.

###  4.4 Estimation of Weights in one round communication

To formulate the WD estimator, the optimal weights
$W_{k}^{*}=(\sum_{s=1}^{K}n_{s}H_{s}^{-1})^{-1}n_{k}H_{k}^{-1}$ have to be
estimated. By the structure of $W_{k}^{*}$, we only need to separately
estimate $H_{k}$, the leading principal submatrix of order $p_{1}$ of the
asymptotic covariance matrix $\tilde{H}_{k}$ of $\hat{\theta}_{k}$. It is
noted that

$\tilde{H}_{k}=(\nabla\Psi_{\theta}(\theta_{k}^{*}))^{-1}\mathbb{E}\\{\psi_{\theta_{k}}(X_{k,1};\theta_{k}^{*})\psi_{\theta_{k}}(X_{k,1};\theta_{k}^{*})^{T}\\}(\nabla\Psi_{\theta}(\theta_{k}^{*}))^{-1}=\begin{pmatrix}H_{k}&*\\\
*&*\\\ \end{pmatrix},$

where
$\Psi_{\theta}(\theta_{k})=\mathbb{E}\psi_{\theta_{k}}(X_{k,1};\theta_{k})$.
We can construct the sandwich estimator [25] to estimate $\tilde{H}_{k}$ and
then $H_{k}$. The distributive procedure to attain the WD estimator is
summarized in the Algorithm 1.

Input: $\\{X_{k,i},k=1,...,K;i=1,...,n_{k}\\}$

Output: $\hat{\phi}^{WD},\hat{\lambda}_{k}$

1 Obtain the initial estimates
$\hat{\theta}_{k}=(\hat{\phi}_{k},\hat{\lambda}_{k})$ based on data block $k$
;

2 Calculate $\hat{H}_{k}(\hat{\theta}_{k})$ in each block, which is the
leading principal sub-matrix of order $p_{1}$ of
$(\nabla_{\theta_{k}}\hat{\Psi}_{\theta_{k}})^{-1}(n_{k}^{-1}\sum_{i=1}^{n_{k}}\psi_{\theta_{k}}(X_{k,i};\hat{\theta}_{k})\psi_{\theta_{k}}(X_{k,i};\hat{\theta}_{k})^{T})(\nabla_{\theta_{k}}\hat{\Psi}_{\theta_{k}})^{-T}$
where
$\hat{\Psi}_{\theta_{k}}=n_{k}^{-1}\sum_{i=1}^{n_{k}}\psi_{\theta_{k}}(X_{k,i};\hat{\theta}_{k})$;

3 Send $(\hat{\phi}_{k},\hat{H}_{k}(\hat{\theta}_{k})^{-1})$ to a central
server and construct
$\tilde{\hat{\phi}}^{WD}:=\\{\sum_{k=1}^{K}n_{k}\hat{H}_{k}(\hat{\theta}_{k})^{-1}\\}^{-1}\sum_{k=1}^{K}n_{k}(\hat{H}_{k}(\hat{\theta}_{k}))^{-1}\hat{\phi}_{k}$;

$\hat{\phi}^{WD}:=\tilde{\hat{\phi}}^{WD}I(\tilde{\hat{\phi}}^{WD}\in\Phi)+\hat{\phi}^{SaC}I(\tilde{\hat{\phi}}^{WD}\not\in\Phi)$,
where $\hat{\phi}^{SaC}:=N^{-1}\sum_{k=1}^{K}n_{k}\hat{\phi}_{k}$.

Algorithm 1 Weighted Distributed M-estimator

The Step 4 in the algorithm is necessary since there is no guarantee that
after weighting the estimator $\tilde{\hat{\phi}}^{WD}$ still belongs to the
set $\Phi$ as required in Assumption 3. However the event
$\\{\tilde{\hat{\phi}}^{WD}\in\Phi\\}$ should happen with probability
approaching one. Hence, the
$\hat{\phi}^{SAC}I(\tilde{\hat{\phi}}^{WD}\not\in\Phi)$ term is negligible
compared with that of $\hat{\phi}^{WD}I(\tilde{\hat{\phi}}^{WD}\in\Phi)$. We
need the following assumption in order to establish the MSE bound and
asymptotic properties of the proposed WD estimator.

######  Assumption 6.

(Boundedness) There exists constants $\rho_{\sigma},c>0$ such that for $k\geq
1$,

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\Sigma_{S,k}(\theta_{k}^{*})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}\leq\rho_{\sigma},H_{k}\succeq
cI_{p_{1}\times p_{1}},$

where
$\Sigma_{S,k}(\theta_{k})=\mathbb{E}\psi_{\theta_{k}}(X_{k,1};\theta_{k})\psi_{\theta_{k}}(X_{k,1};\theta_{k})^{T}$.

By the definition of $H_{k}(\theta_{k})$, we have that

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|H_{k}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\Psi_{\theta}^{\theta}(\theta_{k}^{*})^{-1}\Sigma_{S,k}(\theta_{k}^{*}){\Psi_{\theta}^{\theta}(\theta_{k}^{*})^{-1}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|{\Psi_{\theta}^{\theta}(\theta_{k}^{*})^{-1}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}^{2}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\Sigma_{S,k}(\theta_{k}^{*})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}\leq\frac{\rho_{\sigma}}{\rho_{-}^{2}},$

which implies
$H_{k}(\theta_{k}^{*})^{-1}\succeq\frac{\rho_{-}^{2}}{\rho_{\sigma}}I_{p_{1}\times
p_{1}}$. On the other hand, the above inequality leads to
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|{\Psi_{\theta}^{\theta}(\theta_{k}^{*})^{-1}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}\geq\sqrt{\frac{c}{\rho_{\sigma}}}$,
and this indicate a finite upper bound for the norm of the Hessian matrix,
just as that assumed in Jordan et al. (2019) [11] and Duan et al. (2021) [6].

######  Theorem 2.

Under Assumptions 1 \- 4 and 6, and Assumption 5 with $v,v_{1}\geq 2$ , the
mean-squared error of the WD estimator $\hat{\phi}^{WD}$ satisfies

$\mathbb{E}\|\hat{\phi}^{WD}-\phi^{*}\|_{2}^{2}\leq\frac{C_{1}}{nK}+\frac{C_{2}}{n^{2}}+\frac{C_{3}}{n^{2}K}+\frac{C_{4}}{n^{3}}+\frac{C_{5}K}{n^{\bar{v}}},$

for $n=NK^{-1}$ and $\bar{v}=min\\{v,\frac{v_{1}}{2}\\}$.

The $v$ and $v_{1}$ appeared in Assumption 5 quantify the moments of the first
two orders of the derivatives of the $M$ function and their corresponding
Lipschitz functions. When the number of data blocks
$K=O(n^{min\\{1,\frac{\bar{v}-1}{2}\\}})$, the convergence rate of MSE of
$\hat{\phi}^{WD}$ is $\mathcal{O}((nK)^{-1})$, which is the same as the
standard full sample estimator. However, when there are too many data blocks
such that $K>>n$, the convergence rate is reduced to $\mathcal{O}(n^{-2})$.
Furthermore, if the derivatives of the $M$ function and their corresponding
Lipschitz functions are heavy-tailed, say $\bar{v}<3$, the convergence rate is
further reduced to $\mathcal{O}(Kn^{-\bar{v}})$.

######  Theorem 3.

Under Assumptions 1 \- 4 and 6, and Assumption 5 with $v,v_{1}\geq 2$, if
$K=o(n)$,

$\big{(}\hat{\phi}^{WD}-\phi^{*})^{T}\\{\sum_{k=1}^{K}n_{k}H_{k}(\theta_{k}^{*})^{-1}\\}\big{(}\hat{\phi}^{WD}-\phi^{*})\overset{d}{\rightarrow}\chi_{p_{1}}^{2}.$

Although $\\{H_{k}(\theta_{k}^{*})\\}_{k=1}^{K}$ have bounded spectral norms,
$\sum_{k=1}^{K}\frac{n_{k}}{N}H_{k}(\theta_{k}^{*})^{-1}$ may not converge to
a fixed matrix in presence of heterogeneity. Thus, we can only obtain the
asymptotic normality of the standardized
$\sqrt{N}\\{\sum_{k=1}^{K}\frac{n_{k}}{N}H_{k}(\theta_{k}^{*})^{-1}\\}^{1/2}(\hat{\phi}^{DW}-\phi^{*})$.
This is why Theorem 3 is formulated in a limiting chi-squared distribution
form.

The asymptotic normality implies that we can construct confidence regions for
$\phi$ with confidence level $1-\alpha$ as

$\\{\phi\,|\big{(}\hat{\phi}^{WD}-\phi)^{T}\\{\sum_{k=1}^{K}n_{k}\hat{H}_{k}(\hat{\theta}_{k})^{-1}\\}\big{(}\hat{\phi}^{WD}-\phi)\leq\chi^{2}_{p_{1},\alpha}\\}$
(24)

after replacing $\sum_{k=1}^{K}n_{k}H_{k}(\theta_{k}^{*})^{-1}$ with its
sample counterpart $\sum_{k=1}^{K}n_{k}\hat{H}_{k}(\hat{\theta}_{k})^{-1}$,
where $\chi^{2}_{p_{1},\alpha}$ is the upper $\alpha$ quantile of the
$\chi_{p_{1}}^{2}$ distribution. The block-specific parameter $\lambda_{k}$
can also be of interest. Then given the WD estimator of the common parameter
$\phi^{*}$, a question is that whether a more efficient estimator of
$\lambda_{k}^{*}$ can be obtained. Specifically, we plug in the WD estimator
to each data block and re-estimate $\lambda_{k}$. The corresponding updated
estimator is denoted as $\hat{\lambda}_{k}^{(2)}$. Actually, the answer is
that $\hat{\lambda}_{k}^{(2)}$ is not necessarily more efficient than
$\hat{\lambda}_{k}$. Due to space limit, more discussions on this aspect are
available in Section 1.8 in SM.

##  5 Debiased Estimator for diverging K

It is noted that $K=o(\sqrt{N})$ is required in both Theorems 2 and 3 to
validate the $\mathcal{O}(N^{-1})$ leading order MSE and limiting chi-squared
distribution of the WD estimator. The reason is that the bias of the local
estimator $\hat{\theta}_{k}$ is at order $O_{p}(n_{k}^{-1})$, which can
accumulate across the data blocks by the weighted averaging. This leads to the
bias of $\sqrt{N}(\hat{\phi}^{WD}-\phi^{*})$ being at order
$O_{p}(KN^{-1/2})$, which is not necessarily diminishing to zero unless
$K=o(\sqrt{N})$. It is worth mentioning that Duan et al. (2021) [6] needed the
same $K=o(\sqrt{N})$ order in their MLE framework to obtain the
$\sqrt{N}$-convergence since Li et al. (2003) [15] showed that the MLE is
asymptotically biased when $K/n\rightarrow C\in(0,+\infty)$. This calls for a
debias step for the local estimators before aggregation to allow for larger
$K$, which is needed especially in the Federated Learning scenario where the
number of users (data blocks) can be much larger than the size of local data.

To facilitate the bias correction operation, we have to simplify the
notations. Suppose $F(\theta)$ is a $p\times 1$ vector function, $\nabla
F(\theta)$ is the usual Jacobian whose $l$-th row contains the partial
derivatives of the $l$-th element of $F(\theta)$. Then the matrices of higher
derivatives are defined recursively so that the $j$-th element of the $l$-th
row of $\nabla^{s}L(\theta)$ (a $p\times p^{s}$ matrix) is the $1\times p$
vector $f_{lj}^{v}(\theta)=\partial f_{lj}^{v-1}(\theta)/\partial\theta^{T}$,
where $f_{lj}^{v-1}$ is the $l-$th row and $j$-th element of
$\nabla^{v-1}F(\theta)$. We use $\otimes$ to denote a usual Kronecker product.
Using Kronecker product we can express
$\nabla^{v}F(\theta)=\frac{\partial^{v}F(\theta)}{\partial\theta^{T}\otimes\partial\theta^{T}\otimes\cdots\otimes\partial\theta^{T}}$.
Besides, let
$M_{n,k}(\theta_{k})=n_{k}^{-1}\sum_{i=1}^{n_{k}}M(X_{k,i};\theta_{k})$,

$\displaystyle H_{3,k}(\theta_{k})$ $\displaystyle=$
$\displaystyle\mathbb{E}\nabla_{\theta_{k}}^{2}\psi_{\theta_{k}}(X_{k,1};\theta_{k}),\quad
Q_{k}(\theta_{k})=\\{-\mathbb{E}\nabla_{\theta_{k}}\psi_{\theta_{k}}(X_{k,1};\theta_{k})\\}^{-1},$
$\displaystyle\quad d_{i,k}(\theta_{k})$ $\displaystyle=$ $\displaystyle
Q_{k}(\theta_{k})\psi_{\theta_{k}}(X_{k,i};\theta_{k})\quad\hbox{ and }\quad
v_{i,k}(\theta_{k})=\nabla_{\theta_{k}}\psi_{\theta_{k}}(X_{k,i},\theta_{k})-\nabla_{\theta_{k}}\Psi_{\theta}(\theta_{k}).$

According to [21], the leading order bias of $\hat{\theta}_{k}$ is

$Bias(\hat{\theta}_{k})=n_{k}^{-1}Q_{k}(\theta_{k}^{*})\big{(}\mathbb{E}v_{i,k}(\theta_{k}^{*})d_{i,k}(\theta_{k}^{*})+\frac{1}{2}H_{3,k}(\theta_{k}^{*})\mathbb{E}\\{d_{i,k}(\theta_{k}^{*})\otimes
d_{i,k}(\theta_{k}^{*})\\}\big{)}.$ (25)

Let
$B_{k}(\theta_{k})=Q_{k}(\theta_{k})\big{(}\mathbb{E}v_{i,k}(\theta_{k})d_{i,k}(\theta_{k})+\frac{1}{2}H_{3,k}(\theta_{k})\mathbb{E}\\{d_{i,k}(\theta_{k})\otimes
d_{i,k}(\theta_{k})\\}\big{)}$, whose the first $p_{1}$ dimension associated
with $\phi$ are denoted as $B^{1}_{k}(\theta_{k})$. The empirical estimator of
$B_{k}(\theta_{k})$ is

$\hat{B}_{k}(\theta_{k})=\hat{Q}_{k}(\theta_{k})\big{(}n_{k}^{-1}\sum_{i=1}^{n_{k}}\hat{v}_{i,k}(\theta_{k})\hat{d}_{i,k}(\theta_{k})+\frac{1}{2}\hat{H}_{3,k}(\theta_{k})n_{k}^{-1}\sum_{i=1}^{n_{k}}(\hat{d}_{i,k}(\theta_{k})\otimes\hat{d}_{i,k}(\theta_{k}))\big{)}$
(26)

where
$\hat{H}_{3,k}(\theta_{k})=n_{k}^{-1}\sum_{i=1}^{n_{k}}\nabla_{\theta_{k}}^{2}\psi_{\theta_{k}}(X_{k,i};\theta_{k})$,
$\hat{Q}_{k}(\theta_{k})=\\{-n_{k}^{-1}\sum_{i=1}^{n_{k}}\nabla_{\theta_{k}}\psi_{\theta_{k}}(X_{k,i};\theta_{k})\\}^{-1},$
$\hat{d}_{i,k}(\theta_{k})=\hat{Q}_{k}(\theta_{k})\psi_{\theta_{k}}(X_{k,i};\theta_{k})$
and
$\hat{v}_{i,k}(\theta_{k})=\nabla_{\theta_{k}}\psi_{\theta_{k}}(X_{k,i};\theta_{k})$.
Applying bias correction to each data block, we have the bias-corrected local
estimator

$\hat{\theta}_{k,bc}:=\hat{\theta}_{k}-n_{k}^{-1}\hat{B}_{k}(\hat{\theta}_{k})1_{\mathcal{E}_{k,bc}},$
(27)

where
$\mathcal{E}_{k,bc}=\\{\hat{\theta}_{k}-n_{k}^{-1}\hat{B}_{k}(\hat{\theta}_{k})\in\Theta_{k}\\}$.
The indicator function here is to ensure that $\hat{\theta}_{k,bc}$ is within
the parameter space.

After the local debiased estimators are obtained, we need to aggregate them
with estimated weights. However, a direct aggregation will invalidate the bias
correction procedure due to the correlation between the estimated weights and
the local debiased estimator as they are constructed with the same dataset.
The accumulation of the dependence over a large number of data blocks can make
the bias correction fail.

To remove such correlation between the local estimators and the corresponding
estimated local weights
$\hat{W}_{k}=\\{\sum_{s=1}^{K}\hat{H}_{s}(\hat{\theta}_{s})^{-1}\\}^{-1}\hat{H}_{k}(\hat{\theta}_{k})^{-1}$,
we first divide each local dataset $\\{X_{k,i}\\}_{i=1}^{n_{k}}$ into two
equal-sized subsets $D_{k}^{s}=\\{X_{k,i}^{(s)}\\}_{i=1}^{n_{k}/2},s=1,2$.
Then, for $s=1,2$ we calculate the local M-estimators $\hat{\theta}_{k,s}$ and
obtain $\hat{H}_{k,s}(\hat{\theta}_{k,s})$, which is the leading principal
sub-matrix of order $p_{1}$ of

$(\nabla_{\theta_{k}}\hat{\Psi}_{\theta_{k}})^{-1}(\frac{1}{n_{k}/2}\sum_{i=1}^{n_{k}/2}\psi_{\theta_{k}}(X_{k,i}^{(s)};\hat{\theta}_{k,s})\psi_{\theta_{k}}(X_{k,i}^{(s)};\hat{\theta}_{k,s})^{T})(\nabla_{\theta_{k}}\hat{\Psi}_{\theta_{k}})^{-T},$

where
$\hat{\Psi}_{\theta_{k}}=\frac{1}{n_{k}/2}\sum_{i=1}^{n_{k}/2}\psi_{\theta_{k}}(X_{k,i}^{(s)};\hat{\theta}_{k,s})$.
We then perform the local bias correction to $\\{\hat{\theta}_{k,s}\\}$ based
on data in subset $D_{k}^{s}$ to attain the debiased estimators
$\\{\hat{\theta}_{k,s}^{bc}\\}$. At last, two debiased weighted distributed
estimators of the form

$\tilde{\hat{\phi}}^{dWD}_{s}:=\\{\sum_{k=1}^{K}n_{k}\hat{H}_{k,s}(\hat{\theta}_{k,s})^{-1}\\}^{-1}\sum_{k=1}^{K}n_{k}(\hat{H}_{k,s}(\hat{\theta}_{k,s}))^{-1}\hat{\phi}_{k,2-|s-1|}^{bc}$

for $s=1,2$ are averaged to obtain the final debiased WD (dWD) estimator,
whose procedure is summarized in Algorithm 2. That the weight estimation and
the debiasing are conducted on different data splits remove the correlation,
and realize the gain of bias-correction procedure.

Input: $\\{X_{k,i},k=1,...,K;i=1,...,n_{k}\\}$

Output: $\hat{\phi}^{dWD}$

1 For each data block, split the data set into two non-overlapping equal-sized
subsets and denote those subsets as
$D_{k}^{s}=\\{X_{k,i}^{(s)}\\}_{i=1}^{n_{k}/2},s=1,2$ ;

2 Obtain the initial estimates
$\hat{\theta}_{k,s}=(\hat{\phi}_{k,s},\hat{\lambda}_{k,s})$ based on data from
$D_{k}^{s},s=1,2$ ;

3 Calculate $\hat{H}_{k,s}(\hat{\theta}_{k,s})$ in each block $(s=1,2)$, which
is the leading principal sub-matrix of order $p_{1}$ of
$(\nabla_{\theta_{k}}\hat{\Psi}_{\theta_{k}}^{(s)})^{-1}(2n_{k}^{-1}\sum_{i=1}^{n_{k}/2}\psi_{\theta_{k}}(X_{k,i}^{(s)};\hat{\theta}_{k,s})\psi_{\theta_{k}}(X_{k,i}^{(s)};\hat{\theta}_{k,s})^{T})(\nabla_{\theta_{k}}\hat{\Psi}_{\theta_{k}}^{(s)})^{-T}$
where
$\hat{\Psi}_{\theta_{k}}=2n_{k}^{-1}\sum_{i=1}^{n_{k}/2}\psi_{\theta_{k}}(X_{k,i}^{(s)};\hat{\theta}_{k,s})$;

4 Calculate the bias corrected estimators in each block
$(k=1,2,\cdots,K;s=1,2)$:
$\hat{\theta}_{k,s}^{bc}:=\hat{\theta}_{k,s}-2n_{k}^{-1}\hat{B}_{k,s}(\hat{\theta}_{k,s})1_{\mathcal{E}_{k,bc,s}}$
where
$\mathcal{E}_{k,bc,s}:=\\{\hat{\theta}_{k,s}-2n_{k}^{-1}\hat{B}_{k,s}(\hat{\theta}_{k,s})\in\Theta_{k}\\}$.
Denote the first $p_{1}$ dimensions of $\hat{\theta}_{k,s}^{bc}$ as
$\hat{\phi}_{k,s}^{bc}$;

5

6Send
$\\{\hat{\phi}_{k,s}^{bc},\hat{H}_{k,1}(\hat{\theta}_{k,s})^{-1},s=1,2\\}$ to
a central server and construct
$\tilde{\hat{\phi}}^{dWD}_{s}:=\\{\sum_{k=1}^{K}n_{k}\hat{H}_{k,s}(\hat{\theta}_{k,s})^{-1}\\}^{-1}\sum_{k=1}^{K}n_{k}(\hat{H}_{k,s}(\hat{\theta}_{k,s}))^{-1}\hat{\phi}_{k,2-|s-1|}^{bc}$;

7
$\hat{\phi}^{dWD}_{s}:=\tilde{\hat{\phi}}^{dWD}_{s}I(\tilde{\hat{\phi}}^{dWD}_{s}\in\Phi)+K^{-1}\sum_{k=1}^{K}n_{k}\hat{\phi}_{k,2-|s-1|}^{bc}I(\tilde{\hat{\phi}}^{dWD}_{s}\not\in\Phi)$
for $s=1,2$.;

8 $\hat{\phi}^{dWD}=\frac{1}{2}\sum_{s=1}^{2}\hat{\phi}_{s}^{dWD}.$

Algorithm 2 debiased Weighted Distributed (dWD) Estimator

To provide theoretical guarantee on the bias correction, we need an assumption
on the third derivative of the M-function $M$ (see [30]), which strengthens
part of Assumption 5.

######  Assumption 7.

(Strong smoothness) For each $x\in\mathbb{R}^{p}$, the third order derivatives
of $M(x;\theta_{k})$ with respect to $\theta_{k}$ exist and are $A(x)-$
Lipschitz continuous, i.e.

$\|(\nabla_{\theta_{k}}^{2}\psi_{\theta_{k}}(x;\theta_{k})-\nabla_{\theta_{k}}^{2}\psi_{\theta_{k}}(x;\theta_{k}^{{}^{\prime}}))(u\otimes
u)\|_{2}\leq A(x)\|\theta_{k}-\theta_{k}^{{}^{\prime}}\|_{2}\|u\|_{2}^{2},$

for all $\theta_{k},\theta_{k}^{{}^{\prime}}\in U_{k}$ defined in Assumption 5
and $u\in\mathbb{R}^{p}$, where $\mathbb{E}A(X_{k,i})^{2v}\leq A^{2v}$ for
some $v>0$ and $A<\infty$.

######  Theorem 4.

Under Assumptions 1 \- 4 and 6 \- 7, and Assumption 5 with $v,v_{1}\geq 4$ ,

$\mathbb{E}\|\hat{\phi}^{dWD}-\phi^{*}\|_{2}^{2}\leq\frac{C_{1}}{nK}+\frac{C_{2}}{n^{2}K}+\frac{C_{3}}{n^{3}}+\frac{C_{4}K}{n^{\bar{v}}},$

where $\bar{v}=min\\{v,\frac{v_{1}}{2}\\}$.

The main difference between the upper bounds in Theorem 4 and that of Theorem
2 for the WD estimator is the disappearance of the $\mathcal{O}(n^{-2})$ term
for the WD estimator, which has been dissolved and absorbed into the
$\mathcal{O}((n^{2}K)^{-1})$ and $\mathcal{O}(n^{-3})$ terms for the dWD
estimator. As shown next, this translates to more relaxed $K=o(n^{2})$ as
compared with $K=o(n)$ for the WD estimator in Theorem 3.

######  Theorem 5.

Under the conditions required by Theorem 4, if $K=o(n^{2})$,

$(\hat{\phi}^{dWD}-\phi^{*})^{T}\\{\sum_{k=1}^{K}n_{k}H_{k}(\theta_{k}^{*})^{-1}\\}(\hat{\phi}^{dWD}-\phi^{*})\overset{d}{\rightarrow}\chi_{p_{1}}^{2}.$

Note that the reason why Theorem 5 is formulated in the chi-squared
distribution form is the same as that when we formulate Theorem 3, and similar
confidence region with confidence level $1-\alpha$ can be constructed as

$\\{\phi|\big{(}\hat{\phi}^{dWD}-\phi)^{T}\\{\sum_{k=1}^{K}n_{k}H_{k}(\hat{\theta}_{k})^{-1}\\}\big{(}\hat{\phi}^{dWD}-\phi)\leq\chi_{p_{1},\alpha}^{2}\\}.$
(28)

The fact that the confidence regions of dWD and WD estimators use the same
standardizing matrix $\sum_{k=1}^{K}n_{k}H_{k}(\hat{\theta}_{k})^{-1}$
reflects that the dWD and WD estimators have the same estimation efficiency.
However, the debiased version has more relaxed constraint on $K=O(n^{2})$
(which is equivalent to $K=o(N^{2/3})$) than that of the WD estimator at
$K=o(n)$ ($K=o(\sqrt{N})$).

A more communication-efficient estimator of the common parameter can be
defined as the following debiased SaC (dSaC) estimator:

$\hat{\phi}^{dSaC}=N^{-1}\sum_{k=1}^{K}n_{k}(\hat{\phi}_{k}-n_{k}^{-1}\hat{B}_{k}^{1}(\hat{\theta}_{k})1_{\mathcal{E}_{k,bc}}),$
(29)

which only performs bias correction and may be preferable when the
heterogeneity is not large. The asymptotic property of the dSaC estimator is
summarized in the following proposition.

######  Theorem 6.

Under the conditions required by Theorem 4, if $K=o(n^{2})$,

$\displaystyle\mathbb{E}\|\hat{\phi}^{dSaC}-\phi^{*}\|_{2}^{2}\leq\frac{C_{1}}{nK}+\frac{C_{2}}{n^{2}K}+\frac{C_{3}}{n^{3}}\quad\hbox{
and }$ $\displaystyle
N^{2}(\hat{\phi}^{dSaC}-\phi^{*})^{T}\\{\sum_{k=1}^{K}n_{k}H_{k}(\theta_{k}^{*})\\}^{-1}(\hat{\phi}^{dSaC}-\phi^{*})\overset{d}{\rightarrow}\chi_{p_{1}}^{2}.$

The corresponding confidence region with confidence level $1-\alpha$ can be
constructed as

$\\{\phi|N^{2}\big{(}\hat{\phi}^{dSaC}-\phi)^{T}\\{\sum_{k=1}^{K}n_{k}H_{k}(\hat{\theta}_{k})\\}^{-1}\big{(}\hat{\phi}^{dSaC}-\phi)\leq\chi_{p_{1},\alpha}^{2}\\}.$
(30)

It is noted that the dSaC and SaC estimators have the same asymptotic
distribution. Hence, the confidence regions based on the SaC estimator can be
constructed as (30) with $\hat{\phi}^{dSaC}$ replaced by $\hat{\phi}^{SaC}$.

To compare with the subsampled average mixture method (SAVGM) estimator
proposed in [30] which also performs local bias correction but under the
homogeneous setting, we have the following corollary to Theorem 6.

######  Corollary 1.

Under the homogeneous case such that $\\{X_{k,i},k=1,...,K,i=1,...,n;\\}$ are
IID distributed, and the assumptions required by Theorem 4,

$\mathbb{E}\|\hat{\theta}^{dSaC}-\theta_{1}^{*}\|_{2}^{2}\leq\frac{2\mathbb{E}\|\nabla_{\theta_{1}}\Psi_{\theta}(\theta_{1}^{*})^{-1}\psi_{\theta_{1}}(X_{1,1};\theta_{1}^{*})\|_{2}^{2}}{nK}+\frac{C_{1}}{n^{2}K}+\frac{C_{2}}{n^{3}},$
(31)

where $\theta_{1}^{*}$ is the true parameter for all the $K$ data blocks.

The SAVGM estimator resamples $\lfloor rn_{k}\rfloor$ data points from each
data block $k$ for a $r\in(0,1)$ to obtain a local estimator
$\hat{\theta}^{SaC}_{k,r}$ based on the sub-samples. Then, the SAVGM estimator
is

$\bar{\theta}_{SAVGM}=\frac{\hat{\theta}_{k}^{SaC}-r\hat{\theta}^{SaC}_{k,r}}{1-r},$
(32)

whose MSE bound as given in Theorem 4 of [30] is

$\mathbb{E}\|\bar{\theta}_{SAVGM}-\theta_{1}^{*}\|_{2}^{2}\leq\frac{2+3r}{(1-r)^{2}}\frac{\mathbb{E}\|\nabla_{\theta_{1}}\Psi_{\theta}(\theta_{1}^{*})^{-1}\psi_{\theta_{1}}(X_{1,1};\theta_{1}^{*})\|_{2}^{2}}{nK}+\frac{C_{1}}{n^{2}K}+\frac{C_{2}}{n^{3}}.$
(33)

Thus, the MSE bound (33) of the SAVGM estimator has an inflated factor
$\frac{2+3r}{2(1-r)^{2}}>1$ for $r\in(0,1)$, when compared with that of the
dSaC estimator, although it is computationally more efficient than the dSaC
and dWD estimators as it only draws one subsample in its resampling. For more
comparisons between the dSaC estimator and one-step estimators proposed by
Huang and Huo (2019) [10], see Section 1.10 in SM.

##  6 Simulation Results

We report results from simulation experiments designed to verify two sets
theoretical findings made in the previous sections. One was to confirm the
finding in Section 2 that the full sample estimator $\hat{\phi}_{full}$ is not
necessarily more efficient than the SaC estimator $\hat{\phi}^{SaC}$. The
other was to evaluate the numerical performance of the newly proposed weighted
distributive (WD) , debaised SaC (dSaC) and debiased WD (dWD) estimators of
the common parameter and compare them with the existing SaC and subsampled
average mixture method (SAVGM) (with subsampling rate $r=0.05$) estimators.
Although the SAVGM estimator [30] was proposed under the homogeneous setting,
but since its main bias correction is performed locally on each data block $k$
as shown in (32), similar theoretical bounds as formula (33) can be derived
without much modifications on the original proof. Throughout the simulation
experiments, the results of each simulation setting were based on $B=500$
number of replications and were conducted in R paralleled with a single
10-core Intel(R) Core(TM) i9-10900K @3.7 GHz processor.

In the first simulation experiment, we simulated the errors-in-variables Model
(13) with the objective function (14) to compare the performance of the full
sample, the SaC and the WD estimators: $\hat{\phi}_{full}$, $\hat{\phi}^{SaC}$
and $\hat{\phi}^{WD}$. The simulation was carried out by first generating IID
$\\{Z_{i,k}\\}$ from $\mathcal{N}(\mu_{Z},\sigma^{2}_{Z})$, and then upon
given a $Z_{i,k}$, $(X_{k,i},Y_{i,k})^{T}$ were independently drawn from
$\mathcal{N}\big{(}(Z_{i,k},\phi^{*}+\lambda_{k}^{*}Z_{i,k})^{T},\sigma^{2}I_{2\times
2}\big{)}$. We chose $\phi^{*}=1,K=2,\sigma^{2}=1$ and $n_{1}=n_{2}=5\times
10^{4}=N/2$, and $\lambda_{1}^{*},\lambda_{2}^{*},\mu_{Z}$ and
$\sigma^{2}_{Z}$ were those reported in Table 1 under four scenarios.

As discussed in Section 2, the relative efficiency of $\hat{\phi}_{full}$ to
$\hat{\phi}^{SaC}$ depends on the ratio
$\sigma^{2}(\mathbb{E}Z)^{2}/(var(Z)\mathbb{E}Z^{2})$ as shown in (15). We
designed four scenarios according to the above ratio under
$\lambda_{1}^{*}\not=\lambda_{2}^{*}$ and $\mathbb{E}Z\not=0$, respectively,
which represented the settings where the full sample estimator
$\hat{\phi}_{full}$ would be less (Scenario 1) or more (Scenario 2) efficient
than the SaC estimator as predicted by the ratio, but not as efficient as the
weighted distributed estimator $\hat{\phi}^{WD}$. Scenario 3
($\lambda_{1}^{*}\not=\lambda_{2}^{*},\mathbb{E}Z=0$) was the case when
$\hat{\phi}_{full}$ and $\hat{\phi}^{WD}$ would be asymptotically equivalent,
and both estimators would be more efficient than $\hat{\phi}^{SaC}$. Scenario
4 was the homogeneous case with $\lambda_{1}^{*}=\lambda_{2}^{*}$ in which all
the three estimators would have the same asymptotic efficiency. For all the
four scenarios, the ARE column of the Table 1 confirmed the relative
efficiency as predicted by the asymptotic variances in (15), and was well
reflected in the comparison of the RMSEs, as the bias is of smaller order as
compared with that of the SD and thus negligible.

Table 1: Average root mean squared error (RMSE) and the standard deviation (SD), multiplied by $10^{2}$, of the full sample estimator $\hat{\phi}_{full}$, the SaC estimator $\hat{\phi}^{SaC}$ and the WD estimator $\hat{\phi}^{WD}$ under four scenarios for the errors-in-variables model (12) for $N=10^{5},K=2$ and $n_{1}=n_{2}$. AREs (asymptotic relative efficiency) of $\hat{\phi}_{full}$ to $\hat{\phi}^{SaC}$ are calculated from (15).  |  |  | $\hat{\phi}_{full}$ | $\hat{\phi}^{SaC}$ | $\hat{\phi}^{WD}$
---|---|---|---|---|---
Scenario | $(\lambda_{1}^{*},\lambda_{2}^{*})$ | ARE | RMSE | SD | RMSE | SD | RMSE | SD
Scenario 1 | (0.25,3.25) | 0.89 | 4.55 | 4.51 | 4.12 | 4.09 | 3.91 | 3.89
$(\mu_{Z}=1,\sigma_{Z}^{2}=0.1)$ | (0.5,3.5) | 0.93 | 4.65 | 4.65 | 4.35 | 4.35 | 4.08 | 4.08
| (0.75,3.75) | 0.97 | 4.52 | 4.52 | 4.40 | 4.38 | 4.13 | 4.13
Scenario 2 | (0.25,2.25) | 1.18 | 2.95 | 2.95 | 3.24 | 3.24 | 2.89 | 2.89
$(\mu_{Z}=3,\sigma_{Z}^{2}=0.5)$ | (0.75,2.75) | 1.28 | 3.28 | 3.26 | 3.65 | 3.64 | 3.17 | 3.16
| (1.25,3.25) | 1.31 | 3.71 | 3.71 | 4.16 | 4.07 | 3.64 | 3.61
Scenario 3 | (0.25,2.25) | 1.97 | 0.41 | 0.41 | 0.61 | 0.61 | 0.41 | 0.41
$(\mu_{Z}=0,\sigma_{Z}^{2}=0.5)$ | (0.75,2.75) | 1.92 | 0.51 | 0.51 | 0.70 | 0.70 | 0.51 | 0.51
| (1.25,3.25) | 1.68 | 0.64 | 0.64 | 0.82 | 0.82 | 0.64 | 0.64
Scenario 4 | (0.5,0.5) | 1 | 3.25 | 3.24 | 3.31 | 3.28 | 3.30 | 3.26
$(\mu_{Z}=4,\sigma_{Z}^{2}=0.5)$ | (1.0,1.0) | 1 | 3.53 | 3.53 | 3.59 | 3.59 | 3.59 | 3.59
| (1.5,1.5) | 1 | 4.06 | 4.03 | 4.08 | 4.07 | 4.06 | 4.06

In the second simulation experiment, we evaluated the numerical performance of
the five estimators for the common parameter $\phi$ under a logistic
regression model. For each of $K$ data block with
$K\in\\{10,50,100,250,500,1000,2000\\}$,
$\\{(X_{k,i};Y_{k,i})\\}_{i=1}^{n}\subset\mathbb{R}^{p}\times\\{0,1\\}$ were
independently sampled from the following model:

$X_{k,i}\overset{i.i.d}{\sim}\mathcal{N}(\mathbf{0}_{p\times
1},0.75^{2}I_{p\times p})\quad\hbox{and}\quad
P(Y_{k,i}=1|X_{k,i})=\frac{exp(X_{k,i}^{T}\theta_{k}^{*})}{1+exp(X_{k,i}^{T}\theta_{k}^{*})},$

where $\theta_{k}^{*}=(\phi^{*},\lambda_{k}^{*T})^{T}$, $\phi^{*}=1$,
$\lambda_{k}^{*}=(\lambda_{k,1}^{*},\lambda_{k,2}^{*},\cdots,\lambda_{k,p_{2}}^{*})^{T}$
and $\lambda_{k,j}^{*}=(-1)^{j}10(1-\frac{2(k-1)}{K-1})$. The sample sizes of
the data blocks were equal at $n=NK^{-1}$ with $N=2\times 10^{6}$. Two levels
of the dimension $p_{2}=4$ and $10$ of the nuisance parameter $\lambda_{k}$
were considered. A derivation of the bias correction formula for the logistic
model is given in Section 1.9 of the SM.

Figure 1 reports the root mean square errors (RMSEs) and absolute bias of the
estimators. It is observed that the weighted distributed estimator WD, and the
two debiased estimators dSaC and dWD had smaller RMSE than those of the SaC
and SAVGM for almost all the simulation settings. Between the SaC and SAVGM,
the SAVGM fared better in the lower dimensional case of $p_{2}=4$, but was
another way around for $p_{2}=10$. It was evidence that the WD estimator had
much smaller RMSEs than the SaC and SAVGM estimators for all the block number
$K$, realizing its theoretical promises. In most cases the WD estimator had
smaller bias than the SaC estimator although it was not debiased. It also had
smaller RMSEs than the debiased SaC estimator dSaC for almost all cases of the
block numbers for $p_{2}=4$, while in the higher dimensional $p_{2}=10$ the WD
estimator was advantageous for $K\leq 250$. The latter indicated the need for
conducting the bias correction to the WD estimator. Both bias corrected dWD
and dSaC were very effective in reducing the bias of the WD and SaC
estimators, respectively, especially for larger $K$ when the bias was more
severe. The debiased WD attained the smallest RMSEs and the bias in all
settings, suggesting the need for conducting both weighting and the bias
correction in the distributed inference especially for large $K$. These
empirical results were consistent with Theorems 2 and 4, namely the leading
RMSE term of the WD estimator changes from $\mathcal{O}((Kn)^{-1})$ to
$\mathcal{O}(n^{-2})$ when $K$ surpasses the local sample size $n$, while the
leading RMSEs of the dWD is still $\mathcal{O}((nK)^{-1})$ until $K>>n^{2}$.

We also evaluated the coverage probabilities and widths of the $1-\alpha$
($\alpha=0.01,0.05,0.1$) confidence intervals (CIs) of the common parameter
based on the asymptotic normality as given after Theorems 3 and 5. The SAVGM
estimator was not included as its asymptotic distribution was not made
available in [30]. Table 2 reports the empirical coverage and the average
width of the CIs. It is observed that for the lower dimensional nuisance
parameter case of $p_{2}=4$ the four types of the CIs all had quite adequate
coverage levels when $K\leq 100$. However, for $K\geq 250$, the SaC CIs first
started to lose coverage, followed by those of the WD, while the CIs of the
debiased SaC (dSaC) and debiased weighted distributed (dWD) estimators can
hold up to the promised coverage for all cases of $K$. The outstanding
performance of the dSaC and dWD CIs was largely replicated for the higher
dimensional nuisance parameter case of $p_{2}=10$, while the other two non-
debiased estimator based CIs had their coverage quickly slipped below the
nominal coverage levels. Although the dSaC CIs had comparable coverages with
the dWD CIs, their widths were much wider than those of the dWD. This was
largely due to the fact that the weighted averaging conducted in the weighted
distributed estimation reduced the variation and hence the width of the CIs.
The widths of the WD CIs were largely the same with those of the dWD, and yet
the coverage levels of the dWD CIs were much more accurate indicating the
importance of the bias correction as it shifted the CIs without inflating the
width.

In addition to the simulation experiments on the statistical properties of the
estimators, the computation efficiency of the estimators was also evaluated.
Table 3 reports the average CPU time per simulation run based on 500
replications of the five estimators for a range of $K$ and dimension $p_{2}$
of the nuisance parameter for the logistic regression model with the total
sample size $N=2\times 10^{6}$.

The computation speed of the dSaC and dWD estimators were relatively slower
than those of the SaC, WD and SAVGM estimators. The WD estimator was quite
fast, which means that the re-weighting used less computing time than the
bias-reduction. In comparison, the dWD estimator was the slowest as a cost for
attaining the best RMSE among the five estimators in all settings. It is
observed in Table 3 that the overall computation time for each estimator first
decreased and then increased as $K$ became larger. The decrease in time was
because the benefit of the distributed computation, while the increase was due
to the increase in the number of optimization associated with the M-estimation
performed as $K$ got larger. However, it is worth mentioning that these
results did not account for the potential time expenditure in data
communication among different data blocks.

##  7 Discussion

This paper investigates several distributed M-estimators in the presence of
heterogeneous distributions among the data blocks. The weighted distributed
(WD) estimator is able to improve the estimation efficiency of the ”Split-And-
Conquer” (SaC) estimator for the common parameter. Two debiased estimators (
dWD and dSaC) are proposed to allow for larger numbers of data blocks $K$. The
statistical properties of these three estimators are shown to be advantageous
over the SaC and SAVGM estimators. In particular, the WD estimator has good
performance for smaller $K$ relative to $n$, and the debiased WD estimator
that conducted both bias correction and weighting offers good estimation
accuracy for large $K$.

An important issue for the distributed estimation is the size of $K$ relative
to the local average sample size $n$. This is especially true in Federated
Learning setting where the number of clients (data blocks) are usually very
large. Both SaC and WD estimators require $K=o(\sqrt{N})$ to preserve the
$\mathcal{O}(N^{-1})$ convergence rate for its MSE and the $\sqrt{N}$ rate for
the asymptotic variance. The debiased dWD and dSaC relax the restriction to
$K=o(N^{2/3})$ without compromising the convergence rate. The dSaC may be used
as a computationally cheaper version of the dWD at the cost of larger
variations and wider confidence regions when compared with dWD.

## References

  * Bartlett, [1953] Bartlett, M. (1953).  Approximate confidence intervals.  Biometrika, 40:12–19.
  * Battey et al., [2018] Battey, H., Fan, J., Liu, H., Lu, J., and Zhu, Z. (2018).  Distributed testing and estimation under sparse high dimensional models.  The Annals of Statistics, 46:1352–1382.
  * Chen and Peng, [2021] Chen, S. X. and Peng, L. (2021).  Distributed statistical inference for massive data.  The Annals of Statistics, 49:2851–2869.
  * Chen et al., [2019] Chen, X., Liu, W., and Zhang, Y. (2019).  Quantile regression under memory constraint.  The Annals of Statistics, 47:3244–3273.
  * Chen and Xie, [2014] Chen, X. and Xie, M. (2014).  A split-and-conquer approach for analysis of extraordinarily large data.  Statistica Sinica, 24:1655–1684.
  * Duan et al., [2021] Duan, R., Ning, Y., and Chen, Y. (2021).  Heterogeneity-aware and communication-efficient distributed statistical inference.  Biometrika, to appear.
  * Evgeniou and Pontil, [2004] Evgeniou, T. and Pontil, M. (2004).  Regularized multi–task learning.  KDD-2004 - Proceedings of the Tenth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 109–117.
  * Fuller, [1987] Fuller, W. (1987).  Measurement Error Models.  Wiley.
  * Hansen, [1982] Hansen, L. (1982).  Large sample properties generalized method of moments estimators.  Econometrica, 50:1029–1054.
  * Huang and Huo, [2019] Huang, C. and Huo, X. (2019).  A distributed one-step estimator.  Mathematical Programming, 174:41–76.
  * Jordan et al., [2019] Jordan, M., Lee, J., and Yang, Y. (2019).  Communication-efficient distributed statistical learning.  Journal of the American Statistical Association, 114:668–681.
  * Kairouz et al., [2021] Kairouz, P., McMahan, H., Avent, B., Bellet, A., Bennis, M., Bhagoji, A., Bonawitz, K., Charles, Z., Cormode, G., Cummings, R., D’Oliveira, R., Eichner, H., El Rouayheb, S., Evans, D., Gardner, J., Garrett, Z., Gascón, A., Ghazi, B., Gibbons, P., and Zhao, S. (2021).  Advances and open problems in federated learning.  Foundations and Trends® in Machine Learning, 14:1–210.
  * Kleiner et al., [2011] Kleiner, A., Talwalkar, A., Sarkar, P., and Jordan, M. (2011).  A scalable bootstrap for massive data.  Journal of the Royal Statistical Society Series B (Statistical Methodology), 76:795–816.
  * Lai and Wang, [1993] Lai, T. and Wang, J. (1993).  Edgeworth expansions for symmetric statistics with applications to bootstrap methods.  Statistica Sinica, 3:517–542.
  * Li et al., [2003] Li, H., Lindsay, B., and Waterman, R. (2003).  Efficiency of projected score methods in rectangular array asymptotics.  Journal of the Royal Statistical Society Series B, 65:191–208.
  * Li et al., [2020] Li, T., Sahu, A., Talwalkar, A., and Smith, V. (2020).  Federated learning: Challenges, methods, and future directions.  IEEE Signal Processing Magazine, 37:50–60.
  * Lin and Xi, [2010] Lin, N. and Xi, R. (2010).  Fast surrogates of U-statistics.  Computational Statistics & Data Analysis, 54:16–24.
  * McCullagh, [1983] McCullagh, P. (1983).  Quasi-likelihood functions.  The Annals of Statistics, 11:59–67.
  * McMahan et al., [2017] McMahan, B., Moore, E., Ramage, D., Hampson, S., and Arcas, B. A. y. (2017).  Communication-Efficient Learning of Deep Networks from Decentralized Data.  Proceedings of Machine Learning Research, 54:1273–1282.
  * Reiersol, [1950] Reiersol, O. (1950).  Identifiability of a linear relation between variables which are subject to error.  Econometrica, 18:375–389.
  * Rilstone et al., [1996] Rilstone, P., Srivastava, V., and Ullah, A. (1996).  The second-order bias and mean squared error of nonlinear estimators.  Journal of Econometrics, 124:369–395.
  * Schafer and Purdy, [1996] Schafer, D. and Purdy, K. (1996).  Likelihood analysis for error-in-variables regression with replicate measurements.  Biometrika, 83:813–824.
  * Sengupta et al., [2015] Sengupta, S., Volgushev, S., and Shao, X. (2015).  A subsampled double bootstrap for massive data.  Journal of the American Statistical Association, 111:1222–1232.
  * Smith et al., [2017] Smith, V., Chiang, C.-K., Sanjabi, M., and Talwalkar, A. (2017).  Federated Multi-Task Learning.  Advances in Neural Information Processing Systems.
  * Stefanski and Boos, [2002] Stefanski, L. and Boos, D. (2002).  The Calculus of M-Estimation.  The American Statistician, 56:29–38.
  * van der Vaart, [1999] van der Vaart, A. (1999).  Asymptotic Statistics.  Cambridge University Press.
  * Volgushev et al., [2017] Volgushev, S., Chao, S.-K., and Cheng, G. (2017).  Distributed inference for quantile regression processes.  Annals of Statistics, 47.
  * Yang et al., [2019] Yang, Q., Liu, Y., Chen, T., and Tong, Y. (2019).  Federated machine learning: Concept and applications.  ACM Transactions on Intelligent Systems and Technology, 10:1–19.
  * Yaron et al., [1996] Yaron, A., Hansen, L., and Heaton, J. (1996).  Finite-Sample Properties of Some Alternative GMM Estimators.  Journal of Business & Economic Statistics, 14:262–80.
  * Zhang et al., [2013] Zhang, Y., Duchi, J., and Wainwright, M. (2013).  Comunication-efficient algorithms for statistical optimization.  Journal of Machine Learning Research, 14:3321–3363.
  * Zhang and Yang, [2021] Zhang, Y. and Yang, Q. (2021).  A Survey on Multi-Task Learning.  IEEE Transactions on Knowledge and Data Engineering.
  * Zhao et al., [2014] Zhao, T., Cheng, G., and Liu, H. (2014).  A partially linear framework for massive heterogeneous data.  The Annals of Statistics, 44:1400–1437.

(a) Absolute Bias ($p_{2}=4$) | (b) RMSE ($p_{2}=4$)
---|---
|
(c) Absolute Bias ($p_{2}=10$) | (d) RMSE ($p_{2}=10$)
|
Figure 1: Average simulated bias (a, c) and the root mean square errors (RMSE)
(b,d) of the weighted distributed (WD) (red circle), the SaC (blue triangle),
the debiased SaC (dSaC) (green square), the debiased WD (dWD) (purple cross),
the subsampled average mixture SAVGM (pink square cross) estimators, with
respect to the number of data block $K$ for the logistic regression model with
the dimension $p_{2}$ of the nuisance parameter $\lambda_{k}$ being $4$ and
$10$, respectively with the full sample size $N=2\times 10^{6}$.  Table 2:
Coverage probabilities and widths (in parentheses, multiplied by $100$) of the
$1-\alpha$ confidence intervals for the common parameter $\phi$ in the
logistic regression model based on the asymptotic normality of the SaC, the
WD, the debiased SaC and the debiased WD estimators with respect to the number
of data blocks $K$ for two dimensions of the heterogeneous parameter $p_{2}$
with the full sample size $N=2\times 10^{6}$.

(a) $p_{2}=4$

K |  | SaC | WD | dSaC | dWD
---|---|---|---|---|---
| $1-\alpha$ | 0.99 | 0.95 | 0.90 | 0.99 | 0.95 | 0.90 | 0.99 | 0.95 | 0.90 | 0.99 | 0.95 | 0.90
10 |  | 0.99 | 0.96 | 0.92 | 0.99 | 0.97 | 0.91 | 0.99 | 0.96 | 0.92 | 0.99 | 0.96 | 0.91
|  | (2.45) | (1.87) | (1.57) | (2.03) | (1.55) | (1.30) | (2.45) | (1.87) | (1.57) | (2.03) | (1.55) | (1.30)
50 |  | 0.99 | 0.95 | 0.91 | 0.98 | 0.93 | 0.89 | 0.99 | 0.95 | 0.91 | 0.99 | 0.93 | 0.88
|  | (2.36) | (1.80) | (1.51) | (1.97) | (1.50) | (1.26) | (2.36) | (1.80) | (1.51) | (1.97) | (1.50) | (1.26)
100 |  | 0.98 | 0.94 | 0.91 | 0.99 | 0.95 | 0.91 | 0.99 | 0.95 | 0.91 | 0.99 | 0.95 | 0.91
|  | (2.36) | (1.79) | (1.51) | (1.96) | (1.49) | (1.25) | (2.36) | (1.79) | (1.51) | (1.96) | (1.49) | (1.25)
250 |  | 0.99 | 0.93 | 0.85 | 0.99 | 0.95 | 0.90 | 0.99 | 0.96 | 0.91 | 0.99 | 0.95 | 0.90
|  | (2.36) | (1.79) | (1.50) | (1.96) | (1.49) | (1.25) | (2.36) | (1.79) | (1.50) | (1.96) | (1.49) | (1.25)
500 |  | 0.91 | 0.77 | 0.66 | 0.99 | 0.95 | 0.88 | 0.99 | 0.96 | 0.90 | 0.99 | 0.95 | 0.89
|  | (2.36) | (1.80) | (1.51) | (1.96) | (1.49) | (1.25) | (2.36) | (1.80) | (1.51) | (1.96) | (1.49) | (1.25)
1000 |  | 0.65 | 0.41 | 0.28 | 0.99 | 0.94 | 0.88 | 0.99 | 0.94 | 0.88 | 0.99 | 0.93 | 0.88
|  | (2.38) | (1.81) | (1.52) | (1.96) | (1.49) | (1.25) | (2.38) | (1.81) | (1.52) | (1.97) | (1.50) | (1.25)
2000 |  | 0.01 | 0.01 | 0.00 | 0.99 | 0.91 | 0.81 | 0.98 | 0.94 | 0.88 | 0.99 | 0.94 | 0.90
|  | (2.42) | (1.84) | (1.55) | (1.96) | (1.50) | (1.25) | (2.42) | (1.84) | (1.55) | (1.98) | (1.50) | (1.26)

(b) $p_{2}=10$

K |  | SaC | WD | dSaC | dWD
---|---|---|---|---|---
| $1-\alpha$ | 0.99 | 0.95 | 0.90 | 0.99 | 0.95 | 0.90 | 0.99 | 0.95 | 0.90 | 0.99 | 0.95 | 0.90
10 |  | 0.99 | 0.94 | 0.88 | 1.00 | 0.96 | 0.92 | 1.00 | 0.94 | 0.88 | 1.00 | 0.96 | 0.92
|  | (3.05) | (2.32) | (1.95) | (2.41) | (1.84) | (1.54) | (3.05) | (2.32) | (1.95) | (2.42) | (1.84) | (1.54)
50 |  | 0.99 | 0.93 | 0.87 | 0.99 | 0.95 | 0.88 | 0.98 | 0.94 | 0.88 | 0.99 | 0.96 | 0.88
|  | (2.94) | (2.24) | (1.88) | (2.29) | (1.74) | (1.46) | (2.94) | (2.24) | (1.88) | (2.29) | (1.74) | (1.46)
100 |  | 0.97 | 0.89 | 0.84 | 0.97 | 0.93 | 0.87 | 0.98 | 0.95 | 0.90 | 0.98 | 0.94 | 0.89
|  | (2.93) | (2.23) | (1.87) | (2.28) | (1.74) | (1.46) | (2.93) | (2.23) | (1.87) | (2.29) | (1.74) | (1.46)
250 |  | 0.89 | 0.72 | 0.63 | 0.98 | 0.92 | 0.87 | 1.00 | 0.97 | 0.90 | 1.00 | 0.96 | 0.90
|  | (2.94) | (2.24) | (1.88) | (2.28) | (1.74) | (1.46) | (2.94) | (2.24) | (1.88) | (2.29) | (1.74) | (1.46)
500 |  | 0.51 | 0.28 | 0.18 | 0.93 | 0.81 | 0.70 | 0.99 | 0.94 | 0.90 | 0.98 | 0.94 | 0.88
|  | (2.97) | (2.26) | (1.90) | (2.29) | (1.74) | (1.46) | (2.97) | (2.26) | (1.90) | (2.30) | (1.75) | (1.47)
1000 |  | 0.00 | 0.00 | 0.00 | 0.66 | 0.37 | 0.28 | 0.99 | 0.95 | 0.90 | 0.99 | 0.96 | 0.89
|  | (3.04) | (2.31) | (1.94) | (2.30) | (1.75) | (1.47) | (3.04) | (2.31) | (1.94) | (2.34) | (1.78) | (1.49)
2000 |  | 0.00 | 0.00 | 0.00 | 0.02 | 0.00 | 0.00 | 0.99 | 0.96 | 0.90 | 0.99 | 0.93 | 0.87
|  | (3.22) | (2.45) | (2.06) | (2.34) | (1.78) | (1.49) | (3.22) | (2.45) | (2.06) | (2.40) | (1.82) | (1.53)

Table 3: Average CPU time for each replication based on $B=500$ replications for the SaC, the SAVGM, the WD, the debiased SaC and the debiased WD estimators for the logistic regression model with respect to $K$ and the dimension $p_{2}$ of the nuisance parameter. Total sample size $N=2\times 10^{6}$.  K | $SaC$ | $SAVGM$ | $WD$ | $dSaC$ | $dWD$
---|---|---|---|---|---
|  | $p_{2}=4$ |
10 | 15.65 | 15.97 | 18.50 | 20.00 | 21.95
50 | 9.63 | 9.95 | 10.66 | 12.37 | 14.59
100 | 8.09 | 8.63 | 8.76 | 10.50 | 12.05
250 | 8.49 | 9.69 | 9.07 | 10.84 | 12.82
500 | 9.68 | 11.58 | 10.25 | 11.97 | 14.84
1000 | 11.67 | 13.81 | 12.32 | 13.93 | 19.08
2000 | 15.78 | 19.68 | 16.57 | 18.11 | 28.55
|  | $p_{2}=10$ |
10 | 34.60 | 35.19 | 43.84 | 50.47 | 55.35
50 | 20.13 | 20.18 | 24.16 | 29.99 | 33.69
100 | 15.60 | 16.20 | 17.74 | 23.63 | 24.47
250 | 10.77 | 12.61 | 11.88 | 18.22 | 20.39
500 | 11.55 | 14.50 | 12.56 | 18.80 | 23.73
1000 | 15.23 | 18.27 | 16.28 | 22.38 | 32.24
2000 | 23.42 | 27.99 | 24.62 | 30.43 | 48.05