Byzantine-Robust Distributed Learning: Towards Optimal Statistical Rates

Introduction

Many tasks in computer vision, natural language processing and recommendation systems require learning complex prediction rules from large datasets. As the scale of the datasets in these learning tasks continues to grow, it is crucial to utilize the power of distributed computing and storage. In such large-scale distributed systems, robustness and security issues have become a major concern. In particular, individual computing units—known as worker machines—may exhibit abnormal behavior due to crashes, faulty hardware, stalled computation or unreliable communication channels. Security issues are only exacerbated in the so-called Federated Learning setting, a modern distributed learning paradigm that is more decentralized, and that uses the data owners’ devices (such as mobile phones and personal computers) as worker machines (McMahan and Ramage, 2017, Konečnỳ et al., 2016). Such machines are often more unpredictable, and in particular may be susceptible to malicious and coordinated attacks.

Due to the inherent unpredictability of this abnormal (sometimes adversarial) behavior, it is typically modeled as Byzantine failure (Lamport et al., 1982), meaning that some worker machines may behave completely arbitrarily and can send any message to the master machine that maintains and updates an estimate of the parameter vector to be learned. Byzantine failures can incur major degradation in learning performance. It is well-known that standard learning algorithms based on naive aggregation of the workers’ messages can be arbitrarily skewed by a single Byzantine-faulty machine. Even when the messages from Byzantine machines take only moderate values—and hence are difficult to detect—and when the number of such machines is small, the performance loss can still be significant. We demonstrate such an example in our experiments in Section 7.

In this paper, we aim to develop distributed statistical learning algorithms that are provably robust against Byzantine failures. While this objective is considered in a few recent works (Feng et al., 2014, Blanchard et al., 2017, Chen et al., 2017), a fundamental problem remains poorly understood, namely the optimal statistical performance of a robust learning algorithm. A learning scheme in which the master machine always outputs zero regardless of the workers’ messages is certainly not affected by Byzantine failures, but it will not return anything statistically useful either. On the other hand, many standard distributed algorithms that achieve good statistical performance in the absence of Byzantine failures, become completely unreliable otherwise. Therefore, a main goal of this work is to understand the following questions: what is the best achievable statistical performance while being Byzantine-robust, and what algorithms achieve this performance?

To formalize this question, we consider a standard statistical setting of empirical risk minimization (ERM). Here $nm$ data points are sampled independently from some distribution and distributed evenly among $m$ machines, $\alpha m$ of which are Byzantine. The goal is to learn a parametric model by minimizing some loss function defined by the data. In this statistical setting, one expects that the error in learning the parameter, measured in an appropriate metric, should decrease when the amount of data $nm$ becomes larger and the fraction of Byzantine machines $\alpha$ becomes smaller. In fact, we can show that, at least for strongly convex problems, no algorithm can achieve an error lower than

regardless of communication costs;Throughout the paper, unless otherwise stated, $\Omega(\cdot)$ and $\mathcal{O}(\cdot)$ hide universal multiplicative constants; $\widetilde{\Omega}(\cdot)$ and $\widetilde{\mathcal{O}}(\cdot)$ further hide terms that are independent of $\alpha,n,m$ or logarithmic in $n,m.$ see Observation 1 in Section 6. Intuitively, the above error rate is the optimal rate that one should target, as $\frac{1}{\sqrt{n}}$ is the effective standard deviation for each machine with $n$ data points, $\alpha$ is the bias effect of Byzantine machines, and $\frac{1}{\sqrt{m}}$ is the averaging effect of $m$ normal machines. When there are no or few Byzantine machines, we see the usual scaling $\frac{1}{\sqrt{mn}}$ with the total number of data points; when some machines are Byzantine, their influence remains bounded, and moreover is proportional to $\alpha$. If an algorithm is guaranteed to attain this bound, we are assured that we do not sacrifice the quality of learning when trying to guard against Byzantine failures—we pay a price that is unavoidable, but otherwise we achieve the best possible statistical accuracy in the presence of Byzantine failures.

Another important consideration for us is communication efficiency. As communication between machines is costly, one cannot simply send all data to the master machine. This constraint precludes direct application of standard robust learning algorithms (such as M-estimators (Huber, 2011)), which assume access to all data. Instead, a desirable algorithm should involve a small number of communication rounds as well as a small amount of data communicated per round. We consider a setting where in each round a worker or master machine can only communicate a vector of size $\mathcal{O}(d)$, where $d$ is the dimension of the parameter to be learned. In this case, the total communication cost is proportional to the number of communication rounds.

To summarize, we aim to develop distributed learning algorithms that simultaneously achieve two objectives:

Statistical optimality: attain an $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}})$ rate.

Communication efficiency: $\mathcal{O}(d)$ communication per round, with as few rounds as possible.

To the best of our knowledge, no existing algorithm achieves these two goals simultaneously. In particular, previous robust algorithms either have unclear or sub-optimal statistical guarantees, or incur a high communication cost and hence are not applicable in a distributed setting—we discuss related work in more detail in Section 2.

We propose two robust distributed gradient descent (GD) algorithms, one based on coordinate-wise median, and the other on coordinate-wise trimmed mean. We establish their statistical error rates for strongly convex, non-strongly convex, and non-convex population loss functions. For strongly convex losses, we show that these algorithms achieve order-optimal statistical rates under mild conditions. We further propose a median-based robust algorithm that only requires one communication round, and show that it also achieves the optimal rate for strongly convex quadratic losses. The statistical error rates of these three algorithms are summarized as follows.

Median-based GD: $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}+\frac{1}{n})$, order-optimal for strongly convex loss if $n\gtrsim m$.

Trimmed-mean-based GD: $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}})$, order-optimal for strongly convex loss.

Median-based one-round algorithm: $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}+\frac{1}{n})$, order-optimal for strongly convex quadratic loss if $n\gtrsim m$.

A major technical challenge in our statistical setting here is as follows: the $nm$ data points are sampled once and fixed, and each worker machine has access to a fixed set of data throughout the learning process. This creates complicated probabilistic dependency across the iterations of the algorithms. Worse yet, the Byzantine machines, which have complete knowledge of the data and the learning algorithm used, may create further unspecified probabilistic dependency. We overcome this difficulty by proving certain uniform bounds via careful covering arguments. Furthermore, for the analysis of median-based algorithms, we cannot simply adapt standard techniques (such as those in Minsker et al. (2015)), which can only show that the output of the master machine is as accurate as that of one normal machine, leading to a sub-optimal $\mathcal{O}(\frac{1}{\sqrt{n}})$ rate even without Byzantine failures ($\alpha=0$). Instead, we make use of a more delicate argument based on normal approximation and Berry-Esseen-type inequalities, which allows us to achieve the better $\mathcal{O}(\frac{1}{mn})$ rates when $\alpha$ is small while being robust for a nonzero $\alpha$.

Above we have omitted the dependence on the parameter dimension $d.$; see our main theorems for the precise results. In some settings the rates in these results may not have the optimal dependence on $d$. Understanding the fundamental limits of robust distributed learning in high dimensions, as well as developing algorithms with optimal dimension dependence, is an interesting and important future direction.

2 Notation

Related Work

