# Variance Reduced Median-of-Means Estimator for Byzantine-Robust Distributed
Inference

Jiyuan Tu 111Shanghai Jiao Tong University, Shanghai, China, Email:
<EMAIL_ADDRESS>Weidong Liu 222Shanghai Jiao Tong University, Shanghai,
China, Email<EMAIL_ADDRESS>Xiaojun Mao 333Fundan University, Shanghai,
China, Email<EMAIL_ADDRESS>Xi Chen 444New York University, New York,
USA, Email<EMAIL_ADDRESS>Jiyuan Tu 555Shanghai Jiao Tong University,
Shanghai, China, Email<EMAIL_ADDRESS>Weidong Liu 666Shanghai Jiao Tong
University, Shanghai, China, Email<EMAIL_ADDRESS>Xiaojun Mao 777Fundan
University, Shanghai, China, Email<EMAIL_ADDRESS>Xi Chen 888New York
University, New York, USA, Email<EMAIL_ADDRESS>
###### Abstract

This paper develops an efficient distributed inference algorithm, which is
robust against a moderate fraction of Byzantine nodes, namely arbitrary and
possibly adversarial machines in a distributed learning system. In robust
statistics, the median-of-means (MOM) has been a popular approach to hedge
against Byzantine failures due to its ease of implementation and computational
efficiency. However, the MOM estimator has the shortcoming in terms of
statistical efficiency. The first main contribution of the paper is to propose
a variance reduced median-of-means (VRMOM) estimator, which improves the
statistical efficiency over the vanilla MOM estimator and is computationally
as efficient as the MOM. Based on the proposed VRMOM estimator, we develop a
general distributed inference algorithm that is robust against Byzantine
failures. Theoretically, our distributed algorithm achieves a fast convergence
rate with only a constant number of rounds of communications. We also provide
the asymptotic normality result for the purpose of statistical inference. To
the best of our knowledge, this is the first normality result in the setting
of Byzantine-robust distributed learning. The simulation results are also
presented to illustrate the effectiveness of our method.

keywords: Byzantine robustness, distributed inference, median-of-means,
statistical efficiency

##  1 Introduction

Due to the rapid increase of the scale of data, modern datasets are usually
too large to fit in a single device, and thus have to be stored and processed
in a distributed manner. In a common distributed computing environment, data
are stored across multiple machines/nodes. A single master node is in charge
of maintaining and updating target parameters, and a large number of worker
machines perform local computations and communicate the computed information
with the master node (see Figure 2 in Section 3 for an illustration). As
compared to the traditional single machine setting, where the entire data can
be loaded into the memory for the centralized computation, the distributed
setting poses two major challenges.

The first challenge comes from the tradeoff between communication cost and
statistical accuracy. For example, one-shot communication (e.g., taking
average of local estimators), though incurs low communication cost, has a poor
performance for nonlinear estimation when the number of machines is large
(see, e.g., Li et al. (2013); Zhang et al. (2013, 2015); Zhao et al. (2016);
Rosenblatt and Nadler (2016); Shang and Cheng (2017); Lee et al. (2017)).
Therefore, iterative approaches are adopted in literature (see, e.g., Shamir
et al. (2014); Jordan et al. (2019); Chen et al. (2019); Fan et al. (2019);
Wang et al. (2019); Chen et al. (2020)). For iterative algorithms, since each
iteration of communication requires synchronization, a communicationally
efficient algorithm should run with a small number of iterations. Our goal is
to develop algorithms that achieve communication efficiency without losing
statistical accuracy.

The second challenge comes from the vulnerability of worker machines and
communication channels. In particular, the information sent from a worker
machine can be arbitrarily erroneous due to hardware or software breakdowns,
data crashes, or communication failures. Such an error is usually referred to
as Byzantine failures (Lamport et al., 1982). In other words, a subset of
workers called Byzantine machines, may send arbitrary and even adversarial
messages to the master. Distributed learning under the Byzantine setting has
attracted a lot of research attentions in recent years (see, e.g., Feng et al.
(2014); Chen et al. (2017); Blanchard et al. (2017); Xie et al. (2018);
Alistarh et al. (2018); Yin et al. (2018, 2019); Su and Xu (2019)). However,
as we will survey later, some of these methods suffer from a larger number of
iterations of communications and existing analysis only focuses on the
convergence rate. The statistical inference with Byzantine failures, which
plays an important role in uncertainty quantification, is still largely open.

The goal of this paper is to propose a communication-efficient statistical
inference method, which is robustly against Byzantine failures. We consider a
general risk minimization problem,

$\boldsymbol{\theta}^{*}=\mathop{\mathrm{argmin}}_{\boldsymbol{\theta}\in{\mathbb{R}}^{p}}{\mathbb{E}}_{X\sim{\mathfrak{X}}}\left\\{f(X,\boldsymbol{\theta})\right\\},$
(1)

where $f$ is the convex loss function, $X$ denotes the random sample from a
probability distribution ${\mathfrak{X}}$, and
$\boldsymbol{\theta}^{*}\in\mathbb{R}^{p}$ is the target parameter vector of
interest. To infer the underlying parameter $\boldsymbol{\theta}^{*}$, assume
that $N$ i.i.d. observations $\\{X_{1},...,X_{N}\\}$ are collected and evenly
distributed over $(m+1)$-machines
$\\{{\mathcal{H}}_{0},\dots\mathcal{H}_{m}\\}$, where each machine contains
$n$ observations. We allow diverging $N$, $n$, and $p$ under certain rate
constraints.

In this paper, we consider the Byzantine distributed framework, which allows
for Byzantine failures described as follows. In particular, we assume there
exists an $\alpha_{n}$ fraction of worker machines (a.k.a. Byzantine
machines), whose indices form a subset ${\mathcal{B}}\subseteq\\{1,\dots,m\\}$
with ${\mathrm{C}ard}({\mathcal{B}})=\lfloor\alpha_{n}m\rfloor$. The Byzantine
machines are subject to the following Byzantine failures when communicating
information:

######  Definition 1 (Byzantine Failures).

In each round of communication, assume the information produced by each
machine is $\boldsymbol{v}_{j}$, then the actual information $\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{v}$\kern-1.00006pt}}}\,_{j}$
received from each worker machine ${\mathcal{H}}_{j}$ is as follows

$\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{v}$\kern-1.00006pt}}}\,_{j}=\begin{cases}\boldsymbol{v}_{j}\quad&j\notin{\mathcal{B}},\\\
*\quad&j\in{\mathcal{B}},\end{cases}$

where $*$ denotes an arbitrary value.

Let us start with the most fundamental setting where $\boldsymbol{\theta}^{*}$
is the population mean and the goal is to infer the population mean with the
presence of Byzantine failures. A widely used robust estimator is the median-
of-means (MOM) estimator (Nemirovsky and Yudin, 1983; Jerrum et al., 1986;
Alon et al., 1999), which first computes the local sample mean on each machine
and then aggregate them by taking the median. Due its ease of implementation
and computational efficiency, the MOM estimator has attracted a lot of
attentions (Minsker, 2015; Hsu and Sabato, 2016; Lecué and Lerasle, 2020;
Lugosi and Mendelson, 2019; Minsker, 2019) and served as an important building
block in distributed learning with Byzantine failures (Yin et al. (2018)).
However, despite its popularity, the MOM estimator suffers from low asymptotic
statistical efficiency. More precisely, the asymptotic efficiency of the MOM
estimator is only $2/\pi\approx 0.637$, which is far from 1 for normal mean
problem.

The main contribution of the paper is to propose a computationally efficient
robust mean estimator, which greatly improves the statistical efficiency of
the MOM. Our estimator is called _variance reduced median-of-means (VRMOM)_
estimator. Instead of using the median in MOM, we use multiple quantile levels
to improve the statistical efficiency. By formulating a carefully designed
stochastic optimization problem and leveraging the idea of one-step Newton
iteration, our VRMOM estimator achieves the same order of computational
complexity as the MOM estimator, but improves the asymptotic efficiency from
$2/\pi\approx 0.637$ of MOM to $3/\pi\approx 0.955$ (see Theorem 1). The
proposed VRMOM estimator naturally serves as a more efficient substitute of
MOM in all robust statistical applications that benefit from the MOM
estimator.

As an application of our VRMOM estimator, we describe a communication-
efficient algorithm for the general risk minimization problem in (1) based on
the VRMOM estimator. In a standard distributed gradient descent (GD) approach,
each local machine computes the gradient information, which takes the form of
the _mean of gradients of each local data point_. Then, the master receives
the transmitted gradient information and aggregates the local gradients by
taking the average. However, the averaged gradient is highly sensitive to
Byzantine failure, whose value can be completely skewed by a single Byzantine
worker. To hedge against Byzantine failures, the work by Yin et al. (2018)
proposed to take the coordinate-wise _median_ of the transmitted gradients,
which is essentially an MOM estimator based on gradients of local data. Our
method improves this result from two aspects. First, instead of using the
median, our VRMOM serves as a new _gradient aggregator_ , which is
statistically more efficient in a large class of distributed robust inference
problems. Second, the distributed gradient method would take a large number of
iterations (i.e., $O(\log(N/p))$) to converge, which is communicationally
expensive. To address this issue, we leverage the surrogate loss function in
the Communication-efficient Surrogate Likelihood (CSL) framework (Jordan et
al., 2019) and develop the _robust CSL (RCSL)_ method. In a wide range of
choices of $N$, $m$, and $p$, our RCSL only requires a constant order of
iterations to achieve a fast convergence rate, which greatly saves the total
communication cost. Theoretically, we establish the convergence rates of our
estimator (see Theorem 5 and Theorem 6 in Appendix C.1.2) and provide the
asymptotic normality result (see Theorems 4).

###  1.1 Contributions and Related Works

The median-of-means (MOM), which was introduced by Nemirovsky and Yudin
(1983), has been a popular estimator in robust statistics due to its ease of
implementation and convergence guarantees (Minsker, 2015; Hsu and Sabato,
2016; Lecué and Lerasle, 2020; Lugosi and Mendelson, 2019; Minsker, 2019). The
MOM estimator finds a wide range of applications, including robust PCA
(Minsker, 2015), linear regression (Hsu and Sabato, 2016), sparse linear
regression (Minsker, 2015; Lecué and Lerasle, 2020), robust empirical risk
minimization (Lecué and Lerasle, 2020; Lugosi and Mendelson, 2019; Minsker,
2019). This paper improves the MOM estimator by proposing a variance reduction
scheme, which significantly boosts the statistical efficiency. Our VRMOM
estimator is motivated by the idea that composite quantiles can improve the
efficiency (Zou and Yuan, 2008). However, directly taking multiple sample
quantiles would incur a higher computational cost and our VRMOM is carefully
designed to be computationally efficient and admit a simple closed-form for
the ease of theoretical analysis. The proposed VRMOM estimator can be a
natural substitute for the classical MOM estimator for all aforementioned
applications.

In recent years, statistical learning and optimization with the presence of
Byzantine failures have attracted a lot of attentions (Feng et al., 2014; Chen
et al., 2017; Blanchard et al., 2017; Xie et al., 2018; Alistarh et al., 2018;
Yin et al., 2018, 2019; Su and Xu, 2019). The key idea behind these work is to
let each worker machine compute the gradient (or stochastic gradient)
information, and the gradients from workers are aggregated using some robust
mean estimators instead of the vanilla gradient mean. There are many
applicable estimators like median, trimmed mean (Yin et al., 2018, 2019),
geometric median (Feng et al., 2014; Chen et al., 2017), Krum (Blanchard et
al., 2017), marginal median, mean-around-median (Xie et al., 2018), and
iterative filtering (Su and Xu, 2019; Yin et al., 2019). However, most
existing methods are only based on gradient information, without utilizing any
second order properties. In this paper, we propose the robust CSL method,
which combines the new gradient aggregator — VRMOM estimator, and the
approximate-Newton framework (See, e.g. Shamir et al. (2014); Jordan et al.
(2019), and Fan et al. (2019)). The combination of the VRMOM and approximate-
Newton greatly facilitates communication efficient estimation by reducing the
total number of communication rounds. From a theoretical perspective, we
establish the asymptotic normality result, which has not been well explored in
previous robust distributed learning literature. A more detailed comparison of
the convergence rates with the existing approaches is presented after Theorem
5, after the formal description of our convergence result.

###  1.2 Paper Organization and Notations

The rest of the paper is organized as follows. Section 2 describes the
proposed VRMOM estimator and its theoretical results. In Section 3, we
introduce the Robust CSL (RCSL) method for Byzantine robust machine learning
problem as an important application of the VRMOM estimator. Simulation
experiments are provided in Section 4, which demonstrate the superiority of
our method over some existing methods. Finally, we conclude our work in
Section 5. The proofs of the theories of the VRMOM estimator and the theories
of the RCSL method are relegated to Appendices.

For every vector $\boldsymbol{v}=(v_{1},...,v_{p})^{\rm T}$, denote
$|\boldsymbol{v}|_{2}=\sqrt{\sum_{l=1}^{p}v_{l}^{2}}$. For every matrix
$\boldsymbol{A}$, define
$\left\|\boldsymbol{A}\right\|=\sup_{|\boldsymbol{v}|_{2}=1}|\boldsymbol{A}\boldsymbol{v}|_{2}$
as the operator norm, $\Lambda_{\max}(\boldsymbol{A})$ and
$\Lambda_{\min}(\boldsymbol{A})$ as the largest and smallest eigenvalues of
$\boldsymbol{A}$ respectively. Suppose there is another matrix
$\boldsymbol{B}$, and we denote $\boldsymbol{B}\preceq\boldsymbol{A}$ if and
only if $\boldsymbol{A}-\boldsymbol{B}$ is positive definite. Let
${\mathcal{N}}(0,1)$ be the standard normal distribution. We denote
$\Phi(x)={\mathbb{P}}({\mathcal{N}}(0,1)\leq x)$ and
$\psi(x)=e^{-x^{2}/2}/\sqrt{2\pi}$ to be its cumulative distribution function
and probability density function, respectively. Denote
${\mathbb{S}}^{p-1}(\boldsymbol{\theta})$ and
${\mathbb{B}}^{p}(\boldsymbol{\theta})$ as the unit sphere and the unit ball
centered at $\boldsymbol{\theta}\in{\mathbb{R}}^{p}$ respectively. For
simplicity, we denote ${\mathbb{S}}^{p-1}$ and ${\mathbb{B}}^{p}$ as unit
sphere and unit ball centered at $\boldsymbol{0}$. We will use
${\mathbb{I}}(\cdot)$ as the indicator function. The symbols $\lfloor
x\rfloor$ ($\lceil x\rceil$) denotes the greatest integer (the smallest
integer) not larger than (not less than) $x$. Summation symbol will be heavily
used throughout this article. For the convenience of reading, in each summand,
we will use the subscripts $i(1\leq i\leq n)$ for each data point, $j(0\leq
j\leq m)$ for each machine, $k(1\leq k\leq K)$ for each quantile level and
$l(1\leq l\leq p)$ for each entry of a vector, respectively. Lastly, the
generic constants are assumed to be independent of $m,n,$ and $p$.

##  2 Proposed Methods

In this section, we will firstly introduce the construction of our VRMOM
estimator. Then we provide theoretical guarantees for it.

###  2.1 Variance Reduced Median-of-Means Estimator

To motivate our estimator, let us provide a brief review of the standard MOM
estimator. Let $X_{1},...,X_{N}$ be _i.i.d._ copies of $X$ with
${\mathbb{E}}(X)=\mu$ and $\mathrm{Var}(X)=\sigma^{2}$. For the ease of
presentation, we assume that $N$ observations are evenly partitioned into
$(m+1)$-batches $\\{{\mathcal{H}}_{0},\dots\mathcal{H}_{m}\\}$, where each
${\mathcal{H}}_{j}$ denotes the indices of the samples within the $j$-th
batch. Let $n=N/(m+1)$ be the sample size of each batch and $\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}=\sum_{i\in{\mathcal{H}}_{j}}X_{i}/n$
be the sample mean of the observations in the $j$-th batch. To estimate the
population mean $\mu$, the MOM estimator is defined as

$\widehat{\mu}=\mathrm{med}(\,\hbox{\kern 0.09995pt\vbox{\hrule
height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{0},...,\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{m}),$ (2)

where $\mathrm{med}(\cdot)$ denotes the sample median. The MOM estimator is
computationally efficient and robust against Byzantine failures. Moreover, as
shown in Minsker (2019), when $m\rightarrow\infty$ and $m=o(\sqrt{N})$, under
some mild moment conditions (e.g., ${\mathbb{E}}|X-\mu|^{3}<\infty$), the MOM
estimator admits the following limiting distribution as $N\rightarrow\infty$,

$\sqrt{N}(\widehat{\mu}-\mu)\xrightarrow{d}{\mathcal{N}}(0,\frac{\pi}{2}\sigma^{2}).$

In addition to the robustness, the statistical efficiency is another important
issue. In a classical statistical estimation setting without Byzantine
failures, we can see the relative efficiency of $\widehat{\mu}$ with respect
to the vanilla sample mean is
$(\sigma^{2})/(\frac{\pi}{2}\sigma^{2})=2/\pi\approx 0.637$, which is far from
the optimal efficiency 1. Therefore, a natural question is:

_Is it possible to construct a computationally efficient robust estimator that
achieves a nearly-optimal efficiency?_

##### The key idea behind our VRMOM estimator

To address this challenge, we first note that by the central limit theorem,
for each sample mean $\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}$,
$\sqrt{n}\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}$ asymptotically
obeys the normal distribution ${\mathcal{N}}(\mu,\sigma^{2})$. Moreover, for
the ease of notation, for a fixed $n$, we define

$\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,=\mu+\epsilon,\qquad\epsilon\sim{\mathcal{N}}(0,\sigma^{2}/n),$
(3)

where $\mu$ and $\sigma$ are unknown. Note that for every quantile level
$\tau$, the $\tau$-th population quantile of the normal distribution
${\mathcal{N}}(\mu,\sigma^{2}/n)$ is exactly
$\mu^{\tau}:=\mu+\sigma\Phi^{-1}(\tau)/\sqrt{n}$. To see this,

${\mathbb{P}}(\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,\leq\mu^{\tau})={\mathbb{P}}((\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,-\mu)/(\sigma/\sqrt{n})\leq\Phi^{-1}(\tau))={\mathbb{P}}(N(0,1)\leq\Phi^{-1}(\tau))=\tau.$

Additionally, by symmetry of normal distribution, we have

$\frac{1}{2}(\mu^{\tau}+\mu^{1-\tau})=\frac{1}{2}\Big{[}2\mu+\frac{\sigma}{\sqrt{n}}\left\\{\Phi^{-1}(\tau)+\Phi^{-1}(1-\tau)\right\\}\Big{]}=\mu.$

