Communication-efficient Byzantine-robust distributed learning with statistical guarantee

When the Byzantine failures occur, the aggregation rule (2.2) will be sensitive to the bad gradient values. More precisely, although the master machine communicates with the worker machines via some predefined protocol, the Byzantine machines do not have to obey this protocol and may send arbitrary messages to the master machine. At this time, the gradients $\{\nabla f_{k}(\cdot):k=2,\cdots,m+1\}$ received by the master machine are not always reliable, since the information from Byzantine machine may be completely out of its local data. To state it clearly, we assume that the Byzantine workers can provide arbitrary values written by the symbol “$*$” to the master machine. In this situation, several robust operations should be implemented to substitute the simply average of local gradients as in (2.2).

Inspired by robust techniques developed recently in Yin et al. (2018), we apply for the coordinate-wise median and coordinate-wise trimmed mean to formulate our Byzantine-robust CSL distributed learning algorithm.

See Algorithm 1 below for details, and we call it Algorithm BCSL. In each parallel iteration of Algorithm 1, the 1st machine (the normal master machine) broadcasts the current model parameter to all worker machines. The normal worker machines calculate their own gradients of loss functions based on their local data and then send them to the 1st machine. Considering that the Byzantine machines may send any messages due to their abnormal or adversarial behavior, we implement the coordinate-wise median or trimmed mean operation for these received gradients at the 1st machine. Then the aggregation algorithm in (3.1) is conducted to update the global parameter.

In order to provide an optimization error and statistical error of Algorithm 1, we need to introduce some basic conditions for our theoretical analysis.

For Option I in Algorithm 1: median-based algorithm, some moment conditions of the gradient of $\ell()$ is introduced to control stochastic behaviors.

Assumption 3.2 is standard in the literature and is satisfied in many learning problems. See Proposition 1 in Yin et al. (2018) for a specific linear regression problem.

In most existing studies, it is usually assumed that the population risk $F$ is smooth and strongly convex. The empirical risk $f$ also enjoys such good properties, as long as $\{{\mbox{\boldmath${z}$}}_{i}:i=1,\cdots,N\}$ are i.i.d. and the total sample size $N$ is sufficiently large relative to $p$ (Fan et al., 2019). From (3.1) in Algorithm 1, we know the local data at the 1st machine are used to optimize. So we need to control the gap between $\nabla^{2}f_{1}({\mbox{\boldmath${\theta}$}})$ and $\nabla^{2}F({\mbox{\boldmath${\theta}$}})$ to contract optimization rate of the algorithm. The similarity between $f_{1}$ and $F$ is depicted by Assumption 3.4. Indeed, the empirical risk $f_{1}$ should not be too far away from their population risk $F$ as long as the sample size $n$ of the 1st machine is not too small. Mei et al. (2018) showed it holds with reasonably small $\delta$ and large $R$ with high probability under general conditions. Specially, a large $n$ implies a small homogeneity index $\delta$. Obviously, Assumption 3.4 always holds when taking $\delta=\sup_{{\mbox{\boldmath${\theta}$}}\in B(\hat{{\mbox{\boldmath${\theta}$}}},R)}\left(\|\nabla^{2}f_{1}({\mbox{\boldmath${\theta}$}})\|_{2}+\|\nabla^{2}F({\mbox{\boldmath${\theta}$}})\|_{2}\right)$.

The following theorem establishes the global convergence of the proposed algorithm BCSL, involving a trade off between the optimization error and statistical error.

Theorem 3.1 shows the linear convergence of ${\mbox{\boldmath${\theta}$}}_{t}$, which depends explicitly on the homogeneity index $\delta$, the strong convex index $\rho$, and the fraction of Byzantine machines $\alpha$. The result is viewed as an extension of that in Jordan et al. (2019) and Fan et al. (2019). A significant difference from theirs is that, we allow the initial estimator to be inaccurate, and with high probability we have more explicit rates of convergence on optimization errors under the Byzantine failures.

Specially, the factor $C_{\varepsilon}$ in Theorem 3.1 is a function of $\varepsilon$, for example, $C_{\varepsilon}\approx 4$ if setting $\varepsilon=1/6$. After running $T$ for Algorithm 1, with high probability, we have

Notice that $\log(1-x)\leq-x$. Theorem 3.1 guarantees that after parallel iterating $T\geq\frac{\rho}{\rho-\delta}\log\frac{(\rho-\delta)\|{\mbox{\boldmath${\theta}$}}_{0}-\hat{{\mbox{\boldmath${\theta}$}}}\|_{2}}{2\Delta_{nm\alpha}}$, with high probability we can obtain a solution $\tilde{{\mbox{\boldmath${\theta}$}}}={\mbox{\boldmath${\theta}$}}_{T}$ with an error