Outlier-robust estimation in non-distributed settings is a classical topic in statistics (Huber, 2011). Particularly relevant to us is the so-called median-of-means method, in which one partitions the data into $m$ subsets, computes an estimate from each subset, and finally takes the median of these $m$ estimates. This idea is studied in Nemirovskii et al. (1983), Jerrum et al. (1986), Alon et al. (1999), Lerasle and Oliveira (2011), Minsker et al. (2015), and has been applied to bandit and least square regression problems (Bubeck et al., 2013, Lugosi and Mendelson, 2016, Kogler and Traxler, 2016) as well as problems involving heavy-tailed distributions (Hsu and Sabato, 2016, Lugosi and Mendelson, 2017). In a very recent work, Minsker and Strawn (2017) provide a new analysis of median-of-means using a normal approximation. We borrow some techniques from this paper, but need to address a significant harder problem: 1) we deal with the Byzantine setting with arbitrary/adversarial outliers, which is not considered in their paper; 2) we study iterative algorithms for general multi-dimensional problems with convex and non-convex losses, while they mainly focus on one-shot algorithms for mean-estimation-type problems.

The median-of-means method is used in the context of Byzantine-robust distributed learning in two recent papers. In particular, the work of Feng et al. (2014) considers a simple one-shot application of median-of-means, and only proves a sub-optimal $\widetilde{\mathcal{O}}(\frac{1}{\sqrt{n}})$ error rate as mentioned. The work of Chen et al. (2017) considers only strongly convex losses, and seeks to circumvent the above issue by grouping the worker machines into mini-batches; however, their rate $\widetilde{\mathcal{O}}(\frac{\sqrt{\alpha}}{\sqrt{n}}+\frac{1}{\sqrt{nm}})$ still falls short of being optimal, and in particular their algorithm fails even when there is only one Byzantine machine in each mini-batch.

Other methods have been proposed for Byzantine-robust distributed learning and optimization; e.g., Su and Vaidya (2016a, b). These works consider optimizing fixed functions and do not provide guarantees on statistical error rates. Most relevant is the work by Blanchard et al. (2017), who propose to aggregate the gradients from worker machines using a robust procedure. Their optimization setting—which is at the level of stochastic gradient descent and assumes unlimited, independent access to a strong stochastic gradient oracle—is fundamentally different from ours; in particular, they do not provide a characterization of the statistical errors given a fixed number of data points.

Communication efficiency has been studied extensively in non-Byzantine distributed settings (McMahan et al., 2016, Yin et al., 2017). An important class of algorithms are based on one-round aggregation methods (Zhang et al., 2012, 2015, Rosenblatt and Nadler, 2016). More sophisticated algorithms have been proposed in order to achieve better accuracy than the one-round approach while maintaining lower communication costs; examples include DANE (Shamir et al., 2014), Disco (Zhang and Lin, 2015), distributed SVRG (Lee et al., 2015) and their variants (Reddi et al., 2016, Wang et al., 2017). Developing Byzantine-robust versions of these algorithms is an interesting future direction.

For outlier-robust estimation in non-distributed settings, much progress has been made recently in terms of improved performance in high-dimensional problems (Diakonikolas et al., 2016, Lai et al., 2016, Bhatia et al., 2015) as well as developing list-decodable and semi-verified learning schemes when a majority of the data points are adversarial (Charikar et al., 2017). These results are not directly applicable to our distributed setting with general loss functions, but it is nevertheless an interesting future problem to investigate their potential extension for our problem.

Problem Setup

The parameter space $\mathcal{W}$ is assumed to be convex and compact with diameter $D$, i.e., $\|\mathbf{w}-\mathbf{w}^{\prime}\|_{2}\leq D,\forall\mathbf{w},\mathbf{w}^{\prime}\in\mathcal{W}$. We consider a distributed computation model with one master machine and $m$ worker machines. Each worker machine stores $n$ data points, each of which is sampled independently from $\mathcal{D}$. Denote by $\mathbf{z}^{i,j}$ the $j$-th data on the $i$-th worker machine, and $F_{i}(\mathbf{w}):=\frac{1}{n}\sum_{j=1}^{n}f(\mathbf{w};\mathbf{z}^{i,j})$ the empirical risk function for the $i$-th worker. We assume that an $\alpha$ fraction of the $m$ worker machines are Byzantine, and the remaining $1-\alpha$ fraction are normal. With the notation $[N]:=\{1,2,\ldots,N\}$, we index the set of worker machines by $[m]$, and denote the set of Byzantine machines by $\mathcal{B}\subset[m]$ (thus $|\mathcal{B}|=\alpha m$). The master machine communicates with the worker machines using some predefined protocol. The Byzantine machines need not obey this protocol and can send arbitrary messages to the master; in particular, they may have complete knowledge of the system and learning algorithms, and can collude with each other.

We introduce the coordinate-wise median and trimmed mean operations, which serve as building blocks for our algorithm.

For the analysis, we need several standard definitions concerning random variables/vectors.

Robust Distributed Gradient Descent

We describe two robust distributed gradient descent algorithms, one based on coordinate-wise median and the other on trimmed mean. These two algorithms are formally given in Algorithm 1 as Option I and Option II, respectively, where the symbol $*$ represents an arbitrary vector.

In each parallel iteration of the algorithms, the master machine broadcasts the current model parameter to all worker machines. The normal worker machines compute the gradients of their local loss functions and then send the gradients back to the master machine. The Byzantine machines may send any messages of their choices. The master machine then performs a gradient descent update on the model parameter with step-size $\eta$, using either the coordinate-wise median or trimmed mean of the received gradients. The Euclidean projection $\Pi_{\mathcal{W}}(\cdot)$ ensures that the model parameter stays in the parameter space $\mathcal{W}$.

Below we provide statistical guarantees on the error rates of these algorithms, and compare their performance. Throughout we assume that each loss function $f(\mathbf{w};\mathbf{z})$ and the population loss function $F(\mathbf{w})$ are smooth:

For any $\mathbf{z}\in\mathcal{Z}$, the partial derivative of $f(\cdot;\mathbf{z})$ with respect to the $k$-th coordinate of its first argument, denoted by $\partial_{k}f(\cdot;\mathbf{z})$, is $L_{k}$-Lipschitz for each $k\in[d]$, and the function $f(\cdot;\mathbf{z})$ is $L$-smooth. Let $\widehat{L}:=\sqrt{\sum_{k=1}^{d}L_{k}^{2}}$. Also assume that the population loss function $F(\cdot)$ is $L_{F}$-smooth.

It is easy to see hat $L_{F}\leq L\leq\widehat{L}$. When the dimension of $\mathbf{w}$ is high, the quantity $\widehat{L}$ may be large. However, we will soon see that $\widehat{L}$ only appears in the logarithmic factors in our bounds and thus does not have a significant impact.

We first consider our median-based algorithm, namely Algorithm 1 with Option I. We impose the assumptions that the gradient of the loss function $f$ has bounded variance, and each coordinate of the gradient has coordinate-wise bounded absolute skewness:

For any $\mathbf{w}\in\mathcal{W}$, $\operatorname{Var}(\nabla f(\mathbf{w};\mathbf{z}))\leq V^{2}$.

For any $\mathbf{w}\in\mathcal{W}$, $\|\gamma(\nabla f(\mathbf{w};\mathbf{z}))\|_{\infty}\leq S$.

These assumptions are satisfied in many learning problems with small values of $V^{2}$ and $S$. Below we provide a concrete example in terms of a linear regression problem.

We prove Proposition 1 in Appendix A.1. In this example, the upper bound $V^{2}$ on $\operatorname{Var}(\nabla f(\mathbf{w};\mathbf{x},y))$ depends on dimension $d$ and the diameter of the parameter space. If the diameter is a constant, we have $V=\mathcal{O}(\sqrt{d})$. Moreover, the gradient skewness is bounded by a universal constant $S$ regardless of the size of the parameter space. In Appendix A.2, we provide another example showing that when the features in $\mathbf{x}$ are i.i.d. Gaussian distributed, the coordinate-wise skewness can be upper bounded by $429$.