Therefore, to improve the statistical efficiency, a natural idea is to
approximate $\mu$ by averaging many pairs of estimators for the
$(\tau,1-\tau)$-th quantiles of ${\mathcal{N}}(\mu,\sigma^{2}/n)$, instead of
using a single quantity (i.e., median). More precisely, let $K$ be a pre-fixed
integer. For any $1\leq k\leq K$, let $\tau_{k}:=k/(K+1)$ and $\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{k}}$ be the
$\tau_{k}$-th sample quantile of $\\{\,\hbox{\kern 0.09995pt\vbox{\hrule
height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{0},\dots,\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{m}\\}$. Since
$\tau_{K+1-k}=1-\tau_{k}$, $\,\hbox{\kern 0.09995pt\vbox{\hrule
height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{k}}$ and
$\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{K+1-k}}$ are
symmetrical about $\mu$ and their average is a natural estimator of $\mu$.
Based on this idea, we can take weighted average of $\\{\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{k}}\\}_{k=1}^{K}$
as an estimator of $\mu$, which improves the statistical efficiency. We also
illustrate the main idea of the weighted averaged estimator in Figure 1. Next,
we introduce a computationally more efficient estimator for implementing this
idea. Moreover, since it is a closed-form estimator, which also facilitates
the theoretical analysis.

Figure 1: Let $K=6$. For $1\leq k\leq K$, $\,\hbox{\kern 0.09995pt\vbox{\hrule
height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{k}}$ is defined
as the $k/7$-th sample quantile of $\\{\,\hbox{\kern 0.09995pt\vbox{\hrule
height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{0},\dots,\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{m}\\}$. We can see that
the pairs $(\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{1}},\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{6}}),(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{2}},\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{5}}),(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{3}},\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{4}})$ are
nearly symmetrical about the targeting parameter $\mu$. Thus the weighted
average of $\\{\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{k}}\\}_{k=1}^{K}$
serves as an estimator of $\mu$.

##### Computationally efficient VRMOM estimator

Denote the quantile loss function as $\rho_{\tau}(z)=z(\tau-{\mathbb{I}}(z\leq
0))$, we consider the following stochastic optimization problem:

$\mathop{\mathrm{argmin}}_{x\in{\mathbb{R}}}\left[{\mathbb{E}}\left\\{G(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,,x)\right\\}\right]:=\mathop{\mathrm{argmin}}_{x\in{\mathbb{R}}}\left[{\mathbb{E}}\left\\{\sum_{k=1}^{K}\rho_{\tau_{k}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,-\frac{\sigma\Delta_{k}}{\sqrt{n}}-x\right)\right\\}\right],$
(4)

where $\Delta_{k}:=\Phi^{-1}(\tau_{k})$, and the expectation is taken over
$X$ in (3). We can easily see that $\mu$ is the solution of (4). To
approximate $\mu$ from (4), we adopt the idea of one-step estimator as
follows. Define

$\displaystyle g(x):=$
$\displaystyle\frac{{\mathrm{d}}}{{\mathrm{d}}x}{\mathbb{E}}\left\\{G(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,,x)\right\\}={\mathbb{E}}\left\\{\sum_{k=1}^{K}{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,\leq
x+\frac{\sigma\Delta_{k}}{\sqrt{n}}\right)-\tau_{k}\right\\},$ $\displaystyle
H(x):=$
$\displaystyle\frac{{\mathrm{d}}}{{\mathrm{d}}x}g(x)=\sum_{k=1}^{K}\mathfrak{p}\left(x-\mu+\frac{\sigma\Delta_{k}}{\sqrt{n}}\right),$

as the gradient and Hessian of the loss function in (4) respectively. Here,
$\mathfrak{p}(x)=\sqrt{n}\psi(\sqrt{n}x/\sigma)/\sigma$ denotes the
probability density function of the noise $\epsilon\sim N(0,\sigma^{2}/n)$ in
(3), and $\psi(\cdot)$ denotes the density function of ${\mathcal{N}}(0,1)$.
Given an initial crude estimator $\mu_{0}$ of $\mu$, the one-step estimator
essentially takes the following Newton-Raphson step:

$\widetilde{\mu}_{1}:=\mu_{0}-g(\mu_{0})/H(\mu_{0})=\mu_{0}-\frac{{\mathbb{E}}\\{\sum_{k=1}^{K}{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,\leq\mu_{0}+\sigma\Delta_{k}/\sqrt{n}\right)-\tau_{k}\\}}{\sum_{k=1}^{K}\mathfrak{p}(\mu_{0}-\mu+\sigma\Delta_{k}/\sqrt{n})}.$
(5)

Next, we use the MOM estimator $\widehat{\mu}$ in (2) as the initial estimator
of (4). The unknown parameter $\sigma^{2}$ can be estimated by
$\widehat{\sigma}^{2}:=\sum_{i\in{\mathcal{H}}_{0}}(X_{i}-\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{0})^{2}/n$, the sample
variance of the first batch of observations ${\mathcal{H}}_{0}$. With the
initial estimator $\mu_{0}$ in place, replacing $g(\mu_{0})$ in (5) with its
empirical counterpart and approximating $H(\mu_{0})$ by
$\sum_{k=1}^{K}\mathfrak{p}(\widehat{\sigma}\Delta_{k}/\sqrt{n})\approx\sqrt{n}\sum_{k=1}^{K}\psi(\Delta_{k})/\widehat{\sigma}$,
we derive the following one-step estimator of (4) from (5)

$\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,=\widehat{\mu}-\frac{\widehat{\sigma}}{(m+1)\sqrt{n}\sum_{k=1}^{K}\psi(\Delta_{k})}\sum_{j=0}^{m}\sum_{k=1}^{K}\left\\{{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\leq\widehat{\mu}+\frac{\widehat{\sigma}\Delta_{k}}{\sqrt{n}}\right)-\frac{k}{K+1}\right\\}.$
(6)

To further alleviate the burden of computation, we choose one summand in (6)
and simplify it as follows

$\displaystyle\sum_{k=1}^{K}\left\\{{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\leq\widehat{\mu}+\frac{\widehat{\sigma}\Delta_{k}}{\sqrt{n}}\right)-\frac{k}{K+1}\right\\}$
$\displaystyle=$ $\displaystyle\sum_{k=1}^{K}{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\leq\widehat{\mu}+\frac{\widehat{\sigma}\Delta_{k}}{\sqrt{n}}\right)-\frac{K}{2}$
$\displaystyle=$
$\displaystyle{\mathrm{C}ard}\left\\{k\;:\;\frac{\sqrt{n}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}-\widehat{\mu})}{\widehat{\sigma}}\leq\Delta_{k},\;1\leq
k\leq K\right\\}-\frac{K}{2}$ $\displaystyle=$
$\displaystyle{\mathrm{C}ard}\left\\{k\;:\;\Phi\left(\frac{\sqrt{n}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}-\widehat{\mu})}{\widehat{\sigma}}\right)\leq\frac{k}{K+1},\;1\leq
k\leq K\right\\}-\frac{K}{2}$ $\displaystyle=$
$\displaystyle\frac{K}{2}+1-\left\lceil(K+1)\Phi\left(\frac{\sqrt{n}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}-\widehat{\mu})}{\widehat{\sigma}}\right)\right\rceil.$

In fact, our theoretical result will show that a larger $K$ leads to a better
statistical efficiency. This derivation shows that it is possible to enhance
the efficiency by taking a larger $K$ without incurring additional
computational cost. Our VRMOM estimator in (6) can be rewritten as

$\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,=\widehat{\mu}-\frac{\widehat{\sigma}}{(m+1)\sqrt{n}\sum_{k=1}^{K}\psi(\Delta_{k})}\sum_{j=0}^{m}\left\\{\frac{K}{2}+1-\left\lceil(K+1)\Phi\left(\frac{\sqrt{n}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}-\widehat{\mu})}{\widehat{\sigma}}\right)\right\rceil\right\\}.$
(7)

Although the expression of VRMOM  $\mu$ in (7) seems more complicated than the
MOM estimator $\widehat{\mu}$, the time complexity of  $\mu$ is at the same
order as the MOM estimator. In particular, at each iteration, each worker
machine computes a local sample mean in parallel, which takes $O(n)$ time
complexity. In MOM estimator, it takes another $O(m)$ operations to find the
median $\widehat{\mu}$ (see Paterson (1996)). While in VRMOM estimator, we
only need an extra $O(m+n+K)$ time complexity ($O(n)$ for sample variance
computed in ${\mathcal{H}}_{0}$, $O(m+K)$ for the variance reduction term in
(7)). Therefore the time complexity of both methods is $O(m+n)$ (here $K$ is
fixed). While keeping the same order of computational complexity, the VRMOM
greatly improves the statistical efficiency. As we can see from Theorem 1 in
the following section, the asymptotic efficiency of  $\mu$ approaches
$3/\pi\approx 0.955$ as $K$ grows to infinity, which is nearly optimal. In
fact, by taking $K=5$, the efficiency has already been more than $0.9$ as
compared to $0.637$ of $\widehat{\mu}$.

######  Remark 1.

As illustrated in Figure 1, we can also find the $\tau_{k}$-th sample quantile
$\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{k}}$ and take a
weighted average of $\\{\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,^{\tau_{k}}\\}_{k=1}^{K}$
as an estimator of $\mu$. However, to give the sample quantiles at $K$
different levels, we need to perform a sorting algorithm among the set
$\\{\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\\}_{j=0}^{m}$, which
takes $O(m\log m)$ operations. Therefore the total complexity would be
$O(n+m\log m)$, which it is more costly compared with $O(m+n)$ complexity of
our VRMOM estimator in (7). The inferior in complexity is exacerbated in
multivariate case. When the dimension $p$ is very large, our coordinate-wise
VRMOM estimator has complexity $O(p(m+n))$, while the average of sample
quantiles would take complexity $O(p(m\log m+n))$.

######  Remark 2.

Although we utilize averaging sample quantiles to motivate our method, the
direct average among sample quantiles cannot tolerate even a small fraction of
Byzantine machines. For example, we assume $K=6$, and there are $1/6$ of
machines are Byzantine. In this case, the $1/7$-th and $6/7$-th sample
quantiles can be completely ruined by the Byzantine machines, and further foil
the weighted average. In contrast, our VRMOM estimator tolerates $1/2-\delta$
(where $\delta\in(0,1/2)$ can be arbitrarily small) fraction of Byzantine
machines, which is better than the direct weighted average of sample quantiles
(see Theorem 2 for more details). To see this in a more intuitive way, we can
take a closer look at (6). Each summand of the correction term is bounded in
the interval $[-1,1]$. Noticing that there is a factor of order
$O(1/(m\sqrt{n}))$ multiplying the summation, the overall magnitude of the
correction term is only of the order $O(1/\sqrt{n})$. In consequence, as long
as the initial estimator (e.g., the MOM estimator) is robust, our proposed
VRMOM estimator is Byzantine robust.

######  Remark 3.

To approximate $\mu$ from (4), we can also directly solve the following
optimization problem

$\mathop{\mathrm{argmin}}_{x\in{\mathbb{R}}}\left[\sum_{k=1}^{K}\sum_{j=0}^{m}\rho_{\tau_{k}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}-\frac{\widehat{\sigma}\Delta_{k}}{\sqrt{n}}-x\right)\right],$

which is the empirical version of (4). This formula is similar as the
univariate composite quantile regression in Zou and Yuan (2008). However, it
is much costly to solve such non-smooth optimization problem in a direct way.
Instead, we leverage the idea of Newton-Raphson step, which greatly improves
computation efficiency.

###  2.2 Theories for VRMOM Estimator

Now, we present several theoretical results for the proposed VRMOM estimator.
We firstly provide the asymptotic normality and convergence result of the
VRMOM estimator in one-dimensional case. Then we extend these results to its
multi-dimensional variant. The proofs of results in this section are all
relegated in Appendix A.

######  Theorem 1 (Asymptotic normality of VRMOM).

Let $N=(m+1)n$ i.i.d. random variables $X_{1},...,X_{N}$ be evenly distributed
in $m+1$ subsets ${\mathcal{H}}_{0},...,{\mathcal{H}}_{m}$. There is a subset
of Byzantine machine indices ${\mathcal{B}}\subseteq\\{1,\dots,m\\}$ with
${\mathrm{C}ard}({\mathcal{B}})=\lfloor\alpha_{n}m\rfloor$, where
$\alpha_{n}=o(m^{-1/2})$. Let

$\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}=\begin{cases}\frac{1}{n}\sum_{i\in{\mathcal{H}}_{j}}X_{i}\quad&j\notin{\mathcal{B}},\\\
\quad*\quad&j\in{\mathcal{B}},\end{cases}$ (8)

and  $\mu$ be defined as in (6). Suppose $X$ satisfies
${\mathbb{E}}(X)=\mu,\mathrm{Var}(X)=\sigma^{2}$. Assume there exists some
$\kappa>0$ such that ${\mathbb{E}}[|X-\mu|^{2+\kappa}]<\infty$, and
$m=o(\min\\{n,n^{2\kappa/(2+\kappa)}\\}),\log^{3}n=o(m)$. Then we have

$\sqrt{N}(\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,-\mu)\xrightarrow{d}{\mathcal{N}}(0,\sigma^{2}_{K}),$

where

$\sigma_{K}^{2}=\frac{\sum_{k_{1},k_{2}=1}^{K}\min(\tau_{k_{1}},\tau_{k_{2}})\\{1-\max(\tau_{k_{1}},\tau_{k_{2}})\\}}{\\{\sum_{k=1}^{K}\psi(\Delta_{k})\\}^{2}}\sigma^{2},$
(9)

with $\tau_{k}=k/(K+1)$. Moreover,
$\lim_{K\rightarrow\infty}\sigma_{K}^{2}=\pi\sigma^{2}/3$.

This theorem provides the asymptotic normality result of our VRMOM estimator
$\mu$ and characterizes the asymptotic variance. In particular, it shows that
$\mu$ is a consistent estimator of $\mu$. Comparing with the MOM estimator, we
_improve the efficiency_ of the estimator by reducing the variance from
$\pi\sigma^{2}/2$ (Minsker, 2019) to $\pi\sigma^{2}/3$ when $K$ goes to
infinity. It should be noted that, we impose the rate constraints
$\alpha_{n}=o(m^{-1/2}),m=o(\min\\{n,n^{2\kappa/(2+\kappa)}\\})$, and
$\log^{3}n=o(m)$ in order to obtain asymptotic normality. In the following
theorem, we drop out these conditions and investigate the convergence rate of
the VRMOM estimator.

######  Theorem 2 (Convergence rate of VRMOM).

Let $N=(m+1)n$ i.i.d. random variables $X_{1},...,X_{N}$ be evenly distributed
in $m+1$ subsets ${\mathcal{H}}_{0},...,{\mathcal{H}}_{m}$. There is a subset
of Byzantine machine indices ${\mathcal{B}}\subseteq\\{1,\dots,m\\}$ with
${\mathrm{C}ard}({\mathcal{B}})=\lfloor\alpha_{n}m\rfloor$, where
$\alpha_{n}\leq 1/2-\delta$ for some fixed $\delta\in(0,1/2)$. Let

$\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}=\begin{cases}\frac{1}{n}\sum_{i\in{\mathcal{H}}_{j}}X_{i}\quad&j\notin{\mathcal{B}},\\\
\quad*\quad&j\in{\mathcal{B}},\end{cases}$

and  $\mu$ be defined as in (6). Suppose $X$ satisfies
${\mathbb{E}}(X)=\mu,\text{Var}(X)=\sigma^{2}$. Assume there exists some
$\kappa>0$ such that ${\mathbb{E}}[|X-\mu|^{2+\kappa}]<\infty$. Then we have

$|\,\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,-\mu|=O_{{\mathbb{P}}}\left(\frac{\alpha_{n}}{\sqrt{n}}+\frac{1}{\sqrt{mn}}+\frac{1}{n^{(3\kappa_{2}+2)/(2\kappa_{2}+4)}}+\frac{\log^{3/4}n}{n^{1/2}m^{3/4}}\right),$
(10)

where $\kappa_{2}=\min(\kappa,2)$.

This convergence result shows that our VRMOM estimator is consistent as long
as $\alpha_{n}$ is strictly smaller than $1/2$. The condition on $\alpha_{n}$
is also necessary because clearly the sample median can be ruined when there
are more than $\lceil m/2\rceil$ corruptions. From (10), when
$m=O(\min\\{n,n^{2\kappa/(2+\kappa)}\\}),\log^{3}n=O(m)$, the rate matches the
optimal rate $O(\alpha_{n}/\sqrt{n}+1/\sqrt{mn})$ (See Observation 1 in Yin et
al. (2018)). Further assume that $\alpha_{n}=O(1/\sqrt{m})$, the VRMOM
achieves square root-$N$ consistency.

Next we extend our VRMOM estimator to the multi-dimensional extension setting.
Let $\boldsymbol{X}_{1},\dots,\boldsymbol{X}_{N}$ be i.i.d. copies of the
$p$-dimensional random vectors
$\boldsymbol{X}=(X^{(1)},\dots,X^{(p)})^{{\mathrm{T}}}$ with
${\mathbb{E}}(\boldsymbol{X})=\boldsymbol{\mu}=(\mu^{(1)},\dots,\mu^{(p)})^{{\mathrm{T}}}$
and
${\mathrm{C}ov}(\boldsymbol{X})=\boldsymbol{\Sigma}=(\sigma_{l_{1},l_{2}})_{l_{1},l_{2}=1}^{p}$.
Then the multi-dimensional VRMOM estimator  $\boldsymbol{\mu}$ is defined by
applying (7) on each coordinate $l$, where $1\leq l\leq p$. We first obtain
the convergence rate of the multi-dimensional VRMOM estimator in terms of
$\ell_{2}$-norm.

######  Theorem 3 (Convergence rate of multi-dimensional VRMOM).

Let $N=(m+1)n$ i.i.d. random vectors
$\boldsymbol{X}_{1},...,\boldsymbol{X}_{N}$ be evenly distributed in $m+1$
subsets ${\mathcal{H}}_{0},...,{\mathcal{H}}_{m}$. There is a subset of
Byzantine machine indices ${\mathcal{B}}\subseteq\\{1,\dots,m\\}$ with
${\mathrm{C}ard}({\mathcal{B}})=\lfloor\alpha_{n}m\rfloor$, where
$\alpha_{n}\leq 1/2-\delta$ for some fixed $\delta\in(0,1/2)$. Let
$\boldsymbol{\mu}$ be the multi-dimensional VRMOM estimator defined in (7).
Suppose $\boldsymbol{X}$ satisfies
${\mathbb{E}}(\boldsymbol{X})=\boldsymbol{\mu},{\mathrm{C}ov}(\boldsymbol{X})=\boldsymbol{\Sigma}$.
Moreover, for each coordinate $l\in\\{1,\dots,p\\}$, we assume there exists
some $\kappa>0$ such that
${\mathbb{E}}[|X^{(l)}-\mu^{(l)}|^{2+\kappa}]<\infty$. Then we have

$|\,\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{\mu}$\kern-1.00006pt}}}\,-\boldsymbol{\mu}|_{2}=O_{{\mathbb{P}}}\left(\frac{\alpha_{n}\sqrt{p}}{\sqrt{n}}+\sqrt{\frac{p}{mn}}+\frac{\sqrt{p}}{n^{(3\kappa_{2}+2)/(2\kappa_{2}+4)}}+\frac{p^{1/2}\log^{3/4}n}{n^{1/2}m^{3/4}}\right),$
(11)

where $\kappa_{2}=\min(\kappa,2)$.

As we can see from the theorem, the convergence rate of multi-dimensional
VRMOM estimator is simply the rate in Theorem 2 multiplied with $\sqrt{p}$,
which is not surprising because the VRMOM estimator is applied coordinate-
wisely. Moreover, to guarantee consistency of the proposed estimator, we
require the rate on the right hand side of (11) to be the order of
$o_{{\mathbb{P}}}(1)$, which implies that