In this case, the derived rate of the statistical error between $\tilde{{\mbox{\boldmath${\theta}$}}}$ and the centralized empirical risk minimizer $\hat{{\mbox{\boldmath${\theta}$}}}$ is of the order $\mathcal{O}\left(\frac{\alpha}{\sqrt{n}}+\sqrt{\frac{p\log(nm)}{nm}}+\frac{1}{n}\right)$ up to some constants, alternatively, $\tilde{\mathcal{O}}\left(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}+\frac{1}{n}\right)$ up to the logarithmic factor. Note that

Intuitively, the above error rate is a near 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, $\frac{1}{\sqrt{m}}$ is the averaging effect of $m$ normal machines, and $\frac{1}{n}$ is the effect of the dependence of median on skewness of the gradients. If $n\gtrsim m$, then $\tilde{\mathcal{O}}\left(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}+\frac{1}{n}\right)=\tilde{\mathcal{O}}\left(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}\right)$ is the order-optimal rate (Yin et al., 2018). When $\alpha=0$ (no Byzantine machine), one sees the usual scaling $\frac{1}{\sqrt{nm}}$ with the global sample size; when $\alpha\neq 0$ (some machines are Byzantine), their influence remains bounded and is proportional to $\alpha$. So we do not sacrifice the quality of learning to guard against Byzantine failures, provided that the Byzentine failure proportion satisfies $\alpha=\mathcal{O}(1/\sqrt{m})$ .

Remark also that, our results for convex problems follow for any finite-bounded $R$ of the initial radius.

We next turn to an analysis for Option II in Algorithm 1: The robust distributed learning based on coordinate-wise trimmed mean. Compared to Option I, a stronger assumption on the tail behavior of the partial derivatives of the loss functions is needed as follows.

The sub-exponential assumption implies that all the moments of the derivatives are bounded. Hence, this condition is stronger a little than the bounded absolute skewness (Assumption 3.2(ii)). Fortunately, Assumption 3.5 can be satisfied in some learning problems, and see Proposition 2 in Yin et al. (2018).

Theorem 3.2 also shows the linear convergence of ${\mbox{\boldmath${\theta}$}}_{t}$, which depends explicitly on the homogeneity index $\delta$, the strong convex index $\rho$, and the trimmed mean index $\beta$ with choosing the index to satisfy $\beta\geq\alpha$, where $\alpha$ is a fraction of Byzantine machines. Note that, the hyperparameter $\upsilon$ in Assumption (3.5) only affect the statistical error appearing in $\Delta_{nm\beta}$.

Similar to (3.3), we also have

after $T$-step iterations. By running $T\geq\frac{\rho}{\rho-\delta}\log\frac{(\rho-\delta)\|{\mbox{\boldmath${\theta}$}}_{0}-\hat{{\mbox{\boldmath${\theta}$}}}\|_{2}}{2\Delta_{nm\beta}}$ parallel computations, we can obtain a solution $\tilde{\tilde{{\mbox{\boldmath${\theta}$}}}}={\mbox{\boldmath${\theta}$}}_{T}$ satisfying $\|\tilde{\tilde{{\mbox{\boldmath${\theta}$}}}}-\hat{{\mbox{\boldmath${\theta}$}}}\|_{2}=\mathcal{\tilde{O}}\left(\frac{\beta}{\sqrt{n}}+\frac{1}{\sqrt{nm}}\right)$, since the term $\Delta_{nm\beta}$ can be reduced to be $\Delta_{nm\beta}=\mathcal{O}\left(\frac{\upsilon p}{\varepsilon}\left[\frac{\beta}{\sqrt{n}}+\frac{1}{\sqrt{nm}}\right]\sqrt{\log(nm)}\right)$.

It should be pointed out that, the trimmed mean index $\beta$ is strictly controlled by the fraction of Byzantine machines $\alpha$, that is $\frac{1}{2}-\varepsilon\geq\beta\geq\alpha$. By choosing $\beta=c\alpha$ with $c\geq 1$, we still achieve the optimization error rate $\mathcal{\tilde{O}}\left(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}\right)$, which is also order-optimal in the statistical literature.

We now take comparable analysis for the above two Byzantine-robust CSL distributed learning in Algorithm 1 (Options I and II). The trimmed-mean-based algorithm (Option II) has an order-optimal optimization error rate $\mathcal{\tilde{O}}\left(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}\right)$. By contrast, the median-based algorithm (Option I) has the rate $\tilde{\mathcal{O}}\left(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}+\frac{1}{n}\right)$ involving an additional term $\frac{1}{n}$, and thus the optimality is achieved only for $n\succeq m$. Note that Option I algorithm needs milder moment conditions (bounded skewness) on the tail of the loss derivatives than the Option II algorithm (sub-exponentiality). In other words, this provides a profound insight into the underlying relation between the tail decay of the loss derivatives and the block number of local machines ($m$).