We first consider the case where the population loss function $F(\cdot)$ is strongly convex. Note that we do not require strong convexity of the individual loss functions $f(\cdot;\mathbf{z})$.

Consider Option I in Algorithm 1. Suppose that Assumptions 1, 2, and 3 hold, $F(\cdot)$ is $\lambda_{F}$-strongly convex, and the fraction $\alpha$ of Byzantine machines satisfies

for some $\epsilon>0$. Choose step-size $\eta=1/L_{F}$. Then, with probability at least $1-\frac{4d}{(1+nm\widehat{L}D)^{d}}$, after $T$ parallel iterations, we have

with $\Phi^{-1}(\cdot)$ being the inverse of the cumulative distribution function of the standard Gaussian distribution $\Phi(\cdot)$.

We prove Theorem 1 in Appendix B. In (3), we hide universal constants and a higher order term that scales as $\frac{1}{nm}$, and the factor $C_{\epsilon}$ is a function of $\epsilon$; as a concrete example, $C_{\epsilon}\approx 4$ when $\epsilon=\frac{1}{6}$. Theorem 1 together with the inequality $\log(1-x)\leq-x$, guarantees that after running $T\geq\frac{L_{F}+\lambda_{F}}{\lambda_{F}}\log(\frac{\lambda_{F}}{2\Delta}\|\mathbf{w}^{0}-\mathbf{w}^{*}\|_{2})$ parallel iterations, with high probability we can obtain a solution $\widehat{\mathbf{w}}=\mathbf{w}^{T}$ with error $\|\widehat{\mathbf{w}}-\mathbf{w}^{*}\|_{2}\leq\frac{4}{\lambda_{F}}\Delta$.

Here we achieve an error rate (defined as the distance between $\widehat{\mathbf{w}}$ and the optimal solution $\mathbf{w}^{*}$) of the form $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}+\frac{1}{n})$. In Section 6, we provide a lower bound showing that the error rate of any algorithm is $\widetilde{\Omega}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}})$. Therefore the first two terms in the upper bound cannot be improved. The third term $\frac{1}{n}$ is due to the dependence of the median on the skewness of the gradients. When each worker machine has a sufficient amount of data, more specifically $n\gtrsim m$, we achieve an order-optimal error rate up to logarithmic factors.

Non-strongly Convex Losses:

We next consider the case where the population risk function $F(\cdot)$ is convex, but not necessarily strongly convex. In this case, we need a mild technical assumption on the size of the parameter space $\mathcal{W}$.

We then have the following result on the convergence rate in terms of the value of the population risk function.

Consider Option I in Algorithm 1. Suppose that Assumptions 1, 2, 3 and 4 hold, and that the population loss $F(\cdot)$ is convex, and $\alpha$ satisfies (2) for some $\epsilon>0$. Define $\Delta$ as in (3), and choose step-size $\eta=1/L_{F}$. Then, with probability at least $1-\frac{4d}{(1+nm\widehat{L}D)^{d}}$, after $T=\frac{L_{F}}{\Delta}\|\mathbf{w}^{0}-\mathbf{w}^{*}\|_{2}$ parallel iterations, we have

We prove Theorem 2 in Appendix C. We observe that the error rate, defined as the excess risk $F(\mathbf{w}^{T})-F(\mathbf{w}^{*})$, again has the form \widetilde{\mathcal{O}}\big{(}\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}+\frac{1}{n}\big{)}.

Non-convex Losses:

When $F(\cdot)$ is non-convex but smooth, we need a somewhat different technical assumption on the size of $\mathcal{W}$.

We have the following guarantees on the rate of convergence to a critical point of the population loss $F(\cdot)$.

Consider Option I in Algorithm 1. Suppose that Assumptions 1 2, 3 and 5 hold, and $\alpha$ satisfies (2) for some $\epsilon>0$. Define $\Delta$ as in (3), and choose step-size $\eta=1/L_{F}$. With probability at least $1-\frac{4d}{(1+nm\widehat{L}D)^{d}}$, after $T=\frac{2L_{F}}{\Delta^{2}}(F(\mathbf{w}^{0})-F(\mathbf{w}^{*}))$ parallel iterations, we have

We prove Theorem 3 in Appendix D. We again obtain an $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}+\frac{1}{n})$ error rate in terms of the gap to a critical point of $F(\mathbf{w})$.

2 Guarantees for Trimmed-mean-based Gradient Descent

We next analyze the robust distributed gradient descent algorithm based on coordinate-wise trimmed mean, namely Option II in Algorithm 1. Here we need stronger assumptions on the tail behavior of the partial derivatives of the loss functions—in particular, sub-exponentiality.

We assume that for all $k\in[d]$ and $\mathbf{w}\in\mathcal{W}$, the partial derivative of $f(\mathbf{w};\mathbf{z})$ with respect to the $k$-th coordinate of $\mathbf{w}$, $\partial_{k}f(\mathbf{w};\mathbf{z})$, is $v$-sub-exponential.

The sub-exponential property implies that all the moments of the derivatives are bounded. This is a stronger assumption than the bounded absolute skewness (hence bounded third moments) required by the median-based GD algorithm.

We use the same example as in Proposition 1 and show that the derivatives of the loss are indeed sub-exponential.

Consider the regression problem in Proposition 1. For all $k\in[d]$ and $\mathbf{w}\in\mathcal{W}$, the partial derivative $\partial_{k}f(\mathbf{w};\mathbf{z})$ is $\sqrt{\sigma^{2}+\|\mathbf{w}-\mathbf{w}^{*}\|_{2}^{2}}$-sub-exponential.

Consider Option II in Algorithm 1. Suppose that Assumptions 1 and 6 hold, $F(\cdot)$ is $\lambda_{F}$-strongly convex, and $\alpha\leq\beta\leq\frac{1}{2}-\epsilon$ for some $\epsilon>0$. Choose step-size $\eta=1/L_{F}$. Then, with probability at least $1-\frac{4d}{(1+nm\widehat{L}D)^{d}}$, after $T$ parallel iterations, we have

We prove Theorem 4 in Appendix E. In (5), we hide universal constants and higher order terms that scale as $\frac{\beta}{n}$ or $\frac{1}{nm}$. By running $T\geq\frac{L_{F}+\lambda_{F}}{\lambda_{F}}\log(\frac{\lambda_{F}}{2\Delta^{\prime}}\|\mathbf{w}^{0}-\mathbf{w}^{*}\|_{2})$ parallel iterations, we can obtain a solution $\widehat{\mathbf{w}}=\mathbf{w}^{T}$ satisfying $\|\widehat{\mathbf{w}}-\mathbf{w}^{*}\|_{2}\leq\widetilde{\mathcal{O}}(\frac{\beta}{\sqrt{n}}+\frac{1}{\sqrt{nm}})$. Note that one needs to choose the parameter for trimmed mean to satisfy $\beta\geq\alpha$. If we set $\beta=c\alpha$ for some universal constant $c\geq 1$, we can achieve an order-optimal error rate $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}})$.

Non-strongly Convex Losses:

Again imposing Assumption 4 on the size of $\mathcal{W}$, we have the following guarantee.

Consider Option II in Algorithm 1. Suppose that Assumptions 1, 4 and 6 hold, $F(\cdot)$ is convex, and $\alpha\leq\beta\leq\frac{1}{2}-\epsilon$ for some $\epsilon>0$. Choose step-size $\eta=1/L_{F}$, and define $\Delta^{\prime}$ as in (5). Then, with probability at least $1-\frac{4d}{(1+nm\widehat{L}D)^{d}}$, after $T=\frac{L_{F}}{\Delta^{\prime}}\|\mathbf{w}^{0}-\mathbf{w}^{*}\|_{2}$ parallel iterations, we have