$p=o\Big{(}\min\Big{\\{}\frac{n}{\alpha^{2}_{n}},mn,n^{\frac{3\kappa_{2}+2}{\kappa_{2}+2}},\frac{nm^{3/2}}{\log^{3/2}n}\Big{\\}}\Big{)}.$
(12)

In particular, we are more interested in the asymptotic normality of the
multi-dimensional VRMOM estimator, which will be presented in the next
theorem.

By definition we know
$\sigma_{l_{1},l_{2}}={\mathrm{C}ov}\\{X^{(l_{1})},X^{(l_{2})}\\}$ is the
$(l_{1},l_{2})$-entry of covariance matrix of $\boldsymbol{X}$. Let
$(Z_{l_{1}},Z_{l_{2}})$ admit the following bivariate normal distribution

${\mathcal{N}}\left(\boldsymbol{0},\bm{\Sigma}_{l_{1},l_{2}}\right),\quad\text{where}\quad\bm{\Sigma}_{l_{1},l_{2}}=\begin{pmatrix}1&\frac{\sigma_{l_{1},l_{2}}}{\sqrt{\sigma_{l_{1},l_{1}}\sigma_{l_{2},l_{2}}}}\\\
\frac{\sigma_{l_{1},l_{2}}}{\sqrt{\sigma_{l_{1},l_{1}}\sigma_{l_{2},l_{2}}}}&1\end{pmatrix}.$
(13)

Next we define the $p\times p$ matrix $\boldsymbol{\mathcal{C}}$ with its
$(l_{1},l_{2})$-entry given by the following formula:

$\mathcal{C}_{l_{1},l_{2}}=\frac{\sum_{k_{1},k_{2}=1}^{K}(\tau^{l_{1},l_{2}}_{k_{1},k_{2}}-\tau_{k_{1}}\tau_{k_{2}})}{\\{\sum_{k=1}^{K}\psi(\Delta_{k})\\}^{2}}\sqrt{\sigma_{l_{1},l_{1}}\sigma_{l_{2},l_{2}}},$
(14)

where $\tau_{k}=k/(K+1)$, and
$\tau^{l_{1},l_{2}}_{k_{1},k_{2}}={\mathbb{P}}(Z_{l_{1}}\leq\Delta_{k_{1}},Z_{l_{2}}\leq\Delta_{k_{2}})$.
Then we can prove the following asymptotic normality result:

######  Theorem 4 (Asymptotic normality of multi-dimensional VRMOM estimator).

Under the same assumption as in Theorem 3, and additionally, we assume the
rate constraints
$p=o(\min\\{\frac{m^{1/2}}{\log^{3/2}n},\frac{n^{2\kappa_{2}/(\kappa_{2}+2)}}{m}\\})$,
and $\alpha_{n}=o(1/\sqrt{mp})$. Then for any vector
$\boldsymbol{v}\in{\mathbb{R}}^{p}$ with $|\boldsymbol{v}|_{2}=1$, we have
that

$\frac{\sqrt{N}}{\sigma_{\boldsymbol{v}}}\left\langle\boldsymbol{v},\,\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{\mu}$\kern-1.00006pt}}}\,-\boldsymbol{\mu}\right\rangle\xrightarrow{d}{\mathcal{N}}(0,1),$
(15)

as $n\to\infty$, where
$\sigma^{2}_{\boldsymbol{v}}=\boldsymbol{v}^{{\mathrm{T}}}\boldsymbol{\mathcal{C}}\boldsymbol{v}$.

To prove asymptotic normality result for the multi-dimensional VRMOM
estimator, we need a more restrictive constraint on the dimension $p$ than the
one in (12). Moreover, we require the number of Byzantine machines is
$o(\sqrt{m/p})$, i.e., the fraction $\alpha_{n}=o(1/\sqrt{mp})$. As compared
to the condition $\alpha_{n}=o(1/\sqrt{m})$ in Theorem 1, there is an extra
$1/\sqrt{p}$ in the condition since we are dealing with a $p$-dimensional
multivariate inference problem.

In order to illustrate the efficiency of our VRMOM estimator in multi-
dimensional case, let us consider the multi-dimensional median-of-means (MOM)
estimator $\widehat{\boldsymbol{\mu}}_{\mathrm{MOM}}$. More specifically, we
also apply the MOM estimator at each coordinate and establish the following
parallel asymptotic normality result for the multi-dimensional MOM estimator.

######  Proposition 1 (Asymptotic normality of multi-dimensional MOM
estimator).

Under the same assumption as in Theorem 3, and additionally, we assume the
rate constraints
$p=o(\min\\{\frac{m^{1/2}}{\log^{3/2}n},\frac{n^{2\kappa_{2}/(\kappa_{2}+2)}}{m}\\})$,
and $\alpha_{n}=o(1/\sqrt{mp})$. Then for any vector
$\boldsymbol{v}\in{\mathbb{R}}^{p}$ with $|\boldsymbol{v}|_{2}=1$, we have
that

$\frac{\sqrt{N}}{\sigma_{\mathrm{MOM},\boldsymbol{v}}}\left\langle\boldsymbol{v},\,\widehat{\boldsymbol{\mu}}_{\rm
MOM}-\boldsymbol{\mu}\right\rangle\xrightarrow{d}{\mathcal{N}}(0,1),$ (16)

as $n\to\infty$, where
$\sigma^{2}_{\mathrm{MOM},\boldsymbol{v}}=\boldsymbol{v}^{{\mathrm{T}}}\boldsymbol{\mathcal{C}}_{\rm
MOM}\boldsymbol{v}$, and $\boldsymbol{\mathcal{C}}_{\rm MOM}$ is a $p\times p$
matrix with each $(l_{1},l_{2})$-entry taking the following form,

$\mathcal{C}_{\mathrm{MOM},l_{1},l_{2}}=\left(2\pi\tau^{l_{1},l_{2}}_{(K+1)/2,(K+1)/2}-\frac{\pi}{2}\right)\sqrt{\sigma_{l_{1},l_{1}}\sigma_{l_{2},l_{2}}}.$
(17)

When each coordinate of random vector $\boldsymbol{X}$ is independent, the
off-diagonal entries of the matrices $\boldsymbol{\mathcal{C}}$ and
$\boldsymbol{\mathcal{C}}_{\rm MOM}$ are all zero (in this case
$\tau_{k_{1},k_{2}}^{l_{1},l_{2}}=\tau_{k_{1}}\tau_{k_{2}}$ for $l_{1}\neq
l_{2}$). For diagonal entries, we can readily compute that

$\mathcal{C}_{l,l}=\frac{\sum_{k_{1},k_{2}=1}^{K}\min(\tau_{k_{1}},\tau_{k_{2}})\\{1-\max(\tau_{k_{1}},\tau_{k_{2}})\\}}{\\{\sum_{k=1}^{K}\psi(\Delta_{k})\\}^{2}}\sigma_{l,l},\quad\mathcal{C}_{\mathrm{MOM},l,l}=\frac{\pi}{2}\sigma_{l,l}.$

According to Theorem 1, when $K\rightarrow\infty$, we have
$\mathcal{C}_{l,l}\rightarrow\frac{\pi}{3}\sigma_{l,l}$, which suggests that
our multi-dimensional VRMOM estimator has a higher statistical efficiency than
the corresponding MOM estimator.

######  Remark 4.

In two-dimensional case, the covariance matrix of $\boldsymbol{X}$ can be
written as the following form

$\boldsymbol{\Sigma}=\begin{pmatrix}\sigma_{1,1}&\sin\phi\sqrt{\sigma_{1,1}\sigma_{2,2}}\\\
\sin\phi\sqrt{\sigma_{1,1}\sigma_{2,2}}&\sigma_{2,2}\end{pmatrix},$

where $\phi\in[-\frac{\pi}{2},\frac{\pi}{2}]$. Then we have that our
2-dimensional $\mathrm{VRMOM}$ estimator  $\boldsymbol{\mu}$ has higher
statistical efficiency than the 2-dimensional $\mathrm{MOM}$ estimator
$\widehat{\boldsymbol{\mu}}$ as $K$ tends to infinity. The detailed argument
is relegated to Appendix B. In the higher dimension case when $p>2$, we
believe that the superiority in efficiency of our $\mathrm{VRMOM}$ estimator
still holds. We leave the theoretical investigation as a future work.

##  3 Application for Byzantine Distributed Statistical Optimization

As an important application of the proposed VRMOM estimator, in this section,
we consider the general distributed statistical optimization problem in (1)
under Byzantine setup. In particular, we propose a Byzantine robust
distributed approximate newton method, called Robust Communication-efficient
Surrogate Likelihood (RCSL) Method.

MasterMachine $\mathcal{H}_{0}$ptWorkerMachine $\mathcal{H}_{1}$WorkerMachine
$\mathcal{H}_{2}$$\bullet\quad\bullet\quad\bullet$WorkerMachine
$\mathcal{H}_{m}$$\widehat{\boldsymbol{\theta}}^{(t-1)}$$\widehat{\boldsymbol{\theta}}^{(t-1)}$$\widehat{\boldsymbol{\theta}}^{(t-1)}$$\boldsymbol{g}^{(t-1)}_{1}$$\boldsymbol{g}^{(t-1)}_{2}$$\boldsymbol{g}^{(t-1)}_{m}$pt$\widehat{\boldsymbol{\theta}}^{(t)}$
Figure 2: Communication protocol of the robust CSL (RCSL) method. In the
$t$-th iteration, the master machine ${\mathcal{H}}_{0}$ distributes the
parameter $\widehat{\boldsymbol{\theta}}^{(t-1)}$ to each worker machine. The
$j$-th worker machine computes the local gradient $\boldsymbol{g}_{j}^{(t-1)}$
and sends it back to the master machine. Then ${\mathcal{H}}_{0}$ updates the
new parameter $\widehat{\boldsymbol{\theta}}^{(t)}$ and repeats the procedure.

For the ease of presentation, we adopt the master/worker setting in Jordan et
al. (2019), where ${\mathcal{H}}_{0}$ denotes the master machine and the rest
are worker machines. We assume that the master machine ${\mathcal{H}}_{0}$
stores $n$ observations as each local worker and the data on the master
machine will not be corrupted. In practice, it is easier to use one powerful
machine as the master machine that is robust. Let
$\widehat{\boldsymbol{\theta}}^{(0)}$ be an initial estimator of
$\boldsymbol{\theta}^{*}$. At the beginning, the master machine
${\mathcal{H}}_{0}$ broadcasts the parameter
$\widehat{\boldsymbol{\theta}}^{(0)}$ to all worker machines
${\mathcal{H}}_{1},\dots,{\mathcal{H}}_{m}$. The $j$-th worker computes a
local gradient
$\boldsymbol{g}_{j}^{(0)}=n^{-1}\sum_{i\in{\mathcal{H}}_{j}}\nabla|_{\boldsymbol{\theta}}f(X_{i},\widehat{\boldsymbol{\theta}}^{(0)})$
and sends it back to master. Then the master machine applies the VRMOM
estimator to every coordinate. More precisely, for each of the $l$-th
coordinate (we will use the subscript $l$ to represent the entry of a vector),
master machine computes the VRMOM estimator

$\displaystyle\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$g$\kern-1.00006pt}}}\,_{l}^{(0)}=$
$\displaystyle\widehat{g}_{l}^{(0)}-\frac{\widehat{\sigma}_{l}^{(0)}}{(m+1)\sqrt{n}\sum_{k=1}^{K}\psi(\Delta_{k})}\times$
(18)
$\displaystyle\sum_{j=0}^{m}\left\\{\frac{K}{2}+1-\left\lceil(K+1)\Phi\left(\frac{\sqrt{n}(g_{j,l}^{(0)}-\widehat{g}_{l}^{(0)})}{\widehat{\sigma}^{(0)}_{l}}\right)\right\rceil\right\\},$

where $\widehat{g}_{l}^{(0)}=\text{med}\\{g_{0,l}^{(0)},...,g_{m,l}^{(0)}\\}$
is the median and

$(\widehat{\sigma}_{l}^{(0)})^{2}=\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}\left\\{\nabla|_{\boldsymbol{\theta}}f_{l}(X_{i},\widehat{\boldsymbol{\theta}}^{(0)})-g_{0,l}^{(0)}\right\\}^{2},$

is the local sample variance. Here,
$\nabla|_{\boldsymbol{\theta}}f_{l}(X_{i},\widehat{\boldsymbol{\theta}}^{(0)})$
denotes the $l$-th coordinate of
$\nabla|_{\boldsymbol{\theta}}f(X_{i},\widehat{\boldsymbol{\theta}}^{(0)})$.
Let $\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{g}$\kern-1.00006pt}}}\,^{(0)}=(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$g$\kern-1.00006pt}}}\,_{1}^{(0)},\dots\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$g$\kern-1.00006pt}}}\,_{p}^{(0)})^{\rm T}$
denote the VRMOM aggregated gradient. The master machine solves the following
surrogate loss introduced by Jordan et al. (2019) to update the parameter,

$\widehat{\boldsymbol{\theta}}^{(1)}=\mathop{\mathrm{argmin}}_{\boldsymbol{\theta}\in{\mathbb{R}}^{p}}\left\\{\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}f(X_{i},\boldsymbol{\theta})-\left\langle\boldsymbol{g}_{0}^{(0)}-\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{g}$\kern-1.00006pt}}}\,^{(0)},\boldsymbol{\theta}\right\rangle\right\\},$
(19)

where
$\boldsymbol{g}_{0}^{(0)}=n^{-1}\sum_{i\in{\mathcal{H}}_{0}}\nabla|_{\boldsymbol{\theta}}f(X_{i},\widehat{\boldsymbol{\theta}}^{(0)})$
is the local gradient computed on the master machine. As shown in Jordan et
al. (2019) and our experiments, for a wide range of statistical learning
problems, the surrogate loss in (19) can be easily minimized by existing
optimization solvers. Moreover, this surrogate loss minimization is done on
the master machine, and thus does not involve any communication.

Algorithm 1 Robust CSL (RCSL) Method

Input: The data $\\{X_{1},...,X_{N}\\}$ is evenly distributed on $m+1$
machines $\\{{\mathcal{H}}_{0},\dots,{\mathcal{H}}_{m}\\}$. Let
${\mathcal{H}}_{0}$ be the master and the rest be workers.

1: Compute an initial estimator $\widehat{\boldsymbol{\theta}}^{(0)}$ on the
master machine ${\mathcal{H}}_{0}$.

2: for $t=1,\dots,T$ do

3: Distribute $\widehat{\boldsymbol{\theta}}^{(t-1)}$ to each local machine
$j=1,2,\dots,m$.

4: for $j=0,\dots,m$ do

5: The $j$-th worker machine computes the local gradient
$\boldsymbol{g}_{j}^{(t-1)}=\begin{cases}n^{-1}\sum_{i\in{\mathcal{H}}_{j}}\nabla|_{\boldsymbol{\theta}}f(X_{i},\widehat{\boldsymbol{\theta}}^{(t-1)})\quad&\text{if
${\mathcal{H}}_{j}$ is normal},\\\ *\quad&\text{if ${\mathcal{H}}_{j}$ is
Byzantine},\end{cases}$ where $*$ denotes arbitrary values. Then the $j$-th
worker sends $\boldsymbol{g}_{j}^{(t-1)}$ back to master machine.

6: end for

7: Master machine constructs the VRMOM aggregated gradient $\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{g}$\kern-1.00006pt}}}\,^{(t-1)}=(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$g$\kern-1.00006pt}}}\,^{(t-1)}_{1},\dots\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$g$\kern-1.00006pt}}}\,^{(t-1)}_{p})^{\rm T}$,
where each $l$-th coordinate takes the following form,
$\displaystyle\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$g$\kern-1.00006pt}}}\,_{l}^{(t-1)}=$
$\displaystyle\widehat{g}_{l}^{(t-1)}-\frac{\widehat{\sigma}_{l}^{(t-1)}}{(m+1)\sqrt{n}\sum_{k=1}^{K}\psi(\Delta_{k})}\times$
(20)
$\displaystyle\sum_{j=0}^{m}\left\\{\frac{K}{2}+1-\left\lceil(K+1)\Phi\left(\frac{\sqrt{n}(g_{j,l}^{(t-1)}-\widehat{g}_{l}^{(t-1)})}{\widehat{\sigma}^{(t-1)}_{l}}\right)\right\rceil\right\\},$
where
$\widehat{g}_{l}^{(t-1)}=\text{med}\\{g_{0,l}^{(t-1)},...,g_{m,l}^{(t-1)}\\}$,
and
$\left(\widehat{\sigma}^{(t-1)}_{l}\right)^{2}=\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}\left\\{\nabla|_{\boldsymbol{\theta}}f_{l}(X_{i},\widehat{\boldsymbol{\theta}}^{(t-1)})-g_{0,l}^{(t-1)}\right\\}^{2}.$

8: Master machine solves the following surrogate loss minimization problem
$\widehat{\boldsymbol{\theta}}^{(t)}=\mathop{\mathrm{argmin}}_{\boldsymbol{\theta}\in{\mathbb{R}}^{p}}\left\\{\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}f(X_{i},\boldsymbol{\theta})-\left\langle\,\boldsymbol{g}_{0}^{(t-1)}-\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{g}$\kern-1.00006pt}}}\,^{(t-1)},\,\boldsymbol{\theta}\,\right\rangle\right\\}.$
(21)

9: end for

Output: The final estimator $\widehat{\boldsymbol{\theta}}^{(T)}$.

Repeating the above procedure, we develop a multi-round algorithm named Robust
CSL (RCSL), which is presented in Algorithm 1. We note that in the Byzantine
setting, there is a subset of workers ${\mathcal{B}}\subseteq\\{1,\dots,m\\}$
, which return arbitrary values in each iteration. For $j\in{\mathcal{B}}$, we
will use $\boldsymbol{g}_{j}^{(t-1)}=*$ to represent these nuisance values. To
guarantee the consistency of the initial estimator, in Step 1 of Algorithm 1,
we can compute $\widehat{\boldsymbol{\theta}}^{(0)}$ by the local empirical
risk minimization on the master machine ${\mathcal{H}}_{0}$, i.e.,

$\widehat{\boldsymbol{\theta}}^{(0)}=\mathop{\mathrm{argmin}}_{\boldsymbol{\theta}\in{\mathbb{R}}^{p}}\left\\{\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}f(X_{i},\boldsymbol{\theta})\right\\}.$
(22)

We note that our theoretical result only requires the consistency of the
initial estimator, and thus other consistent estimators could also be used as
the initial estimator.

Now we briefly comment on the communication cost of Algorithm 1. In each
round, the communication cost is $O(mp)$, which is at the same order as other
gradient descent algorithms in the literature (e.g., Yin et al. (2018, 2019)).
Moreover, from Theorem 5 below, our RCSL only takes a constant number of
rounds of communication, as opposed to the order of $O(\log(N/p))$ rounds in
other gradient-based methods. Therefore, the total communication cost of RCSL
is only $O(mp)$. For the sake of clarity, we only present the result of multi-
round convergence rate of the RCSL method in the following. More detailed
technical conditions and theoretical results are relegated to Appendix C.

######  Theorem 5 (Multi-round convergence rate of RCSL method).