On the other hand, Algorithm 1 based on Option II has an additional parameter $\beta$ such that $1/2>\beta\geq\alpha$, which requires that the fraction of Byzantine machines is absolutely dominated by the normal ones for guaranteeing robustness. In contrast to this, Algorithm 1 based on Option I has a weaker restriction on $\alpha$.

We now unpack Theorems 3.1 and 3.2 in generalized linear models, taking into consideration the effects of iterations in the proposed estimator, the roles of the initial estimator and the Byzantine failures. We will find an explicit rate of convergence of $\delta$ in Assumption 3.4 in the setting of generalized linear models. Theorem 3.1 and 3.2 guarantee that after running $T$ parallel iterations, with high probability we achieve an optimization error with a linear rate. Moreover, through finite steps, the optimization errors are eventually negligible in comparison with the statistical errors. We will give a specific analysis below.

For GLMs, the loss function is the negative partial log-likelihood of an exponential-type variable of the response given any input feature ${x}$. Suppose that the i.i.d. pairs $\{{\mbox{\boldmath${z}$}}_{i}=({\mbox{\boldmath${x}$}}_{i}^{T},y_{i})^{T}\}_{i=1}^{N}$ are drawn from a generalized linear model. Recall that the conditional density of $y_{i}$ given ${\mbox{\boldmath${x}$}}_{i}$ takes the form

where $b(\cdot)$ is some known convex function, and $c(\cdot)$ is a known function such that $h(\cdot)$ is a valid probability density function. The loss function corresponding to the negative log likelihood of the whole data is given by $f({\mbox{\boldmath${\theta}$}})=\frac{1}{m+1}\sum_{k=1}^{m+1}f_{k}({\mbox{\boldmath${\theta}$}})$ with

We further impose the following standard regularity technical conditions.

Assumptions (3.6)(i)(ii)(iv)(v)-(3.7) have been proposed by Fan et al. (2019) to establish the rate of optimization errors for a distributed statistical inference. In our paper, these assumptions mainly are used to obtain asymptotic properties of ${\mbox{\boldmath${\theta}$}}_{t}$ in Option I algorithm. For establishing the rate of optimization errors of Option II algorithm, we need Assumption (3.6) (iii) (sub-exponentiality) instead of Assumption (3.6)(ii) (sub-Gaussian), which is similar to Theorem 3.2.

From (3.4), we easily obtain

By Lemma A.5 in Fan et al. (2019), we immediately get

as long as $n>cp$ for a given positive constant $c$. So, $\delta=O\left(\|{\mbox{\boldmath${\Sigma}$}}\|_{2}\sqrt{p(\log p+\log n)/n}\right)$ can be chosen with high probability. From Theorems 3.1 and 3.2, we have the contraction factor $O\left(\kappa\sqrt{p(\log p+\log n)/n}\right)$ with $\kappa=\|{\mbox{\boldmath${\Sigma}$}}\|_{2}/\rho$. The explicit parameter $\kappa$ is comparable to the condition number in Jordan et al. (2019), where finite $p$ and $\kappa$ were imposed.

Equipped with these above facts, we have the following corollary.

Theorems 3.3 and 3.4 clearly present how Algorithm 1 depend on structural parameters of the problem. When $\kappa$ is bounded, $\eta=\mathcal{O}(p(\log p+\log n)/n)$, through finite steps, $\eta^{t/2}$ can be much smaller than $\Delta_{nm\alpha}(\Delta_{nm\beta})$. Thus,

So, By contrast to that results in Jordan et al. (2019), we allow an inaccurate initial value ${\mbox{\boldmath${\theta}$}}_{0}$ and give more explicit rates of optimization error even when $p$ diverges.

We know that the statistical error of the estimator ${\mbox{\boldmath${\theta}$}}_{t}$ can be controlled by the optimization error of ${\mbox{\boldmath${\theta}$}}_{t}$ and statistical error of $\hat{{\mbox{\boldmath${\theta}$}}}_{t}$, that is,

Note that the first term is not a deterministic optimization error, it holds in probability. Therefore, it is an optimization error in a statistical sense. We call it statistical optimization error. The second term is of order $\mathcal{O}_{p}(\sqrt{p/(nm)})$ under mild conditions, which has been well studied in statistics. Thus, the statistical error of ${\mbox{\boldmath${\theta}$}}_{t}$ is controlled by the magnitude of the first term. If we adopt two-step iteration, and $\|{\mbox{\boldmath${\theta}$}}_{0}-{\mbox{\boldmath${\theta}$}}^{*}\|_{2}=\mathcal{O}_{p}(\sqrt{p/n})$ when ${\mbox{\boldmath${\theta}$}}_{0}$ is obtained in 1st machine, one gets