The proof of Theorem 5 is similar to that of Theorem 2, and we refer readers to Remark 1 in Appendix E. Again, by choosing $\beta=c\alpha$ $(c\geq 1)$, we obtain the $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}})$ error rate in the function value of $F(\mathbf{w})$.

Non-convex Losses:

In this case, imposing a version of Assumption 5 on the size of $\mathcal{W}$, we have the following.

Consider Option II in Algorithm 1, and define $\Delta^{\prime}$ as in (5). Suppose that Assumptions 1 and 6 hold, Assumption 5 holds with $\Delta$ replaced by $\Delta^{\prime}$, and $\alpha\leq\beta\leq\frac{1}{2}-\epsilon$ for some $\epsilon>0$. Choose step-size $\eta=1/L_{F}$. Then, with probability at least $1-\frac{4d}{(1+nm\widehat{L}D)^{d}}$, after $T=\frac{2L_{F}}{\Delta^{\prime 2}}(F(\mathbf{w}^{0})-F(\mathbf{w}^{*}))$ parallel iterations, we have

The proof of Theorem 6 is similar to that of Theorem 3; see Remark 1 in Appendix E. By choosing $\beta=c\alpha$ with $c\geq 1$, we again achieve the statistical rate $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}})$.

3 Comparisons

We compare the performance guarantees of the above two robust distribute GD algorithms. The trimmed-mean-based algorithm achieves the statistical error rate $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}})$, which is order-optimal for strongly convex loss. In comparison, the rate of the median-based algorithm is $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}+\frac{1}{n})$, which has an additional $\frac{1}{n}$ term and is only optimal when $n\gtrsim m$. In particular, the trimmed-mean-based algorithm has better rates when each worker machine has small local sample size—the rates are meaningful even in the extreme case $n=\mathcal{O}(1)$. On the other hand, the median-based algorithm requires milder tail/moment assumptions on the loss derivatives (bounded skewness) than its trimmed-mean counterpart (sub-exponentiality). Finally, the trimmed-mean operation requires an additional parameter $\beta$, which can be any upper bound on the fraction $\alpha$ of Byzantine machines in order to guarantee robustness. Using an overly large $\beta$ may lead to a looser bound and sub-optimal performance. In contrast, median-based GD does not require knowledge of $\alpha$. We summarize these observations in Table 1. We see that the two algorithms are complementary to each other, and our experiment results corroborate this point.

Robust One-round Algorithm

As mentioned, in our distributed computing framework, the communication cost is proportional to the number of parallel iterations. The above two GD algorithms both require a number iterations depending on the desired accuracy. Can we further reduce the communication cost while keeping the algorithm Byzantine-robust and statistically optimal?

A natural candidate is the so-called one-round algorithm. Previous work has considered a standard one-round scheme where each local machine computes the empirical risk minimizer (ERM) using its local data and the master machine receives all workers’ ERMs and computes their average (Zhang et al., 2012). Clearly, a single Byzantine machine can arbitrary skew the output of this algorithm. We instead consider a Byzantine-robust one-round algorithm. As detailed in Algorithm 2, we employ the coordinate-wise median operation to aggregate all the ERMs.

The loss function $f(\mathbf{w};\mathbf{z})$ is quadratic if it can be written as

where $\mathbf{z}=(\mathbf{H},\mathbf{p},c)$, $\mathbf{H}$, and $\mathbf{p}$, and $c$ are drawn from the distributions $\mathcal{D}_{H}$, $\mathcal{D}_{p}$, and $\mathcal{D}_{c}$, respectively.

Denote by $\mathbf{H}_{F}$, $\mathbf{p}_{F}$, and $c_{F}$ the expectations of $\mathbf{H}$, $\mathbf{p}$, and $c$, respectively. Thus the population risk function takes the form $F(\mathbf{w})=\frac{1}{2}\mathbf{w}^{\rm T}\mathbf{H}_{F}\mathbf{w}+\mathbf{p}_{F}^{\rm T}\mathbf{w}+c_{F}$.

We need a technical assumption which guarantees that each normal worker machine has unique ERM.

With probability $1$, the empirical risk minimization function $F_{i}(\cdot)$ on each normal machine is strongly convex.

Note that this assumption is imposed on $F_{i}(\mathbf{w})$, rather than on the individual loss $f(\mathbf{w};\mathbf{z})$ associated with a single data point. This assumption is satisfied, for example, when all $f(\cdot;\mathbf{z})$’s are strongly convex, or in the linear regression problems with the features $\mathbf{x}$ drawn from some continuous distribution (e.g. isotropic Gaussian) and $n\geq d$. We have the following guarantee for the robust one-round algorithm.

for some $\epsilon>0$, where $\widetilde{C}$ is a quantity that depends on $\mathcal{D}_{H}$, $\mathcal{D}_{p}$, $\lambda_{F}$ and is monotonically decreasing in $n$. Then, with probability at least $1-\frac{4}{nm}$, the output $\widehat{\mathbf{w}}$ of the robust one-round algorithm satisfies

where $C_{\epsilon}$ is defined as in (4) and

with $\mathbf{H}$ and $\mathbf{p}$ drawn from $\mathcal{D}_{H}$ and $\mathcal{D}_{p}$, respectively.

We prove Theorem 7 and provide an explicit expression of $\widetilde{C}$ in Appendix F. In terms of the dependence on $\alpha$, $n$, and $m$, the robust one-round algorithm achieves the same error rate as the robust gradient descent algorithm based on coordinate-wise median, i.e., $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}+\frac{1}{n})$, for quadratic problems. Again, this rate is optimal when $n\gtrsim m$. Therefore, at least for quadratic loss functions, the robust one-round algorithm has similar theoretical performance as the robust gradient descent algorithm with significantly less communication cost. Our experiments show that the one-round algorithm has good empirical performance for other losses as well.

Lower Bound

In this section, we provide a lower bound on the error rate for strongly convex losses, which implies that the $\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}$ term is unimprovable. This lower bound is derived using a mean estimation problem, and is an extension of the lower bounds in the robust mean estimation literature such as Chen et al. (2015), Lai et al. (2016).

We consider the problem of estimating the mean $\bm{\mu}$ of some random variable $\mathbf{z}\sim\mathcal{Z}$, which is equivalent to solving the following minimization problem:

Note that this is a special case of the general learning problem (1). We consider the same distributed setting as in Section 4, with a minor technical difference regarding the Byzantine machines. We assume that each of the $m$ worker machines is Byzantine with probability $\alpha$, independently of each other. The parameter $\alpha$ is therefore the expected fraction of Byzantine machines. This setting makes the analysis slightly easier, and we believe the result can be extended to the original setting.

In this setting we have the following lower bound.

Consider the distributed mean estimation problem in (6) with Byzantine failure probability $\alpha$, and suppose that $\mathcal{Z}$ is Gaussian distribution with mean $\bm{\mu}$ and covariance matrix $\sigma^{2}\mathbf{I}$ $(\sigma=\mathcal{O}(1))$. Then, any algorithm that computes an estimation $\widehat{\bm{\mu}}$ of the mean from the data has a constant probability of error $\|\widehat{\bm{\mu}}-\bm{\mu}\|_{2}=\Omega(\frac{\alpha}{\sqrt{n}}+\sqrt{\frac{d}{nm}})$.