Suppose Assumption A-G (see Appendix C.1.1) hold and the initial estimator
$\widehat{\boldsymbol{\theta}}^{(0)}$ satisfies
$|\widehat{\boldsymbol{\theta}}^{(0)}-\boldsymbol{\theta}^{*}|_{2}=O_{{\mathbb{P}}}(r_{n})$.
Further assume the fraction $\alpha_{n}$ of Byzantine machines satisfies
$\alpha_{n}\leq 1/2-\delta$ for some fixed $\delta\in(0,1/2)$. The RCSL
estimator in the $t$-th iteration $\widehat{\boldsymbol{\theta}}^{(t)}$
defined in (21) satisfies

$|\widehat{\boldsymbol{\theta}}^{(t)}-\boldsymbol{\theta}^{*}|_{2}=O_{{\mathbb{P}}}\left(\frac{\alpha_{n}\sqrt{p}}{\sqrt{n}}+\sqrt{\frac{p\log
n}{mn}}+\frac{p^{1/2}\log^{3/4}n}{n^{1/2}m^{3/4}}+r_{n}p^{t}\left(\frac{\log
n}{n}\right)^{t/2}\right).$ (23)

The proof of this theorem can be found in Appendix C.1.2. We note that the
second term $\sqrt{p\log n/(mn)}$ in the right hand side of (23) matches the
optimal rate $\sqrt{p/(mn)}$ up to a logarithmic factor. The third term is
inherited from the last term in (10) of Theorem 2. It becomes $O(\sqrt{p\log
n/(mn)})$ when $\log n=O(m)$. The first term is the price paid for the
Byzantine failures. As we can see, after each iteration, the fourth term in
(23) is improved by a factor $p\sqrt{\log n/n}$, which is of the order $o(1)$
by the rate constraints in Assumption G (see Appendix C.1.1). In stark
contrast, for the Byzantine robust gradient descent, the convergence rate is
only improved by a constant factor $c<1$ after each round of communication
(See, e.g. Theorem 3 in Yin et al. (2019) and Theorem 3.4 in Alistarh et al.
(2018)). Therefore, Theorem 5 suggests that the proposed RCSL method enjoys
faster convergence rate than vanilla Byzantine robust gradient descent,
therefore is more communication-efficient. We can also demonstrate the
communication-efficiency of our RCSL method in another way. From the
expression of (23), we can see that when the number of $t$ is sufficiently
large, i.e.,

$t\geq\frac{\log m+\log p+2\log r_{n}}{\log n-2\log p-\log\log n}+1,$ (24)

the last term will become the order of $O(\sqrt{p\log n/(mn)})$. Moreover,
with $p=O(n^{1/3}/\sqrt{\log n})$ and $r_{n}=o(1)$ in Assumption G, we have

$\frac{\log m+\log p+2\log r_{n}}{\log n-2\log p-\log\log n}+1\leq\frac{3\log
m}{\log n}+c_{0},\quad\text{for some constant }c_{0}>0.$

Therefore, when $t\geq c_{0}+3\log m/\log n$, the rate of $t$-th iteration is
dominated by the first three terms. It would save a lot of communication cost
as compared to gradient-based algorithms, where at least $O(\log(N/p))$
iterations of communications are necessary.

Assume $\log n=O(m)$ and the number of iterations $t$ is sufficiently large,
our rate in (23) will become
$O_{{\mathbb{P}}}(\alpha_{n}\sqrt{p}/\sqrt{n}+\sqrt{p\log n/(mn)})$. It is
interesting to compare this rate with contemporary results of the Byzantine
perturbed gradient method with different aggregators (see Theorem 2 in Yin et
al. (2019)). For example, median aggregator leads to the rate
$\alpha_{n}\sqrt{p}/\sqrt{n}+p/\sqrt{mn}$ (we save a $\sqrt{p}$ in the second
term), the trimmed-mean aggregator to the rate
$\alpha_{n}p/\sqrt{n}+p/\sqrt{mn}$ (saving a $\sqrt{p}$ in both terms), and
the iterative filtering to the rate $\sqrt{\alpha_{n}/n}+\sqrt{p/mn}$ (saving
a $\sqrt{\alpha_{n}}$ but losing a $\sqrt{p}$ in the first term). However,
their filtering estimator involves solving convex programs iteratively.
Moreover, in the theory of the filtering estimator, the Byzantine fraction
$\alpha_{n}$ is required to be not larger than $1/4$ (see, e.g. Theorem 1 in
Su and Xu (2019) and Theorem 5 in Yin et al. (2019)), which is more
restrictive than ours ($\alpha_{n}\leq 1/2-\delta$). We would also like to
note that the lower bound on the convergence rate is known to be
$\Omega(\alpha_{n}/\sqrt{n}+\sqrt{p/mn})$. As compared to the lower bound, our
upper bound is missing a $\sqrt{p}$ factor in the first term. When there is no
Byzantine worker (i.e., $\alpha_{n}=0$) or the dimensionality $p$ is a
constant, our rate is optimal upto logarithmic factors. We believe this extra
$\sqrt{p}$ in the first term comes out because the gradients are aggregated
coordinate-wisely. It would be interesting as a future direction to develop a
new multi-variate aggregator based on VRMOM, which is both statistically and
computationally efficient and achieves the optimal convergence rate.

It is also worthwhile noting that we can extend our algorithm to a general
scheme by replacing (19) with the following surrogate loss minimization,

$\mathop{\mathrm{argmin}}_{\boldsymbol{\theta}\in{\mathbb{R}}^{p}}\left\\{\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}f(X_{i},\boldsymbol{\theta})-\left\langle\boldsymbol{g}^{(0)}_{0}-\mathrm{Aggr}(\boldsymbol{g}^{(0)}_{0},\dots,\boldsymbol{g}^{(0)}_{m}),\boldsymbol{\theta}\right\rangle\right\\},$
(25)

where $\mathrm{Aggr}(\boldsymbol{g}^{(0)}_{0},\dots,\boldsymbol{g}^{(0)}_{m})$
can be any consistent estimator of the population gradient
${\mathbb{E}}\\{\nabla|_{\boldsymbol{\theta}}f(X,\widehat{\boldsymbol{\theta}}^{(0)})\\}$
given $\widehat{\boldsymbol{\theta}}^{(0)}$. Any robust aggregator in the
literature can be adopted in (25), e.g., median-of-means, trimmed mean (Yin et
al., 2018), geometric median (Feng et al., 2014; Chen et al., 2017), Krum
(Blanchard et al., 2017), marginal median (Xie et al., 2018). In the original
CSL framework (Jordan et al., 2019), the aggregator is chosen as the
coordinate-wise average, which is sensitive to corruptions.

######  Remark 5.

It is worthwhile noting that when the target parameter admits some specific
structures, e.g., sparsity structure, it is straightforward to extend the
proposed framework (25) to the following regularized problem

$\mathop{\mathrm{argmin}}_{\boldsymbol{\theta}\in{\mathbb{R}}^{p}}\left\\{\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}f(X_{i},\boldsymbol{\theta})-\left\langle\boldsymbol{g}^{(0)}_{0}-\mathrm{Aggr}(\boldsymbol{g}^{(0)}_{0},\dots,\boldsymbol{g}^{(0)}_{m}),\boldsymbol{\theta}\right\rangle+\lambda_{n}{\mathcal{R}}(\boldsymbol{\theta})\right\\},$
(26)

where ${\mathcal{R}}(\boldsymbol{\theta})$ is some regularizer, for example,
the $\ell_{1}$-penalty (Tibshirani, 1996), smooth clipped absolute deviation
(SCAD) (Fan and Li, 2001), and minimax concave penalty (MCP) (Zhang, 2010).
With the above formulation, we are able to address the sparse learning problem
in Byzantine robust setup. We leave more rigorous theoretical investigation of
the penalized Byzantine robust estimation to future research.

##  4 Simulation Studies

In our simulation studies, we conduct several experiments to show the
effectiveness of the VRMOM and the robust CSL (RCSL) Method.

###  4.1 Results for VRMOM

In this section, we show the performance of the proposed VRMOM estimator for
robust mean estimation problem. We first demonstrate how the number of
quantile levels $K$ in (7) affects the estimation accuracy of the VRMOM
estimator, and then compare the statistical efficiency of our VRMOM estimator
and the MOM estimator defined in (2). We generate the random vectors $X_{i}$s
from the normal distribution
$\mathcal{N}(\boldsymbol{\mu}^{\ast},\boldsymbol{\Sigma}_{X})$ where
$\boldsymbol{\mu}^{\ast}=p^{-1/2}(1,(p-2)/(p-1),(p-3)/(p-1),\dots,0)$ and
$\boldsymbol{\Sigma}_{{\mathbf{X}}}=\text{diag}(1,\dots,1)$. We choose $p=1$
and $p=30$ to consider both univariate and multivariate cases. The entire
sample size is $N=1000\times(100+1)$. By dividing the data into one master
machine ${\mathcal{H}}_{0}$ and $100$ worker machines
$\\{{\mathcal{H}}_{1},\dots,{\mathcal{H}}_{100}\\}$. Note that the master
machine ${\mathcal{H}}_{0}$ can never be corrupted in our setup. Therefore,
each local sample size is $n=1000$. We consider the following settings: (1)
$\alpha_{n}=0$ which denotes no Byzantine machine, (2) $\alpha_{n}>0$, which
denotes the existence of Byzantine machine case. We vary the fraction of
Byzantine machines $\alpha_{n}=0.05,0.1,0.15$. When $\alpha_{n}>0$, we replace
the sample means in each Byzantine machine by a random vector whose entries
are generated from $\mathcal{N}(0,200{\mathbb{I}})$ independently. For each
experiment, we repeat 500 independent simulations and report averaged root
mean square estimation errors and standard deviations.

####  4.1.1 Effect of $K$

In the first experiment, we vary the number of quantile levels $K$ from
$\\{10,20,50,100\\}$ and investigate the estimation accuracy of the VRMOM
estimator for different dimensions $p\in\\{1,30\\}$ and different fractions
$\alpha_{n}\in\\{0,0.05,0.1,0.15\\}$. The results of the root mean square
errors and the standard errors are presented in Table 1. As we can see from
the table, for each fraction of Byzantine machine $\alpha_{n}$, the root mean
square errors of the VRMOM estimator for different $K$s are almost the same.
Based on this observation, in the following experiment, we fix $K$ to be $10$
for the ease of computation.

Table 1: The root mean square errors (RMSEs) and their standard errors (in
parentheses) of the VRMOM under sample size $N=1000\times(100+1)$, local
sample size $n=1000$ and number of quantile levels $K=10,20,50,100$. The
sample means sent from Byzantine machines are generated from Gaussian
$\mathcal{N}(0,200{\mathbb{I}})$.
$p$ | $K$ | $\alpha_{n}=0$ | $\alpha_{n}=0.05$ | $\alpha_{n}=0.1$ | $\alpha_{n}=0.15$
---|---|---|---|---|---
1 | 10 | 0.0025 (0.0018) | 0.0027 (0.0020) | 0.0030 (0.0022) | 0.0032 (0.0024)
20 | 0.0025 (0.0018) | 0.0027 (0.0021) | 0.0030 (0.0022) | 0.0034 (0.0026)
| 50 | 0.0025 (0.0018) | 0.0028 (0.0021) | 0.0030 (0.0022) | 0.0033 (0.0024)
| 100 | 0.0025 (0.0018) | 0.0028 (0.0020) | 0.0030 (0.0022) | 0.0032 (0.0024)
30 | 10 | 0.0175 (0.0022) | 0.0192 (0.0024) | 0.0209 (0.0026) | 0.0227 (0.0028)
20 | 0.0174 (0.0022) | 0.0192 (0.0024) | 0.0209 (0.0026) | 0.0228 (0.0030)
| 50 | 0.0174 (0.0022) | 0.0192 (0.0025) | 0.0208 (0.0026) | 0.0230 (0.0029)
| 100 | 0.0174 (0.0022) | 0.0192 (0.0024) | 0.0208 (0.0027) | 0.0230 (0.0030)

####  4.1.2 Comparison between VRMOM and MOM

Table 2: The root mean square errors (RMSEs) and their standard errors (in
parentheses) of the VRMOM and MOM, and the ratios of RMSEs between VRMOM and
MOM under sample size $N=1000\times(100+1)$, local sample size $n=1000$ and
integer $K=10$. The sample means sent from Byzantine machines are generated
from Gaussian $\mathcal{N}(0,200{\mathbb{I}})$.
$p$ |  | $\alpha_{n}=0$ | $\alpha_{n}=0.05$ | $\alpha_{n}=0.1$ | $\alpha_{n}=0.15$
---|---|---|---|---|---
1 | VRMOM | 0.0025 (0.0018) | 0.0027 (0.0020) | 0.0030 (0.0022) | 0.0032 (0.0024)
MOM | 0.0030 (0.0022) | 0.0031 (0.0023) | 0.0033 (0.0025) | 0.0035 (0.0026)
| Ratio | 0.8613 | 0.8901 | 0.9044 | 0.9129
30 | VRMOM | 0.0175 (0.0022) | 0.0192 (0.0024) | 0.0209 (0.0026) | 0.0227 (0.0028)
MOM | 0.0211 (0.0028) | 0.0223 (0.0028) | 0.0234 (0.0030) | 0.0249 (0.0034)
| Ratio | 0.8285 | 0.8601 | 0.8921 | 0.9108

In the second experiment, we compare the performance of the VRMOM estimator
and the MOM estimator in terms of the root mean square errors and their
standard errors. We fixed the total sample size as $N=1000\times(100+1)$,
local sample size $n=1000$. We let the dimension $p\in\\{1,30\\}$ and the
fraction of Byzantine machines varies from
$\alpha_{n}\in\\{0,0.05,0.1,0.15\\}$. Throughout our experiment, we fix the
number of quantiles $K$ in (7) to be $K=10$.

From Table 2, we observe that VRMOM has smaller root mean square errors than
MOM as all the ratios are greater than 1. With the increase of the fraction of
Byzantine machines, both methods has larger root mean square errors. The
difference between VRMOM and MOM tends to be smaller with more Byzantine
machines. It is interesting to note that, when the dimension $p$ is $30$, the
ratio of the root mean square errors between VRMOM and MOM is smaller than
that when $p=1$. It suggests that the variance reduction effect of our VRMOM
estimator becomes better for higher dimensions, although we have only proved
the superiority when $p=1$ and $2$ in this paper (see Remark 4).

###  4.2 Results for Robust CSL Method

In this section, we consider the linear model and logistic regression model to
demonstrate our robust CSL method.

##### Settings for the linear regression model

For the linear model experiment, the data are generated as follows:

$Y_{i}=\boldsymbol{X}_{i}^{\rm
T}\boldsymbol{\theta}^{\ast}+\epsilon_{i},\qquad i=1,2,\dots,n,$

where each $\boldsymbol{X}_{i}=(X_{i,1},\dots,X_{i,p})^{\rm T}$ is a
$p$-dimensional covariate vector and $(X_{i,1},\dots,X_{i,p})$s are drawn
$i.i.d.$ from a multivariate normal distribution
$\mathcal{N}(0,\bm{\Sigma}_{X})$. The covariance $\bm{\Sigma}_{X}$ is a
symmetric Toeplitz matrix with $\Sigma_{ij}=0.5^{|i-j|}$ for $1\leq i,j\leq
p$, which enforces correlation structure among covariates. We fix the
dimension $p=30$ and generate the entries of the true coefficient vector
$\boldsymbol{\theta}^{\ast}$ to be
$p^{-1/2}(1,(p-2)/(p-1),(p-3)/(p-1),\dots,0)$. Similar to the previous
experiment, the entire sample size is $N=1000\times(100+1)$ and we divide the
data into one master machine ${\mathcal{H}}_{0}$ and $100$ worker machines
$\\{{\mathcal{H}}_{1},\dots,{\mathcal{H}}_{100}\\}$. so that each local sample
size $n=1000$. We consider the standard normal noise distribution where the
noise $\epsilon_{i}\sim\mathcal{N}(0,1)$. For the initial estimator
$\widehat{\boldsymbol{\theta}}^{(0)}$, according to (22), we use the least
square estimator with the data only on the master machine ${\mathcal{H}}_{0}$,
i.e.,
$\widehat{\boldsymbol{\theta}}^{(0)}=(\sum_{i\in{\mathcal{H}}_{0}}\boldsymbol{X}_{i}\boldsymbol{X}_{i}^{\rm
T})^{-1}(\sum_{i\in{\mathcal{H}}_{0}}\boldsymbol{X}_{i}Y_{i}).$ The
computation of the initial estimator is very efficient since it has closed
form and does not require any communication. We also note that for the
surrogate loss minimization problem (21) at $t-$th iteration, we directly
obtain the closed-form solution
$\widehat{\boldsymbol{\theta}}^{(t)}=(2n^{-1}\sum_{i\in{\mathcal{H}}_{0}}\boldsymbol{X}_{i}\boldsymbol{X}_{i}^{\rm
T})^{-1}(2n^{-1}\sum_{i\in{\mathcal{H}}_{0}}\boldsymbol{X}_{i}Y_{i}+\boldsymbol{g}_{0}^{(t-1)}-\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{g}$\kern-1.00006pt}}}\,^{(t-1)})$,
which is also computationally efficient. Moreover, we generate corrupted
gradients sent from Byzantine machines from the following attack models,

  * (a)

Gaussian attack: We replace the gradient vectors in the Byzantine machines by
random vectors in which all the entries are generated from
$\mathcal{N}(0,200{\mathbb{I}})$ independently.

  * (b)

Omniscient attack: For the Byzantine machines, we replace the true gradient
vectors by the scaled negative gradients where the scale constant is extremely
large ($1e10$ in our experiment).

  * (c)

Bit-flip attack: For the Byzantine machines, we replace the true gradient
vectors by flipping the sign fo the first five dimensions.

##### Settings for the logistic regression model

For the logistic regression model experiment, the data are generated from the
following:

$Y_{i}=\begin{cases}\quad 1\quad&\text{with probability
}{\mathcal{L}}(\boldsymbol{X}_{i}^{\rm T}\boldsymbol{\theta}^{\ast}),\\\ \quad
0\quad&\text{with probability }1-{\mathcal{L}}(\boldsymbol{X}_{i}^{\rm
T}\boldsymbol{\theta}^{\ast}),\end{cases}\qquad i=1,2,\dots,n$

where the link function ${\mathcal{L}}(x)={e^{x}}/{(1+e^{x})}$ and each
$\boldsymbol{X}_{i}=(X_{i,1},\dots,X_{i,p})^{\rm T}$ is a $p$-dimensional
covariate vector which is drawn $i.i.d.$ from a multivariate normal
distribution $\mathcal{N}(\bm{\mu}_{x},\bm{\Sigma}_{X})$, which is coincident
with the setting in the linear regression model. We choose
$\bm{\mu}_{x}=(\mu_{x},\dots,\mu_{x})^{\rm T}$ and adopt two settings that
$\mu_{x}=0$ and $\mu_{x}=0.5$. Here $\mu_{x}=0$ corresponds to the balanced
response case that $50\%$ $Y_{i}$s are $1$ and $50\%$ $Y_{i}$s are $0$. And
$\mu_{x}=0.5$ corresponds to the imbalanced response case where $76\%$
$Y_{i}$s are $1$ and $24\%$ $Y_{i}$s are $0$. We also fix the dimension $p=30$
and adopt the same true coefficient vector $\boldsymbol{\theta}^{\ast}$ as
before. For each setting, we repeat 500 independent simulations and report
averaged root mean square estimation errors and standard deviations. For the
initial estimator $\widehat{\boldsymbol{\theta}}^{(0)}$, we use the logistic
regression estimator with the data only on master machine ${\mathcal{H}}_{0}$,
i.e.,
$\widehat{\boldsymbol{\theta}}^{(0)}=\mathop{\mathrm{argmin}}_{\boldsymbol{\theta}\in{\mathbb{R}}^{p}}\left[\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}\left\\{\log\left(1+e^{\boldsymbol{X}_{i}^{\rm
T}\boldsymbol{\theta}}\right)-Y_{i}\boldsymbol{X}_{i}^{\rm
T}\boldsymbol{\theta}\right\\}\right].$ The surrogate loss minimization
problem (21) at $t-$th iteration, i.e.,