by (3.5)-(3.6), where $E_{n}$ is defined in (3.5), provided that for Option I

or for Option II

It implies that Algorithm 1 has the order-optimal statistical error rate, for Option I which needs $n\gtrsim m$.

First, recall that the proximal operator $\mathrm{prox}_{h}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$ is defined by

By the proximal operator of the function $\lambda^{-1}h$ with $\lambda>0$, the proximal algorithm for minimizing $h$ iteratively computes

starting from some initial value ${\mbox{\boldmath${x}$}}_{0}$. Rockafellar (1976) showed the $\{{\mbox{\boldmath${x}$}}_{t}\}_{t=0}^{\infty}$ converges linearly to some $\hat{{\mbox{\boldmath${x}$}}}\in\mathrm{argmin}_{\mathbb{R}^{p}}h({\mbox{\boldmath${x}$}})$.

For our problem (2.1), the proximal iteration algorithm is

In our setting, we adopt the distributed learning, and optimization is mainly on the 1st worker machine. The $g(\cdot)$ is a penalty function, which is used to the optimization step, and keep the local data of the 1st worker machine in $f_{1}(\cdot)$. So, the penalty function in our proximal algorithm becomes $g({\mbox{\boldmath${\theta}$}})+\frac{\lambda}{2}\|{\mbox{\boldmath${\theta}$}}-{\mbox{\boldmath${\theta}$}}_{t}\|_{2}^{2}$. Further, the optimization (3.1) in Algorithm 1 is replaced by

where $h({\mbox{\boldmath${\theta}$}}_{t})$ is the gradient information at ${\mbox{\boldmath${\theta}$}}_{t}$ from the other worker machines. This optimization (4.1) can make the Byzantine-robust CSL-proximal distributed learning converges rapidly. See the following Algorithm 2. We call it Algorithm BCSLp.

The above Algorithm 2 is a communication-efficient Byzantine-robust accurate statistical learning, which adopts coordinate-wise median and coordinate-wise trimmed mean cope with Byzantine fails, and use the proximal algorithm as the backbone. In each iteration, it has one round of communication and one optimization step similar to Algorithm 1. It is regularized version of Algorithm 1 by adding a strict convex quadratic term in the objective function. The technique has been used in the distributed stochastic optimization such as accelerating first-order algorithm (Lee et al., 2017) and regularizing sizes of updates (Wang et al., 2017b), and in the communication-efficient accurate distributed statistical estimation (Fan et al., 2019).

Now, we give contraction guarantees for Algorithm 2.

Theorems 4.1 and 4.2 present the linear convergence of Algorithm 2 Options I and II, respectively. Obviously, the contraction factor consists of two parts: $\frac{\frac{\delta}{\rho+\lambda}\sqrt{\rho^{2}+2\lambda\rho}}{\rho+\lambda}$ which comes from the error of the inexact proximal update $\|{\mbox{\boldmath${\theta}$}}_{t+1}-\mathrm{prof}_{\lambda^{-1}(f+g)}({\mbox{\boldmath${\theta}$}}_{t})\|_{2}$, and $\frac{\lambda}{\rho+\lambda}$ which comes from the residual of proximal point $\|\mathrm{prof}_{\lambda^{-1}(f+g)}({\mbox{\boldmath${\theta}$}}_{t})-\hat{{\mbox{\boldmath${\theta}$}}}\|_{2}$. Remarking similar to Theorems 3.1 and 3.2, within a finite $T$-step , we have $\|{\mbox{\boldmath${\theta}$}}_{T}-\hat{{\mbox{\boldmath${\theta}$}}}\|_{2}=\tilde{\mathcal{O}}_{p}\left(\frac{\alpha}{\sqrt{n}}+\frac{1}{\sqrt{nm}}+\frac{1}{n}\right)$ for Option I, which is a order-optimal if $n\gtrsim m$; and $\|{\mbox{\boldmath${\theta}$}}_{T}-\hat{{\mbox{\boldmath${\theta}$}}}\|_{2}=\mathcal{\tilde{O}}_{p}\left(\frac{\beta}{\sqrt{n}}+\frac{1}{\sqrt{nm}}\right)$ for Option II, which also is a order optimal. These results just need that $\{f_{k}\}_{1}^{m+1}$ are convex and smooth while the penalty $g$ is allowed to be non-smooth, for example, $\ell_{1}$ norm. However, Most distributed statistical learning algorithms are only designed for smooth problems, and don’t consider Byzantine problems, for instance, Shamir et al. (2014), Wang et al. (2017a), Jordan et al. (2019), and so on.