We prove Observation 1 in Appendix G. According to this observation, we see that the $\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}$ dependence cannot be avoided, which in turn implies the order-optimality of the results in Theorem 1 (when $n\gtrsim m$) and Theorem 4.

Experiments

We conduct experiments to show the effectiveness of the median and trimmed mean operations. Our experiments are implemented with Tensorflow (Abadi et al., 2016) on Microsoft Azure system. We use the MNIST (LeCun et al., 1998) dataset and randomly partition the 60,000 training data into $m$ subsamples with equal sizes. We use these subsamples to represent the data on $m$ machines.

In the first experiment, we compare the performance of distributed gradient descent algorithms in the following four settings: 1) $\alpha=0$ (no Byzantine machines), using vanilla distributed gradient descent (aggregating the gradients by taking the mean), 2) $\alpha>0$, using vanilla distributed gradient descent, 3) $\alpha>0$, using median-based algorithm, and 4) $\alpha>0$, using trimmed-mean-based algorithm. We generate the Byzantine machines in the following way: we replace every training label $y$ on these machines with $9-y$, e.g., is replaced with $9$, $1$ is replaced with $8$, etc, and the Byzantine machines simply compute gradients based on these data. We also note that when generating the Byzantine machines, we do not simply add extreme values in the features or gradients; instead, the Byzantine machines send messages to the master machine with moderate values.

We train a multi-class logistic regression model and a convolutional neural network model using distributed gradient descent, and for each model, we compare the test accuracies in the aforementioned four settings. For the convolutional neural network model, we use the stochastic version of the distributed gradient descent algorithm; more specifically, in every iteration, each worker machine computes the gradient using $10\%$ of its local data. We periodically check the test errors, and the convergence performances are shown in Figure 1. The final test accuracies are presented in Tables 2 and 3.

As we can see, in the adversarial settings, the vanilla distributed gradient descent algorithm suffers from severe performance loss, and using the median and trimmed mean operations, we observe significant improvement in test accuracy. This shows these two operations can indeed defend against Byzantine failures.

In the second experiment, we compare the performance of distributed one-round algorithms in the following three settings: 1) $\alpha=0$, mean aggregation, 2) $\alpha>0$, mean aggregation, and 3) $\alpha>0$, median aggregation. In this experiment, the training labels on the Byzantine machines are i.i.d. uniformly sampled from $\{0,\ldots,9\}$, and these machines train models using the faulty data. We choose the multi-class logistic regression model, and the test accuracies are presented in Table 4.

As we can see, for the one-round algorithm, although the theoretical guarantee is only proved for quadratic loss, in practice, the median-based one-round algorithm still improves the test accuracy in problems with other loss functions, such as the logistic loss here.

Conclusions

In this paper, we study Byzantine-robust distributed statistical learning algorithms with a focus on statistical optimality. We analyze two robust distributed gradient descent algorithms — one is based on coordinate-wise median and the other is based on coordinate-wise trimmed mean. We show that the trimmed-mean-based algorithm can achieve order-optimal $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}})$ error rate, whereas the median-based algorithm can achieve $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}+\frac{1}{n})$ under weaker assumptions. We further study learning algorithms that have better communication efficiency. We propose a simple one-round algorithm that aggregates local solutions using coordinate-wise median. We show that for strongly convex quadratic problems, this algorithm can achieve $\widetilde{\mathcal{O}}(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}+\frac{1}{n})$ error rate, similar to the median-based gradient descent algorithm. Our experiments validates the effectiveness of the median and trimmed mean operations in the adversarial setting.

D. Yin is partially supported by Berkeley DeepDrive Industry Consortium. Y. Chen is partially supported by NSF CRII award 1657420 and grant 1704828. K. Ramchandran is partially supported by NSF CIF award 1703678 and Gift award from Huawei. P. Bartlett is partially supported by NSF grant IIS-1619362. Cloud computing resources are provided by a Microsoft Azure for Research award.

References

Appendix

Appendix A Variance, Skewness, and Sub-exponential Property

We use the simplified notation $f(\mathbf{w}):=f(\mathbf{w};\mathbf{x},y)$. One can directly compute the gradients:

Define $\Delta(\mathbf{w}):=\nabla f(\mathbf{w})-\nabla F(\mathbf{w})$ with its $k$-th element being $\Delta_{k}(\mathbf{w})$. We now compute the variance and absolute skewness of $\Delta_{k}(\mathbf{w})$.

Then we proceed to bound $\gamma(\Delta_{k}(\mathbf{w}))$. By Jensen’s inequality, we know that

We first find a lower bound for $\operatorname{Var}(\Delta_{k}(\mathbf{w}))^{3}$. According to (8), we know that

where in the last inequality we use the moments of Gaussian random variables. Then, we compute the first term in (15). By algebra, one can obtain

A.2 Example of Regression with Gaussian Features

with some $\mathbf{w}^{*}\in\mathcal{W}$. Assume that the elements of $\mathbf{x}$ are i.i.d. samples of standard Gaussian distribution, and that the noise $\xi$ is independent of $\mathbf{x}$ and drawn from Gaussian distribution $\mathcal{N}(0,\sigma^{2})$. Define the quadratic loss function $f(\mathbf{w};\mathbf{x},y)=\frac{1}{2}(y-\mathbf{x}^{\rm T}\mathbf{w})^{2}$. Then, we have

We use the same simplified notation as in Appendix A.1. One can also see that (7) still holds for in the Gaussian setting. Thus,

Then we proceed to bound $\gamma(\Delta_{k}(\mathbf{w}))$. By Jensen’s inequality, we know that

We first find a lower bound for $\operatorname{Var}(\Delta_{k}(\mathbf{w}))^{3}$. According to (18), we know that

Define the $W_{1}$, $W_{2}$, and $W_{3}$ as in (10), (11), and (12). We can also see that (13) still holds, and thus

where in the last inequality we use the moments of Gaussian random variables. Then, we compute the first term in (22). By algebra, one can obtain

A.3 Proof of Proposition 2

We use the same notation as in Appendix A.1. We have

Appendix B Proof of Theorem 1

The proof of Theorem 1 consists of two parts: 1) the analysis of coordinate-wise median estimator of the population gradients, and 2) the convergence analysis of the robustified gradient descent algorithm.

Recall that at iteration $t$, the master machine sends $\mathbf{w}^{t}$ to all the worker machines. For any normal worker machine, say machine $i\in[m]\setminus\mathcal{B}$, the gradient of the local empirical loss function $\mathbf{g}^{i}(\mathbf{w}^{t})=\nabla F_{i}(\mathbf{w}^{t})$ is computed and returned to the center machine, while the Byzantine machines, say machine $i\in\mathcal{B}$, the returned message $\mathbf{g}^{i}(\mathbf{w}^{t})$ can be arbitrary or even adversarial. The master machine then compute the coordinate-wise median, i.e.,

The following theorem provides a uniform bound on the distance between $\mathbf{g}(\mathbf{w}^{t})$ and $\nabla F(\mathbf{w}^{t})$.

and the coordinate-wise median of $\mathbf{g}^{i}(\mathbf{w})$:

Suppose that Assumptions 1, 2, and 3 hold, and inequality (2) is satisfied with some $\epsilon>0$. Then, we have with probability at least $1-\frac{4d}{(1+nm\widehat{L}D)^{d}}$,

for all $\mathbf{w}\in\mathcal{W}$, where $C_{\epsilon}$ is defined as in (4).

Then, we proceed to analyze the convergence of the robust distributed gradient descent algorithm. We condition on the event that the bound in (27) is satisfied for all $\mathbf{w}\in\mathcal{W}$. Then, in the $t$-th iteration, we define