$\displaystyle\widehat{\boldsymbol{\theta}}^{(t)}=$
$\displaystyle\mathop{\mathrm{argmin}}_{\boldsymbol{\theta}\in{\mathbb{R}}^{p}}\Big{[}\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}\Big{\\{}\log\left(1+e^{\boldsymbol{X}_{i}^{\rm
T}\boldsymbol{\theta}}\right)\Big{\\}}$
$\displaystyle-\Big{\\{}\Big{(}\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}Y_{i}\boldsymbol{X}_{i}^{\rm
T}\Big{)}+\boldsymbol{g}_{0}^{(t-1)}-\,\hbox{\kern 0.09995pt\vbox{\hrule
height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{g}$\kern-1.00006pt}}}\,^{(t-1)}\Big{\\}}\boldsymbol{\theta}\Big{]},$

can be efficiently solved by standard gradient descent or quasi-Newton
approaches (Nocedal and Wright, 2006) on the center machine without any
communication.

As for the attack mode of Byzantine machines in the logistic regression model,
we simulate the transmitted message in the following way. We replace every
response $Y$ by $1-Y$ and compute the gradients based on these transformed
local data on each Byzantine machine.

Table 3: The root mean square errors (RMSEs) and their standard errors (in
parentheses) of the RCSL and MOM-RCSL, and the ratios of RMSEs between RCSL
and MOM-RCSL under sample size $N=1000\times(100+1)$, local sample size
$n=1000$ and integer $K=10$. The corrupted gradients sent from Byzantine
machines are generated from Gaussian, Omniscient and Bit-flip attacks. The
tolerance parameter for the stopping rule is set to $e_{r}=10^{-4}$.
Attack | None
---|---
$\alpha_{n}$ | $0$
RCSL | 0.0231 (0.0036)
MOM-RCSL | 0.0319 (0.0050)
Ratio | 0.7243
Attack | Gaussian
$\alpha_{n}$ | $0.05$ | $0.1$ | $0.15$
RCSL | 0.0270 (0.0044) | 0.0312 (0.0049) | 0.0351 (0.0060)
MOM-RCSL | 0.0343 (0.0054) | 0.0369 (0.0058) | 0.0398 (0.0063)
Ratio | 0.7863 | 0.8434 | 0.8817
Attack | Omniscient
$\alpha_{n}$ | $0.05$ | $0.1$ | $0.15$
RCSL | 0.0276 (0.0042) | 0.0328 (0.0051) | 0.0396 (0.0060)
MOM-RCSL | 0.0355 (0.0057) | 0.0395 (0.0061) | 0.0449 (0.0069)
Ratio | 0.7774 | 0.8296 | 0.8815
Attack | Bit-flip
$\alpha_{n}$ | $0.05$ | $0.1$ | $0.15$
RCSL | 0.0236 (0.0037) | 0.0242 (0.0039) | 0.0250 (0.0041)
MOM-RCSL | 0.0325 (0.0051) | 0.0334 (0.0053) | 0.0343 (0.0058)
Ratio | 0.7276 | 0.7248 | 0.7291

##### More settings in the simulation

For the fraction of Byzantine machines, we consider the following settings:
(1) $\alpha_{n}=0$ which denotes no Byzantine machine, (2) $\alpha_{n}>0$,
which denotes the existence of Byzantine machine case. We vary the fraction of
Byzantine machines $\alpha_{n}=0.05,0.1,0.15$. We also discuss the stopping
criteria for Algorithm 1. Throughout the experiments, we use the tolerance
parameter $e_{r}=10^{-4}$ as the stopping criterion. In particular, at the
$t$-th iteration of Algorithm 1, we compute
$e=|\widehat{\boldsymbol{\theta}}^{(t)}-\widehat{\boldsymbol{\theta}}^{(t-1)}|_{2}^{2}/|\widehat{\boldsymbol{\theta}}^{(t-1)}|_{2}^{2}$
and stop the algorithm once $e\leq e_{r}$. In our experiments, it only
requires $4$ to $8$ iterations to stop. We also provide the results with
simple fixed number of iterations with $T=5$ and $T=10$.

Since the main focus of this paper is the variance reduction effect of the
proposed VRMOM estimator, in the simulation study, we mainly compare the
performance of our Robust CSL algorithm (RCSL) with the MOM-based Robust CSL
algorithm (MOM-RCSL). More specifically, let
$\widehat{\boldsymbol{\theta}}^{(0)}$ be the initial estimator obtained in the
master machine ${\mathcal{H}}_{0}$, the refined estimator
$\widehat{\boldsymbol{\theta}}^{(1)}_{\rm{MOM}}$ is defined as the solution of
(25) with
$\mathrm{Aggr}(\boldsymbol{g}^{(0)}_{0},\dots,\boldsymbol{g}^{(0)}_{m})$ being
the MOM aggregator of the gradients. Repeat the procedure, and we denote the
$t$-th round MOM-RCSL estimator by
$\widehat{\boldsymbol{\theta}}^{(t)}_{\rm{MOM}}$.

####  4.2.1 Results for linear regression model

The results for linear regression model are presented in Table 3 (for adaptive
stopping criterion) and 4 (for fixed number of iterations). From Table 3, RCSL
has smaller root mean square errors than MOM-RCSL and all the ratios are
greater than 1 for different kinds of attacks. With the increase of Byzantine
fractions $\alpha_{n}$, both methods have larger root mean square errors. The
difference between RCSL and MOM-RCSL tends to be smaller with more Byzantine
machines, which is coincident with the phenomenon of the mean estimation
problem (See Table 2). Among these three attack models, it seems that the
omniscient attacker has the largest root mean square error in general. It is
quite natural since this attacker makes the parameter vector go into the
opposite direction by negative gradients with an extremely large scale factor
(i.e., $1e10$), which exerts the most negative impact on the gradient
aggregation step. On the other hand, even for such a strong attack model, our
RCSL method still performs quite well. It is also interesting to compare Table
2 with the Gaussian attacker part in Table 3 and 4. All simulations share the
same attack mode and the same dimension. It seems that the variance reduction
effect of VRMOM is more significant when it is performed as an iterative
gradient aggregator than as a mean estimator.

Table 4: The root mean square errors (RMSEs) and their standard errors (in
parentheses) of the RCSL and MOM-RCSL, and the ratios of RMSEs between RCSL
and MOM-RCSL under sample size $N=1000\times(100+1)$, local sample size
$n=1000$ and integer $K=10$. The iteration numbers $T$ are fixed to be $5$ and
$10$. The corrupted gradients sent from Byzantine machines are generated from
Gaussian, omniscient and bit-flip attacks.
$T$ | Attack | None
---|---|---
$\alpha_{n}$ | $0$
$5$ | RCSL | 0.0231 (0.0036)
MOM-RCSL | 0.0319 (0.0051)
Ratio | 0.7236
$10$ | RCSL | 0.0231 (0.0036)
MOM-RCSL | 0.0319 (0.0051)
Ratio | 0.7233
$T$ | Attack | Gaussian
$\alpha_{n}$ | $0.05$ | $0.1$ | $0.15$
$5$ | RCSL | 0.0271 (0.0043) | 0.0310 (0.005) | 0.0354 (0.0057)
MOM-RCSL | 0.0343 (0.0056) | 0.0373 (0.006) | 0.0400 (0.0063)
Ratio | 0.7897 | 0.8305 | 0.8859
$10$ | RCSL | 0.0272 (0.0042) | 0.0313 (0.0051) | 0.0348 (0.0058)
MOM-RCSL | 0.0344 (0.0053) | 0.0368 (0.0058) | 0.0398 (0.0065)
Ratio | 0.7905 | 0.8483 | 0.8750
$T$ | Attack | Omniscient
$\alpha_{n}$ | $0.05$ | $0.1$ | $0.15$
$5$ | RCSL | 0.0276 (0.0042) | 0.0328 (0.0052) | 0.0396 (0.0061)
MOM-RCSL | 0.0355 (0.0057) | 0.0395 (0.0061) | 0.045 (0.0069)
Ratio | 0.7768 | 0.8304 | 0.8811
$10$ | RCSL | 0.0276 (0.0042) | 0.0329 (0.0052) | 0.0398 (0.0061)
MOM-RCSL | 0.0355 (0.0057) | 0.0396 (0.0061) | 0.0451 (0.0069)
Ratio | 0.7769 | 0.8311 | 0.8820
$T$ | Attack | Bit-flip
$\alpha_{n}$ | $0.05$ | $0.1$ | $0.15$
$5$ | RCSL | 0.0236 (0.0037) | 0.0242 (0.0039) | 0.0250 (0.0041)
MOM-RCSL | 0.0325 (0.0051) | 0.0334 (0.0053) | 0.0343 (0.0058)
Ratio | 0.7268 | 0.7240 | 0.7281
$10$ | RCSL | 0.0236 (0.0037) | 0.0242 (0.0039) | 0.0250 (0.0041)
MOM-RCSL | 0.0325 (0.0051) | 0.0335 (0.0053) | 0.0344 (0.0058)
Ratio | 0.7259 | 0.7235 | 0.7260

In Table 4, we fix the iteration number to be $5$ and $10$. Table 4 shows
similar patterns as Table 3. There are almost no difference between $T=5$ and
$T=10$ iterations. In other words, it shows that only using a very small
number of iterations (i.e., $T=5$), our RCSL estimator has already converged.
This experiment suggests that the RCSL is communicationally efficient.

####  4.2.2 Results for logistic regression model

Table 5: The root mean square errors (RMSEs) and their standard errors (in
parentheses) of the RCSL and MOM-RCSL, and the ratios of RMSEs between RCSL
and MOM-RCSL under sample size $N=1000\times(100+1)$, local sample size
$n=1000$ and integer $K=10$. The tolerance parameter for the stopping rule is
set to $e_{r}=10^{-4}$. The parameter $\mu_{x}$ controls the class balance.
$\mu_{x}$ | 0
---|---
$\alpha_{n}$ | $0$ | $0.05$ | $0.1$ | $0.15$
RCSL | 0.0504 (0.0075) | 0.0531 (0.0075) | 0.0600 (0.0072) | 0.0701 (0.0075)
MOM-RCSL | 0.0699 (0.0109) | 0.0716 (0.0107) | 0.0765 (0.0108) | 0.0830 (0.0112)
Ratio | 0.7215 | 0.7418 | 0.7844 | 0.8452
$\mu_{x}$ | 0.5
$\alpha_{n}$ | $0$ | $0.05$ | $0.1$ | $0.15$
RCSL | 0.0583 (0.0087) | 0.0601 (0.0087) | 0.0632 (0.0091) | 0.0669 (0.0096)
MOM-RCSL | 0.0830 (0.0135) | 0.0841 (0.0132) | 0.0868 (0.0134) | 0.0905 (0.0141)
Ratio | 0.7024 | 0.7142 | 0.7281 | 0.7395

The results for logistic regression model are presented in Table 5 (for
adaptive stopping criterion) and 6 (for fixed number of iterations). From
Table 5, RCSL has smaller root mean square errors than MOM-RCSL since all the
ratios are greater than 1. With the increase of Byzantine fractions
$\alpha_{n}$, both methods has larger root mean square errors. Moreover, the
RCSL leads to a more significant improvement over the MOM-RCSL for imbalanced
class case. Similar observations can be made from Table 6, where the number of
iterations have been pre-determined. Table 6 also shows that our RCSL is
communicationally efficient since $T=5$ iterations have been sufficient for
the convergence.

Table 6: The root mean square errors (RMSEs) and their standard errors (in
parentheses) of the RCSL and MOM-RCSL, and the ratios of RMSEs between RCSL
and MOM-RCSL under sample size $N=1000\times(100+1)$, local sample size
$n=1000$ and integer $K=10$. The iteration numbers $T$ are fixed to be $5$ and
$10$. The parameter $\mu_{x}$ controls the class balance.
$T$ | $\mu_{x}$ | $0$
---|---|---
$\alpha_{n}$ | $0$ | $0.05$ | $0.1$ | $0.15$
$5$ | RCSL | 0.0505 (0.0075) | 0.0531 (0.0076) | 0.0600 (0.0072) | 0.0702 (0.0075)
MOM-RCSL | 0.0699 (0.0109) | 0.0716 (0.0107) | 0.0765 (0.0108) | 0.0830 (0.0113)
Ratio | 0.7112 | 0.7371 | 0.7901 | 0.8646
$10$ | RCSL | 0.0504 (0.0075) | 0.0531 (0.0075) | 0.0600 (0.0072) | 0.0701 (0.0075)
MOM-RCSL | 0.0699 (0.0109) | 0.0716 (0.0107) | 0.0765 (0.0108) | 0.0830 (0.0112)
Ratio | 0.7218 | 0.7418 | 0.7845 | 0.8452
$T$ | $\mu_{x}$ | $0.5$
$\alpha_{n}$ | $0$ | $0.05$ | $0.1$ | $0.15$
$5$ | RCSL | 0.0583 (0.0087) | 0.0601 (0.0087) | 0.0632 (0.0091) | 0.0670 (0.0096)
MOM-RCSL | 0.0829 (0.0135) | 0.0841 (0.0132) | 0.0868 (0.0134) | 0.0906 (0.0140)
Ratio | 0.7028 | 0.7145 | 0.7281 | 0.7395
$10$ | RCSL | 0.0583 (0.0087) | 0.0601 (0.0087) | 0.0632 (0.0091) | 0.0670 (0.0096)
MOM-RCSL | 0.0830 (0.0135) | 0.0841 (0.0132) | 0.0868 (0.0134) | 0.0905 (0.0141)
Ratio | 0.7025 | 0.7145 | 0.7281 | 0.7395

##  5 Conclusions and Future Work

In this paper, we design a Byzantine tolerant algorithm to address a general
class of estimation and inference problems in a distributed setting. The first
contribution is to improve the statistical efficiency of the widely used
Median-of-Means (MOM) estimator by proposing a new Variance Reduced MOM
(VRMOM) estimator. It achieves a nearly optimal convergence rate upto a
logarithmic factor and has the same order of computational complexity as the
MOM estimator. Inspired by the VRMOM estimator, we further develop the Robust
CSL (RCSL) method for general statistical inference problem. The convergence
rate improves the previous results in Yin et al. (2018) using either median or
trimmed-mean as the gradient aggregator. Moreover, we establish the asymptotic
normality result for our RCSL method.

To highlight our main idea behind VRMOM, we choose to focus on the one-
dimensional case and extend to multi-variate case in a coordinate-wise way.
Therefore, the convergence rate in Theorem 5 has an extra $\sqrt{p}$ in the
term related to Byzantine failures (i.e., $\alpha_{n}\sqrt{p}/\sqrt{n}$). How
to remove this extra $\sqrt{p}$ in a computationally efficient way while still
achieving the same statistical efficiency is still an open problem. As an
important future direction, we will address this challenge by developing new
multivariate efficient aggregators based on the VRMOM estimator. On the other
hand, by adding additional regularizers as in (26), we can address a wider
class of problems in Byzantine setup. We leave these problems in the future
study.

## Acknowledgement

Weidong Liu and Xiaojun Mao are the co-corresponding authors. Weidong Liu’s
research is supported by National Program on Key Basic Research Project (973
Program, 2018AAA0100704), NSFC Grant No. 11825104 and 11690013, Youth Talent
Support Program, and a grant from Australian Research Council. Xiaojun Mao’s
research is supported by NSFC Grant No. 12001109, Shanghai Sailing Program
19YF1402800, Major Research Plan of NSFC Grant No. 92046021, and the Science
and Technology Commission of Shanghai Municipality grant 20dz1200600. The
authors would like to thank the action editor and two anonymous referees for
their constructive comments, which greatly improves the quality of the paper.

## Appendix

The appendix consists of four parts. In Appendix A, we provide detailed proof
for the theoretical results of VRMOM estimator presented in Section 2.2. In
Appendix B, we prove the positive definiteness of
$\boldsymbol{\mathcal{C}}_{\mathrm{MOM}}-\boldsymbol{\mathcal{C}}$ in two-
dimensional case, which have been mentioned in Remark 4 of the main paper. In
Appendix C, we present the theoretical results and some technical assumptions
for the RCSL method. Lastly, in Appendix D, we will show that a large class of
generalized linear models and $M$-estimators suffice the proposed assumptions
in Appendix C.1.1.

##  Appendix A Proof of Theories for VRMOM Estimator

In this appendix, we mainly prove the theoretical results of the proposed
VRMOM estimator in Section 2.2. In Appendix A.1, we introduce several lemmas
which are useful for proofing the main theorems. Next, we present the main
proofs in Appendix A.2.

###  A.1 Technical Lemmas

######  Lemma 1.

(Berry-Esseen Theorem, Theorem 9.1.3 in Chow and Teicher (2012)) If
$\\{X_{i},\;i\geq 1\\}$ are i.i.d. mean-zero random variables with
${\mathbb{E}}|X_{1}|^{2}=\sigma^{2},{\mathbb{E}}|X_{1}|^{2+\kappa}<\infty$,
where $\kappa\in(0,1]$ . Then there exists a constant $C_{\kappa}>0$ such that

$\sup_{-\infty<x<\infty}\left|{\mathbb{P}}\left\\{\sum_{i=1}^{n}X_{i}<x\sigma
n^{1/2}\right\\}-\Phi(x)\right|\leq
C_{\kappa}\frac{{\mathbb{E}}|X_{1}|^{2+\kappa}}{\sigma^{2+\kappa}n^{\kappa/2}}.$

######  Lemma 2.

(Exponential Inequality, Lemma 1 in Cai and Liu (2011)) Let $X_{1},...,X_{n}$
be i.i.d. random variables with zero mean. Suppose that there exist some
$\eta>0$ and $C>0$ such that ${\mathbb{E}}(X_{1}^{2}e^{\eta|X_{1}|})\leq C$.
Then uniformly for $0<x\leq C$ and $n\geq 1$, there is

${\mathbb{P}}\left\\{\frac{1}{n}\sum_{i=1}^{n}X_{i}\geq(\eta+\eta^{-1})x\right\\}\leq\exp\left(-\frac{nx^{2}}{C}\right).$

This lemma will be the workhorse throughout all proofs in this article. For
ease of notations, we will use

${\mathfrak{E}}(\eta,X):={\mathbb{E}}\left(X^{2}e^{\eta|X|}\right),$ (27)

when the expression of $X$ is too complicated.

######  Lemma 3.

Let $X_{1},...,X_{m}$ be i.i.d. random variables with cumulative distribution
function $G(x)={\mathbb{P}}(X_{1}\leq x)$. And there exists some constants
$C>0$ and $\kappa\in(0,1]$ such that there is

$|G(x_{1})-G(x_{2})|\leq C|x_{1}-x_{2}|+\frac{C}{n^{\kappa/2}},$ (28)

holds for any $x_{1},x_{2}\in{\mathbb{R}}$. Further define the function