Thus, we have $\mathbf{w}^{t+1}=\Pi_{\mathcal{W}}(\widehat{\mathbf{w}}^{t+1})$. By the property of Euclidean projection, we know that

Since $F(\mathbf{w})$ is $\lambda_{F}$-strongly convex, by the co-coercivity of strongly convex functions (see Lemma 3.11 in Bubeck et al. (2015) for more details), we obtain

where in the second inequality we use the fact that $\lambda_{F}\leq L_{F}$. Using the fact $\sqrt{1-x}\leq 1-\frac{x}{2}$, we get

Then we can complete the proof by iterating (31).

The proof of Theorem 8 relies on careful analysis of the median of means estimator in the presence of adversarial data and a covering net argument.

for some $\epsilon>0$. Then, with probability at least $1-4e^{-2t}$, we have

where $C_{\epsilon}$ is defined as in (4).

We further define the distribution function of all the $m$ machines as $\widehat{p}(z):=\frac{1}{m}\sum_{i\in[m]}\mathds{1}(\bar{x}^{i}\leq z)$. We have the following direct corollary on $\widehat{p}(z)$ and the median of means estimator ${\sf med}\{\bar{x}^{i}:i\in[m]\}$.

Suppose that condition (32) is satisfied. Then, with probability at least $1-4e^{-2t}$, we have,

Thus, we have with probability at least $1-4e^{-2t}$,

We use the symbol $\partial_{k}$ to denote the partial derivative of any function with respect to its $k$-th argument. We also use the simplified notation $\sigma_{k}^{2}(\mathbf{w}):=\operatorname{Var}(\partial_{k}f(\mathbf{w};\mathbf{z}))$, and $\gamma_{k}(\mathbf{w}):=\gamma(\partial_{k}f(\mathbf{w};\mathbf{z}))$. Then, according to Lemma 1, when (32) is satisfied, for any fixed $\mathbf{w}\in\mathcal{W}$ and $k\in[d]$, we have with probability at least $1-4e^{-2t}$,

Further, according to Corollary 1, we know that with probability $1-4e^{-2t}$,

Here, the inequality (42) gives a bound on the accuracy of the median of means estimator for the gradient at any fixed $\mathbf{w}$ and any coordinate $k\in[d]$. To extend this result to all $\mathbf{w}\in\mathcal{W}$ and all the $d$ coordinates, we need to use union bound and a covering net argument.

Again, by gathering all the $k$ coordinates we get

where we use the fact that $(a+b)^{2}\leq 2(a^{2}+b^{2})$. Then, by Assumption 2 and 3, we further obtain

where we use the fact that $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$. Combining (43) and (46), we conclude that for any $\delta>0$, with probability at least $1-4dN_{\delta}e^{-2t}$, (46) holds for all $\mathbf{w}\in\mathcal{W}$. We simply choose $\delta=\frac{1}{nm\widehat{L}}$, and $t=d\log(1+nm\widehat{L}D)$. Then, we know that with probability at least $1-\frac{4d}{(1+nm\widehat{L}D)^{d}}$, we have

B.2 Proof of Lemma 1

We recall the Berry-Esseen Theorem (Berry, 1941, Esseen, 1942, Shevtsova, 2014) and the bounded difference inequality, which are useful in this proof.

where $\bar{Y}=\frac{1}{n}\sum_{i=1}^{n}Y_{i}$ and $\Phi(s)$ is the cumulative distribution function of the standard normal random variable.

Let $X_{1},\ldots,X_{n}$ be i.i.d. random variables, and assume that $Z=g(X_{1},\ldots,X_{n})$, where $g$ satisfies that for all $j\in[n]$ and all $x_{1},x_{2},\ldots,x_{j},x^{\prime}_{j},\ldots,x_{n}$,

on the draw of $W_{i}$, $i\in[m]\setminus\mathcal{B}$ with probability at least $1-2e^{-2t}$. Let $z_{1}\geq z_{2}$ be such that $\Phi_{n}(z_{1})\geq\frac{1}{2}+\alpha+\sqrt{\frac{t}{m(1-\alpha)}}$, and $\Phi_{n}(z_{2})\leq\frac{1}{2}-\alpha-\sqrt{\frac{t}{m(1-\alpha)}}$. Then, by union bound, we know that with probability at least $1-4e^{-2t}$, $\widetilde{\Phi}_{n}(z_{1})\geq 1/2+\alpha$ and $\widetilde{\Phi}_{n}(z_{2})\leq 1/2-\alpha$. The next step is to choose $z_{1}$ and $z_{2}$. According to Theorem 9, we know that

and thus, it suffices to find $z_{1}$ such that

By mean value theorem, we know that there exists $\xi\in[0,z_{1}]$ such that

Suppose that for some fix constant $\epsilon\in(0,1/2)$, we have

Then, we know that $z_{1}\leq\Phi^{-1}(1-\epsilon)$, and thus we have

For simplicity, let $C_{\epsilon}:=\sqrt{2\pi}\exp(\frac{1}{2}(\Phi^{-1}(1-\epsilon))^{2})$. We conclude that with probability $1-4e^{-2t}$, we have

Appendix C Proof of Theorem 2

Since Theorem 8 holds without assuming the convexity of $F(\mathbf{w})$, when $F(\mathbf{w})$ is non-strongly convex, the event that (27) holds for all $\mathbf{w}\in\mathcal{W}$ still happens with probability at least $1-\frac{4d}{(1+nm\widehat{L}D)^{d}}$. We condition on this event. We first show that when Assumption 4 is satisfied and we choose $\eta=\frac{1}{L_{F}}$, the iterates $\mathbf{w}^{t}$ stays in $\mathcal{W}$ without using projection. Namely, define

for $T=0,1,\ldots,T-1$, then $\mathbf{w}^{t}\in\mathcal{W}$ for all $t=0,1,\ldots,T$. To see this, we have

where the inequality is due to the co-coercivity of convex functions. Thus, we get

and since $T=\frac{L_{F}D_{0}}{\Delta}$, according to Assumption 4 we know that $\mathbf{w}^{t}\in\mathcal{W}$ for all $t=0,1,\ldots,T$. Then, we proceed to study the algorithm without projection. Here, we define $D_{t}:=\|\mathbf{w}^{0}-\mathbf{w}^{*}\|_{2}+\frac{t\Delta}{L_{F}}$ for $t=0,1,\ldots,T$.

Using the smoothness of $F(\mathbf{w})$, we have

Since $\eta=\frac{1}{L_{F}}$ and $\|\mathbf{g}(\mathbf{w}^{t})-\nabla F(\mathbf{w}^{t})\|_{2}\leq\Delta$, by simple algebra, we obtain

Condition on the event that (27) holds for all $\mathbf{w}\in\mathcal{W}$. When $F(\mathbf{w})$ is convex, by running $T=\frac{L_{F}D_{0}}{\Delta}$ parallel iterations, there exists $t\in\{0,1,2,\ldots,T\}$ such that

We first notice that since $T=\frac{L_{F}D_{0}}{\Delta}$, we have $D_{t}\leq 2D_{0}$ for all $t=0,1,\ldots,T$. According to the first order optimality of convex functions, for any $\mathbf{w}$,

Suppose that there exists $t\in\{0,1,\ldots,T-1\}$ such that $\|\nabla F(\mathbf{w}^{t})\|_{2}<\sqrt{2}\Delta$. Then we have

Otherwise, for all $t\in\{0,1,\ldots,T-1\}$, $\|\nabla F(\mathbf{w}^{t})\|_{2}\geq\sqrt{2}\Delta$. Then, according to (49) and (50), we have for all $t<T$,