$\displaystyle Z_{i}(x)={\mathbb{I}}\\{X_{i}\leq x\\}-G(x).$

Let $\delta_{n}=O(1)$ be some rate, then there exists $\widetilde{C}$ large
enough such that

${\mathbb{P}}\left[\frac{1}{m}\sup_{|x|\leq\delta_{n}}\left|\sum_{i=1}^{m}\left\\{Z_{i}(x)-Z_{i}(0)\right\\}\right|\geq\widetilde{C}\left(\sqrt{\frac{\delta_{n}\log
n}{m}}+\frac{1}{n^{\kappa/2}}\right)\right]=O(n^{-\gamma}),$

provided the rates constraints $\log n=O(m)$ and $\max\\{m^{-1}\log
n,n^{-\kappa/2}\\}=O(\delta_{n})$.

###### Proof.

Evenly divide the interval $[-\delta_{n},\delta_{n}]$ into $2n$ pieces and
denote the set
${\mathfrak{N}}=\\{-\delta_{n},-\frac{n-1}{n}\delta_{n},...,\delta_{n}\\}$.
Then

$\displaystyle\frac{1}{m}\sup_{|x|\leq\delta_{n}}\left|\sum_{i=1}^{m}\\{Z_{i}(x)-Z_{i}(0)\\}\right|\leq$
$\displaystyle\max_{\widetilde{x}\in{\mathfrak{N}}}\left|\frac{1}{m}\sum_{i=1}^{m}\\{Z_{i}(\widetilde{x})-Z_{i}(0)\\}\right|$
(29)
$\displaystyle+\max_{\widetilde{x}\in{\mathfrak{N}}}\sup_{\\{x:|\widetilde{x}-x|\leq\delta_{n}/n\\}}\left|\frac{1}{m}\sum_{i=1}^{m}\\{Z_{i}(\widetilde{x})-Z_{i}(x)\\}\right|.$

For the second term, notice that

$\displaystyle\sup_{\\{x:|\widetilde{x}-x|\leq\delta_{n}/n\\}}\left|\frac{1}{m}\sum_{i=1}^{m}\\{Z_{i}(\widetilde{x})-Z_{i}(x)\\}\right|\leq\frac{1}{m}\sum_{i=1}^{m}{\mathbb{I}}\left\\{|X_{i}-\widetilde{x}|\leq\frac{\delta_{n}}{n}\right\\}+C\left(\frac{2\delta_{n}}{n}+\frac{1}{n^{\kappa/2}}\right).$

For every fixed $\widetilde{x}\in{\mathfrak{N}}$, we know

$\displaystyle{\mathbb{E}}\left\\{{\mathbb{I}}\left(|X_{i}-\widetilde{x}|\leq\frac{\delta_{n}}{n}\right)\right\\}={\mathbb{P}}\left(|X_{i}-\widetilde{x}|\leq\frac{\delta_{n}}{n}\right)\leq
C\left(\frac{2\delta_{n}}{n}+\frac{1}{n^{\kappa/2}}\right),$
$\displaystyle\mathfrak{E}\left\\{1,{\mathbb{I}}\left(|X_{i}-\widetilde{x}|\leq\frac{\delta_{n}}{n}\right)-{\mathbb{P}}\left(|X_{i}-\widetilde{x}|\leq\frac{\delta_{n}}{n}\right)\right\\}$
$\displaystyle\leq$ $\displaystyle
e{\mathbb{P}}\left(|X_{i}-\widetilde{x}|\leq\frac{\delta_{n}}{n}\right)\left\\{{\mathbb{P}}\left(|X_{i}-\widetilde{x}|\leq\frac{\delta_{n}}{n}\right)+1\right\\}=O\left(\frac{\log
n}{m}+\frac{1}{n^{\kappa/2}}\right),$

since $\delta_{n}/n=O(n^{-\kappa/2})$ from the rates constraints. Then
applying Lemma 2 we have

$\displaystyle{\mathbb{P}}\left[\max_{\widetilde{x}\in{\mathfrak{N}}}\sup_{\\{x:|\widetilde{x}-x|\leq\delta_{n}/n\\}}\frac{1}{m}\left|\sum_{i=1}^{m}\left\\{Z_{i}(\widetilde{x})-Z_{i}(x)\right\\}\right|\geq
C_{1}\left(\frac{1}{n^{\kappa/2}}+\frac{\log
n}{m}\right)+C\left(\frac{2\delta_{n}}{n}+\frac{1}{n^{\kappa/2}}\right)\right]$
(30) $\displaystyle\leq$ $\displaystyle
2n\max_{\widetilde{x}\in{\mathfrak{N}}}{\mathbb{P}}\left\\{\frac{1}{m}\sum_{i=1}^{m}{\mathbb{I}}\left(|X_{i}-\widetilde{x}|\leq\frac{\delta_{n}}{n}\right)-{\mathbb{P}}\left(|X_{1}-\widetilde{x}|\leq\frac{\delta_{n}}{n}\right)\geq
C_{1}\left(\frac{\log n}{m}+\frac{1}{n^{\kappa/2}}\right)\right\\}$
$\displaystyle=$ $\displaystyle O(n^{-\gamma}),$

by letting $C_{1}$ large enough. Now for the first term in (29), again we
apply Lemma 2 with

$\displaystyle{\mathbb{E}}[Z_{i}(x)-Z_{i}(0)]=0,$
$\displaystyle\sup_{-\delta_{n}\leq x\leq\delta_{n}}$
$\displaystyle{\mathbb{E}}\left[\\{Z_{i}(x)-Z_{i}(0)\\}^{2}e^{|Z_{i}(x)-Z_{i}(0)|}\right]=O(\delta_{n}),$

there is

$\displaystyle{\mathbb{P}}\left[\max_{\widetilde{x}\in{\mathfrak{N}}}\left|\frac{1}{m}\sum_{i=1}^{m}\\{Z_{i}(\widetilde{x})-Z_{i}(0)\\}\right|\geq
C_{2}\sqrt{\frac{\delta_{n}\log n}{m}}\right]$ (31) $\displaystyle\leq$
$\displaystyle
2n\max_{\widetilde{x}\in{\mathfrak{N}}}{\mathbb{P}}\left[\frac{1}{m}\left|\sum_{i=1}^{m}\\{Z_{i}(\widetilde{x})-Z_{i}(0)\\}\right|\geq
C_{2}\sqrt{\frac{\delta_{n}\log n}{m}}\right]=O(n^{-\gamma}),$

for some $C_{2}$ large enough. Then the lemma is proved by combining (30) and
(31). ∎

######  Lemma 4.

(Concentration of median-of-means with Byzantine machines) Let $N$($=(m+1)n$)
i.i.d. random variables $X_{1},...,X_{N}$ evenly distributed in $m+1$ subsets
${\mathcal{H}}_{0},...,{\mathcal{H}}_{m}$. There is a subset of index
${\mathcal{B}}\subseteq\\{1,\dots,m\\}$ with
${\mathrm{C}ard}({\mathcal{B}})=\lfloor\alpha_{n}m\rfloor$. Define

$\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}=\begin{cases}\frac{1}{n}\sum_{i\in{\mathcal{H}}_{j}}X_{i}\quad&j\notin{\mathcal{B}},\\\
\quad*\quad&j\in{\mathcal{B}},\end{cases}\qquad\widehat{X}=\text{med}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\mid 0\leq j\leq
m\right).$

Suppose
${\mathbb{E}}(X_{1})=0,\mathrm{Var}(X_{1})=\sigma^{2},{\mathbb{E}}|X_{1}|^{2+\kappa}<\infty$,
where $\kappa\in(0,1]$. The fraction $\alpha_{n}$ satisfies $\alpha_{n}\leq
1/2-\delta$ for some $\delta>0$. Then for every $\gamma>1$, there exists
$\widetilde{C}>0$ large enough, such that

${\mathbb{P}}\left\\{\left|\widehat{X}\right|\geq\frac{\widetilde{C}}{\sqrt{n}}\left(\alpha_{n}+\frac{1}{n^{\kappa/2}}+\sqrt{\frac{\log
n}{m}}\right)\right\\}=O(n^{-\gamma}).$

###### Proof.

By definition of medians, for any $x>0$, there is

$\displaystyle{\mathbb{P}}\left(\widehat{X}\geq
x\right)={\mathbb{P}}\left\\{\sum_{j=0}^{m}{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}<x\right)\leq\frac{m+1}{2}\right\\}\leq{\mathbb{P}}\left\\{\sum_{j\notin{\mathcal{B}}}{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}<x\right)\leq\frac{m+1}{2}\right\\}$
$\displaystyle=$
$\displaystyle{\mathbb{P}}\left\\{\frac{1}{(1-\alpha_{n})(m+1)}\sum_{j\notin{\mathcal{B}}}{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}<x\right)-{\mathbb{P}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{1}<x\right)\leq\frac{1}{2}+\frac{\alpha_{n}}{2(1-\alpha_{n})}-{\mathbb{P}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{1}<x\right)\right\\}$
$\displaystyle\leq$
$\displaystyle{\mathbb{P}}\left\\{\frac{1}{(1-\alpha_{n})(m+1)}\sum_{j\notin{\mathcal{B}}}{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}<x\right)-{\mathbb{P}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{1}<x\right)\leq\frac{1}{2}+\frac{\alpha_{n}}{2(1-\alpha_{n})}+C_{\kappa}\frac{{\mathbb{E}}|X_{1}|^{2+\kappa}}{\sigma^{2+\kappa}n^{\kappa/2}}-\Phi\left(\frac{\sqrt{n}x}{\sigma}\right)\right\\}$
$\displaystyle=$
$\displaystyle{\mathbb{P}}\left\\{\frac{1}{(1-\alpha_{n})(m+1)}\sum_{j\notin{\mathcal{B}}}{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}<x\right)-{\mathbb{P}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{1}<x\right)\leq-\sqrt{\frac{c\log
n}{(1-\alpha_{n})(m+1)}}\right\\},$

where the last line uses Berry-Esseen Theorem (Lemma 1), and $x$ is given by

$x=\frac{\sigma}{\sqrt{n}}\Phi^{-1}\left(\frac{1}{2}+\frac{\alpha_{n}}{2(1-\alpha_{n})}+C_{\kappa}\frac{{\mathbb{E}}|X_{1}|^{2+\kappa}}{\sigma^{2+\kappa}n^{\kappa/2}}+\sqrt{\frac{c\log
n}{(1-\alpha_{n})(m+1)}}\right).$

Now apply Lemma 2 for the i.i.d. sequence ${\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}<x\right)-{\mathbb{P}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{1}<x\right)$ with

$\mathfrak{E}\left\\{1,{\mathbb{I}}\left(\,\hbox{\kern 0.09995pt\vbox{\hrule
height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}<x\right)-{\mathbb{P}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{1}<x\right)\right\\}\leq
e,$

we have

${\mathbb{P}}\left\\{\frac{1}{(1-\alpha_{n})(m+1)}\sum_{j\notin{\mathcal{B}}}{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}<x\right)-{\mathbb{P}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{1}<x\right)\leq-\sqrt{\frac{c\log
n}{(1-\alpha_{n})(m+1)}}\right\\}=O(n^{-\gamma}),$

with $c=4e\gamma$. Moreover, we have the following elementary facts

$\Phi^{-1}(x_{0})=\Phi^{-1}(x_{0})-\Phi^{-1}\left(1/2\right)\leq\left(x_{0}-1/2\right)\left(\Phi^{-1}\right)^{\prime}(x_{0})\leq\frac{x_{0}-1/2}{\psi\left\\{\Phi^{-1}(1-\delta/2)\right\\}},$

holds for any $1/2\leq x_{0}<1-\delta/2$. Also we know
$1/2\leq\Phi(\sqrt{n}x/\sigma)<1-\delta/2$ holds for $m,n$ sufficiently large.
Denote $C_{\delta}=1/\psi\\{\Phi^{-1}(1-\delta/2)\\}$, then there is

$\displaystyle{\mathbb{P}}\left\\{\widehat{X}\geq\frac{\sigma
C_{\delta}}{\sqrt{n}}\left(\frac{\alpha_{n}}{2(1-\alpha_{n})}+C_{\kappa}\frac{{\mathbb{E}}|X_{1}|^{2+\kappa}}{\sigma^{2+\kappa}n^{\kappa/2}}+\sqrt{\frac{c\log
n}{(1-\alpha_{n})(m+1)}}\right)\right\\}\leq{\mathbb{P}}\left(\widehat{X}\geq
x\right)$ $\displaystyle\leq$
$\displaystyle{\mathbb{P}}\left\\{\frac{1}{(1-\alpha_{n})(m+1)}\sum_{j\notin{\mathcal{B}}}{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}<x\right)-{\mathbb{P}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{1}<x\right)\leq-\sqrt{\frac{c\log
n}{(1-\alpha_{n})(m+1)}}\right\\}\leq n^{-\gamma}.$

Now do the same thing for ${\mathbb{P}}(\widehat{X}\leq-x)$, we finally get
the desired result. ∎

###  A.2 Proofs of the results in Section 2.2

We firstly provide the proofs related to the univariate VRMOM estimator.

###### Proof of Theorem 1, Theorem 2, and Theorem 3.

Denote $\kappa_{1}=\min\\{1,\kappa\\}$ and $\kappa_{2}=\min\\{2,\kappa\\}$. To
prove this result, let’s firstly give a convergence rate for sample variance
$\widehat{\sigma}-\sigma$. From the moment bound
${\mathbb{E}}|X_{1}-\mu|^{2+\kappa}<\infty$, by Marcinkiewicz-Zygmund theorem
(Theorem 5.2.2 in Chow and Teicher (2012)), we have

$\displaystyle\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}(X_{i}-\mu)^{2}-\sigma^{2}=o_{{\mathbb{P}}}\left(\frac{1}{n^{\kappa_{2}/(2+\kappa_{2})}}\right),\quad\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}(X_{i}-\mu)=o_{{\mathbb{P}}}\left(\frac{1}{n^{\kappa_{2}/(2+\kappa_{2})}}\right).$

Therefore we have

$\displaystyle\left|\widehat{\sigma}-\sigma\right|=$
$\displaystyle\frac{1}{\widehat{\sigma}+\sigma}\left|\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}\left(X_{i}-\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{0}\right)^{2}-\sigma^{2}\right|$
(32) $\displaystyle\leq$
$\displaystyle\frac{1}{\sigma}\left|\frac{1}{n}\sum_{i\in{\mathcal{H}}_{0}}(X_{i}-\mu)^{2}-\sigma^{2}-\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{0}-\mu\right)^{2}\right|=o_{{\mathbb{P}}}\left(\frac{1}{n^{\kappa_{2}/(2+\kappa_{2})}}\right).$

Next pick out one summand in (6). Denote

$\displaystyle G_{j}(x):={\mathbb{P}}\left\\{\frac{\sqrt{n}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}-\mu)}{\sigma}\leq
x\right\\},\quad I_{j}(x):={\mathbb{I}}\left\\{\frac{\sqrt{n}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}-\mu)}{\sigma}\leq
x\right\\},$

then

$\displaystyle{\mathbb{I}}\left\\{\,\hbox{\kern 0.09995pt\vbox{\hrule
height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\leq\widehat{\mu}+\frac{\widehat{\sigma}\Delta_{k}}{\sqrt{n}}\right\\}-\frac{k}{K+1}$
$\displaystyle=$ $\displaystyle
I_{j}\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}}{\sigma}\Delta_{k}\right\\}-G_{j}\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}}{\sigma}\Delta_{k}\right\\}+\underbrace{G_{j}\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}}{\sigma}\Delta_{k}\right\\}-\Phi(\Delta_{k})}_{T}.$

For the term $T$,

$\displaystyle T=$ $\displaystyle
G_{j}\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}}{\sigma}\Delta_{k}\right\\}-\Phi\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}}{\sigma}\Delta_{k}\right\\}+\Phi\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}}{\sigma}\Delta_{k}\right\\}-\Phi(\Delta_{k})$
(33) $\displaystyle=$ $\displaystyle
O(n^{-\kappa_{1}/2})+\psi(\Delta_{k})\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\psi(\Delta_{k})\frac{\widehat{\sigma}-\sigma}{\sigma}\Delta_{k}+O\left(\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}-\sigma}{\sigma}\Delta_{k}\right)^{2}$
$\displaystyle=$
$\displaystyle\psi(\Delta_{k})\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+O_{{\mathbb{P}}}\left(\alpha_{n}^{2}+\frac{\log
n}{m}+\frac{1}{n^{\kappa_{2}/(2+\kappa_{2})}}\right),$

where line 1 to line 2 uses Berry-Esseen Inequality (Lemma 1) and Taylor
expansion, line 2 to line 3 uses (32) and concentration inequalities for
$\widehat{\mu}$ (Lemma 4). By using Lemma 4 and (32) again we have that

$\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}-\sigma}{\sigma}\Delta_{k}=O_{{\mathbb{P}}}\left(\alpha_{n}+\sqrt{\frac{\log
n}{m}}+\frac{1}{n^{\kappa_{2}/(2+\kappa_{2})}}\right).$

Also note that Berry-Esseen Inequality guarantees that
$\sqrt{n}(Y_{j}-\mu)/\sigma$ satisfies

$|G_{j}(x_{1})-G_{j}(x_{2})|\leq\psi(0)|x_{1}-x_{2}|+2C_{\kappa_{1}}\frac{{\mathbb{E}}|X_{1}-\mu|^{2+\kappa_{1}}}{\sigma^{2+\kappa_{1}}n^{\kappa_{1}/2}}.$

So we can apply Lemma 3 with $\delta_{n}=O(\alpha_{n}+\sqrt{\frac{\log
n}{m}}+\frac{1}{n^{\kappa_{2}/(2+\kappa_{2})}})$, which yields

$\displaystyle\frac{1}{m+1}\left|\sum_{j\notin{\mathcal{B}}}\left[I_{j}\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}}{\sigma}\Delta_{k}\right\\}-G_{j}\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}}{\sigma}\Delta_{k}\right\\}-I_{j}(\Delta_{k})+G_{j}(\Delta_{k})\right]\right|$
$\displaystyle\leq$
$\displaystyle\frac{1}{m+1}\sup_{|x-\Delta_{k}|\leq\delta_{n}}\left|\sum_{j\notin{\mathcal{B}}}\left\\{I_{j}(x+\Delta_{k})-G_{j}(x+\Delta_{k})-I_{j}(\Delta_{k})+G_{j}(\Delta_{k})\right\\}\right|$
$\displaystyle=$ $\displaystyle
O_{{\mathbb{P}}}\left(\sqrt{\frac{\delta_{n}\log
n}{m}}+\frac{1}{n^{\kappa_{1}/2}}\right)$ $\displaystyle=$ $\displaystyle
O_{{\mathbb{P}}}\left(\sqrt{\frac{\alpha_{n}\log n}{m}}+\Big{(}\frac{\log
n}{m}\Big{)}^{3/4}+\frac{\log^{1/2}n}{m^{1/2}n^{\kappa_{2}/(4+2\kappa_{2})}}+\frac{1}{n^{\kappa_{1}/2}}\right).$

Thus it implies that