Multiplying both sides by $[(F(\mathbf{w}^{t+1})-F(\mathbf{w}^{*}))(F(\mathbf{w}^{t})-F(\mathbf{w}^{*})]^{-1}$ and rearranging the terms, we obtain

Then, we obtain $F(\mathbf{w}^{T})-F(\mathbf{w}^{*})\leq 16D_{0}\Delta$ using the fact that $T=\frac{L_{F}D_{0}}{\Delta}$. ∎

Next, we show that $F(\mathbf{w}^{T})-F(\mathbf{w}^{*})\leq 16D_{0}\Delta+\frac{1}{2L_{F}}\Delta^{2}$. More specifically, let $t=t_{0}$ be the first time that $F(\mathbf{w}^{t})-F(\mathbf{w}^{*})\leq 16D_{0}\Delta$, and we show that for any $t>t_{0}$, $F(\mathbf{w}^{t})-F(\mathbf{w}^{*})\leq 16D_{0}\Delta+\frac{1}{2L_{F}}\Delta^{2}$. If this statement is not true, then we let $t_{1}>t_{0}$ be the first time that $F(\mathbf{w}^{t})-F(\mathbf{w}^{*})>16D_{0}\Delta+\frac{1}{2L_{F}}\Delta^{2}$. Then there must be $F(\mathbf{w}^{t_{1}-1})<F(\mathbf{w}^{t_{1}})$. According to (49), there should also be

Then according to (49), this implies $F(\mathbf{w}^{t_{1}})\leq F(\mathbf{w}^{t_{1}-1})$, which contradicts with the fact that $F(\mathbf{w}^{t_{1}-1})<F(\mathbf{w}^{t_{1}})$.

Appendix D Proof of Theorem 3

Since Theorem 8 holds without assuming the convexity of $F(\mathbf{w})$, when $F(\mathbf{w})$ is non-convex, the event that (27) holds for all $\mathbf{w}\in\mathcal{W}$ still happens with probability at least $1-\frac{4d}{(1+nm\widehat{L}D)^{d}}$. We condition on this event. We first show that when Assumption 5 is satisfied and we choose $\eta=\frac{1}{L_{F}}$, the iterates $\mathbf{w}^{t}$ stays in $\mathcal{W}$ without using projection. Since we have

Then, we know that by running $T=\frac{2L_{F}}{\Delta^{2}}(F(\mathbf{w}^{0})-F(\mathbf{w}^{*}))$ parallel iterations, using Assumption 5, we know that $\mathbf{w}^{t}\in\mathcal{W}$ for $t=0,1,\ldots,T$ without projection.

We proceed to study the convergence rate of the algorithm. By the smoothness of $F(\mathbf{w})$, we know that when choosing $\eta=\frac{1}{L_{F}}$, the inequality (49) still holds. More specifically, for all $t=0,1,\ldots,T-1$,

Sum up (51) for $t=0,1,\ldots,T-1$. Then, we get

Appendix E Proof of Theorem 4

The proof of Theorem 4 consists of two parts: 1) the analysis of coordinate-wise trimmed mean of means estimator of the population gradients, and 2) the convergence analysis of the robustified gradient descent algorithm. Since the second part is essentially the same as the proof of Theorem 1, we mainly focus on the first part here.

and the coordinate-wise trimmed mean of $\mathbf{g}^{i}(\mathbf{w})$:

Suppose that Assumptions 1 and 6 are satisfied, and that $\alpha\leq\beta\leq\frac{1}{2}-\epsilon$. Then, with probability at least $1-\frac{2d(m+1)}{(1+nm\widehat{L}D)^{d}}$,

The rest of the proof is essentially the same as the proof of Theorem 1. In fact, we essentially analyze a gradient descent algorithm with bounded noise in the gradients. In the proof of Theorem 1 in Appendix B. The bound on the noise in the gradients is

while here we replace $\Delta$ with $\Delta^{\prime}$:

and the same analysis can still go through. Therefore, we omit the details of the analysis here.

The same arguments still go through when the population risk function $F(\mathbf{w})$ is non-strongly convex or non-convex. One can simply replace the bound on the noise in the gradients $\Delta$ in Theorems 2 and 3 with $\Delta^{\prime}$ here. Thus we omit the details here.

Suppose that the one dimensional samples on all the normal machines are i.i.d. $v$-sub-exponential with mean $\mu$. Then, we have for any $t\geq 0$,

and when $\beta\geq\alpha$, $|\frac{1}{(1-\alpha)m}\sum_{i\in[m]\setminus\mathcal{B}}\bar{x}^{i}-\mu|\leq t$, and $\max_{i\in[m]\setminus\mathcal{B}}\{|\bar{x}^{i}-\mu|\}\leq s$, we have

Lemma 3 can be directly applied to the $k$-th partial derivative of the loss functions. Since we assume that for any $k\in[d]$ and $\mathbf{w}\in\mathcal{W}$, $\partial_{k}f(\mathbf{w};\mathbf{z})$ is $v$-sub-exponential, we have for any $t\geq 0$, $s\geq 0$,

and consequently with probability at least

hold for all $\mathbf{w}\in\mathcal{W}$. This implies that for all $\mathbf{w}\in\mathcal{W}$ and $k\in[d]$,

The proof is completed by choosing $\delta=\frac{1}{nm\widehat{L}}$,

and using the fact that $\beta\leq\frac{1}{2}-\epsilon$.

E.2 Proof of Lemma 3

We first recall Bernstein’s inequality for sub-exponential random variables.

Suppose that $X_{1},X_{2},\ldots,X_{n}$ are i.i.d. $v$-sub-exponential random variables with mean $\mu$. Then for any $t\geq 0$,

Similarly, for any $i\in[m]\setminus\mathcal{B}$, and any $s\geq 0$

We proceed to analyze the trimmed mean of means estimator. To simplify notation, we define $\mathcal{M}=[m]\setminus\mathcal{B}$ as the set of all normal worker machines, $\mathcal{U}\subseteq[m]$ as the set of all untrimmed machines, and $\mathcal{T}\subseteq[m]$ as the set of all trimmed machines. The trimmed mean of means estimator simply computes

We also know that $|\sum_{i\in\mathcal{M}\cap\mathcal{T}}(\bar{x}^{i}-\mu)|\leq 2\beta m\max_{i\in\mathcal{M}}\{|\bar{x}^{i}-\mu|\}$. In addition, since $\beta\geq\alpha$, without loss of generality, we assume that $\mathcal{M}\cap\mathcal{T}\neq\emptyset$, and then $|\sum_{i\in\mathcal{B}\cap\mathcal{U}}(\bar{x}^{i}-\mu)|\leq\alpha m\max_{i\in\mathcal{M}}\{|\bar{x}^{i}-\mu|\}$. Then we directly obtain the desired result.

Appendix F Proof of Theorem 7

Since the loss functions are quadratic, we denote the loss function $f(\mathbf{w};\mathbf{z}^{i,j})$ by

We further define $\mathbf{H}_{i}:=\frac{1}{n}\sum_{j=1}^{n}\mathbf{H}_{i,j}$, $\mathbf{p}_{i}:=\frac{1}{n}\sum_{j=1}^{n}\mathbf{p}_{i,j}$, and $c_{i}:=\frac{1}{n}\sum_{j=1}^{n}c_{i,j}$. Thus the empirical risk function on the $i$-th machine is