$\displaystyle\frac{1}{m+1}\sum_{j=0}^{m}\left[I_{j}\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}}{\sigma}\Delta_{k}\right\\}-G_{j}\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}}{\sigma}\Delta_{k}\right\\}\right]$
(34) $\displaystyle=$
$\displaystyle\frac{1}{m+1}\sum_{j=0}^{m}\left\\{I_{j}(\Delta_{k})-G_{j}(\Delta_{k})\right\\}$
$\displaystyle+O_{{\mathbb{P}}}\left(\sqrt{\frac{\alpha_{n}\log
n}{m}}+\Big{(}\frac{\log
n}{m}\Big{)}^{3/4}+\frac{\log^{1/2}n}{m^{1/2}n^{\kappa_{2}/(4+2\kappa_{2})}}+\frac{1}{n^{\kappa_{1}/2}}\right).$

Again from (32) we have

$\displaystyle\left|\frac{\widehat{\sigma}-\sigma}{(m+1)\sqrt{n}\sum_{k=1}^{K}\psi(\Delta_{k})}\sum_{k=1}^{K}\sum_{j=0}^{m}\left[{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\leq\widehat{\mu}+\frac{\widehat{\sigma}\Delta_{k}}{\sqrt{n}}\right)-\frac{k}{K+1}\right]\right|$
(35) $\displaystyle\leq$
$\displaystyle\frac{K|\widehat{\sigma}-\sigma|}{\sqrt{n}\sum_{k=1}^{K}\psi(\Delta_{k})}=O_{{\mathbb{P}}}\left(\frac{1}{n^{(3\kappa_{2}+2)/(2\kappa_{2}+4)}}\right).$

Combining equations (33) (34) (35), we have

$\displaystyle\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,-\mu=$
$\displaystyle\widehat{\mu}-\mu-\frac{\widehat{\sigma}}{(m+1)\sqrt{n}\sum_{k=1}^{K}\psi(\Delta_{k})}\sum_{k=1}^{K}\sum_{j=0}^{m}\left\\{{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\leq\widehat{\mu}+\frac{\widehat{\sigma}\Delta_{k}}{\sqrt{n}}\right)-\frac{k}{K+1}\right\\}$
(36) $\displaystyle=$
$\displaystyle\widehat{\mu}-\mu-\frac{\sigma}{(m+1)\sqrt{n}\sum_{k=1}^{K}\psi(\Delta_{k})}\sum_{k=1}^{K}\sum_{j\notin{\mathcal{B}}}\left\\{{\mathbb{I}}\left(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\leq\widehat{\mu}+\frac{\widehat{\sigma}\Delta_{k}}{\sqrt{n}}\right)-\frac{k}{K+1}\right\\}$
$\displaystyle+O_{{\mathbb{P}}}\left(\frac{\alpha_{n}}{\sqrt{n}}+\frac{1}{n^{(3\kappa_{2}+2)/(2\kappa_{2}+4)}}\right)$
$\displaystyle=$
$\displaystyle\frac{\sigma}{\sqrt{n}\sum_{k=1}^{K}\psi(\Delta_{k})}\sum_{k=1}^{K}\frac{1}{m+1}\sum_{j\notin{\mathcal{B}}}\left[G_{j}\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}}{\sigma}\Delta_{k}\right\\}-I_{j}\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}+\frac{\widehat{\sigma}}{\sigma}\Delta_{k}\right\\}\right]$
$\displaystyle+O_{{\mathbb{P}}}\left(\frac{\alpha_{n}}{\sqrt{n}}+\frac{1}{n^{(3\kappa_{2}+2)/(2\kappa_{2}+4)}}+\frac{\log
n}{m\sqrt{n}}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{m+1}\sum_{j\notin{\mathcal{B}}}\frac{\sigma}{\sqrt{n}\sum_{k=1}^{K}\psi(\Delta_{k})}\sum_{k=1}^{K}\\{G_{j}(\Delta_{k})-I_{j}(\Delta_{k})\\}$
$\displaystyle+O_{{\mathbb{P}}}\left(\frac{\alpha_{n}}{\sqrt{n}}+\frac{\log
n}{m\sqrt{n}}+\frac{1}{n^{(3\kappa_{2}+2)/(2\kappa_{2}+4)}}+\frac{\log^{3/4}n}{n^{1/2}m^{3/4}}\right).$

From central limit theorem, we have

$\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,-\mu=O_{{\mathbb{P}}}\left(\frac{\alpha_{n}}{\sqrt{n}}+\frac{1}{\sqrt{mn}}+\frac{1}{n^{(3\kappa_{2}+2)/(2\kappa_{2}+4)}}+\frac{\log^{3/4}n}{n^{1/2}m^{3/4}}\right),$
(37)

which proves Theorem 2. Moreover, when
$m=o(n^{2\kappa_{2}/(\kappa_{2}+2)})=o(\min\\{n^{2\kappa/(\kappa+2)},n\\}),\log^{3}n=o(m)$
and $\alpha_{n}=o(1/\sqrt{m})$, the remainder in (36) will be of the order
$o_{{\mathbb{P}}}(\frac{1}{\sqrt{mn}})$, while the major term becomes a
summation of i.i.d. sequence with

$\displaystyle{\mathbb{E}}\left[\sum_{k=1}^{K}\left\\{G_{j}(\Delta_{k})-I_{j}(\Delta_{k})\right\\}\right]=0,\>\text{
and}$
$\displaystyle\text{Var}\left[\sum_{k=1}^{K}\left\\{G_{j}(\Delta_{k})-I_{j}(\Delta_{k})\right\\}\right]$
$\displaystyle=$
$\displaystyle\sum_{k_{1},k_{2}=1}^{K}G_{j}\left\\{\min(\Delta_{k_{1}},\Delta_{k_{2}})\right\\}-G_{j}(\Delta_{k_{1}})G_{j}(\Delta_{k_{2}})$
$\displaystyle=$
$\displaystyle\sum_{k_{1},k_{2}=1}^{K}\min(\tau_{k_{1}},\tau_{k_{2}})\\{1-\max(\tau_{k_{1}},\tau_{k_{2}})\\}+O(n^{-1/2}),$

where the last line again follows from Berry-Esseen theorem (1). Applying
standard central limit theorem we have

$\sqrt{N}(\,\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,-\mu)\xrightarrow{d}\mathcal{N}(0,\sigma_{K}^{2}),$

where $\sigma^{2}_{K}$, as defined in (9), tends to $\pi\sigma^{2}/3$
according to Lemma 6 in Appendix B. Therefore Theorem 1 is proved.

To prove Theorem 3, we just apply (37) to each coordinate and obtain that

$|\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{\mu}$\kern-1.00006pt}}}\,-\boldsymbol{\mu}|_{2}=\sqrt{\sum_{l=1}^{p}|\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,_{l}-\mu_{l}|^{2}}=O_{{\mathbb{P}}}\left(\frac{\alpha_{n}\sqrt{p}}{\sqrt{n}}+\sqrt{\frac{p}{mn}}+\frac{\sqrt{p}}{n^{(3\kappa_{2}+2)/(2\kappa_{2}+4)}}+\frac{p^{1/2}\log^{3/4}n}{n^{1/2}m^{3/4}}\right).$

∎

###### Proof of Theorem 4.

To prove asymptotic normality of multi-dimensional VRMOM estimator, using (36)
and the rate constraints, for each coordinate $l$ (where $1\leq l\leq p$),
there is

$\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\mu$\kern-1.00006pt}}}\,_{l}-\mu_{l}=\frac{1}{m+1}\sum_{j\notin{\mathcal{B}}}\frac{\sqrt{\sigma_{l,l}}}{\sqrt{n}\sum_{k=1}^{K}\psi(\Delta_{k})}\sum_{k=1}^{K}\\{G_{j,l}(\Delta_{k})-I_{j,l}(\Delta_{k})\\}+o_{{\mathbb{P}}}\left(\frac{1}{\sqrt{pmn}}\right),$

where

$\displaystyle G_{j,l}(x):={\mathbb{P}}\left\\{\frac{\sqrt{n}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j,l}-\mu_{l})}{\sqrt{\sigma_{l,l}}}\leq
x\right\\},\quad I_{j,l}(x):={\mathbb{I}}\left\\{\frac{\sqrt{n}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j,l}-\mu_{l})}{\sqrt{\sigma_{l,l}}}\leq
x\right\\}.$

For any vector $|\boldsymbol{v}|_{2}=1$, there is

$\displaystyle\left\langle\,\hbox{\kern 0.09995pt\vbox{\hrule
height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{\mu}$\kern-1.00006pt}}}\,-\boldsymbol{\mu},\boldsymbol{v}\right\rangle=\frac{1}{m+1}\sum_{j\notin{\mathcal{B}}}\frac{1}{\sqrt{n}\sum_{k=1}^{K}\psi(\Delta_{k})}\sum_{k=1}^{K}\sum_{l=1}^{p}\sqrt{\sigma_{l,l}}\big{\\{}I_{j,l}(\Delta_{k})-G_{j,l}(\Delta_{k})\big{\\}}v_{l}+o_{{\mathbb{P}}}\left(\frac{1}{\sqrt{mn}}\right).$

Now we apply central limit theorem and yield

$\displaystyle\frac{\sqrt{(m+1)n}}{\widetilde{\sigma}_{\boldsymbol{v}}}\big{\langle}\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{\mu}$\kern-1.00006pt}}}\,-\boldsymbol{\mu},\boldsymbol{v}\big{\rangle}\xrightarrow{d}{\mathcal{N}}(0,1),$
$\displaystyle\text{where
}\quad\widetilde{\sigma}_{\boldsymbol{v}}^{2}=\boldsymbol{v}^{\rm
T}\widetilde{\boldsymbol{\mathcal{C}}}\boldsymbol{v}.$

Here $\widetilde{\boldsymbol{\mathcal{C}}}\in{\mathbb{R}}^{p\times p}$ has its
$(l_{1},l_{2})$-entry defined as

$\displaystyle\widetilde{\mathcal{C}}_{l_{1},l_{2}}=$
$\displaystyle\frac{\sqrt{\sigma_{l_{1},l_{1}}\sigma_{l_{2},l_{2}}}}{\\{\sum_{k=1}^{K}\psi(\Delta_{k})\\}^{2}}{\mathbb{E}}\left[\sum_{k=1}^{K}\\{I_{0,l_{1}}(\Delta_{k})-G_{0,l_{1}}(\Delta_{k})\\}\sum_{k=1}^{K}\\{I_{0,l_{2}}(\Delta_{k})-G_{0,l_{2}}(\Delta_{k})\\}\right]$
$\displaystyle=$
$\displaystyle\frac{\sqrt{\sigma_{l_{1},l_{1}}\sigma_{l_{2},l_{2}}}}{\\{\sum_{k=1}^{K}\psi(\Delta_{k})\\}^{2}}\sum_{k_{1},k_{2}}\left\\{{\mathbb{P}}\left(\frac{\sqrt{n}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j,l_{1}}-\mu_{l_{1}})}{\sqrt{\sigma_{l_{1},l_{1}}}}\leq\Delta_{k_{1}},\frac{\sqrt{n}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j,l_{2}}-\mu_{l_{2}})}{\sqrt{\sigma_{l_{2},l_{2}}}}\leq\Delta_{k_{2}}\right)-G_{0,l_{1}}(\Delta_{k_{1}})G_{0,l_{2}}(\Delta_{k_{2}})\right\\},$

Moreover, we can apply multivariate Berry-Esseen theorem (See Theorem 1.3 in
Götze (1991)) and give

$\displaystyle\widetilde{\mathcal{C}}_{l_{1},l_{2}}$
$\displaystyle=\mathcal{C}_{l_{1},l_{2}}+O(n^{-1/2}),$

where $\mathcal{C}_{l_{1},l_{2}}$ is defined in (14). Thus the theorem is
proved. ∎

###### Proof of Proposition 1.

For the asymptotic normality of multi-dimensional MOM estimator, together with
Lemma 5 below and the rate constraint, we can show that

$\widehat{\mu}_{\mathrm{MOM},l}-\mu_{l}=\frac{\sqrt{2\pi\sigma_{l,l}}}{(m+1)\sqrt{n}}\sum_{j\notin{\mathcal{B}}}\left\\{G_{j,l}(0)-I_{j,l}(0)\right\\}+o_{{\mathbb{P}}}\left(\frac{1}{\sqrt{pmn}}\right).$

For any vector $|\boldsymbol{v}|_{2}=1$, there is

$\displaystyle\left\langle\widehat{\boldsymbol{\mu}}_{\mathrm{MOM}}-\boldsymbol{\mu},\boldsymbol{v}\right\rangle=\frac{\sqrt{2\pi}}{(m+1)\sqrt{n}}\sum_{j\notin{\mathcal{B}}}\sum_{l=1}^{p}\sqrt{\sigma_{l,l}}\big{\\{}I_{j,l}(\Delta_{k})-G_{j,l}(\Delta_{k})\big{\\}}v_{l}+o_{{\mathbb{P}}}\left(\frac{1}{\sqrt{mn}}\right).$

Now we apply central limit theorem and yield

$\displaystyle\frac{\sqrt{(m+1)n}}{\widetilde{\sigma}_{\boldsymbol{v}}}\big{\langle}\widehat{\boldsymbol{\mu}}_{\mathrm{MOM}}-\boldsymbol{\mu},\boldsymbol{v}\big{\rangle}\xrightarrow{d}{\mathcal{N}}(0,1),$
$\displaystyle\text{where
}\quad\widetilde{\sigma}_{\boldsymbol{v}}^{2}=\boldsymbol{v}^{\rm
T}\widetilde{\boldsymbol{\mathcal{C}}}_{\mathrm{MOM}}\boldsymbol{v}.$

Here
$\widetilde{\boldsymbol{\mathcal{C}}}_{\mathrm{MOM}}\in{\mathbb{R}}^{p\times
p}$ has its $(l_{1},l_{2})$-entry defined as

$\displaystyle\widetilde{\mathcal{C}}_{\mathrm{MOM},l_{1},l_{2}}=$
$\displaystyle
2\pi\sqrt{\sigma_{l_{1},l_{1}}\sigma_{l_{2},l_{2}}}\;{\mathbb{E}}\left[\\{I_{0,l_{1}}(0)-G_{0,l_{1}}(0)\\}\\{I_{0,l_{2}}(0)-G_{0,l_{2}}(0)\\}\right]$
$\displaystyle=$ $\displaystyle
2\pi\sqrt{\sigma_{l_{1},l_{1}}\sigma_{l_{2},l_{2}}}\left\\{{\mathbb{P}}\left(\frac{\sqrt{n}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j,l_{1}}-\mu_{l_{1}})}{\sqrt{\sigma_{l_{1},l_{1}}}}\leq
0,\frac{\sqrt{n}(\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j,l_{2}}-\mu_{l_{2}})}{\sqrt{\sigma_{l_{2},l_{2}}}}\leq
0\right)-G_{0,l_{1}}(0)G_{0,l_{2}}(0)\right\\},$

Moreover, we can apply multivariate Berry-Esseen theorem (See Theorem 1.3 in
Götze (1991)) and give

$\displaystyle\widetilde{\mathcal{C}}_{\mathrm{MOM},l_{1},l_{2}}$
$\displaystyle=\mathcal{C}_{\mathrm{MOM},l_{1},l_{2}}+O(n^{-1/2}),$

where $\mathcal{C}_{\mathrm{MOM},l_{1},l_{2}}$ is defined in (17). Thus the
proposition is proved. ∎

######  Lemma 5.

Let $N=(m+1)n$ i.i.d. random variables $X_{1},...,X_{N}$ evenly distributed in
$(m+1)$ subsets ${\mathcal{H}}_{0},\dots,{\mathcal{H}}_{m}$. There is a subset
of index ${\mathcal{B}}\subset\\{1,\dots,m\\}$ with
${\mathrm{C}ard}({\mathcal{B}})=\lfloor\alpha_{n}m\rfloor$. Define

$\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}=\begin{cases}\frac{1}{n}\sum_{i\in{\mathcal{H}}_{j}}X_{i}\quad&j\notin{\mathcal{B}},\\\
*\quad&j\in{\mathcal{B}},\end{cases}\quad\widehat{\mu}=\mathrm{med}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\mid 0\leq j\leq m).$

Suppose $X_{1}$ satisfies
${\mathbb{E}}(X_{1})=\mu,\mathrm{Var}(X_{1})=\sigma^{2}$, and
${\mathbb{E}}|X_{1}-\mu|^{3}<\infty$. The fraction $\alpha_{n}$ satisfies
$\alpha_{n}\leq 1/2-\delta$ for some $\delta>0$. Then $\widehat{\mu}$ admits
the following representation:

$\widehat{\mu}=\mu-\frac{\sqrt{2\pi}\sigma}{(m+1)\sqrt{n}}\sum_{j\notin{\mathcal{B}}}\left\\{{\mathbb{I}}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\leq\mu)-{\mathbb{P}}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\leq\mu)\right\\}+O_{{\mathbb{P}}}\left(\frac{\alpha_{n}}{\sqrt{n}}+\frac{\log
n}{m\sqrt{n}}+\frac{\log^{3/4}n}{m^{3/4}n^{1/2}}+\frac{1}{n}\right).$

###### Proof.

Using Taylor expansion we have

$\displaystyle\Phi\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}\right\\}-\Phi(0)=$
$\displaystyle\frac{\sqrt{n}}{\sqrt{2\pi}\sigma}(\widehat{\mu}-\mu)+O(n(\widehat{\mu}-\mu)^{2})$
$\displaystyle=$
$\displaystyle\frac{\sqrt{n}}{\sqrt{2\pi}\sigma}(\widehat{\mu}-\mu)+O_{{\mathbb{P}}}\left(\alpha_{n}^{2}+\frac{\log
n}{m}\right),$ (38)

since
$\widehat{\mu}-\mu=O_{{\mathbb{P}}}(\frac{\alpha_{n}}{\sqrt{n}}+\sqrt{\frac{\log
n}{mn}})$ by Lemma 4. On the other hand, denote

$G(x)={\mathbb{P}}\left\\{\frac{\sqrt{n}(\,\hbox{\kern 0.09995pt\vbox{\hrule
height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{1}-\mu)}{\sigma}\leq
x\right\\},\quad I_{j}(x)={\mathbb{I}}\left\\{\frac{\sqrt{n}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}-\mu)}{\sigma}\leq
x\right\\}.$

From Berry-Essen theorem (Lemma 1) we know Lemma 3 is applicable and yields

$\displaystyle\frac{1}{m+1}\sum_{j=0}^{m}I_{j}\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}\right\\}-G\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}\right\\}-\frac{1}{m+1}\sum_{j=0}^{m}I_{j}(0)+G(0)$
$\displaystyle=$ $\displaystyle
O_{{\mathbb{P}}}\left(\alpha_{n}+\frac{\log^{3/4}n}{m^{3/4}}+\frac{1}{\sqrt{n}}\right).$

By definition of median we know

$\Phi(0)=\frac{1}{2}=\frac{1}{m+1}\sum_{j=0}^{m}I_{j}\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}\right\\}+O\left(\frac{1}{m}\right).$

Thus

$\displaystyle\Phi\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}\right\\}-\Phi(0)=$
$\displaystyle
G\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}\right\\}-\frac{1}{m+1}\sum_{j=0}^{m}I_{j}\left\\{\frac{\sqrt{n}(\widehat{\mu}-\mu)}{\sigma}\right\\}+O\left(\frac{1}{m}+\frac{1}{\sqrt{n}}\right)$
$\displaystyle=$
$\displaystyle-\frac{1}{m+1}\sum_{j\notin{\mathcal{B}}}\\{I_{j}(0)-G(0)\\}+O_{{\mathbb{P}}}\left(\alpha_{n}+\frac{\log^{3/4}n}{m^{3/4}}+\frac{1}{\sqrt{n}}\right).$
(39)

Combining (38) and (39) we have

$\displaystyle\widehat{\mu}=$
$\displaystyle\mu-\frac{\sqrt{2\pi}\sigma}{(m+1)\sqrt{n}}\sum_{j\notin{\mathcal{B}}}\left\\{I_{j}(0)-G(0)\right\\}+O_{{\mathbb{P}}}\left(\frac{\alpha_{n}}{\sqrt{n}}+\frac{\log
n}{m\sqrt{n}}+\frac{\log^{3/4}n}{m^{3/4}n^{1/2}}+\frac{1}{n}\right)$
$\displaystyle=$
$\displaystyle\mu-\frac{\sqrt{2\pi}\sigma}{(m+1)\sqrt{n}}\sum_{j\notin{\mathcal{B}}}\left\\{{\mathbb{I}}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\leq\mu)-{\mathbb{P}}(\,\hbox{\kern
0.09995pt\vbox{\hrule height=0.5pt\kern
1.42082pt\hbox{\kern-1.00006pt$X$\kern-1.00006pt}}}\,_{j}\leq\mu)\right\\}+O_{{\mathbb{P}}}\left(\frac{\alpha_{n}}{\sqrt{n}}+\frac{\log
n}{m\sqrt{n}}+\frac{\log^{3/4}n}{m^{3/4}n^{1/2}}+\frac{1}{n}\right),$

which is exactly what we want to prove. ∎

##  Appendix B Positive Definiteness of
$\mathcal{C}_{\mathrm{MOM}}-\mathcal{C}$

From Theorem 4 and Proposition 1, in order to show that our proposed multi-
dimensional VRMOM estimator  $\boldsymbol{\mu}$ has higher statistical
efficiency than the corresponding MOM estimator $\widehat{\boldsymbol{\mu}}$,
it is equivalent to prove that
$\boldsymbol{\mathcal{C}}\preceq\boldsymbol{\mathcal{C}}_{\mathrm{MOM}}$ holds
true. In this appendix, we will verify that the covariance difference
$\boldsymbol{\mathcal{C}}_{\mathrm{MOM}}-\boldsymbol{\mathcal{C}}$ is positive
definite in dimension 2 as $K$ tends to infinity. First of all, we shall give
a general formulation for the entry $\mathcal{C}_{l_{1},l_{2}}$ as
$K\rightarrow\infty$.

######  Lemma 6.

Denote $\mathcal{C}_{l_{1},l_{2}}^{K}$ as the $(l_{1},l_{2})$-entry of the
matrix $\boldsymbol{\mathcal{C}}$ defined in (14). Then we have

$\lim_{K\rightarrow\infty}\mathcal{C}_{l_{1},l_{2}}^{K}=\Big{\\{}4\pi\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\psi(y_{1})\psi(y_{2})F_{l_{1},l_{2}}(y_{1},y_{2}){\mathrm{d}}y_{1}{\mathrm{d}}y_{2}-\pi\Big{\\}}\sqrt{\sigma_{l_{1},l_{1}}\sigma_{l_{2},l_{2}}}.$
(40)

In particular, when $l_{1}=l_{2}=l$, we have

$\lim_{K\rightarrow\infty}\mathcal{C}_{l,l}^{K}=\frac{\pi}{3}\sigma_{l,l}.$
(41)

###### Proof.

For the denominator in (14), we compute that

$\displaystyle\lim_{K\rightarrow\infty}\frac{1}{K}\sum_{k=1}^{K}\psi(\Delta_{k})=$
$\displaystyle\lim_{K\rightarrow\infty}\frac{K+1}{K}\sum_{k=1}^{K}\psi\Big{(}\Phi^{-1}\Big{(}\frac{k}{K+1}\Big{)}\Big{)}\frac{1}{K+1}$
(42) $\displaystyle=$
$\displaystyle\int_{0}^{1}\psi\big{(}\Phi^{-1}(x)\big{)}{\mathrm{d}}x$
$\displaystyle(\text{change variable }x=\Phi(y))=$
$\displaystyle\int_{-\infty}^{\infty}\psi^{2}(y){\mathrm{d}}y$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-x^{2}}{\mathrm{d}}x=\frac{1}{2\sqrt{\pi}}.$

For the numerator, on the one hand, we have

$\displaystyle\lim_{K\rightarrow\infty}\frac{1}{K^{2}}\sum_{k_{1},k_{2}=1}^{K}\tau_{k_{1}}\tau_{k_{2}}=$
$\displaystyle\lim_{K\rightarrow\infty}\Big{(}\frac{1}{K}\sum_{k=1}^{K}\tau_{k}\Big{)}^{2}$
(43) $\displaystyle=$
$\displaystyle\lim_{K\rightarrow\infty}\Big{(}\frac{1}{K}\sum_{k=1}^{K}\frac{k}{K+1}\Big{)}^{2}$
$\displaystyle=$
$\displaystyle\Big{(}\int_{0}^{1}x{\mathrm{d}}x\Big{)}^{2}=\frac{1}{4}.$

On the other hand, we have that

$\displaystyle\lim_{K\rightarrow\infty}\frac{1}{K^{2}}\sum_{k_{1},k_{2}=1}^{K}\tau^{l_{1},l_{2}}_{k_{1},k_{2}}=$
$\displaystyle\int_{0}^{1}\int_{0}^{1}F_{l_{1},l_{2}}(\Phi^{-1}(x_{1}),\Phi^{-1}(x_{2})){\mathrm{d}}x_{1}{\mathrm{d}}x_{2}$
(44) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\psi(y_{1})\psi(y_{2})F_{l_{1},l_{2}}(y_{1},y_{2}){\mathrm{d}}y_{1}{\mathrm{d}}y_{2},$

where $F_{l_{1},l_{2}}(y_{1},y_{2})={\mathbb{P}}(Z_{l_{1}}\leq
y_{1},Z_{l_{2}}\leq y_{2})$. Combining (42), (43) and (44) we have

$\displaystyle\lim_{K\rightarrow\infty}\mathcal{C}_{l_{1},l_{2}}^{K}=\Big{\\{}4\pi\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\psi(y_{1})\psi(y_{2})F_{l_{1},l_{2}}(y_{1},y_{2}){\mathrm{d}}y_{1}{\mathrm{d}}y_{2}-\pi\Big{\\}}\sqrt{\sigma_{l_{1},l_{1}}\sigma_{l_{2},l_{2}}}.$
(45)

In particular, when $l_{1}=l_{2}=l$, we have

$F_{l,l}(y_{1},y_{2})=\Phi(\min(y_{1},y_{2})).$

Substitute it in (44), we can obtain

$\displaystyle\lim_{K\rightarrow\infty}\frac{1}{K^{2}}\sum_{k_{1},k_{2}=1}^{K}\tau^{l,l}_{k_{1},k_{2}}=$
$\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\psi(y_{1})\psi(y_{2})\Phi(\min(y_{1},y_{2})){\mathrm{d}}y_{1}{\mathrm{d}}y_{2}$
$\displaystyle=$ $\displaystyle
2\int_{-\infty}^{\infty}\psi(y_{1})\int_{-\infty}^{y_{1}}\psi(y_{2})\Phi(y_{2}){\mathrm{d}}y_{2}{\mathrm{d}}y_{1}$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\psi(y_{1})\Phi^{2}(y_{1}){\mathrm{d}}y_{1}=\frac{1}{3}.$

Therefore we have

$\lim_{K\rightarrow\infty}\mathcal{C}_{l,l}=\frac{1/3-1/4}{1/(4\pi)}\sigma_{l,l}=\frac{\pi}{3}\sigma_{l,l},$

which completes the proof. ∎

###### Verification of
$\boldsymbol{\mathcal{C}}\preceq\boldsymbol{\mathcal{C}}_{\mathrm{MOM}}$.

In the case of dimension $2$, we assume the gradient $\nabla
f(X,\boldsymbol{\theta}^{*})=(\nabla_{1}f(X,\boldsymbol{\theta}^{*}),$
$\nabla_{2}f(X,\boldsymbol{\theta}^{*}))^{{\mathrm{T}}}$ has covariance matrix

$\boldsymbol{\Sigma}=\begin{pmatrix}\sigma_{1,1}&\sin\phi\sqrt{\sigma_{1,1}\sigma_{2,2}}\\\
\sin\phi\sqrt{\sigma_{1,1}\sigma_{2,2}}&\sigma_{2,2}\end{pmatrix},\quad\text{therefore,
}\quad\boldsymbol{\Sigma}_{1,2}=\begin{pmatrix}1&\sin\phi\\\
\sin\phi&1\end{pmatrix}.$

From Lemma 6, we have that
$\mathcal{C}_{\mathrm{MOM},l,l}-\mathcal{C}_{l,l}=\pi\sigma_{l,l}/6$ as
$K\rightarrow\infty$, and

$\displaystyle\mathcal{C}_{1,2}=$ $\displaystyle
4\pi\left\\{\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\psi(y_{1})\psi(y_{2})F_{1,2}(y_{1},y_{2}){\mathrm{d}}y_{1}{\mathrm{d}}y_{2}-1/4\right\\}\sqrt{\sigma_{1,1}\sigma_{2,2}}$
$\displaystyle=$ $\displaystyle
4\pi\left\\{\int_{-\infty}^{\infty}\psi(y_{2}){\mathrm{d}}y_{2}\int_{-\infty}^{y_{2}}{\mathrm{d}}x_{2}\int_{-\infty}^{\infty}\int_{-\infty}^{y_{1}}\psi(y_{1})\psi_{\phi}(x_{1},x_{2}){\mathrm{d}}x_{1}{\mathrm{d}}y_{1}-1/4\right\\}\sqrt{\sigma_{1,1}\sigma_{2,2}}$
$\displaystyle=$ $\displaystyle
4\pi\left\\{\int_{-\infty}^{\infty}(1-\Phi(x_{1})){\mathrm{d}}x_{1}\int_{-\infty}^{\infty}\int_{-\infty}^{y_{2}}\psi(y_{2})\psi_{\phi}(x_{1},x_{2}){\mathrm{d}}x_{2}{\mathrm{d}}y_{2}-1/4\right\\}\sqrt{\sigma_{1,1}\sigma_{2,2}}$
$\displaystyle=$
$\displaystyle\left\\{4\pi\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\Phi(-x_{1})\Phi(-x_{2})\psi_{\phi}(x_{1},x_{2}){\mathrm{d}}x_{1}{\mathrm{d}}x_{2}-\pi\right\\}\sqrt{\sigma_{1,1}\sigma_{2,2}},$

where $\psi_{\phi}(\cdot,\cdot)$ denotes the probability density function of
the multivariate normal distribution with covariance matrix $\Sigma_{1,2}$,
more precisely, we have

$\psi_{\phi}(x_{1},x_{2})=\frac{1}{2\pi\cos\phi}\exp\Big{\\{}-\frac{x_{1}^{2}-2\sin\phi\cdot
x_{1}x_{2}+x_{2}^{2}}{2\cos^{2}\phi}\Big{\\}}.$

By symmetry we clearly see that

$\psi_{\phi}(x_{1},x_{2})=\psi_{\phi}(-x_{1},-x_{2}),\quad\psi_{\phi}(x_{1},-x_{2})=\psi_{\phi}(-x_{1},x_{2})=\psi_{-\phi}(x_{1},x_{2}),\quad\Phi(x)+\Phi(-x)=1.$

Then $\mathcal{C}_{1,2}$ can be further simplified as follows

$\displaystyle\mathcal{C}_{1,2}/\sqrt{\sigma_{1,1}\sigma_{2,2}}$
$\displaystyle=$ $\displaystyle
4\pi\int_{-\infty}^{0}\int_{-\infty}^{0}\big{\\{}\Phi(-x_{1})\Phi(-x_{2})+\Phi(-x_{1})\Phi(x_{2})+\Phi(x_{1})\Phi(-x_{2})+\Phi(x_{1})\Phi(x_{2})\big{\\}}\psi_{\phi}(x_{1},x_{2}){\mathrm{d}}x_{1}{\mathrm{d}}x_{2}$
$\displaystyle+4\pi\int_{-\infty}^{0}\int_{-\infty}^{0}\big{\\{}\Phi(-x_{1})\Phi(x_{2})+\Phi(x_{1})\Phi(-x_{2})\big{\\}}\big{\\{}\psi_{\phi}(x_{1},-x_{2})-\psi_{\phi}(x_{1},x_{2})\big{\\}}{\mathrm{d}}x_{1}{\mathrm{d}}x_{2}-\pi$
$\displaystyle=$ $\displaystyle
2\pi\int_{-\infty}^{0}\int_{-\infty}^{0}\big{\\{}1-2\Phi(x_{2})\big{\\}}\big{\\{}1-2\Phi(x_{1})\big{\\}}\big{\\{}\psi_{\phi}(x_{1},x_{2})-\psi_{-\phi}(x_{1},x_{2})\big{\\}}{\mathrm{d}}x_{1}{\mathrm{d}}x_{2}.$

Similarly, $\mathcal{C}_{\mathrm{MOM},1,2}$ can be written in a similar form

$\mathcal{C}_{\mathrm{MOM},1,2}/\sqrt{\sigma_{1,1}\sigma_{2,2}}=\pi\int_{-\infty}^{0}\int_{-\infty}^{0}\big{\\{}\psi_{\phi}(x_{1},x_{2})-\psi_{-\phi}(x_{1},x_{2})\big{\\}}{\mathrm{d}}x_{1}{\mathrm{d}}x_{2}=\phi.$

Therefore, to prove the positive definiteness of the matrix
$\boldsymbol{\mathcal{C}}_{\mathrm{MOM}}-\boldsymbol{\mathcal{C}}$, it left to
prove that

$|\mathcal{C}_{\mathrm{MOM},1,2}-\mathcal{C}_{1,2}|/\sqrt{\sigma_{1,1}\sigma_{2,2}}\leq|\mathcal{C}_{\mathrm{MOM},1,1}-\mathcal{C}_{1,1}|/\sigma_{1,1}=\frac{\pi}{6}.$

It is equivalent to the following inequality

$\displaystyle|h(\phi)|:=$
$\displaystyle\left|\frac{\phi}{\pi}-2\int_{-\infty}^{0}\int_{-\infty}^{0}\big{\\{}1-2\Phi(x_{2})\big{\\}}\big{\\{}1-2\Phi(x_{1})\big{\\}}\big{\\{}\psi_{\phi}(x_{1},x_{2})-\psi_{-\phi}(x_{1},x_{2})\big{\\}}{\mathrm{d}}x_{1}{\mathrm{d}}x_{2}\right|$
(46) $\displaystyle\leq$ $\displaystyle\frac{1}{6},$

holds for all $\phi\in[-\frac{\pi}{2},\frac{\pi}{2}]$. In order to show this
bound, we can draw the graph of $h(\phi)$ numerically.

Figure 3: The graph of function $h(\phi)$ on the interval $[-\pi/2,\pi/2]$. We
can see that $h$ is increasing and has absolute value uniformly bounded by
$1/6$.

As shown in Figure 3, we can see that (46) holds true, which implies that
$\boldsymbol{\mathcal{C}}\preceq\boldsymbol{\mathcal{C}}_{\mathrm{MOM}}$. ∎

##  Appendix C Theories and Proofs for RCSL Estimator

This appendix consists of the theoretical results and proofs for the robust
CSL estimator. In Appendix C.1, we present the technical assumptions and the
main theories for the robust CSL estimator. The proofs of the results will be
given in Appendix C.2.

###  C.1 Theoretical Results for Robust CSL Estimator

In this part, we present the theoretical results for the robust CSL estimator.
Before that, we first introduce some notations and our technical assumptions.

####  C.1.1 Notations and technical assumptions

As we will demonstrate in Appendix D, all the technical assumptions hold on
common statistical models, which suggests the wide applicability of these
assumptions. For the loss function $f(x,\boldsymbol{\theta})$, we assume that
$f(x,\boldsymbol{\theta})$ is differentiable with respect to
$\boldsymbol{\theta}$ and denote $\nabla
f(x,\boldsymbol{\theta}):=\nabla|_{\boldsymbol{\theta}}f(x,\boldsymbol{\theta})$
as the gradient of $f(x,\boldsymbol{\theta})$ at $\boldsymbol{\theta}$. For a
given $\boldsymbol{\theta}$, we define the expected gradient
$\boldsymbol{\mu}(\boldsymbol{\theta})$ and population standard deviation of
the gradient for each coordinate $l$, $\sigma_{l}(\boldsymbol{\theta})$, as
follows,

$\displaystyle\boldsymbol{\mu}(\boldsymbol{\theta})=(\mu_{1}(\boldsymbol{\theta}),\dots\mu_{p}(\boldsymbol{\theta}))^{{\mathrm{T}}}:={\mathbb{E}}_{X\sim{\mathfrak{X}}}\\{\nabla
f(X,\boldsymbol{\theta})\\},$ (47)
$\displaystyle\sigma_{l}^{2}(\boldsymbol{\theta}):={\mathbb{E}}_{X\sim{\mathfrak{X}}}\left[\\{\nabla
f_{l}(X,\boldsymbol{\theta})-\mu_{l}(\boldsymbol{\theta})\\}^{2}\right],\;\;\text{for}\;1\leq
l\leq p.$

For notational simplicity, we will denote the expectation taken over the
randomness of $X$ by ${\mathbb{E}}:={\mathbb{E}}_{X\sim{\mathfrak{X}}}$.
Recall that $\boldsymbol{\theta}^{*}$ is the true parameter that minimizes the
loss function ${\mathbb{E}}f(X,\boldsymbol{\theta})$. Throughout the paper, we
assume that $\boldsymbol{\mu}(\boldsymbol{\theta}^{*})={\mathbb{E}}\\{\nabla
f(X,\boldsymbol{\theta}^{*})\\}=\boldsymbol{0}$, which holds as long as the
expectation and $\nabla$ can be interchanged. Finally, for the ease of
presentation, some parameter-independent constants such as $\rho,\eta,C_{M}$
will be used across different assumptions when there is no confusion.

######  Assumption A.

For every fixed $x\in{\mathfrak{X}}$, $f(x,\boldsymbol{\theta})$ is a convex
function of $\boldsymbol{\theta}$ on ${\mathbb{R}}^{p}$.

######  Assumption B.

There exists a constant $\rho>0$ such that

$\rho\leq\Lambda_{\min}\left[\nabla\boldsymbol{\mu}(\boldsymbol{\theta}^{*})\right]\leq\Lambda_{\max}\left[\nabla\boldsymbol{\mu}(\boldsymbol{\theta}^{*})\right]\leq\rho^{-1}.$

######  Assumption C.

The Hessian of the population loss $\nabla\boldsymbol{\mu}$ is Lipschitz
continuous. In particular, there exists a constant $C_{H}>0$ such that

$\left\|\,\nabla\boldsymbol{\mu}(\boldsymbol{\theta}_{1})-\nabla\boldsymbol{\mu}(\boldsymbol{\theta}_{2})\,\right\|\leq
C_{H}|\boldsymbol{\theta}_{1}-\boldsymbol{\theta}_{2}|_{2},$

holds for arbitrary
$\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2}\in{\mathbb{R}}^{p}$.

######  Assumption D.

For every $\boldsymbol{v}\in{\mathbb{S}}^{p-1}$ and $x\in{\mathfrak{X}}$,
define

$\displaystyle