Then, for any worker machine $i\in[m]\setminus\mathcal{B}$, $\widehat{\mathbf{w}}^{i}=-\mathbf{H}_{i}^{-1}\mathbf{p}_{i}$. In addition, the population risk minimizer is $\mathbf{w}^{*}=-\mathbf{H}_{F}^{-1}\mathbf{p}_{F}$. We further define $\mathbf{U}_{i,j}:=\mathbf{H}_{i,j}-\mathbf{H}_{F}$, $\mathbf{U}_{i}=\mathbf{H}_{i}-\mathbf{H}_{F}$, $\mathbf{v}_{i,j}=\mathbf{p}_{i,j}-\mathbf{p}_{F}$, and $\mathbf{v}_{i}=\mathbf{p}_{i}-\mathbf{p}_{F}$. Then

Let $\mathbf{e}_{k}$ be the $k$-th vector in the standard basis, i.e., the $k$-th column of the $d\times d$ identity matrix. We proceed to study the distribution of the $k$-th coordinate of $\widehat{\mathbf{w}}^{i}-\mathbf{w}^{*}$, $i\in[m]\setminus\mathcal{B}$, i.e.,

Similar to the proof of Theorem 1, we need to obtain a Berry-Esseen type bound for $\widehat{w}^{i}_{k}-w^{*}_{k}$. However, here, $\widehat{w}^{i}_{k}$ is not a sample mean of $n$ i.i.d. random variables, and thus we cannot directly apply the vanilla Berry-Esseen bound. Instead, we apply the following bound in Pinelis and Molzon (2016) on functions of sample means.

where $C=C_{0}+C_{1}\varsigma^{3}+(C_{20}+C_{21}\varsigma)\nu_{2}^{2}+(C_{30}+C_{31}\varsigma)\nu_{3}^{2}+C_{4}$, with

where $\|\cdot\|_{F}$ denotes the Frobenius norm of matrices. We then provide the following lemma on $\psi_{k}(\mathbf{U},\mathbf{v})$.

Lemma 4 tells us that the condition (61) is satisfied with $\theta=\frac{\lambda_{F}}{2}$ and $M_{\theta}=\frac{2\lambda_{F}+4\|\mathbf{p}_{F}\|_{2}}{\lambda_{F}^{3}}$. For all normal worker machine $i\in[m]\setminus\mathcal{B}$, denote the distribution of $\mathbf{U}_{i,j}$ and $\mathbf{v}_{i,j}$ by $\mathcal{D}_{U}$ and $\mathcal{D}_{v}$, respectively. Since $\widehat{w}^{i}_{k}-w^{*}_{k}=\psi_{k}(\frac{1}{n}\sum_{j=1}^{n}\mathbf{U}_{i,j},\frac{1}{n}\sum_{j=1}^{n}\mathbf{v}_{i,j})$, Theorem 13 directly gives us the following lemma.

We have the following lemma on $\widetilde{p}(z;k)$.

for some $\epsilon>0$. Then, with probability at least $1-4e^{-2t}$, we have

where $C_{\epsilon}$ is defined as in (4).

The proof is essentially the same as the proof of Lemma 1. One can simply replace $\sigma$ in Lemma 1 with $\widetilde{\sigma}_{k}$ and $0.4748\gamma(x)$ in Lemma 1 with $C_{k}$. Then the same arguments still apply. Thus, we skip the details of this proof. ∎

Then, define $\widehat{p}(z;k)=\frac{1}{m}\sum_{i\in[m]}\mathds{1}(\widehat{w}^{i}_{k}-w^{*}_{k}\leq z)$. Using the same arguments as in Corollary 1, we know that

which implies that $|{\sf med}\{\widehat{w}^{i}_{k}:i\in[m]\}-w^{*}_{k}|\leq C_{\epsilon}\frac{\widetilde{\sigma}_{k}}{\sqrt{n}}(\alpha+\sqrt{\frac{t}{m(1-\alpha)}}+\frac{C_{k}}{\sqrt{n}})$. Then, let

and $\widetilde{C}=\max_{k\in[d]}C_{k}$, we have with probability at least $1-4de^{-2t}$,

We complete the proof by choosing $t=\frac{1}{2}\log(nmd)$.

Let $\mathbf{H}\sim\mathcal{D}_{H}$ and $\mathbf{p}\sim\mathcal{D}_{p}$ and define

Then, $\widetilde{C}=\max_{k\in[d]}C_{k}$, with where

F.1 Proof of Lemma 4

We use $\|\cdot\|_{2}$ and $\|\cdot\|_{F}$ to denote the operator norm and the Frobenius norm of matrices, respectively. We have the identity

Let us consider the set of matrices such that $\|\mathbf{U}\|_{F}\leq\frac{\lambda_{F}}{2}$. One can check that for any such matrix, we have $\|\mathbf{H}_{F}^{-1}\mathbf{U}\|_{2}\leq\frac{1}{2}$. Let

Denote the operator norm of matrices by $\|\cdot\|_{2}$. We further have for any $r\geq 1$,

where we use the fact $\|\mathbf{U}\|_{2}\leq\|\mathbf{U}\|_{F}$. In addition, for any $r\geq 2$,

Then, we plug (71) and (72) into (70), and obtain

Appendix G Proof of Observation 1

We show that for two $d$ dimensional distributions $P_{1}=\mathcal{N}(\bm{\mu}_{1},\sigma^{2}\mathbf{I})$ and $P_{2}=\mathcal{N}(\bm{\mu}_{2},\sigma^{2}\mathbf{I})$, there exist two $dn$ dimensional distributions $Q_{1}$ and $Q_{2}$ such that

If this happens, then no algorithm can distinguish between $P_{1}$ and $P_{2}$. Let $\phi_{1}$ and $\phi_{2}$ be the PDF of $P_{1}^{n}$ and $P_{2}^{n}$, respectively. Let $\bm{\mu}_{1}$ and $\bm{\mu}_{2}$ be such that the total variation distance between $P_{1}^{n}$ and $P_{2}^{n}$ is

By the results of the total variation distance between Gaussian distributions, we know that

Let $Q_{1}$ be the distribution with PDF $\frac{1-\alpha}{\alpha}(\phi_{2}-\phi_{1})\mathds{1}_{\phi_{2}\geq\phi_{1}}$ and $Q_{2}$ be the distribution with PDF $\frac{1-\alpha}{\alpha}(\phi_{1}-\phi_{2})\mathds{1}_{\phi_{1}\geq\phi_{2}}$. One can verify that (73) is satisfied. In this case, by the lower bound in (74), we get

This implies that for two Gaussian distributions such that $\|\bm{\mu}_{1}-\bm{\mu}_{2}\|_{2}=\Omega(\frac{\alpha}{\sqrt{n}})$, in the worst case it can be impossible to distinguish these two distributions due to the existence of the adversary. Thus, to estimate the mean $\bm{\mu}$ of a Gaussian distribution in the distributed setting with $\alpha$ fraction of Byzantine machines, any algorithm that computes an estimation $\widehat{\bm{\mu}}$ of the mean has a constant probability of error $\|\widehat{\bm{\mu}}-\bm{\mu}\|_{2}=\Omega(\frac{\alpha}{\sqrt{n}})$.

Further, according to the standard results from minimax theory (Wu, 2017), we know that using $\mathcal{O}(nm)$ data, there is a constant probability that $\|\widehat{\bm{\mu}}-\bm{\mu}\|_{2}=\Omega(\sqrt{\frac{d}{nm}})$. Combining these two results, we know that $\|\widehat{\bm{\mu}}-\bm{\mu}\|_{2}=\Omega(\frac{\alpha}{\sqrt{n}}+\sqrt{\frac{d}{nm}})$.