Data Encoding for Byzantine-Resilient Distributed Optimization

TL;DR

This paper introduces a data encoding method with error correction to enable Byzantine-resilient distributed optimization, tolerating up to half the worker nodes being malicious, with efficient schemes for PGD and CD algorithms.

Contribution

It proposes the first secure coordinate descent method against adversarial attacks using a novel data encoding and error correction scheme.

Findings

Tolerates up to .5m corrupt nodes, optimal by information theory.

Achieves constant resource overhead for up to m/3 corrupt nodes.

Extends efficiently to streaming data and stochastic gradient descent.

Abstract

We study distributed optimization in the presence of Byzantine adversaries, where both data and computation are distributed among $m$ worker machines, $t$ of which may be corrupt. The compromised nodes may collaboratively and arbitrarily deviate from their pre-specified programs, and a designated (master) node iteratively computes the model/parameter vector for generalized linear models. In this work, we primarily focus on two iterative algorithms: Proximal Gradient Descent (PGD) and Coordinate Descent (CD). Gradient descent (GD) is a special case of these algorithms. PGD is typically used in the data-parallel setting, where data is partitioned across different samples, whereas, CD is used in the model-parallelism setting, where data is partitioned across the parameter space. In this paper, we propose a method based on data encoding and error correction over real numbers to combat adversarial attacks. We can tolerate up to $t\leq \lfloor\frac{m-1}{2}\rfloor$ corrupt worker nodes, which is information-theoretically optimal. We give deterministic guarantees, and our method does not assume any probability distribution on the data. We develop a {\em sparse} encoding scheme which enables computationally efficient data encoding and decoding. We demonstrate a trade-off between the corruption threshold and the resource requirements (storage, computational, and communication complexity). As an example, for $t\leq\frac{m}{3}$, our scheme incurs only a {\em constant} overhead on these resources, over that required by the plain distributed PGD/CD algorithms which provide no adversarial protection. To the best of our knowledge, ours is the first paper that makes CD secure against adversarial attacks. Our encoding scheme extends efficiently to the data streaming model and for stochastic gradient descent (SGD). We also give experimental results to show the efficacy of our proposed schemes.

Equations

$$\displaystyle\min_{{\bf w}\in\mathbb{R}^{d}}\left(\bigg{(}\frac{1}{n}\sum_{i=1}^{n}f_{i}({\bf w})\bigg{)}+h({\bf w})\right),$$

$${\mathbf{w}}_{t+1}={\text{prox}}_{h,{\alpha }_t}({\mathbf{w}}_t-{\alpha }_t\nabla f({\mathbf{w}}_t)),t=1,2,3,\dots$$

$${\text{prox}}_{h,\alpha }(\mathbf{w})=\arg \underset{\mathbf{z}\in {\mathbb{R}}^d}{\min }\frac{1}{2\alpha }\Vert \mathbf{z}-\mathbf{w}{\Vert }_2^2+h(\mathbf{z}).$$

$$(S_{\lambda \alpha }(\mathbf{z}){)}_j=\{\begin{array}{ll}{z}_j+\lambda \alpha & \text{ if }{z}_j<-\lambda \alpha ,\\ 0& \text{ if }-\lambda \alpha \le z_j\le \lambda \alpha ,\\ z_j-\lambda \alpha & \text{ if }{z}_j>\lambda \alpha .\end{array}$$

$$\underset{\mathbf{w}\in \mathcal{C}}{\min }f(\mathbf{w})⟺\underset{\mathbf{w}\in {\mathbb{R}}^d}{\min }f(\mathbf{w})+I_{\mathcal{C}}(\mathbf{w}).$$

$$f_i(\mathbf{w})=-y_i\log \left(\frac{1}{1+e^{-u_i}}\right)-(1-y_i)\log \left(\frac{e^{-u_i}}{1+e^{-u_i}}\right),$$

$$\nabla f({\mathbf{w}}_t)=\frac{1}{m}\sum _{i=1}^m{\nabla }_{i}f({\mathbf{w}}_t).$$

$$\underset{\mathbf{w}\in {\mathbb{R}}^d}{\min }\left(\varphi (\mathbf{X}\mathbf{w};\mathbf{y}):=\sum _{i=1}^n\ell (⟨{\mathbf{x}}_i,\mathbf{w}⟩;y_i)\right).$$

$${\mathbf{w}}_{i\mathcal{U}}\leftarrow {\mathbf{w}}_{i\mathcal{U}}-\alpha {\nabla }_{i\mathcal{U}}\varphi (\mathbf{X}\mathbf{w};\mathbf{y})$$

$$\nabla f(\mathbf{w})={\mathbf{X}}^T{f}^{\prime }(\mathbf{w}),\forall \mathbf{w}\in {\mathbb{R}}^d.$$

$${\tilde{h}}_i(\mathbf{v})={\tilde{\mathbf{S}}}_i\mathbf{A}\mathbf{v}+{\tilde{\mathbf{e}}}_i,i\in [p]$$

$${\mathbf{f}}_i:=\mathbf{F}{\tilde{h}}_i(\mathbf{v})=\mathbf{F}{\tilde{\mathbf{e}}}_i,$$

$${\mathbf{F}}^{\perp }={\left[\begin{array}{cccc}{b}_{11}& b_{21}& \dots & b_{q1}\\ b_{12}& b_{22}& \dots & b_{q2}\\ ⋮& ⋮& ⋮& ⋮\\ b_{1m}& b_{2m}& \dots & b_{qm}\end{array}\right]}_{m\times q}$$

$${\tilde{\mathbf{S}}}_i=\left[\begin{array}{cccccccccc}0& \dots & 0& b_{11}& b_{21}& \dots & b_{l1}& 0& \dots & 0\\ 0& \dots & 0& b_{12}& b_{22}& \dots & b_{l2}& 0& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& \dots & 0& b_{1m}& b_{2m}& \dots & b_{lm}& 0& \dots & 0\end{array}\right]$$

$${\mathbf{S}}_i={\left[\begin{array}{cccc}{b}_{1i}\dots b_{qi}& & & \\ & \ddots & & \\ & & b_{1i}\dots b_{qi}& \\ & & & b_{1i}\dots b_{li}\end{array}\right]}_{p\times n_r}$$

$${\mathbf{r}}_{\mathcal{T}}={\mathbf{S}}_{\mathcal{T}}\mathbf{A}\mathbf{v},$$

$${\tilde{\mathbf{r}}}_j$$

$$\mathbf{F}={\left[\begin{array}{ccccc}1& 1& 1& \dots & 1\\ z_1& z_2& z_3& \dots & z_m\\ z_1^2& z_2^2& z_3^2& \dots & z_m^2\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ z_1^{k-1}& z_2^{k-1}& z_3^{k-1}& \dots & z_m^{k-1}\end{array}\right]}_{k\times m}$$

$$\displaystyle\big{(}p_{1}(d+1)$$

$$<[(1+ϵ)n+m](d+1)+[(1+ϵ)d+m]n$$

$$=(1+ϵ)n(2d+1)+m(n+d+1).$$

$$\arg \underset{\mathbf{w}\in {\mathbb{R}}^d}{\min }\varphi (\mathbf{X}\mathbf{w}):=\sum _{i=1}^n\ell (⟨{\mathbf{x}}_i,\mathbf{w}⟩).$$

$$\arg \underset{\mathbf{v}\in {\mathbb{R}}^{pm}}{\min }\varphi ({\tilde{\mathbf{X}}}^R\mathbf{v}).$$

$${\mathbf{v}}_{i\mathcal{U}}^{t+1}={\mathbf{v}}_{i\mathcal{U}}^t-{\alpha }_t{\nabla }_{i\mathcal{U}}\varphi ({\tilde{\mathbf{X}}}^R{\mathbf{v}}^t),$$

$${\mathbf{v}}_{i\mathcal{U}}^{t+1}={\mathbf{v}}_{i\mathcal{U}}^t-{\alpha }_t{\mathbf{R}}_{i\mathcal{U}}^T{\mathbf{X}}^T{\varphi }^{\prime }({\tilde{\mathbf{X}}}^R{\mathbf{v}}^t).$$

$${\mathbf{w}}_{f(\mathcal{U})}^{t+1}={\mathbf{w}}_{f(\mathcal{U})}^t-{\alpha }_t{\mathbf{X}}_{f(\mathcal{U})}^T{\varphi }^{\prime }(\mathbf{X}{\mathbf{w}}^t).$$

$${\mathbf{S}}_i={\left[\begin{array}{cccc}{b}_{1i}\dots b_{qi}& & & \\ & \ddots & & \\ & & b_{1i}\dots b_{qi}& \\ & & & b_{1i}\dots b_{li}\end{array}\right]}_{p\times d}$$

$$f(i):=\{\begin{array}{ll}[(i-1)\ast q+1:i\ast q]& \text{ if }1\le i<p,\\ [(p-1)\ast q+1:d]& \text{ if }i=p,\end{array}$$

$${\mathbf{S}}_{i\mathcal{U}}\mathbf{w}={\mathbf{S}}_{i\mathcal{U},f(\mathcal{U})}{\mathbf{w}}_{f(\mathcal{U})},\text{ for every }\mathbf{w}\in {\mathbb{R}}^d,$$

$${\mathbf{S}}_{i\mathcal{U}}{\mathbf{X}}^T={\mathbf{S}}_{i\mathcal{U},f(\mathcal{U})}{\mathbf{X}}_{f(\mathcal{U})}^T.$$

$${\mathbf{v}}_{i\mathcal{U}}^{t+1}={\mathbf{S}}_{i\mathcal{U},f(\mathcal{U})}{\mathbf{w}}_{f(\mathcal{U})}^{t+1},$$

Full text

At the core of our solutions to both these algorithms is a method for Byzantine-resilient matrix-vector (MV) multiplication; and for that, we propose a method based on data encoding and error correction over real numbers to combat adversarial attacks. We can tolerate up to $t\leq\lfloor\frac{m-1}{2}\rfloor$ corrupt worker nodes, which is information-theoretically optimal. We give deterministic guarantees, and our method does not assume any probability distribution on the data. We develop a sparse encoding scheme which enables computationally efficient data encoding and decoding. We demonstrate a trade-off between the corruption threshold and the resource requirements (storage, computational, and communication complexity). As an example, for $t\leq\frac{m}{3}$, our scheme incurs only a constant overhead on these resources, over that required by the plain distributed PGD/CD algorithms which provide no adversarial protection. To the best of our knowledge, ours is the first paper that connects MV multiplication with CD and designs a specific encoding matrix for MV multiplication whose structure we can leverage to make CD secure against adversarial attacks.

Our encoding scheme extends efficiently to (i) the data streaming model, in which data samples come in an online fashion and are encoded as they arrive, and (ii) making stochastic gradient descent (SGD) Byzantine-resilient. In the end, we give experimental results to show the efficacy of our proposed schemes.

1 Introduction

Map-reduce architecture [DG08] is implemented in many distributed learning tasks, where there is one designated machine (called the master) that computes the model iteratively, based on the inputs from the worker machines at each iteration, typically using descent techniques, like (proximal) gradient descent, coordinate descent, stochastic gradient descent, the Newton’s method, etc. The worker nodes perform the required computations using local data, distributed to the nodes [ZWLS10]. Several other architectures, including having no hierarchy among the nodes have been explored [LZZ*+*17].

In several applications of distributed learning, including the Internet of Battlefield Things (IoBT) [A*+*18], federated optimization [Kon17], the recruited worker nodes might be partially trusted with their computation. Therefore, an important question is whether we can reliably perform distributed computation, taking advantage of partially trusted worker nodes. These Byzantine adversaries can collaborate and arbitrarily deviate from their pre-specified programs. The problem of distributed computation with Byzantine adversaries has a long history [LSP82], and there has been recent interest in applying this computational model to large-scale distributed learning [BMGS17, CWCP18, CSX17].

In this paper, we study Byzantine-tolerant distributed optimization to learn a regularized generalized linear model (GLM) (e.g., linear/ridge regression, logistic regression, Lasso, SVM dual, constrained minimization, etc.). We consider two frameworks for distributed optimization: (i) data-parallelism architecture, where data points are distributed across different worker nodes, and in each iteration, they all parallelly compute gradients on their local data and master aggregates them to update the parameter vector using gradient descent (GD) [BT89, Bot10, DCM*+*12]; and (ii) model-parallelism architecture, where data points are partitioned across features, and several worker nodes work in parallel, updating different subsets of coordinates of the model/parameter vector through coordinate descent (CD) [BKBG11, Wri15, RT16]. Note that GD requires full gradients to update the parameter vector; and if full gradients are too costly to compute, we can reduce the per-iteration cost by using CD,111Alternatively, we can also use SGD to reduce the per-iteration cost, and we give a method for making SGD Byzantine-resilient in Section 6.1. which also has been shown to be very effective for solving generalized linear models, and is particularly widely used for sparse logistic regression, SVM, and Lasso [BKBG11]. Given its simplicity and effectiveness, CD can be chosen over GD in such applications [Nes12]. Computing gradients in the presence of Byzantine adversaries has been recently studied [BMGS17, CSX17, CWCP18, YCRB18, AAL18, SX19, XKG19, YCRB19, GV19, RWCP19, LXC*+*19, GHYR19, YLR*+*19, DD20b, DD20a, HKJ20], and we discuss them in detail Section 3 where we also put our work in context. However, as far as we know, making CD robust to Byzantine adversaries has not received much attention, and to the best of our knowledge, ours is the first paper that studies CD against Byzantine attacks and provides an efficient solution for that.

1.1 Our Contributions

We propose Byzantine-resilient distributed optimization algorithms both for PGD and CD based on data encoding and error correction (over real numbers). As mentioned above, there have been several papers that provide different methods for gradient computation in the presence of Byzantine adversaries, however, our proposed algorithm differs from them in one or more of the following aspects: (i) it does not make statistical assumptions on the data or Byzantine attack patterns; (ii) it can tolerate up to a constant fraction ($<1/2$) of the worker nodes being Byzantine, which is information-theoretically optimal; and (iii) it enables a trade-off (in terms of storage and computation/communication overhead at the master and the worker nodes) with Byzantine adversary tolerance, without compromising the efficiency at the master node. We give the same guarantees for CD also.

First we design a coding scheme for distributed matrix-vector (MV) multiplication, specifically, for operating in the presence of Byzantine adversaries, and use that in both our algorithms for PGD and CD to learn GLMs. Note that the connection of MV multiplication with gradient computation is straightforward and has been known for some time (see, for example, [LLP*+*18, DCG16]), however, it is not clear whether we can use MV multiplication methods for CD also. Indeed, since each CD update has a different requirement than that of gradient computation, a general-purpose algorithm for MV multiplication may not be applicable for CD. One distinction is that in gradient computation, we only need to encode the data to compute the MV multiplication, whereas, in CD, in addition to data encoding, since workers update few coordinates of different parts of the parameter vector in parallel, we need to encode the parameter vector as well for master to be able to decode that. In this paper, we design our encoding matrix for MV multiplication in such a way that it is sparse and has a regular structure of non-zero entries (see (11) for the encoding matrix for any worker), which makes it applicable for CD too. This leads to efficient solutions for both PGD and CD, which are our main focus in this paper.

Inspired from the real-error correction (or sparse reconstruction) problem [CT05], we develop efficient encoding/decoding procedures for MV multiplication, where we encode the data matrix and distribute it to the $m$ worker nodes, and to recover the MV product at the master, we reduce the decoding problem to the sparse reconstruction or real-error correction problem [CT05]. Note that in PGD, we only need to encode the data, whereas, in CD, we also need to encode the parameter vector, and our coding scheme should facilitate the requirement that the update on a small fraction of the encoded parameter vector should affect only a small fraction of the original parameter vector. This is a non-trivial requirement, and our coding scheme for MV multiplication is designed in such a way that it supports this requirement in an efficient manner; see Section 2.2 for a description on plain distributed CD, Section 2.5 for our approach to making CD robust to Byzantine attacks, and Section 5 for a complete solution for Byzantine-resilient CD. In the context of PGD/CD, for decoding, the master node processes the inputs from the worker nodes, either to compute the true gradient in the case of PGD or to facilitate the computation at the worker nodes in the case of CD. We take a two-round approach in each iteration of both these algorithms. Our main results are summarized in Theorem 1 (on page 1) for PGD and Theorem 2 (on page 2) for CD, and demonstrate a trade-off between the Byzantine resilience (in terms of the number of adversarial nodes) and the resource requirement (storage, computational, and communication complexity). As an example, for $t\leq\frac{m}{3}$, our scheme incurs only a constant overhead on these resources, over that required by the plain distributed PGD and CD algorithms which provide no adversarial protection. Our coding schemes can handle both Byzantine attacks and missing updates (e.g., caused by delay or asynchrony of worker nodes). Our encoding process is also efficient. Though data encoding is a one-time process, it has to be efficient to harness the advantage of distributed computation. We design a sparse encoding process, based on real-error correction, which enables efficient encoding, and the worker nodes encode data using the sparse structure. This allows encoding with storage redundancy222Storage redundancy is defined as the ratio of the size of the encoded matrix and the size of the raw data matrix. of $\frac{2m}{m-2t}$ (which is a constant, even if $t$ is a constant ($<\frac{1}{2}$) fraction of $m$), and a one-time total computation cost for encoding is $O((1+2t)nd)$. Note that the time for data encoding is a factor of $(1+2t)$ (where $t$ is the corruption threshold) more than the time required for plain data distribution which is $O(nd)$, the size of the data matrix.

We extend our encoding scheme in a couple of important ways: first, to make the stochastic gradient descent (SGD) algorithm Byzantine-resilient without compromising much on the resource requirements; and second, to handle streaming data efficiently, where data points arrives one by one (and we encode them as they arrive), rather than being available at the beginning of the computation; we also give few more applications of our method. For the streaming model, more specifically, our encoding requires the same amount of time, irrespective of whether we encode all the data at once, or we get data points one by one (or in batches) and we encode them as they arrive. This setting encompasses a more realistic scenario, in which we design our coding scheme with the initial set of data points and distribute the encoded data among the workers. Later on, when we get some more samples, we can easily incorporate them into our existing encoded setup. See Section 6 for details on these extensions.

1.2 Paper Organization

We present our problem formulation, description of the plain distributed PGD and CD algorithms, and the high-level ideas of our Byzantine-resilient algorithms for both PGD and CD along-with our main results in Section 2. We give detailed related work in Section 3. We present our full coding schemes for MV multiplication and also for gradient computation for PGD along-with a complete analysis of their resource requirements in Section 4. In Section 5, we provide a complete solution to CD. In Section 6, we show how our method can be extended to SGD and to the data streaming model. We also discuss applicability of our method to a few more important applications in that section. In Section 7, we show numerical results of our method: we show the efficiency of our method for both gradient descent (GD) and coordinate descent (CD) by running them to solve linear regression on two datasets (moderate and large) and plotting the running time with varying number of corrupt worker nodes (up to <1/2 fraction).

1.3 Notation

We denote vectors by bold small letters (e.g., ${\bf x},{\bf y},{\bf z}$, etc.) and matrices by bold capital letters (e.g., ${\bf A},{\bf F},{\bf S},{\bf X}$, etc.). We denote the amount of storage required by a matrix ${\bf X}$ by $|{\bf X}|$. For any positive integer $n\in\mathbb{N}$, we denote the set $\{1,2,\ldots,n\}$ by $[n]$. For $n_{1},n_{2}\in\mathbb{N}$, where $n_{1}\leq n_{2}$, we write $[n_{1}:n_{2}]$ to denote the set $\{n_{1},n_{1}+1,\ldots,n_{2}\}$. For any vector ${\bf u}\in\mathbb{R}^{n}$ and any set $\mathcal{S}\subset[n]$, we write ${\bf u}_{\mathcal{S}}$ to denote the $|\mathcal{S}|$-length vector, which is the restriction of ${\bf u}$ to the coordinates in the set $\mathcal{S}$. The support of a vector ${\bf u}\in\mathbb{R}^{n}$ is defined by $\textsf{supp}({\bf u}):=\{i\in[n]:u_{i}\neq 0\}$. We say that a vector ${\bf u}\in\mathbb{R}^{n}$ is $t$-sparse if $|\textsf{supp}({\bf u})|\leq t$. While stating our results, we assume that performing the basic arithmetic operations (addition, subtraction, multiplication, and division) on real numbers takes unit time.

2 Problem Setting and Our Results

Given a dataset consisting of $n$ labelled data points $({\bf x}_{i},y_{i})\in\mathbb{R}^{d}\times\mathbb{R}$, $i\in[n]$, we want to learn a model/parameter vector ${\bf w}\in\mathbb{R}^{d}$, which is a minimizer of the following empirical risk minimization problem:

[TABLE]

where $f_{i}({\bf w})$, $i=1,2,\ldots,n$, denotes the risk associated with the $i$’th data point with respect to ${\bf w}$ and $h({\bf w})$ denotes a regularizer. We call $f({\bf w}):=\frac{1}{n}\sum_{i=1}^{n}f_{i}({\bf w})$ the average empirical risk associated with the $n$ data points with respect to ${\bf w}$. Our main focus in this paper is on generalized linear models (GLM), where $f_{i}({\bf w})=\ell(\langle{\bf x}_{i},{\bf w}\rangle;y_{i})$ for some differentiable loss function $\ell$. Here, each $f_{i}:\mathbb{R}^{d}\to\mathbb{R}$ is differentiable, $h:\mathbb{R}^{d}\to\mathbb{R}$ is convex but not necessarily differentiable, and $\langle{\bf x}_{i},{\bf w}\rangle$ is the dot product of ${\bf x}_{i}$ and ${\bf w}$. We do not necessarily need each $f_{i}$ to be convex, but we require $f({\bf w})$ to be a convex function. Note that $f({\bf w})+h({\bf w})$ is a convex function. In the following we study different algorithms for solving (1) to learn a GLM.

2.1 Proximal Gradient Descent

We can solve (1) using Proximal Gradient Descent (PGD). This is an iterative algorithm, in which we choose an arbitrary/random initial ${\bf w}_{0}\in\mathbb{R}^{d}$, and then update the parameter vector according to the following update rule:

[TABLE]

where $\alpha_{t}$ is the step size or the learning rate at the $t$’th iteration, determining the convergence behaviour. There are standard choices for it; see, for example, [BV04, Chapter 9]. For any $h$ and $\alpha$, the proximal operator $\textsf{prox}_{h,\alpha}:\mathbb{R}^{d}\to\mathbb{R}$ is defined as

[TABLE]

Observe that if $h=0$, then $\textsf{prox}_{h,\alpha}({\bf w})={\bf w}$ for every ${\bf w}\in\mathbb{R}^{d}$, and PGD reduces to the classical gradient descent (GD). This encompasses several important optimization problems related to learning, for which $prox$ operator has a closed form expression; some of these problems are given below.

• •

Lasso. Here $f_{i}({\bf w})=\frac{1}{2}(\langle{\bf x}_{i},{\bf w}\rangle-y_{i})^{2}$ and $h({\bf w})=\lambda\|{\bf w}\|_{1}$. It turns out that $\textsf{prox}_{h,\alpha}({\bf z})$ for Lasso is equal to the soft-thresholding operator $S_{\lambda\alpha}({\bf z})$ [Tib15], which, for $j\in[d]$, is defined as

[TABLE]

• •

SVM dual. Jaggi [Jag13] showed an equivalence between the dual formulation of Support Vector Machines (SVM) and Lasso. Hence, SVM dual is also a special case of (1).

• •

Constrained optimization. We want to solve a constrained minimization problem $\min_{{\bf w}\in\mathcal{C}}f({\bf w})$, where $\mathcal{C}\subseteq\mathbb{R}^{d}$ is a closed, convex set. Define an indicator function $I_{\mathcal{C}}$ for $\mathcal{C}$ as follows: $I_{\mathcal{C}}({\bf w}):=0$, if ${\bf w}\in\mathcal{C}$; and $I_{\mathcal{C}}({\bf w}):=\infty$, otherwise. Now, observe the following equivalence

[TABLE]

If we solve the RHS using PGD, then it can be easily verified that the corresponding proximal operator is equal to the projection operator onto the set $\mathcal{C}$ [Tib15]. So, the proximal gradient update step is to compute the usual gradient and then project it back onto the set $\mathcal{C}$.

• •

Logistic regression. Here $f_{i}$ is the logistic function, defined as

[TABLE]

where $u_{i}=\langle{\bf x}_{i},{\bf w}\rangle$, and $h=0$. As noted earlier, since $h=0$, PGD reduces to GD for logistic regression.

• •

Ridge regression. Here $f_{i}({\bf w})=\frac{1}{2}(\langle{\bf x}_{i},{\bf w}\rangle-y_{i})^{2}$ and $h({\bf w})=\frac{\lambda}{2}\|{\bf w}\|_{2}^{2}$. Since $f_{i}$’s and $h$ are differentiable, we can alternatively solve this simply using GD.

Let ${\bf X}\in\mathbb{R}^{n\times d}$ denote the data matrix, whose $i$’th row is equal to the $i$’th data point ${\bf x}_{i}$. For simplicity, assume that $m$ divides $n$, and let ${\bf X}_{i}$ denote the $\frac{n}{m}\times d$ matrix, whose $j$’th row is equal to ${\bf x}_{(i-1)\frac{n}{m}+j}$. In a distributed setup, all the data is distributed among $m$ worker machines (worker $i$ has ${\bf X}_{i}$) and master updates the parameter vector using the update rule (2). At the $t$’th iteration, master sends ${\bf w}_{t}$ to all the workers; worker $i$ computes the gradient (denoted by $\nabla_{i}f({\bf w}_{t})$) on its local data and sends it to the master; master aggregates all the received $m$ local gradients to obtain the global gradient

[TABLE]

Now, master updates the parameter vector according to (2) and obtains ${\bf w}_{t+1}$. Repeat the process until convergence.

If full gradients are too costly to compute.

Updating the parameter vector in each iteration of PGD according to (2) requires computing full gradients. This may be prohibitive in large-scale applications, where each machine in a distributed framework has a lot of data, and computing full gradients at local machines may be too expensive and becomes the bottleneck. In such scenarios, there are two alternatives to reduce this per-iteration cost: (i) Coordinate Descent (CD), in which we pick a few coordinates (at random), compute the partial gradient along those, and descent along those coordinates only, and (ii) Stochastic Gradient Descent (SGD), in which we sample a data point at random, compute the gradient on that point, and descent along that direction. These are discussed in Section 2.2 and Section 6.1, respectively.

2.2 Coordinate Descent

For the clear exposition of ideas, we focus on the non-regularized empirical risk minimization from (1) (i.e., taking $h=0$) for learning a generalized linear model (GLM). This can be generalized to objectives with (non-)differentiable regularizers [BKBG11, ST11]. Let ${\bf X}\in\mathbb{R}^{n\times d}$ denote the data matrix and ${\bf y}\in\mathbb{R}^{n}$ the corresponding label vector. To make it distinct from the last section, we denote the objective function by $\phi$ and write it as $\phi({\bf X}{\bf w};{\bf y})$ to emphasize that we want to learn a GLM, where the objective function depends on the data points only through their inner products with the parameter vector. Formally, we want to optimize333Here we are not optimizing the average of loss functions – since $n$ is a fixed number, this does not affect the solution space.

[TABLE]

For $\mathcal{U}\subseteq[d]$, we write $\nabla_{\mathcal{U}}\phi({\bf X}{\bf w};{\bf y})$ to denote the gradient of $\phi({\bf X}{\bf w};{\bf y})$ with respect to ${\bf w}_{\mathcal{U}}$, where ${\bf w}_{\mathcal{U}}$ denotes the $|\mathcal{U}|$-length vector obtained by restricting ${\bf w}$ to the coordinates in $\mathcal{U}$. To make the notation less cluttered, let $\phi^{\prime}({\bf X}{\bf w};{\bf y})$ denote the $n$-length vector, whose $i$’th entry is equal to $\ell^{\prime}(\langle{\bf x}_{i},{\bf w}\rangle;y_{i}):=\frac{\partial}{\partial u}\ell(u;y_{i})|_{u=\langle{\bf x}_{i},{\bf w}\rangle}$. Note that $\nabla\phi({\bf X}{\bf w};{\bf y})={\bf X}^{T}\phi^{\prime}({\bf X}{\bf w};{\bf y})$ and that $\nabla_{\mathcal{U}}\phi({\bf X}{\bf w};{\bf y})={\bf X}_{\mathcal{U}}^{T}\phi^{\prime}({\bf X}{\bf w};{\bf y})$, where ${\bf X}_{\mathcal{U}}$ denotes the $n\times|\mathcal{U}|$ matrix obtained by restricting the column indices of ${\bf X}$ to the elements in $\mathcal{U}$.

Coordinate descent (CD) is an iterative algorithm, where, in each iteration, we choose a set of coordinates and update only those coordinates (while keeping the other coordinates fixed). In distributed CD, we take advantage of the parallel architecture to improve the running time of (centralized) CD. In the distributed setting, we divide the data matrix vertically into $m$ parts and store the $i$’th part at the $i$’th worker node. Concretely, assume, for simplicity, that $m$ divides $d$. Let ${\bf X}=[{\bf X}_{1}\ {\bf X}_{2}\ \ldots\ {\bf X}_{m}]$ and ${\bf w}=[{\bf w}_{1}^{T}\ {\bf w}_{2}^{T}\ \ldots\ {\bf w}_{m}^{T}]^{T}$, where each ${\bf X}_{i}$ is an $n\times\frac{d}{m}$ matrix and each ${\bf w}_{i}$ is a length $\frac{d}{m}$ vector. Each worker $i$ stores ${\bf X}_{i}$ and is responsible for updating (a few coordinates of) ${\bf w}_{i}$ – hence the terminology, model-parallelism. We store the label vector ${\bf y}$ at the master node. In coordinate descent, since we update only a few coordinates in each round, there are a few options on how to update these coordinates in a distributed manner:

Subset of workers:

Master picks a subset $\mathcal{S}\subset[m]$ of workers and asks them to update their ${\bf w}_{i}$’s [RT16]. This may not be good in the adversarial setting, because if only a small subset of workers are updating their parameters, the adversary can corrupt those workers and disrupt the computation.

Subset of coordinates for all workers:

All the worker nodes update only a subset of the coordinates of their local parameter vector ${\bf w}_{i}$’s. Master can (deterministically or randomly) pick a subset $\mathcal{U}$ (which may or may not be different for all workers) of $f\leq d/m$ coordinates and asks each worker to updates only those coordinates. If master picks $\mathcal{U}$ deterministically, it can cycle through and update all coordinates of the parameter vector in $\lceil d/mf\rceil$ iterations.

In Algorithm 1, we give the distributed CD algorithm with the second approach, where all worker nodes update the coordinates of their local parameter vectors for a single subset $\mathcal{U}$. We will adopt this approach in our method to make the distributed CD Byzantine-resilient. Let $r=\frac{d}{m}$. For any $i\in[m]$, let ${\bf w}_{i}=[w_{i1}\ w_{i2}\ldots w_{ir}]^{T}$ and ${\bf X}_{i}=[{\bf X}_{i1}\ {\bf X}_{i2}\ldots{\bf X}_{ir}]$, where ${\bf X}_{ij}$ is the $j$’th column of ${\bf X}_{i}$. For any $i\in[m]$ and $\mathcal{U}\subseteq[r]$, let ${\bf w}_{i\mathcal{U}}$ denote the $|\mathcal{U}|$-length vector that is obtained from ${\bf w}_{i}$ by restricting its entries to the coordinates in $\mathcal{U}$; similarly, let ${\bf X}_{i\mathcal{U}}$ denote the $n\times|\mathcal{U}|$ matrix obtained by restricting the column indices of ${\bf X}_{i}$ to the elements in $\mathcal{U}$.

In Algorithm 1, for each worker $i$ to update ${\bf w}_{i}$ according to (6), where the partial gradient of $\phi$ with respect to ${\bf w}_{i\mathcal{U}}$ is equal to $\nabla_{i\mathcal{U}}\phi({\bf X}{\bf w};{\bf y})$ = ${\bf X}_{i\mathcal{U}}^{T}\phi^{\prime}(\sum_{j=1}^{m}{\bf X}_{j}{\bf w}_{j};{\bf y})$ and worker $i$ has only $({\bf X}_{i},{\bf w}_{i})$, every other worker $j$ sends ${\bf X}_{j}{\bf w}_{j}$ to the master, who computes $\phi^{\prime}(\sum_{j=1}^{m}{\bf X}_{j}{\bf w}_{j};{\bf y})$555Note that even after computing ${\bf X}{\bf w}$, master needs access to the labels $y_{i},i=1,2,\ldots,n$ to compute $\phi^{\prime}({\bf X}{\bf w};{\bf y})$. Since ${\bf y}\in\mathbb{R}^{n}$ is just a vector, we can either store that at master, or, alternatively, we can encode ${\bf y}$ distributedly at the workers and master can recover that using the method developed in Section 4 for Byzantine-resilient distributed matrix-vector multiplication, where the matrix is an identity matrix and vector is equal to ${\bf y}$. and sends it back to all the workers. Observe that, even if one worker is corrupt, it can send an adversarially chosen vector to make the computation at the master deviate arbitrarily from the desired computation, which may adversely affect the update at all the worker nodes subsequently.666Specifically, suppose the $i$’th worker is corrupt and the adversary wants master to compute $\phi^{\prime}({\bf X}{\bf w}+{\bf e};{\bf y})$ for any arbitrary vector ${\bf e}\in\mathbb{R}^{n}$ of its choice, then the $i$’th worker can send ${\bf X}_{i}{\bf w}_{i}+{\bf e}$ to the master. Similarly, corrupt workers can send adversarially chosen information to affect the stopping criterion.

2.3 Adversary Model

We want to perform the distributed computation described in Section 2.1 and Section 2.2 under adversarial attacks, where the corrupt nodes may provide erroneous vectors to the master node. Our adversarial model is described next.

In our adversarial model, the adversary can corrupt at most $t<\frac{m}{2}$ worker nodes777Our results also apply to a slightly different adversarial model, where the adversary can adaptively choose which of the $t$ worker nodes to attack at each iteration. However, in this model, the adversary cannot modify the local stored data of the attacked node, as otherwise, over time, it can corrupt all the data, making any defense impossible., and the compromised nodes may collaborate and arbitrarily deviate from their pre-specified programs. If a worker is corrupt, then instead of sending the true vector, it may send an arbitrary vector to disrupt the computation. We refer to the corrupt nodes as erroneous or under the Byzantine attack. We can also handle asynchronous updates, by dropping the straggling nodes beyond a specified delay, and still compute the correct gradient due to encoding. Therefore we treat updates from these nodes as being “erased”. We refer to these as erasures/stragglers. For every worker $i$ that sends a message to the master, we can assume, without loss of generality, that the master receives ${\bf u}_{i}+{\bf e}_{i}$, where ${\bf u}_{i}$ is the true vector and ${\bf e}_{i}$ is the error vector, where ${\bf e}_{i}={\bf 0}$ if the $i$’th node is honest, otherwise can be arbitrary. We assume that at most $t$ nodes can be adversarially corrupt and at most $s$ nodes can be stragglers, where $s$ and $t$ are some constants less than $\frac{1}{2}$ that we will decide later. Note that the master node does not know which $t$ worker nodes are corrupted (which makes this problem non-trivial to solve), but knows $t$. We propose a method that mitigates the effects of both of these anomalies.

Remark 1**.**

A well-studied problem is that of asynchronous distributed optimization, where the workers can have different delays in updates [DB13]. One mechanism to deal with this is to wait for a subset of responses, before proceeding to the next iteration, treating the others as missing (or erasures) [KSDY17]. Byzantine attacks are quite distinct from such erasures, as the adversary can report wrong local gradients, requiring the master node to create mechanisms to overcome such attacks. If the master node simply aggregates the collected updates as in (4), the computed gradient could be arbitrarily far away from the true one, even with a single adversary [MGR18].

2.4 Our Approach to Gradient Computation

Recall that $f_{i}({\bf w})=\ell(\langle{\bf x}_{i},{\bf w}\rangle;y_{i})$ for some differentiable loss function $\ell$, and the gradient of $f_{i}$ at ${\bf w}$ is equal to $\nabla f_{i}({\bf w})=({\bf x}_{i})^{T}\ell^{\prime}(\langle{\bf x}_{i},{\bf w}\rangle;y_{i})$, where $\ell^{\prime}(\langle{\bf x}_{i},{\bf w}\rangle;y_{i}):=\frac{\partial}{\partial u}\ell(u;y_{i})|_{u=\langle{\bf x}_{i},{\bf w}\rangle}$. Note that $\nabla f_{i}({\bf w})\in\mathbb{R}^{d}$ is a column vector. Let $f^{\prime}({\bf w})$ denote the $n$-length vector whose $i$’th entry is equal to $\ell^{\prime}(\langle{\bf x}_{i},{\bf w}\rangle;y_{i})$. With this notation, since $f({\bf w})=\frac{1}{n}\sum_{i=1}^{n}f_{i}({\bf w})$, we have $\nabla f({\bf w})=\frac{1}{n}{\bf X}^{T}f^{\prime}({\bf w})$. Since $n$ is a constant, it is enough to compute ${\bf X}^{T}f^{\prime}({\bf w})$. So, for simplicity, in the rest of the paper we write

[TABLE]

A natural approach to computing the gradient $\nabla f({\bf w})$ is to compute it in two rounds: (i) compute $f^{\prime}({\bf w})$ in the 1st round by first multiplying ${\bf X}$ with ${\bf w}$ and then master locally computes $f^{\prime}({\bf w})$ from ${\bf X}{\bf w}$ (master can do this locally, because ${\bf X}{\bf w}$ is an $n$-dimensional vector whose $i$’th entry is equal to $\langle{\bf x}_{i},{\bf w}\rangle$ and $(f^{\prime}({\bf w}))_{i}=\ell^{\prime}(\langle{\bf x}_{i},{\bf w}\rangle;y_{i})$);888Note that even after computing ${\bf X}{\bf w}$, master needs access to the labels $y_{i},i=1,2,\ldots,n$ to compute $f^{\prime}({\bf w})$. See Footnote 5 for a discussion on how master can get access to the labels. and then (ii) compute $\nabla f({\bf w})={\bf X}^{T}f^{\prime}({\bf w})$ in the 2nd round by multiplying ${\bf X}^{T}$ with $f^{\prime}({\bf w})$. So, the task of each gradient computation reduces to two matrix-vector (MV) multiplications, where the matrices are fixed and vectors may be different each time. To combat against the adversarial worker nodes, we do both of these MV multiplications using data encoding and real-error correction; see Figure 1 on page 1 for a pictorial description of our approach.

A two-round approach for gradient computation has been proposed for straggler mitigation in [LLP*+*18], but our method for MV multiplication differs from that fundamentally, as we have to provide adversarial protection. Note that in the case of stragglers/erasures we know who the straggling nodes are, but this information is not known in the case of adversarial nodes, and master needs to decode without this information in the context of Byzantine adversaries. This is slightly different from the standard error correcting codes (over finite fields) as the matrix entries in machine learning applications are from reals. In this case, we use ideas from real-error correction (or sparse reconstruction) from the compressive sensing literature [CT05], and using which we develop an efficient decoding at master, which also gives rise to our sparse encoding matrix; see Section 4 for more details. For decoding efficiently, we crucially leverage the block error pattern and design a decoding method at master, which, interestingly, requires just one application of the sparse recovery method on a vector of size $m$, the number of workers, which may be much smaller than the data dimensions $n$ and $d$, thereby making the decoding computationally efficient. Our encoding matrix (given in (11), designed for MV multiplication) is very sparse and has a regular pattern of non-zero entries, which also makes it applicable for making coordinate-descent (CD) Byzantine-resilient. We emphasize that a general-purpose code for MV multiplication may not be applicable for CD, as each CD iteration requires updating only a few coordinates of the parameter vector, which makes it fundamentally different (and arguably more complicated to robustify) than GD iterations; see Section 3.2 and Section 5 for more details. Since iterative algorithms (such as GD and CD) require repeated parameter updates, it is crucial to have a method that has low computational complexity, both at the worker nodes as well as at the master node, and our coding solutions for both GD and CD achieve that, in addition to being highly storage efficient; see Theorem 1 for GD and Theorem 2 for CD.

Coming back to our two-round approach for gradient computations using MV multiplications, for the 1st round, we encode ${\bf X}$ using a sparse encoding matrix ${\bf S}^{(1)}=[({\bf S}_{1}^{(1)})^{T},\ldots,({\bf S}_{m}^{(1)})^{T}]^{T}$ and store ${\bf S}_{i}^{(1)}{\bf X}$ at the $i$’th worker node; and for the 2nd round, we encode ${\bf X}^{T}$ using another sparse encoding matrix ${\bf S}^{(2)}=[({\bf S}_{1}^{(2)})^{T},\ldots,({\bf S}_{m}^{(2)})^{T}]^{T}$, and store ${\bf S}_{i}^{(2)}{\bf X}^{T}$ at the $i$’th worker node. Now, in the 1st round of the gradient computation at ${\bf w}$, the master node broadcasts ${\bf w}$ and the $i$’th worker node replies with ${\bf S}_{i}^{(1)}{\bf X}{\bf w}$ (a corrupt worker may report an arbitrary vector); upon receiving all the vectors, the master node applies error-correction procedure to recover ${\bf X}{\bf w}$ and then locally computes $f^{\prime}({\bf w})$ as described above. In the 2nd round, the master node broadcasts $f^{\prime}({\bf w})$ and similarly can recover ${\bf X}^{T}f^{\prime}({\bf w})$ (which is equal to the gradient) at the end of the 2nd round. So, it suffices to devise a method for multiplying a vector ${\bf v}$ to a fixed matrix ${\bf A}$ in a distributed and adversarial setting. Since this is a linear operation, we can apply error correcting codes over real numbers to perform this task. We describe it briefly below.

A trivial approach.

Take a generator matrix ${\bf G}$ of any real-error correcting linear code. Encode ${\bf A}$ as ${\bf A}^{T}{\bf G}=:{\bf B}$. Divide the columns of ${\bf B}$ into $m$ groups as ${\bf B}=[{\bf B}_{1}\ {\bf B}_{2}\ldots{\bf B}_{m}]$, where worker $i$ stores ${\bf B}_{i}$. Master broadcasts ${\bf v}$ and each worker $i$ responds with ${\bf v}_{T}{\bf B}_{i}+{\bf e}_{i}^{T}$, where ${\bf e}_{i}={\bf 0}$ if the $i$’th worker is honest, otherwise can be arbitrary. Note that at most $t$ of the ${\bf e}_{i}$’s can be non-zero. Responses from the workers can be combined as ${\bf v}^{T}{\bf B}+{\bf e}^{T}$. Since every row of ${\bf B}$ is a codeword, ${\bf v}^{T}{\bf B}={\bf v}^{T}{\bf A}^{T}{\bf G}$ is also a codeword. Therefore, one can take any off-the-shelf decoding algorithm for the code whose generator matrix is ${\bf G}$ and obtain ${\bf v}^{T}{\bf A}^{T}$. For example, we can use the Reed-Solomon codes (over real numbers) for this purpose, which only incurs a constant storage overhead and tolerates optimal number of corruptions (up to $<\frac{1}{2}$). Note that we need fast decoding, as it is performed in every iteration of the gradient computation by the master. As far as we know, any off-the-shelf decoding algorithm “over real numbers” requires at least a quadratic computational complexity, which leads to $\Omega(n^{2}+d^{2})$ decoding complexity per gradient computation, which could be impractical.

The trivial scheme does not exploit the block error pattern which we crucially exploit in our coding scheme to give a $\sim O((n+d)m)$ time decoding per gradient computation, which could be a significant improvement over the trivial scheme, since $m$ typically is much smaller than $n$ and $d$ for large-scale problems. In fact, our coding scheme enables a trade-off (in terms of storage and computation/communication overhead at the master and the worker nodes) with Byzantine adversary tolerance, without compromising the efficiency at the master node. We also want encoding to be efficient (otherwise it defeats the purpose of data encoding) and our sparse encoding matrix achieves that. Our main result for the Byzantine-resilient distributed gradient computation is as follows, which is proved in Section 4:

Theorem 1** (Gradient Computation).**

Let ${\bf X}\in\mathbb{R}^{n\times d}$ denote the data matrix. Let $m$ denote the total number of worker nodes. We can compute the gradient exactly in a distributed manner in the presence of $t$ corrupt worker nodes and $s$ stragglers, with the following guarantees, where $\epsilon>0$ is a free parameter.

• •

$(s+t)\leq\left\lfloor\frac{\epsilon}{1+\epsilon}\cdot\frac{m}{2}\right\rfloor$.

• •

Total storage requirement is roughly $2(1+\epsilon)|{\bf X}|$.

• •

Computational complexity for each gradient computation:

• –

at each worker node is $O((1+\epsilon)\frac{nd}{m})$.

• –

at the master node is $O((1+\epsilon)(n+d)m)$.

• •

Communication complexity for each gradient computation:

• –

each worker sends $\left((1+\epsilon)\frac{n+d}{m}\right)$ real numbers.

• –

master broadcasts $(n+d)$ real numbers.

• •

Total encoding time is $O\left(nd\left(\frac{\epsilon}{1+\epsilon}m+1\right)\right)$.

Remark 2**.**

The statement of Theorem 1 allows for any $s$ and $t$ as long as $(s+t)\leq\left\lfloor\frac{\epsilon}{1+\epsilon}\cdot\frac{m}{2}\right\rfloor$. As we are handling both erasures and errors in the same way999When there are only stragglers, one can design an encoding scheme where both the master and the worker nodes operate oblivious to encoding, while solving a slightly altered optimization problem [KSDY17], in which gradients are computed approximately, leading to more efficient straggler-tolerant GD. the corruption threshold does not have to handle $s$ and $t$ separately. To simplify the discussion, for the rest of the paper, we consider only Byzantine corruption, and denote the corrupted set by $\mathcal{I}\subset[m]$ with $|\mathcal{I}|\leq t$, with the understanding that this can also work with stragglers.

In Theorem 1, $\epsilon$ is a design choice and a free parameter that can take any value in the interval $[0,m-1]$, where $\epsilon=0$ implies no corruption and $\epsilon=m-1$ implies that corruption threshold $t$ can be anything up to $\frac{m-1}{2}$. If we want to tolerate $t$ corrupt workers, then $\epsilon$ must satisfy $\epsilon\geq\frac{2t}{m-2t}$.101010We could have written everything in terms of $t,m,n,d$, but we chose to introduce another variable $\epsilon$ which, in our opinion, clearly brings out the tradeoff between the corruption threshold and the resource requirements without cluttering the expressions.

Remark 3** (Comparison with the plain distributed PGD).**

We compare the resource requirements of our method with the plain distributed PGD (which provides no adversarial protection), where all the data points are evenly distributed among the $m$ workers. In each iteration, master sends the parameter vector ${\bf w}$ to all the workers; upon receiving ${\bf w}$, all workers compute the gradients on their local data in $O(\frac{nd}{m})$ time (per worker) and send them to the master; master aggregates them in $O(md)$ time to obtain the global gradient and then updates the parameter vector using (2).

In our scheme (i) the total storage requirement is $O(1+\epsilon)$ factor more;111111For example, by taking $\epsilon=2$, our method can tolerate $m/3$ corrupt worker nodes. So, we can tolerate linear corruption with a constant overhead in the resource requirement, compared to the plain distributed gradient computation which does not provide any adversarial protection. (see also Remark 4) (ii) the amount of computation at each worker node is $O(1+\epsilon)$ factor more; (iii) the amount of computation at the master node is $O((1+\epsilon)(1+\frac{n}{d}))$ factor more, which is comparable in cases where $n$ is not much bigger than $d$; (iv) master broadcasts $(1+\frac{n}{d})$ factor more data, which is comparable if $n$ is not much bigger than $d$; and (v) each worker sends $O\left((1+\epsilon)\frac{1+\nicefrac{{n}}{{d}}}{m}\right)$ factor more data, which is $O(1+\epsilon)$ – a constant factor – as long as $n=O(dm)$.

Remark 4**.**

Let $m$ be an even number. Note that we can get the corruption threshold $t$ to be any number less than $m/2$, but at the expense of increased storage and computation. For any $\delta>0$, if we want to get $\delta$ close to m/2, i.e., $t=m/2-\delta$, then we must have $(1+\epsilon)\geq m/2\delta$. In particular, at $\epsilon=2$, we can tolerate up to $m/3$ corrupt nodes, with constant overhead in the total storage as well as on the computational complexity.

Note that when $\delta$ is a constant, i.e., $t$ is close to $\frac{m-1}{2}$, then $\epsilon$ grows linearly with $m$; for example, if $t=\frac{m-1}{2}$, then $\epsilon=m-1$. In this case, our storage redundancy factor is $O(m)$. In contrast, the trivial scheme (see “trivial approach” on page 2.4) does better in this regime and has only a constant storage overhead, but at the expense of an increased decoding complexity at the master, which is at least quadratic in the problem dimensions $d$ and $n$, whereas, our decoding complexity at the master always scales linearly with $d$ and $n$. If we always want a constant storage redundancy for all values of the corruption threshold $t$, we can use our coding scheme if $t\leq c\cdot\frac{m-1}{2}$, where $c<1$ is a constant, and use the trivial scheme if $t$ is close to $\frac{m-1}{2}$.

Our encoding is also efficient and requires $O\left(nd\left(\frac{\epsilon}{1+\epsilon}m+1\right)\right)$ time. Note that $O(nd)$ is equal to the time required for distributing the data matrix ${\bf X}$ among $m$ workers (for running the distributed gradient descent algorithms without the adversary); and the encoding time in our scheme (which results in an encoded matrix that provides Byzantine-resiliency) is a factor of $(2t+1)$ more.

Remark 5**.**

Our scheme is not only efficient (both in terms of computational complexity and storage requirement), but it can also tolerate up to $\lfloor\frac{m-1}{2}\rfloor$ corrupt worker nodes (by taking $\epsilon=m-1$ in Theorem 1). It is not hard to prove that this bound is information-theoretically optimal, i.e., no algorithm can tolerate $\lceil\frac{m}{2}\rceil$ corrupt worker nodes, and at the same time correctly computes the gradient.

2.5 Our Approach to Coordinate Descent

We use data encoding and add redundancy to enlarge the parameter space. Specifically, we encode the data matrix ${\bf X}$ using an encoding matrix ${\bf R}=[{\bf R}_{1}\ {\bf R}_{2}\ \ldots\ {\bf R}_{m}]$, where each ${\bf R}_{i}$ is a $d\times p$ matrix (with $pm\geq d$), and store ${\bf X}{\bf R}_{i}$ at the $i$’th worker. Define $\widetilde{\bf X}^{R}:={\bf X}{\bf R}$. Now, instead of solving (5), we solve the encoded problem $\arg\min_{{\bf v}\in\mathbb{R}^{pm}}\phi(\widetilde{\bf X}^{R}{\bf v};{\bf y})$ using Algorithm 1 (together with decoding at the master); see Figure 2 on page 2 for a pictorial description of our algorithm. We design the encoding matrix ${\bf R}$ such that at every iteration of our algorithm, updating any (small) subset of coordinates of ${\bf v}_{i}$’s (let ${\bf v}=[{\bf v}_{1}^{T}\ {\bf v}_{2}^{T}\ \ldots\ {\bf v}_{m}^{T}]$) automatically updates some (small) subset of coordinates of ${\bf w}$; and, furthermore, by updating those coordinates of ${\bf v}_{i}$’s, we can efficiently recover the correspondingly updated coordinates of ${\bf w}$, despite the errors injected by the adversary. In fact, at any iteration $t$, the encoded parameter vector ${\bf v}_{t}$ and the original parameter vector ${\bf w}_{t}$ satisfies ${\bf v}_{t}={\bf R}^{+}{\bf w}_{t}$, where ${\bf R}^{+}:={\bf R}^{T}({\bf R}{\bf R}^{T})^{-1}$ is the Moore-Penrose pseudo-inverse of ${\bf R}$, and ${\bf w}_{t}$ evolves in the same way as if we are running Algorithm 1 on the original problem.

We will be effectively updating the coordinates of the parameter vector ${\bf w}$ in chunks of size $(m-2t)$ or its integer multiples (where $t$ is the number of corrupt workers). In particular, if each worker $i$ updates $k$ coordinates of ${\bf v}_{i}$, then $k(m-2t)$ coordinates of ${\bf w}$ will get updated. For comparison, Algorithm 1 updates $km$ coordinates of the parameter vector ${\bf w}$ in each iteration, if each worker updates $k$ coordinates in that iteration.

As described in Algorithm 1 for the Byzantine-free CD, in order to update its local parameter vector ${\bf w}_{i}$ according to (6), worker $i$ needs access to $\phi^{\prime}({\bf X}{\bf w};{\bf y})$, which master computes after receiving $\{{\bf X}_{j}{\bf w}_{j}\}_{j\in[m]}$ from the workers. In our Byzantine-resilient algorithm for CD also master will need to compute ${\bf X}{\bf w}$ in every CD iteration, and for this purpose, we employ the same encoding-decoding procedure for MV multiplication that we used in the first round of gradient computation, as described in Section 2.4. In particular, to make the notation distinct from gradient computation, in order to compute ${\bf X}{\bf w}$, we encode ${\bf X}$ using an encoding matrix ${\bf L}=[{\bf L}_{1}^{T}\ {\bf L}_{2}^{T}\ \ldots\ {\bf L}_{m}^{T}]^{T}$, where each ${\bf L}_{i}$ is a $p^{\prime}\times n$ matrix (with $p^{\prime}m\geq n$) and worker $i$ stores $\widetilde{\bf X}_{i}^{L}={\bf L}_{i}{\bf X}$.

Note that in order to compute ${\bf X}{\bf w}$, in the first round of gradient computation as described in Section 2.4, master broadcasts ${\bf w}$ to all the workers and each worker $i$ computes $\widetilde{\bf X}_{i}^{L}{\bf w}$ and sends it the the master (corrupt workers may report arbitrary vectors), who then decodes and obtains ${\bf X}{\bf w}$. However, in coordinate descent, though master wants to compute ${\bf X}{\bf w}$ in each CD iteration, we can significantly improve the computation required at each worker: since only a few coordinates of the original parameter vector ${\bf w}$ are updated in each CD iteration, master needs to send only those updated coordinates, and workers need to preform MV multiplication with a much smaller matrix, whose number of columns is equal to the number of updated coordinates of ${\bf w}$ that they receive from master. Thus, the computational complexity in each CD iteration at worker is proportional to the number of coordinates updated in each CD iteration, as desired.

Our main result for the Byzantine-resilient distributed coordinate descent is stated below, which is proved in Section 5.

Theorem 2** (Coordinate Descent).**

Under the setting of Theorem 1, our Byzantine-resilient distributed CD algorithm has the following guarantees, where $\epsilon>0$ is a free parameter.

• •

$(s+t)\leq\left\lfloor\frac{\epsilon}{1+\epsilon}\cdot\frac{m}{2}\right\rfloor$.

• •

Total storage requirement is roughly $2(1+\epsilon)|{\bf X}|$.

• •

If each worker $i$ updates $\tau$ coordinates of ${\bf v}_{i}$, then

• –

$\frac{\tau m}{1+\epsilon}$* coordinates of the corresponding ${\bf w}$ gets updated.*

• –

the computational complexity in each iteration

at each worker node is $O(n\tau)$.

• *

at the master node is $O((1+\epsilon)nm+\tau m^{2})$.

• –

the communication complexity in each iteration

each worker sends $\left(\tau+(1+\epsilon)\frac{n}{m}\right)$ real numbers.

• *

master broadcasts $\left(\frac{\tau m}{1+\epsilon}+n\right)$ real numbers.

• •

Total encoding time is $O\left(nd\left(\frac{\epsilon}{1+\epsilon}m+1\right)\right)$.

Remark 6** (Comparison with the plain distributed CD).**

We compare the resource requirements of our method with the plain distributed CD described in Algorithm 1 that does not provide any adversarial protection. Let $\epsilon$ be any number in the interval $[0,m-1]$ – for illustration, we can take $\epsilon=2$, which means $t\leq\frac{m}{3}$ workers are corrupt. In Algorithm 1, if each worker $i$ updates $\frac{\tau}{1+\epsilon}$ coordinates of ${\bf w}_{i}$ (in total $\frac{\tau m}{1+\epsilon}$ coordinates of ${\bf w}$) in each iteration, then (i) each worker requires $O(\frac{n\tau}{1+\epsilon})$ time to multiply ${\bf X}_{i}$ with the updated part of ${\bf w}_{i}$; (ii) master requires $O(nm)$ time to compute $\sum_{i=1}^{m}{\bf X}_{i}{\bf w}_{i}$ from $\{{\bf X}_{i}{\bf w}_{i}\}_{i\in[m]}$; (iii) each worker sends $n$ real numbers (required for ${\bf X}_{i}{\bf w}_{i}$) to master; and (iv) master broadcasts $n$ real numbers (required for $\phi^{\prime}({\bf X}{\bf w};{\bf y})$).

In our scheme (i) the total storage requirement is $O(1+\epsilon)$ factor more; (ii) the amount of computation at each worker node is $O(1+\epsilon)$ factor more; (iii) the amount of computation at the master node is $O((1+\epsilon)+\frac{\tau m}{n})$ factor more – typically, since $\tau$ is a constant and number of workers is much less than $n$, this again could be $O(1+\epsilon)$; (iv) master broadcasts $\left(1+\frac{\tau m}{(1+\epsilon)n}\right)$ factor more data, which could be a constant if $\tau m$ is smaller than $(1+\epsilon)n$; and (v) each worker sends $\left(\frac{\tau}{n}+\frac{(1+\epsilon)}{nm}\right)$ factor more data, where the 1st term is much smaller than 1 as $\tau$ is typically a constant, and the 2nd term is close to zero as $(1+\epsilon)$ is always upper-bounded by $m$.

Remark 7** (Comparison with the replication-based strategy).**

One simple way to make Algorithm 1 Byzantine-resilient is using repetition code, where we first divide the set of $m$ workers into $\frac{m}{2t+1}$ groups of size $(2t+1)$ each and also divide the data matrix as ${\bf X}=[{\bf X}_{1}\ {\bf X}_{2}\ \ldots\ {\bf X}_{\frac{m}{2t+1}}]$ (assume, for simplicity, that $(2t+1)$ divides $m$). Now, store the $i$’th block ${\bf X}_{i}$ at the $(2t+1)$ workers in the $i$’th group of workers. Let the parameter vector be divided as ${\bf w}=[{\bf w}_{1}^{T}\ {\bf w}_{2}^{T}\ \ldots\ {\bf w}_{\frac{m}{2t+1}}^{T}]^{T}$. In each CD iteration, the local parameter updates in any ${\bf w}_{i}$ is replicated at $(2t+1)$ different workers in the $i$’th group of workers, and since at most $t$ workers are corrupt, master can do a majority vote for decoding. Note that the total storage and the computation at workers in this scheme grow linearly by a factor of $(2t+1)$, where $t$ is the number of corruption, which could be significant. In contrast, the method that we propose can tolerate linear corruption, say, $t=\frac{m}{3}$, with a constant overhead in storage and computational complexity.

The Remarks 2, 4, 5 are also applicable for Theorem 2.

3 Related Work

There has been a significant recent interest in using coding-theoretic techniques to mitigate the well-known straggler problem [DB13], including gradient coding [TLDK17, RTDT18, CP18, HRSH18], encoding computation [LLP*+*18, DCG16, DCG19], and data encoding [KSDY17, KSDY19]. However, one cannot directly apply the methods for straggler mitigation to the Byzantine attacks case, as we do not know which updates are under attack. Distributed computing with Byzantine adversaries is a richly investigated topic since [LSP82], and has received recent attention in the context of large-scale distributed optimization and learning [BMGS17, CSX17, CWCP18, YCRB18, AAL18, SX19, XKG19, YCRB19, GV19, RWCP19, LXC*+*19, GHYR19, YLR*+*19, DD20b, DD20a, HKJ20]. These can be divided into three categories: (i) One which assume explicit statistical models for data across workers (e.g., data drawn i.i.d. from a probability distribution) and analyze gradient descent [CSX17, YCRB18, SX19, YCRB19, GHYR19]. (ii) Other set of works make no probabilistic assumption on data, and optimize through stochastic methods (e.g., stochastic gradient descent) [BMGS17, AAL18, GV19, XKG19, LXC*+*19, RWCP19, DD20a, DD20b, HKJ20] and also with deterministic methods (e.g., gradient descent) [DD20a, DD20b]. Note that none of these two sets of works do data encoding and work with data as it is, and provide Byzantine resilience by applying some robust aggregation procedures (e.g., geometric median, coordinate-wise median, outlier-filtering, etc.) at the master for aggregating gradients. (iii) Another line of work which is most relevant to ours provide Byzantine resiliency using redundant computations, either by encoding the gradients [CWCP18] or by encoding the data itself [YLR*+*19]. Note that [RWCP19] combines both redundant computations and do a hierarchical robust aggregation and not is directly comparable to ours.

Note that the statistical nature of data/analysis in the first two sets of works leads to a statistical approximation error in the convergence rates, which is also intensified by the inaccuracy of the robust gradient aggregation procedure. One of the main focuses in these works is typically on obtaining faster convergence (where the goal is to match the convergence rate of plain SGD/GD) and as good an approximation error as possible. Note that the approximation error in all these works scales at least as $\Omega(\sqrt{d})$, where $d$ is the dimension of the model parameter vector, which may be significant in high-dimensional settings. Moreover, in all these works, since we are not allowed to pre-process the data (such as, doing data encoding, etc.), we need to make some assumptions on the data, and furthermore, master has to apply a non-trivial decoding for gradient aggregation, which requires significantly more time than what our decoding requires. For example, filtering-based decoding [SX19, DD20a, DD20b], median-based decoding [CSX17, YCRB18], and heuristic approaches [BMGS17], all have a super-linear complexity in $m$ – in fact, the filtering-based method as in [SX19, DD20a, DD20b] (which is the most effective in terms of the approximation error) requires $O(m^{3}d)$ time. In contrast, our decoding has a linear dependence on both $m$ and $d$. Note that, unlike the first two categories, the third line of work (to which ours also belongs) gives deterministic guarantees and work with arbitrary datasets, with no probabilistic assumptions; we elaborate on these and do a detailed comparison with ours below. We skip the comparison with the first two categories, as it would not be a fair comparison because the underlying setting is different – results in the first two categories are based on statistical assumptions on data/algorithm and inaccurate gradient recovery, whereas, results in the third category make no assumption on the data/algorithm and allow exact gradient recovery.

We want to emphasize that all these works use gradient descent (GD) or stochastic gradient descent (SGD) as their optimization algorithm, which is a data-parallelization method; in this paper, additionally, we also use coordinate descent (CD) algorithm for optimization, which is a model-parallelization method and is preferred over GD in some applications; see Section 1 for more details on this. As will be evident from Section 5, making CD secure against Byzantine attacks is arguably more intricate than securing GD.

We divide this section into three categories: first we compare the redundancy-based methods for GD in Section 3.1, and then CD in Section 3.2. Since we use matrix-vector (MV) multiplication as a core subroutine for both GD and CD, we also compare related work on this in Section 3.3.

3.1 Gradient Descent (GD)

In this section, we do a detailed comparison with [CWCP18] and [YLR*+*19], which are the closest related works that also combat Byzantine adversaries using redundant computations.

For the sake of comparison, assume that $t\leq\frac{m-1}{2}$ workers are corrupt. The coding scheme of Chen et al. [CWCP18], which they called Draco, requires repetition of each data point $(2t+1)$ times, storing each copy at different workers. This gives the storage redundancy factor of $(2t+1)$ in Draco, whereas, our coding method requires storage redundancy factor of $2(1+\epsilon)=\frac{2m}{m-2t}$, which is a constant even if $t$ is a constant ($<\frac{1}{2}$) fraction of $m$.121212To highlight the storage redundancy gain of our method over that of Draco, consider the following two concrete scenarios, where the data matrix ${\bf X}\in\mathbb{R}^{n\times d}$ consists of $nd$ real numbers: (i) In a large setup with $m=1000$ worker nodes, if we want resiliency against $t=100$ corrupt nodes (1/10 nodes are corrupt), our method requires redundancy of 2.5, whereas Draco requires redundancy of 201 (i.e., we need to store only $2.5\times nd$ real numbers, whereas Draco stores $201\times nd$ real numbers), a multiplicative-factor of $>$ 80 more than ours. (ii) In a moderate setup with $m=150$ and $t=50$ (1/3 nodes are corrupt), the redundancy of our method is 6, whereas Draco requires redundancy of 101, a multiplicative-factor of $\approx 17$ more than ours. Since each worker in Draco is doing $(2t+1)$-factor more computation for each GD iteration (than simply computing the gradients as in plain distributed GD), the computational cost at workers also grows by the same factor, which is a significant downside of their scheme. In contrast, our scheme only requires $O(\frac{m}{m-2t})$ more computation at worker, which is a constant even if $t$ is a constant ($<\frac{1}{2}$) fraction of $m$. This significantly reduces the computation time at the worker nodes in our scheme compared to Draco, without sacrificing much on the computation time required by the master node – the decoding at master in Draco takes $O(md)$ time, whereas, our scheme requires $O(\frac{m}{m-2t}(n+d)m)$ time, which is a factor of $O(\frac{m}{m-2t}(1+\frac{n}{d}))$ more than Draco. In high-dimensional settings, where $n$ is not much bigger than $d$, and $t$ is a constant ($<\frac{1}{2}$) fraction of $m$, this overhead is constant. Overall, for a constant fraction of corruption, say, $t=\frac{m}{3}$, Draco requires $\Omega(t)$ times more storage and computation at workers than our scheme (which could be significant in large-scale settings), and requires $\Omega(1+\frac{n}{d})$ times less computation at master. Note that the computation time at workers scales at least as $\Omega(\frac{nd}{m})$, which dominates the time taken by master (since $n,d$ are typically much larger than $m$), so our scheme will be faster than Draco with respect to the overall running time. Note that the coding in Draco is restricted to data replication redundancy, as they encode the gradient as done in [TLDK17], enabling application to (non)-convex problems; in contrast, we encode the data enabling significantly smaller redundancy, and apply it to learn generalized linear models, and is also applicable to MV multiplication.

Yu et al. [YLR*+*19] (which is a concurrent work131313Yu et al. [YLR*+*19] is concurrent to our conference versions in Allerton 2018 [DSD18] and ISIT 2019 [DSD19, DD19], on which this paper is based.) proposes Lagrange coded computing in a distributed framework to compute any multivariate polynomial of the input data and simultaneously provides resilience against stragglers, security against adversaries, and privacy of the dataset against collusion of workers. They leverage the Lagrange polynomial to create computation redundancy among workers, and using standard Reed-Solomon decoding, they can tolerate both erasures/stragglers and errors/adversaries. Their method provide privacy by adding random elements from the field (which in the case of gradient computation is the field of all matrices of a certain dimension) while doing the polynomial interpolation. This is a standard method in Shamir secret sharing scheme [Sha79] that is widely used in information-theoretically secure MPC protocols [CDN15] to provide privacy of users’ data. For the sake of comparison of the resource requirements of our scheme and the one in [YLR*+*19], consider the task of linear regression (the concrete machine learning application studied in [YLR*+*19]). In the following, we assume that $\frac{m-1}{2}-\delta$ workers are corrupt, which corresponds to $\epsilon=\frac{m}{1+2\delta}-1$ in our setting; here $\delta$ can take any value in $[0:\frac{m-1}{2}]$. (i) The storage overhead of our scheme is $\frac{m}{\delta+\nicefrac{{1}}{{2}}}$, whereas, in [YLR*+*19], it is $\frac{m}{\delta+1}$, which is roughly the same as ours. For example, to tolerate $\frac{m}{3}$ corrupt workers (i.e., $\delta=\frac{m-3}{6}$), the storage overhead of our scheme and of [YLR*+*19] is a multiplicative factor of $6$ and $\frac{6}{1+\nicefrac{{3}}{{m}}}\approx 6$, respectively. (ii) The encoding time complexity of our scheme is $O(nd(m-2\delta))$, whereas, it is $O(m\log^{2}(m)\frac{nd}{\delta+1})$ in [YLR*+*19]. Note that for constant $\delta$ (i.e., corruption close to $1/2$), the encoding time of our scheme is much less (by a factor of $O(m\log^{2}(m))$) than that of [YLR*+*19], whereas, for corruption $cm$, where $c<\frac{1}{2}$, the scheme of [YLR*+*19] takes $O(\frac{m}{\log^{2}(m)})$-factor less time in encoding than ours. (iii) The computation time at each worker per gradient computation in both our scheme and [YLR*+*19] is roughly the same – ours requires $O(\frac{nd}{1+2\delta})$ time and [YLR*+*19] requires $O(\frac{nd}{1+\delta})$ time. (iv) The decoding time complexity per gradient computation in [YLR*+*19] is $O(m\log^{2}(m)d)$, whereas, ours requires $O((1+\epsilon)(n+d)m)$ time. Note that when $n$ is not much bigger than $d$ and we want a constant fraction of corruption, say, $\frac{m}{3}$ corruption, then their decoding complexity is worse than ours by a logarithmic factor. Also note that our decoding algorithm is arguably simpler than theirs. (v) For per gradient computation, each worker respectively sends $\frac{n+d}{1+2\delta}$ and $d$ real numbers in ours and the scheme in [YLR*+*19]. Note that if $n\leq dm$ and to tolerate a constant fraction of corruption, say, $\frac{m}{3}$ corruption, each worker sends roughly $O(m)$ less data in our scheme than that of [YLR*+*19]. Overall, if we want tolerance against $\frac{m}{3}$ corrupt worker nodes, then both our scheme and the one in [YLR*+*19] have similar resource requirements, except for that our scheme has a much better communication complexity (by a factor of $O(m)$) from workers to the master, whereas, the encoding time complexity (which is a one-time process) of [YLR*+*19] is better than ours by a factor of $O(\frac{m}{\log^{2}(m)})$.

3.2 Coordinate Descent (CD)

Even for the straggler problem, we are only aware of one work by Karakus et al. [KSDY19] that, in addition to distributed GD, also studies distributed CD, and that for quadratic problems (e.g., linear/ridge regression) only. It also does data encoding and achieves low redundancy and low complexity, by allowing convergence to an approximate rather than exact solution. As far as we know, ours is the first work that studies distributed CD under Byzantine attacks and provides an efficient solution, much better than the replication-based solution (see Remark 7). At the heart of our solution for CD is the matrix-vector (MV) multiplication procedure that we develop in this paper; and it is the specific regular structure of our encoding matrix (given in (11), designed for the MV multiplication) that allows for partially updating the coordinates of the parameter vector in each CD iteration. Note that a general-purpose encoding matrix for MV multiplication may not be applicable for the CD algorithm.

It has been observed earlier in several works (see, for example, [LLP*+*18, DCG16]) that gradient computation in GD for linear regression can be reduced to MV multiplication, and any general-purpose code for MV multiplication can be used to provide a solution for gradient computation. As far as we know, ours is the first paper that makes the connection of CD and MV multiplication, and provides an efficient solution for CD (which is also resilient to Byzantine attacks) for learning generalized linear models. Note that, unlike GD, not any general-purpose code for MV multiplication can be used for CD: the main challenge in CD comes from the fact that we only update a small number of coordinates of the parameter vector in each CD iteration; when we encode the data and iteratively update some coordinates of the (encoded) parameter vector using the encoded data, we need to make sure that this update in the encoded parameter vector is reconciled with the update in the original parameter vector. This is fundamentally different from GD iterations. See Section 5 for more details.

3.3 Matrix-Vector Multiplication

For the task of a more fundamental problem of matrix-vector (MV) multiplication in the presence of Byzantine adversaries, which is at the core of the optimization algorithms in this paper, we are only aware of two concurrent works [YLR*+*19] (see Footnote 13) and [DCG19]141414The conference version [DCG16] only studies the straggler problem, and the journal version [DCG19] briefly mentions how their results from [DCG16] can be extended to handle adversarial nodes, and we describe that in this section. that provide (coding-theoretic) solutions to this problem. In the following, we do a detailed comparison of our solution with both of these works and also discuss the (dis)similarities.

We have already done a detailed comparison with Yu et al. [YLR*+*19] (concurrent work, see Footnote 13) with respect to gradient descent in Section 3.1. For the problem of MV multiplication, the storage requirement, computation time per worker, and communication complexity to/from workers is the same in both ours and [YLR*+*19]. The comparison of encoding time complexity is same as above; however, for a constant corruption, say, $\frac{m}{3}$ corrupt workers, our method outperforms the one in [YLR*+*19] in terms of the decoding time complexity by a factor of $O(\log^{2}(m))$. Note that, unlike [YLR*+*19], we make a fundamental connection of handling Byzantine errors with the sparse reconstruction (or the real-error correction) problem from the compressive sensing literature [CT05].

Dutta et al. [DCG19] (concurrent work, see Footnote 14) focuses on matrix-vector (MV) multiplication. Though their main focus is on providing resilience against stragglers, they also mention that handling stragglers is very different than handling errors, as it requires to correct errors over real numbers, and, unlike stragglers, we do not know which workers are corrupt. Similar to our observation, they also note that since the matrices and vectors have entries from real numbers, the decoding problem reduces to the sparse reconstruction problem from the compressive sensing literature [CT05] and they also provide such a reduction. Apart from these similarities, our solution for MV multiplication differs from that of [DCG19] in several important ways: (i) [DCG19] provides a detailed solution to the distributed MV multiplication for the straggler problem for the case when the number of rows in the matrix is smaller than the number of workers nodes. As mentioned in [DCG19], this method can be easily generalized to the more general case when the matrix is of arbitrary dimension, in which case, first we can divide the rows of the matrix into several sub-matrices, each having number of rows smaller than the number of workers, and then apply the above method independently to each sub-matrix. This simple extension may work (without losing efficiency) for the straggler/erasure problem, however, leads to a highly inefficient solution for the adversary/error problem. The reason being that, in the presence of Byzantine workers, if we solve the sparse reconstruction problem for each sub-matrix separately, this would be inefficient, as the decoding would then be computationally expensive. To remedy this, we exploit the block error pattern and use a simple idea of linearly combining the response vectors from each worker using coefficients drawn from an absolutely continuous distribution, so that we only need to do just one computation for solving the sparse construction problem. This significantly reduces the decoding complexity; see Section 4.1 for details. (ii) [DCG19] only shows a connection to the sparse recovery problem, whereas, we provide a complete solution, with a concrete sparse recovery (or real-error correction) matrix and resource (encoding/decoding time, storage, communication) requirement analysis. (iii) Our encoding matrix (given in (11)) to encode data matrices of arbitrary dimensions is very sparse and highly structured which allows us to apply that construction to CD algorithm, which, as far we know, has not been connected with MV multiplication before. Also, ours is the first paper that provides a non-trivial and efficient (data encoding) solution to CD in the presence of a Byzantine adversary. (iv) We also want to mention that the focus in [DCG19] is on making the encoded matrix sparse (at the expense of increased computation at workers) so that workers need to compute shorter dot products, whereas, in this paper, we make the encoding matrix sparse (much sparser than the encoded matrix of [DCG19]) to get efficient encoding/decoding.

4 Our Solution to Gradient Computation

In this section, we describe the core technical part of our two-round approach for gradient computation described in Section 2.4 – a method for performing matrix-vector (MV) multiplication in a distributed manner in the presence of a malicious adversary who can corrupt at most $t$ of the $m$ worker nodes. Here, the matrix is fixed and we want to right-multiply a vector with this matrix.

Given a fixed matrix ${\bf A}\in\mathbb{R}^{{n_{r}}\times{n_{c}}}$ and a vector ${\bf v}\in\mathbb{R}^{{n_{c}}}$, we want to compute ${\bf A}{\bf v}$ in a distributed manner in the presence of at most $t$ corrupt worker nodes; see Section 2.3 for details on our adversary model. Our method is based on data encoding and error correction over real numbers, where the matrix ${\bf A}$ is encoded and distributed among all the worker nodes, and the master node recovers the MV product ${\bf A}{\bf v}$ using real-error correction; see Figure 1. We will think of our encoding matrix as ${\bf S}=[{\bf S}_{1}^{T}\ {\bf S}_{2}^{T},\ldots,{\bf S}_{m}^{T}]$, where each ${\bf S}_{i}$ is a $p\times{n_{r}}$ matrix and $pm\geq{n_{r}}$. We will derive the matrix ${\bf S}$ in Section 4.2. For the value of $p$, looking ahead, we will set $p=\lceil\frac{n}{m-2t}\rceil$, which is a constant multiple of $\frac{n}{m}$ even if $t$ is a constant ($<\frac{1}{2}$) fraction of $m$ (e.g., if $t=\frac{m}{3}$, we would have $p=\frac{3n}{m}$). For $i\in[m]$, we store the matrix ${\bf S}_{i}{\bf A}$ at the $i$’th worker node. As described in Section 2, the computation proceeds as follows: The master sends ${\bf v}$ to all the worker nodes and receives $\{{\bf S}_{i}{\bf A}{\bf v}+{\bf e}_{i}\}_{i=1}^{m}$ back from them. Let ${\bf e}_{i}=[e_{i1},e_{i2},\ldots,e_{ip}]^{T}$ for every $i\in[p]$. Note that ${\bf e}_{i}={\bf 0}$ if the $i$’th node is honest, otherwise can be arbitrary. In order to find the set of corrupt worker nodes, master equivalently writes $\{{\bf S}_{i}{\bf A}{\bf v}+{\bf e}_{i}\}_{i=1}^{m}$ as $p$ systems of linear equations.

[TABLE]

where, for every $i\in[p]$, $\tilde{\bf e}_{i}=[e_{1i},e_{2i},\ldots,e_{mi}]^{T}$, and $\tilde{\bf S}_{i}$ is an $m\times{n_{r}}$ matrix whose $j$’th row is equal to the $i$’th row of ${\bf S}_{j}$, for every $j\in[m]$. Note that at most $t$ entries in each $\tilde{\bf e}_{i}$ are non-zero. Observe that $\{{\bf S}_{i}{\bf A}{\bf v}+{\bf e}_{i}\}_{i=1}^{m}$ and $\{\tilde{\bf S}_{i}{\bf A}{\bf v}+\tilde{\bf e}_{i}\}_{i=1}^{p}$ are equivalent systems of linear equations, and we can get one from the other.

Note that $\tilde{\bf S}_{i}$’s constitute the encoding matrix ${\bf S}$, which we have to design. In the following, we will design these matrices $\tilde{\bf S}_{i}$’s (which in turn will determine the encoding matrix ${\bf S}$), with the help of another matrix ${\bf F}$, which will be used to find the error locations, i.e., identities of the compromised worker nodes. We will design the matrix ${\bf F}$ (of dimension $k\times m$, where $k<m$ – here $k$ is determined by the error-correction capability, and we will set $k=2t$; see Section 4.4 for more details) and the matrices $\tilde{\bf S}_{i}$’s such that

1. C.1

${\bf F}\tilde{\bf S}_{i}=0$ for every $i\in[p]$. 2. C.2

For any $t$-sparse ${\bf u}\in\mathbb{R}^{m}$, we can efficiently find all the non-zero locations of ${\bf u}$ from ${\bf F}{\bf u}$. 3. C.3

For any $\mathcal{T}\subset[m]$ such that $|\mathcal{T}|\geq(m-t)$, let ${\bf S}_{\mathcal{T}}$ denote the $|\mathcal{T}|p\times{n_{r}}$ matrix obtained from ${\bf S}$ by restricting it to all the ${\bf S}_{i}$’s for which $i\in\mathcal{T}$. We want ${\bf S}_{\mathcal{T}}$ to be of full column rank.

If we can find such matrices, then we can recover the desired MV multiplication ${\bf A}{\bf v}$ exactly: briefly, C.1 and C.2 will allow us to locate the corrupt worker nodes; once we have found them, we can discard all the information that the master node had received from them. This will yield ${\bf S}_{\mathcal{T}}{\bf A}{\bf v}$, where ${\bf S}_{\mathcal{T}}$ is the $|\mathcal{T}|p\times{n_{r}}$ matrix obtained from ${\bf S}$ by restricting it to ${\bf S}_{i}$’s for all $i\in\mathcal{T}$, where $\mathcal{T}$ is the set of all honest worker nodes. Now, by C.3, since ${\bf S}_{\mathcal{T}}$ is of full column rank, we can recover ${\bf A}{\bf v}$ from ${\bf S}_{\mathcal{T}}{\bf A}{\bf v}$ exactly. Details follow.

Suppose we have matrices ${\bf F}$ and $\tilde{\bf S}_{i}$’s such that C.1 holds. Now, multiplying (8) by ${\bf F}$ yields

[TABLE]

for every $i\in[p]$, where $\|\tilde{\bf e}_{i}\|_{0}\leq t$. In Section 4.1, we give our approach for finding all the corrupt worker nodes with the help of any error locator matrix ${\bf F}$. Then, in Section 4.2, we give a generic construction for designing $\tilde{\bf S}_{i}$’s (and, in turn, our encoding matrix ${\bf S}$) such that C.1 and C.3 hold. In Section 4.3, we show how to compute the desired matrix-vector product ${\bf A}{\bf v}$ efficiently, once we have discarded all the data from the corrupt works nodes. Then, in Section 4.4, we will give details of the error locator matrix ${\bf F}$ that we use in our construction.

Remark 8**.**

As we will see in Section 4.2, the structure of our encoding matrix ${\bf S}$ is independent of our error locator matrix ${\bf F}$. Specifically, the repetitive structure of the non-zero entries of ${\bf S}$ as well as their locations will not change irrespective of what the ${\bf F}$ matrix is. This makes our construction very generic, as we can choose whichever ${\bf F}$ suits our needs the best (in terms of how many erroneous indices it can locate and with what decoding complexity), and it won’t affect the structure of our encoding matrix at all – only the non-zero entries might change, neither their repetitive format, nor their locations!

4.1 Finding The Corrupt Worker Nodes

Observe that $\textsf{supp}(\tilde{\bf e}_{i})$ may not be the same for all $i\in[p]$, but we know, for sure, that the non-zero locations in all these error vectors occur within the same set of $t$ locations. Let $\mathcal{I}=\bigcup_{i=1}^{p}\textsf{supp}(\tilde{\bf e}_{i})$, which is the set of all corrupt worker nodes. Note that $|\mathcal{I}|\leq t$. We want to find this set $\mathcal{I}$ efficiently, and for that we note the following crucial observation. Since the non-zero entries of all the error vectors $\tilde{\bf e}_{i}$’s occur in the same set $\mathcal{I}$, a random linear combination of $\tilde{\bf e}_{i}$’s has support equal to $\mathcal{I}$ with probability one, if the coefficients of the linear combination are chosen from an absolutely continuous probability distribution. This idea has appeared before in [ME08] in the context of compressed sensing for recovering arbitrary sets of jointly sparse signals that have been measured by the same measurement matrix.

Definition 1**.**

A probability distribution is called absolutely continuous, if every event of measure zero occurs with probability zero.

It is well-known that a distribution is absolutely continuous if and only if it can be represented as an integral over an integrable density function [Bil95, Theorem 31.8, Chapter 6]. Since Gaussian and uniform distributions have an explicit integrable density function, both are absolutely continuous. Conversely, discrete distributions are not absolutely continuous. Now we state a lemma from [ME08] that shows that a random linear combination of the error vectors (where coefficients are chosen from an absolutely continuous distribution) preserves the support with probability one.

Lemma 1** ([ME08]).**

Let $\mathcal{I}=\bigcup_{i=1}^{p}\textsf{supp}(\tilde{\bf e}_{i})$, and let $\hat{\bf e}=\sum_{i=1}^{p}\alpha_{i}\tilde{\bf e}_{i}$, where $\alpha_{i}$’s are sampled i.i.d. from an absolutely continuous distribution. Then with probability 1, we have $\textsf{supp}(\hat{\bf e})=\mathcal{I}$.

From (9) we have ${\bf f}_{i}={\bf F}\tilde{\bf e}_{i}$ for every $i\in[p]$. Take a random linear combination of ${\bf f}_{i}$’s with coefficients $\alpha_{i}$’s chosen i.i.d. from an absolutely continuous distribution, for example, the Gaussian distribution. Let $\tilde{\bf f}=\alpha_{i}\left(\sum_{i=1}^{p}{\bf f}_{i}\right)=\alpha_{i}\left(\sum_{i=1}^{p}{\bf F}\tilde{\bf e}_{i}\right)={\bf F}\left(\sum_{i=1}^{p}\alpha_{i}\tilde{\bf e}_{i}\right)={\bf F}\tilde{\bf e}$, where $\tilde{\bf e}=\sum_{i=1}^{p}\alpha_{i}\tilde{\bf e}_{i}$. Note that, with probability 1, $\textsf{supp}(\tilde{\bf e})$ is equal to the set of all corrupt worker nodes, and we want to find this set efficiently. In other words, given ${\bf F}\tilde{\bf e}$, we want to find $\textsf{supp}(\tilde{\bf e})$ efficiently. For this, we need to design a $k\times m$ matrix ${\bf F}$ (where $k<m$) such that for any sparse error vector ${\bf e}\in\mathbb{R}^{m}$, we can efficiently find $\textsf{supp}({\bf e})$ from ${\bf f}={\bf F}{\bf e}$. Many such matrices have been known in the literature that can handle different levels of sparsity with varying decoding complexity. We can choose any of these matrices depending on our need, and this will not affect the design of our encoding matrix ${\bf S}$. In particular, we will use a $k\times m$ Vandermonde matrix along with the Reed-Solomon type decoding, which can correct up to $k/2$ errors and has decoding complexity of $O(m^{2})$; see Section 4.4 for details.

Time required in finding the corrupt worker nodes.

The time taken in finding the corrupt worker nodes is equal to the sum of the time taken in the following 3 tasks. (i) Computing ${\bf F}\tilde{\bf e}_{i}$ for every $i\in[p]$: Note that we can get ${\bf F}\tilde{\bf e}_{i}$ by multiplying (8) with ${\bf F}$. Since ${\bf F}$ is a $k\times m$ matrix, and we compute ${\bf F}\tilde{h}_{i}({\bf v})$ for $p$ systems, this requires $O(pkm)$ time. (ii) Taking a random linear combination of $p$ vectors each of length $m$, which takes $O(pm)$ time. (iii) Applying Lemma 2 (in Section 4.4) once to find the error locations, which takes $O(m^{2})$ time. Since $p$ is much bigger than $m$, the total time complexity is $O(pkm)$.

4.2 Designing The Encoding Matrix ${\bf S}$

Now we give a generic construction for designing $\tilde{\bf S}_{i}$’s such that C.1 and C.3 hold. Fix any $k\times m$ matrix ${\bf F}$ such that we can efficiently find ${\bf e}$ from ${\bf F}{\bf e}$, provided ${\bf e}$ is sufficiently sparse. We can assume, without loss of generality, that ${\bf F}$ has full row-rank; otherwise, there will be redundant observations in ${\bf F}{\bf e}$ that we can discard and make ${\bf F}$ smaller by discarding the redundant rows. Let $\mathcal{N}({\bf F})\subset\mathbb{R}^{m}$ denote the null-space of ${\bf F}$. Since $\textsf{rank}({\bf F})=k$, dimension of $\mathcal{N}({\bf F})$ is $q=(m-k)$. Let $\{{\bf b}_{1},{\bf b}_{2},\ldots,{\bf b}_{q}\}$ be a basis of $\mathcal{N}({\bf F})$, and let ${\bf b}_{i}=[b_{i1}\ b_{i2}\ldots b_{im}]^{T}$, for every $i\in[q]$. We set ${\bf b}_{i}$’s the columns of the following matrix ${\bf F}^{\perp}$:

[TABLE]

The following property of ${\bf F}^{\perp}$ will be used for recovering the MV product in Section 4.3.

Claim 1**.**

For any subset $\mathcal{T}\subset[m]$, such that $|\mathcal{T}|\geq(m-t)$, let ${\bf F}_{\mathcal{T}}^{\perp}$ be the $|\mathcal{T}|\times q$ matrix, which is equal to the restriction of ${\bf F}^{\perp}$ to the rows in $\mathcal{T}$. Then ${\bf F}_{\mathcal{T}}^{\perp}$ is of full column rank.

Proof.

Note that $q=m-k$, where $k=2t$. So, if we show that any $q$ rows of ${\bf F}^{\perp}$ are linearly independent, then, this in turn will imply that for every $\mathcal{T}\subset[m]$ with $|\mathcal{T}|\geq(m-t)$, the sub-matrix ${\bf F}_{\mathcal{T}}^{\perp}$ will have full column rank. In the following we show that any $q$ rows of ${\bf F}^{\perp}$ are linearly independent. To the contrary, suppose not; and let $\mathcal{T}^{\prime}\subset[m]$ with $|\mathcal{T}^{\prime}|=q$ be such that the $q\times q$ matrix ${\bf F}_{\mathcal{T}^{\prime}}^{\perp}$ is not a full rank matrix. This implies that there exists a non-zero ${\bf c}^{\prime}\in\mathbb{R}^{q}$ such that ${\bf F}_{\mathcal{T}^{\prime}}^{\perp}{\bf c}^{\prime}={\bf 0}$. Let ${\bf b}={\bf F}^{\perp}{\bf c}^{\prime}$. Note that ${\bf b}\neq{\bf 0}$ (because columns of ${\bf F}^{\perp}$ are linearly independent) and also that $\|{\bf b}\|_{0}\leq m-q=k$. Now, since ${\bf F}{\bf F}^{\perp}={\bf 0}$, we have ${\bf F}{\bf b}={\bf 0}$, which contradicts the fact that any $k$ columns of ${\bf F}$ are linearly independent. ∎

Now we design $\tilde{\bf S}_{i}$’s. For $i\in[p]$, we set $\tilde{\bf S}_{i}$ as follows:

[TABLE]

where $l=q$ if $i<p$; otherwise $l={n_{r}}-(p-1)q$. The first $(i-1)q$ and the last ${n_{r}}-[(i-1)q+l]$ columns of $\tilde{\bf S}_{i}$ are zero. This also implies that the number of rows in each ${\bf S}_{i}$ is $p=\lceil{n_{r}}/q\rceil$.

Claim 2**.**

For every $i\in[p]$, we have ${\bf F}\tilde{\bf S}_{i}=0$.

Proof.

By construction, the null-space of ${\bf F}$ is $\mathcal{N}({\bf F})=\textsf{span}\{{\bf b}_{1},{\bf b}_{2},\ldots,{\bf b}_{q}\}$, which implies that ${\bf F}{\bf b}_{i}={\bf 0}$, for every $i\in[q]$. Since all the columns of $\tilde{\bf S}_{i}$’s are either ${\bf 0}$ or ${\bf b}_{j}$ for some $j\in[q]$, the claim follows. ∎

The above constructed matrices $\tilde{\bf S}_{i}$’s give the following encoding matrix ${\bf S}_{i}$ for the $i$’th worker node:

[TABLE]

All the unspecified entries of ${\bf S}_{i}$ are zero. The matrix ${\bf S}_{i}$ is for encoding the data for worker $i$. By stacking up the ${\bf S}_{i}$’s on top of each other gives us our desired encoding matrix ${\bf S}$.

To get efficient encoding, we want ${\bf S}$ to be as sparse as possible. Since ${\bf S}$ is completely determined by ${\bf F}^{\perp}$, whose columns are the basis vectors of $\mathcal{N}({\bf F})$, it suffices to find a sparse basis for $\mathcal{N}({\bf F})$. It is known that finding the sparsest basis for the null-space of a matrix is NP-hard [CP86]. Note that we can always find the basis vectors of $\mathcal{N}({\bf F})$ by reducing ${\bf F}$ to its row-reduced-echelon-form (RREF) using the Gaussian elimination [HK71]. This will result in ${\bf F}^{\perp}$ whose last $q$ rows forms a $q\times q$ identity matrix. Note that $q=m-k$, where $k=2t$. So, if the corruption threshold $t$ is very small as compared to $m$, the ${\bf F}^{\perp}$ that we obtain by the RREF will be very sparse – only the first $2t$ rows may be dense. Since computing ${\bf S}$ is equivalent to computing ${\bf F}^{\perp}$, and we can compute ${\bf F}^{\perp}$ in $O(k^{2}m)$ time using the Gaussian elimination, the time complexity of computing ${\bf S}$ is also $O(k^{2}m)$.

Now we prove an important property of the encoding matrix ${\bf S}$ that will be crucial for recovery of the desired matrix-vector product.

Claim 3**.**

For any $\mathcal{T}\subset[m]$ such that $|\mathcal{T}|\geq(m-t)$, let ${\bf S}_{\mathcal{T}}$ denote the $|\mathcal{T}|p\times{n_{r}}$ matrix obtained from ${\bf S}$ by restricting it to all the blocks ${\bf S}_{i}$’s for which $i\in\mathcal{T}$. Then ${\bf S}_{\mathcal{T}}$ is of full column rank.

Proof.

For $i\in[p-1]$, let $\mathcal{B}_{i}=[(i-1)q+1:iq]$ and $\mathcal{B}_{p}=[(p-1)q+1:{n_{r}}-(p-1)q]$, where we see $\mathcal{B}_{i}$’s as a collection of some column indices. Consider any two distinct $i,j\in[p]$. It is clear that for any two vectors ${\bf u}_{1}\in\mathcal{B}_{i},{\bf u}_{2}\in\mathcal{B}_{j}$, we have $\textsf{supp}({\bf u}_{1})\cap\textsf{supp}({\bf u}_{2})=\phi$, which means that all the columns in distinct $\mathcal{B}_{i}$’s are linearly independent. So, to prove the claim, we only need to show that the columns within the same $\mathcal{B}_{i}$’s are linearly independent. Fix any $i\in[p]$, and consider the $|\mathcal{T}|p\times q$ sub-matrix ${\bf S}_{\mathcal{T}}^{(i)}$ of ${\bf S}_{\mathcal{T}}$, which is obtained by restricting ${\bf S}_{\mathcal{T}}$ to the columns in $\mathcal{B}_{i}$. There are precisely $|\mathcal{T}|$ non-zero rows in ${\bf S}_{\mathcal{T}}^{(i)}$, which are equal to the rows of the matrix ${\bf F}_{\mathcal{T}}^{\perp}$ defined in Claim 1. We have already shown in the proof of Claim 1 that ${\bf F}_{\mathcal{T}}^{\perp}$ is of full column rank. Therefore, ${\bf S}_{\mathcal{T}}^{(i)}$ is also of full column rank. This concludes the proof of Claim 3. ∎

Since ${\bf S}_{\mathcal{T}}$ is of full column rank, in principle, we can recover any vector ${\bf u}\in\mathbb{R}^{{n_{r}}}$ from ${\bf S}_{\mathcal{T}}{\bf u}$. In the next section, we show an efficient way for this recovery.

4.3 Recovering The Matrix-Vector Product ${\bf A}{\bf v}$

Once the master has found the set $\mathcal{I}$ of corrupt worker nodes, it discards all the data received from them. Let $\mathcal{T}=[m]\setminus\mathcal{I}=\{i_{1},i_{2},\ldots,i_{f}\}$ be the set of all honest worker nodes, where $f=(m-|\mathcal{I}|)\geq(m-t)$. Let ${\bf r}=[{\bf r}_{1}^{T}{\bf r}_{2}^{T}\ldots{\bf r}_{m}^{T}]$, where ${\bf r}_{i}={\bf S}_{i}{\bf A}{\bf v}+{\bf e}_{i}$. All the ${\bf r}_{i}$’s from the honest worker nodes can be written as

[TABLE]

where ${\bf S}_{\mathcal{T}}$ is as defined in Claim 3, and ${\bf r}_{\mathcal{T}}$ is also defined analogously and equal to the restriction of ${\bf r}$ to all the ${\bf r}_{i}$’s for which $i\in\mathcal{T}$. Since ${\bf S}_{\mathcal{T}}$ has full column rank (by Claim 3), in principle, we can recover ${\bf A}{\bf v}$ from (12). Next we show how to recover ${\bf A}{\bf v}$ efficiently, by exploiting the structure of ${\bf S}$.

Let $\tilde{\bf r}_{j}=[r_{i_{1}j},r_{i_{2}j},\ldots,r_{i_{f}j}]^{T}$, for every $j\in[p]$. The repetitive structure of ${\bf S}_{i}$’s (see (11)) allows us to write (12) equivalently in terms of $p$ smaller systems.

[TABLE]

where, for $j\in[p-1]$, $\mathcal{B}_{i}=[(i-1)q+1:iq]$ and ${\bf F}_{j}={\bf F}_{\mathcal{T}}^{\perp}$, and $\mathcal{B}_{p}=[(p-1)q+1:{n_{r}}-(p-1)q]$ and ${\bf F}_{p}$ is equal to the restriction of ${\bf F}_{\mathcal{T}}^{\perp}$ to its first $({n_{r}}-(p-1)q)$ columns. Since ${\bf F}_{\mathcal{T}}^{\perp}$ has full column rank (by Claim 1), we can compute $({\bf A}{\bf v})_{\mathcal{B}_{i}}$ for all $i\in[p]$, by multiplying (13) by ${\bf F}_{j}^{+}=({\bf F}_{j}^{T}{\bf F}_{j})^{-1}{\bf F}_{j}^{T}$, which it called the Moore-Penrose inverse of ${\bf F}_{j}$. Since ${\bf A}{\bf v}=[({\bf A}{\bf v})_{\mathcal{B}_{1}}^{T},({\bf A}{\bf v})_{\mathcal{B}_{2}}^{T},\ldots,({\bf A}{\bf v})_{\mathcal{B}_{p}}^{T})]^{T}$, we can recover the desired MV product ${\bf A}{\bf v}$.

Time Complexity analysis.

The task of obtaining ${\bf A}{\bf v}$ from ${\bf S}_{\mathcal{T}}{\bf A}{\bf v}$ reduces to (i) computing ${\bf F}_{j}^{+}=({\bf F}_{\mathcal{T}}^{\perp})^{+}$ once, which takes $O(q^{2}|\mathcal{T}|)$ time naïvely; (ii) computing ${\bf F}_{p}^{+}$ once, which takes at most $O(q^{2}|\mathcal{T}|)$ time naïvely; and (iii) computing the MV products ${\bf F}_{j}^{+}\tilde{\bf r}_{j}$ for every $j\in[p]$, which takes $O(pq|\mathcal{T}|)$ time in total. Since $p$ is much bigger than $q$, the total time taken in recovering ${\bf A}{\bf v}$ from ${\bf S}_{\mathcal{T}}{\bf A}{\bf v}$ is $O(pq|\mathcal{T}|)=O(pm^{2})$.

4.4 Designing The Error Locator Matrix ${\bf F}$

In this section, we design a $k\times m$ matrix ${\bf F}$ (where $k<m$) such that for any sparse error vector ${\bf e}\in\mathbb{R}^{m}$, we can uniquely and efficiently recover ${\bf e}$ (and, therefore, $\textsf{supp}({\bf e})$) from the under-determined system of linear equations ${\bf f}={\bf F}{\bf e}\in\mathbb{R}^{k}$. This is related to the sparse representation problem, where one would like to find the sparsest representation of ${\bf f}$ in terms of the linear combination of the columns of ${\bf F}$, i.e., minimizing $\|{\bf e}\|_{0}$ subject to the constraint that ${\bf f}={\bf F}{\bf e}$. This problem is of combinatorial nature and is known to be NP-hard [CT05]. To make this problem computationally tractable, Candes and Tao [CT05] showed that if ${\bf F}$ satisfies a certain regularity condition (which they named the restricted isometry property (RIP)), then the sparsest reconstruction problem can be reduced to minimizing $\|{\bf e}\|_{1}:=\sum_{i=1}^{m}|e_{i}|$ subject to the constraint that ${\bf f}={\bf F}{\bf e}$, which can be efficiently solved using a linear program. They also showed that a random Gaussian matrix satisfies the RIP condition. A common problem with such random constructions is that they may not work with small block-lengths (in our setting, $m$ is the number of workers which may not be a big number), and can only correct a constant fraction of errors, where the constant is very small. We need a deterministic construction that can handle a constant fraction (ideally up to 1/2) of errors and that works with small block-lengths.

Akçakaya and Tarokh [AT08] proposed an efficient solution to the sparse representation problem using Vandermonde matrices. To construct them, take $m$ distinct non-zero elements $z_{1},z_{2},\ldots,z_{m}$ from $\mathbb{R}$, and consider the following $k\times m$ Vandermonde matrix ${\bf F}$.

[TABLE]

For the above ${\bf F}$, it was shown in [AT08] that, if $|\textsf{supp}({\bf e})|\leq k/2$, then the Reed-Solomon type decoding can be used for exact reconstruction of ${\bf e}$ from ${\bf f}={\bf F}{\bf e}$.151515Note that, since any $k$ columns of ${\bf F}$ (which is the Vandermonde matrix) are linearly independent, if there exists a vector ${\bf e}$ such that $|\textsf{supp}({\bf e})|\leq k/2$ and ${\bf e}$ satisfies ${\bf f}={\bf F}{\bf e}$ for a fixed ${\bf f}$, then ${\bf e}$ is unique. Furthermore, their decoding algorithm is efficient and runs in $O(m^{2})$ time. The results in [AT08] are given for complex vector spaces, and they hold over real numbers also. Below we state the sparse recovery result (specialized to reals) from [AT08].

Lemma 2** ([AT08]).**

Let ${\bf F}$ be the $k\times m$ matrix as defined in (14). Let ${\bf e}\in\mathbb{R}^{m}$ be an arbitrary vector with $|\textsf{supp}({\bf e})|\leq k/2$. We can exactly recover the vector ${\bf e}$ from ${\bf f}={\bf F}{\bf e}$ in $O(m^{2})$ time.

Note that ${\bf F}$ is a $k\times m$ matrix, where $k<m$. Choosing $k$ is in our hands, and larger the $k$, more the number of errors we can correct (but at the expense of increased storage and computation); see Section 4.5 for more details.

4.5 Resource Requirement Analysis

In this section, we analyze the total amount of resources (storage, computation, and communication) required by our method for computing gradients in the presence of $t$ (out of $m$) adversarial worker nodes and prove Theorem 1. Fix an $\epsilon>0$. Let the corruption threshold $t$ satisfy $t\leq\lfloor(\epsilon/(1+\epsilon))\cdot(m/2)\rfloor$.

As described earlier in Section 2.4, we compute the gradient $\nabla f({\bf w})={\bf X}^{T}f^{\prime}({\bf w})$ in two-rounds; and in each round we use the Byzantine-tolerant MV multiplication, which we have developed in Section 4, as a subroutine; see Figure 1 for a pictorial representation of our scheme. We encode ${\bf X}$ to compute $f^{\prime}({\bf w})$ in the 1st round: first compute ${\bf X}{\bf w}$ using MV multiplication and then locally compute $f^{\prime}({\bf w})$. To compute ${\bf X}^{T}f^{\prime}({\bf w})$ (which is equal to the gradient) in the 2nd round, we encode ${\bf X}^{T}$ and compute ${\bf X}^{T}f^{\prime}({\bf w})$. Let ${\bf S}^{(1)}$ and ${\bf S}^{(2)}$ be the encoding matrices of dimensions $p_{1}m\times n$ and $p_{2}m\times d$, respectively, to encode ${\bf X}$ and ${\bf X}^{T}$, respectively. Here, $p_{1}=\lceil n/q\rceil$ and $p_{2}=\lceil d/q\rceil$, where $q=m-k$. Since $k=2t$ (by Lemma 2), we have $q=(m-k)\geq m/(1+\epsilon)$.

4.5.1 Storage Requirement

Each worker node $i$ stores two matrices ${\bf S}_{i}^{(1)}{\bf X}$ and ${\bf S}_{i}^{(2)}{\bf X}^{T}$. The first one is a $p_{1}\times(d+1)$ matrix, and the second one is a $p_{2}\times n$ matrix. So, the total amount of storage at all worker nodes is equal to storing $(p_{1}(d+1)+p_{2}n)\times m$ real numbers. Since $p_{1}\leq\lceil(1+\epsilon)\frac{n}{m}\rceil$ and $p_{2}\leq\lceil(1+\epsilon)\frac{d}{m}\rceil$, the total storage is

[TABLE]

where the first term is roughly equal to a $2(1+\epsilon)$ factor more than the size of ${\bf X}$. Note that the second term does not contribute much to the total storage as compared to the first term, because the number of worker nodes $m$ is much smaller than both $n$ and $d$. In fact, if $m-k$ divides both $n$ and $d$, then the second term vanishes. Since $|{\bf X}|$ is an $n\times d$ matrix, the total storage at each worker node is almost equal to $2(1+\epsilon)\frac{|{\bf X}|}{m}$, which is a constant factor of the optimal, that is, $\frac{|{\bf X}|}{m}$, and the total storage is roughly equal to $2(1+\epsilon)|{\bf X}|$.

4.5.2 Computational Complexity

We can divide the computational complexity of our scheme as follows:

• •

Encoding the data matrix. Since, for every $i\leq k$ and $j>k$, the total number of non-zero entries in ${\bf S}_{i}^{(1)}$ and ${\bf S}_{j}^{(1)}$ are at most $n$ and $p_{1}$, respectively (see Section 4.2 for details), the computational complexity for computing ${\bf S}_{i}^{(1)}{\bf X}$ for each $i\leq k$, and ${\bf S}_{j}^{(1)}{\bf X}$ for each $j>k$, is $O(nd)$ and $O(p_{1}d)$, respectively. So, the encoding time for computing ${\bf S}^{(1)}{\bf X}$ is equal to $O\left(k(nd)+(m-k)(p_{1}d)\right)=O\left((\frac{\epsilon}{1+\epsilon}m+1)nd\right)$. Similarly, we can show that the encoding time for computing ${\bf S}^{(2)}{\bf X}^{T}$ is also equal to $O\left((\frac{\epsilon}{1+\epsilon}m+1)nd\right)$. Note that computing ${\bf S}^{(1)}$ and ${\bf S}^{(2)}$ take $O(k^{2}m)$ time each, which is much smaller, as compared to the encoding time. So, the total encoding time is $O\left((\frac{\epsilon}{1+\epsilon}m+1)nd\right)$. Note that this encoding is to be done only once.

• •

Computation at each worker node. In the first round, upon receiving ${\bf w}$ from the master node, each worker $i$ computes $({\bf S}_{i}^{(1)}{\bf X}){\bf w}$, and reports back the resulting vector. Similarly, in the second round, upon receiving $f^{\prime}({\bf w})$ from the master node, each worker $i$ computes $({\bf S}_{i}^{(2)}{\bf X}^{T})f^{\prime}({\bf w})$, and reports back the resulting vector. Since ${\bf S}_{i}^{(1)}{\bf X}$ and ${\bf S}_{i}^{(2)}{\bf X}^{T}$ are $p_{1}\times(d+1)$ and $p_{2}\times n$ matrices, respectively, each worker node $i$ requires $O(p_{1}d+p_{2}n)=O((1+\epsilon)\frac{nd}{m})$ time.

• •

Computation at the master node. The total time taken by the master node in both the rounds is the sum of the time required in (i) finding the corrupt worker nodes in the 1st and 2nd rounds, which requires $O(p_{1}km)$ and $O(p_{2}km)$ time, respectively (see Section 4.1), (ii) recovering ${\bf X}{\bf w}$ from ${\bf S}_{\mathcal{T}}^{(1)}{\bf X}{\bf w}$ in the 1st round, which requires $O(p_{1}m^{2})$ time, (iii) computing $f^{\prime}({\bf w})$ from ${\bf X}{\bf w}$, which takes $O(n)$ time, and (iv) recovering ${\bf X}^{T}f^{\prime}({\bf w})$ from ${\bf S}_{\mathcal{T}}^{(2)}{\bf X}^{T}f^{\prime}({\bf w})$ in the 2nd round, which requires $O(p_{2}m^{2})$ time (see Section 4.3). Since $k<m$, the total time is equal to $O((p_{1}+p_{2})m^{2})=O((1+\epsilon)(n+d)m)$.

4.5.3 Communication Complexity

In each gradient computation, (i) master broadcasts $(n+d)$ real numbers, $d$ in the first round and $n$ in the second round; and (ii) each worker sends $\left((1+\epsilon)\frac{n+d}{m}\right)$ real numbers to master, $(1+\epsilon)\frac{n}{m}$ in the first round and $(1+\epsilon)\frac{d}{m}$ in the second round.

5 Our Solution to Coordinate Descent

In this section, we give a solution to the distributed coordinate descent (CD) under Byzantine attacks and prove Theorem 2. To make our notation simpler, we remove the dependence on the label vector ${\bf y}$ in the problem expression (5) and rewrite it as follows (this is without loss of generality in the light of Footnote 5 and Algorithm 1):

[TABLE]

We want to optimize (15) using distributed CD, described in Section 2.2. As outlined in Section 2.5, we use data encoding and error correction over real numbers for that. To combat the effect of adversary, we add redundancy to enlarge the parameter space. Let $\widetilde{\bf X}^{R}={\bf X}{\bf R}$, where ${\bf R}=[{\bf R}_{1}\ {\bf R}_{2}\ \ldots\ {\bf R}_{m}]\in\mathbb{R}^{d\times pm}$ with $pm\geq d$, and each ${\bf R}_{i}$ is a $p\times d$ matrix. We will determine the encoding matrix ${\bf R}$ later, after describing what properties we want from it. For the value of $p$, looking ahead, when $t$ is the number of corrupt workers, we will choose $p=\frac{d}{m-2t}$, which is a constant multiple of $\frac{d}{m}$ even if $t$ is a constant fraction ($<\frac{1}{2}$) of $m$ (e.g., for $t=\frac{m}{3}$, we have $p=\frac{3d}{m}$). We consider ${\bf R}$’s which are of full row-rank. Let ${\bf R}^{+}:={\bf R}^{T}({\bf R}{\bf R}^{T})^{-1}$ denote its Moore-Penrose inverse such that ${\bf R}{\bf R}^{+}=I_{d}$, where $I_{d}$ is the $d\times d$ identity matrix. Note that ${\bf R}^{+}$ is of full column-rank. Let ${\bf v}={\bf R}^{+}{\bf w}$ be the transformed vector, which lies in a larger (than $d$) dimensional space. Let ${\bf R}^{+}=[({\bf R}_{1}^{+})^{T}\ ({\bf R}_{2}^{+})^{T}\ \ldots\ ({\bf R}_{m}^{+})^{T}]^{T}$, where each ${\bf R}_{i}^{+}:=({\bf R}^{+})_{i}$ is a $p\times d$ matrix. With this, by letting ${\bf v}=[{\bf v}_{1}^{T}\ {\bf v}_{2}^{T}\ \ldots\ {\bf v}_{m}^{T}]^{T}$, we have that ${\bf v}_{i}={\bf R}_{i}^{+}{\bf w}$ for every $i\in[m]$. Now, consider the following modified problem over the encoded data.

[TABLE]

Observe that, since ${\bf R}$ is of full row-rank, $\min_{{\bf w}\in\mathbb{R}^{d}}\phi({\bf X}{\bf w})$ is equal to $\min_{{\bf v}\in\mathbb{R}^{pm}}\phi(\widetilde{\bf X}^{R}{\bf v})$; and from an optimal solution to one problem we can obtain an optimal solution to the other problem. We design an encoding/decoding scheme such that when we optimize the encoded problem (16) using Algorithm 1, the vector ${\bf v}$ that we get in each iteration is of the form ${\bf v}={\bf R}^{+}{\bf w}$ for some vector ${\bf w}\in\mathbb{R}^{d}$.161616If such a ${\bf w}$ exists, then it is unique. This follows from the fact that ${\bf R}^{+}$ is of full column-rank. In fact, our encoding/decoding will ensure that the ${\bf w}$ for which ${\bf v}={\bf R}^{+}{\bf w}$ would be equal to the original parameter vector in that iteration if we had run Algorithm 1 to solve (15). We need this property because in any CD iteration $t$, we need access to the original parameter vector ${\bf w}^{t}$ (such that ${\bf v}^{t}={\bf R}^{+}{\bf w}^{t}$) to facilitate the local parameter updates of ${\bf v}_{1}^{t},\ldots,{\bf v}_{m}^{t}$ at the workers. See the paragraph after (18) for more details.

Now, instead of solving (15), we solve its encoded form (16) using Algorithm 1 (with decoding at the master), where each worker $i$ stores $\widetilde{\bf X}_{i}^{R}={\bf X}{\bf R}_{i}$ and is responsible for updating (some coordinates of) ${\bf v}_{i}$. In the following, let $\mathcal{U}\subseteq[p]$ be a fixed arbitrary subset of $[p]$. Let ${\bf v}^{0}:={\bf R}^{+}{\bf w}^{0}$ for some ${\bf w}^{0}$ at time $t=0$. Suppose, at the beginning of the $t$’th iteration, we have ${\bf v}^{t}={\bf R}^{+}{\bf w}^{t}$ for some ${\bf w}^{t}$, and each worker $i$ updates ${\bf v}_{i\mathcal{U}}^{t}$ according to

[TABLE]

where $\nabla_{i\mathcal{U}}\phi(\widetilde{\bf X}^{R}{\bf v}^{t})=(\widetilde{\bf X}_{i\mathcal{U}}^{R})^{T}\phi^{\prime}(\widetilde{\bf X}^{R}{\bf v}^{t})$. Recall that each ${\bf R}_{i}$ is a $d\times p$ matrix, and each ${\bf R}_{i}^{+}:=({\bf R}^{+})_{i}$ is a $p\times d$ matrix. We denote by ${\bf R}_{i\mathcal{U}}$ the $d\times|\mathcal{U}|$ matrix obtained by restricting the columns of ${\bf R}_{i}$ to the elements of $\mathcal{U}$. Analogously, we denote by ${\bf R}_{i\mathcal{U}}^{+}:=({\bf R}^{+})_{i\mathcal{U}}$ the $|\mathcal{U}|\times d$ matrix obtained by restricting the rows of ${\bf R}_{i}^{+}$ to the elements of $\mathcal{U}$. With this, we can write $\widetilde{\bf X}_{i\mathcal{U}}^{R}={\bf X}{\bf R}_{i\mathcal{U}}$. Now, (17) can be equivalently written as

[TABLE]

In order to update ${\bf v}_{i\mathcal{U}}^{t}$, worker $i$ requires $\phi^{\prime}(\widetilde{\bf X}^{R}{\bf v}^{t})$, where $\widetilde{\bf X}^{R}{\bf v}^{t}=\sum_{j=1}^{m}\widetilde{\bf X}_{j}^{R}{\bf v}_{j}^{t}$ and worker $i$ has only $(\widetilde{\bf X}_{i}^{R},{\bf v}_{i}^{t})$. Since ${\bf v}^{t}={\bf R}^{+}{\bf w}^{t}$, we have $\widetilde{\bf X}^{R}{\bf v}^{t}={\bf X}{\bf R}{\bf v}^{t}={\bf X}{\bf w}^{t}$. So, it suffices to compute ${\bf X}{\bf w}^{t}$ at the master node – once master has ${\bf X}{\bf w}^{t}$, it can locally compute $\phi^{\prime}({\bf X}{\bf w}^{t})$ and send it to all the workers. Computing ${\bf X}{\bf w}^{t}$ is the distributed matrix-vector (MV) multiplication problem, where the matrix ${\bf X}$ is fixed and we want to compute ${\bf X}{\bf w}^{t}$ for any vector ${\bf w}^{t}$ in the presence of an adversary. In Section 4, we give a method for performing distributed MV multiplication in the presence of an adversary. Now we give an overview, together-with an improvement on its computational complexity.

We encode ${\bf X}$ using an encoding matrix ${\bf L}\in\mathbb{R}^{(p^{\prime}m)\times n}$. Let ${\bf L}=[{\bf L}_{1}^{T}\ {\bf L}_{2}^{T}\ \ldots\ {\bf L}_{m}^{T}]^{T}$, where each ${\bf L}_{i}$ is a $p^{\prime}\times n$ matrix with $p^{\prime}=\lceil\frac{n}{m-2t}\rceil$. Each ${\bf L}_{i}$ has $p^{\prime}$ rows and $n$ columns, and has the same structure as that of ${\bf S}_{i}$ from (11). Worker $i$ stores $\widetilde{\bf X}_{i}^{L}={\bf L}_{i}{\bf X}$. To compute ${\bf X}{\bf w}$, master sends ${\bf w}$ to all the workers; worker $i$ responds with ${\bf L}_{i}{\bf X}{\bf w}+{\bf e}_{i}$, where ${\bf e}_{i}={\bf 0}$ if the $i$’th worker is honest, otherwise can be arbitrary; upon receiving $\{{\bf L}_{i}{\bf X}{\bf w}+{\bf e}_{i}\}_{i=1}^{m}$, where at most $t$ of the ${\bf e}_{i}$’s can be non-zero, master applies the decoding procedure and recovers ${\bf X}{\bf w}$ back. We can improve the computational complexity of this method significantly by observing that, in each iteration of our distributed CD algorithm, only a few coordinates of ${\bf w}$ get updated and the rest of the coordinates remain unchanged. (Looking ahead, when each worker updates ${\bf v}_{i\mathcal{U}}$’s according to (17), it automatically updates ${\bf w}_{f(\mathcal{U})}$ according to (6) – for a specific function $f$ as defined in (21) – where ${\bf v}$ and ${\bf w}$ satisfy ${\bf v}={\bf R}^{+}{\bf w}$.) This implies that for computing ${\bf X}{\bf w}$, master only needs to send the updated coordinates to the workers and keeps the result from the previous MV product with itself. This significantly reduces the local computation at the worker nodes, as now they only need to perform a local MV product of a matrix of size $p^{\prime}\times|f(\mathcal{U})|$ and a vector of length $|f(\mathcal{U})|$. See Section 4 for details.

Our goal in each iteration of CD is to update some coordinates of the original parameter vector ${\bf w}$; instead, by solving the encoded problem, we are updating coordinates of the transformed vector ${\bf v}$. We would like to design an algorithm/encoding such that it has exactly the same convergence properties as if we are running the distributed CD on the original problem without any adversary. For this, naturally, we would like our algorithm to satisfy the following property:

Update on any (small) subset of coordinates of ${\bf w}$ should be achieved by updating some (small) subset of coordinates of ${\bf v}_{i}$’s; and, by updating those coordinates of ${\bf v}_{i}$’s, we should be able to efficiently recover the correspondingly updated coordinates of ${\bf w}$. Furthermore, this should be doable despite the errors injected by the adversary in every iteration of the algorithm.

Note that if each coordinate of ${\bf v}$ depends on too many coordinates of ${\bf w}$, then updating a few coordinates of ${\bf v}$ may affect many coordinates of ${\bf w}$, and it becomes information-theoretically impossible to satisfy the above property (even without the presence of an adversary).171717To see this, consider the case when each worker $i$ updates only the first coordinate of ${\bf v}_{i}$ and no worker is corrupt. Master receives $m$ linear equations ${\bf v}_{i1}={\bf R}^{+}_{i1}{\bf w},\ i=1,2,\ldots,m$, where ${\bf R}^{+}_{i1}$ is the first row of ${\bf R}^{+}_{i}$ for every $i\in[m]$. Assume, for simplicity, that these $m$ equations are linearly independent. When $m$ is smaller than $d$ (which is always the case), there are infinite solutions to this system of linear equations, unless at most $m$ elements of ${\bf w}$ are involved in the $m$ linear equations (i.e., the number of unknowns are at most the number of equations), which is equivalent to saying that the rows ${\bf R}^{+}_{i1}$ for $i=1,2,\ldots,m$ are sparse. Our encoding matrix will satisfy this property; see Section 5.1 for more detail. This imposes a restriction that each row of ${\bf R}^{+}$ must have few non-zero entries, in such a way that updating ${\bf v}_{i\mathcal{U}}^{t}$’s, for any choice of $\mathcal{U}\subseteq[p]$, will collectively update only a subset (which may potentially depend on $\mathcal{U}$) of coordinates of the original parameter vector ${\bf w}^{t}$, and we can uniquely and efficiently recover those updated coordinates of ${\bf w}^{t}$, even from the erroneous vectors $\{{\bf v}_{i\mathcal{U}}^{t+1}+{\bf e}_{i\mathcal{U}}\}_{i=1}^{m}$, where at most $t$ out of $m$ error vectors $\{{\bf e}_{i\mathcal{U}}\}_{i=1}^{m}$ are non-zero and may have arbitrary entries. In order to achieve this, we will design a sparse encoding matrix ${\bf R}^{+}$ (which in turn determines ${\bf R}$), that satisfies the following properties:

1. P.1

${\bf R}^{+}$ has structured sparsity, which induces a map $f:[p]\to\mathcal{P}([d])$ (where $\mathcal{P}([d])$ denotes the power set of $[d]$) such that

1. (a)

$\{f(i):i\in[p]\}$ partitions $\{1,2,\ldots,d\}$, i.e., for every $i,j\in[p]$, such that $i\neq j$, we have $f(i)\cap f(j)=\emptyset$ and that $\bigcup_{i=1}^{p}f(i)=[d]$. 2. (b)

$|f(i)|=|f(j)|$ for every $i,j\in[p-1]$, and $|f(p)|\leq|f(i)|$, for any $i\in[p-1]$. 3. (c)

For any $\mathcal{U}\subseteq[p]$, define $f(\mathcal{U}):=\cup_{j\in\mathcal{U}}f(j)$. If we update ${\bf v}_{i\mathcal{U}}^{t}$, $\forall i\in[m]$, according to (18), it automatically updates ${\bf w}_{f(\mathcal{U})}^{t}$ according to

[TABLE]

If we set ${\bf v}_{i\overline{\mathcal{U}}}^{t+1}:={\bf v}_{i\overline{\mathcal{U}}}^{t}$ and ${\bf w}_{\overline{f(\mathcal{U})}}^{t+1}:={\bf w}_{\overline{f(\mathcal{U})}}^{t}$, then ${\bf v}^{t+1}={\bf R}^{+}{\bf w}^{t+1}$, i.e., our invariant holds.

Note that (19) is the same update rule if we run the plain CD algorithm to update ${\bf w}_{f(\mathcal{U})}$. In fact, our encoding matrix satisfies a stronger property, that ${\bf v}_{i\mathcal{U}}^{t+1}={\bf R}_{i\mathcal{U},f(\mathcal{U})}^{+}{\bf w}_{f(\mathcal{U})}^{t+1}$ holds for every $i\in[m]$, $\mathcal{U}\subseteq[p]$, where ${\bf R}_{i\mathcal{U},f(\mathcal{U})}^{+}$ denotes the $|\mathcal{U}|\times|f(\mathcal{U})|$ matrix obtained from ${\bf R}_{i\mathcal{U}}^{+}$ by restricting its column indices to the elements in $f(\mathcal{U})$.

1. P.2

We can efficiently recover ${\bf w}_{f(\mathcal{U})}^{t+1}$ from the erroneous vectors $\{{\bf v}_{i\mathcal{U}}^{t+1}+{\bf e}_{i\mathcal{U}}\}_{i=1}^{m}$, where at most $t$ of ${\bf e}_{i\mathcal{U}}$’s are non-zero and may have arbitrary entries. Since ${\bf v}_{i\mathcal{U}}^{t+1}={\bf R}_{i\mathcal{U},f(\mathcal{U})}^{+}{\bf w}_{f(\mathcal{U})}^{t+1}$, for every $i\in[m]$, $\mathcal{U}\subseteq[p]$, this property requires that not only ${\bf R}^{+}$, but its sub-matrices also have error correcting capabilities.

Remark 9**.**

Note that P.1 implies that for every $i\in[p]$, we have $|f(i)|\leq d/p$. As we see later, this will be equal to $m/(1+\epsilon)$ for some $\epsilon>0$ which is determined by the corruption threshold. This means that in each iteration of the CD algorithm running on the modified encoded problem, we will be effectively updating the coordinates of the parameter vector ${\bf w}$ in chunks of size $m/(1+\epsilon)$ or its integer multiples. In particular, if each worker $i$ updates $k$ coordinates of ${\bf v}_{i}$, then $km/(1+\epsilon)$ coordinates of ${\bf w}$ will get updated. For comparison, Algorithm 1 updates $km$ coordinates of the parameter vector ${\bf w}$ in each iteration, if each worker updates $k$ coordinates in that iteration.

Now we design an encoding matrix ${\bf R}^{+}$ and a decoding method that satisfy P.1 and P.2.

5.1 Encoding and Decoding

In this section, we first design an encoding matrix ${\bf R}^{+}$ that satisfies P.1. ${\bf R}^{+}$ will be such that it has orthonormal rows, so, ${\bf R}$ is easy to compute, ${\bf R}=({\bf R}^{+})^{T}$. For simplicity, we denote ${\bf R}^{+}$ by ${\bf S}$. We show that the encoding matrix that we design for the MV multiplication in Section 4 satisfies all the properties that we want.181818The encoding and decoding of this section is based on the corresponding algorithms from Section 4. In the MV multiplication, we had a fixed matrix ${\bf A}$ and the master node wants to compute ${\bf A}{\bf w}$ for any vector ${\bf w}$ of its choice. In the solution presented in Section 4, we encode ${\bf A}$ and store ${\bf S}_{i}{\bf A}$ at the $i$’th worker node. Now, the master sends ${\bf w}$ to all the worker nodes, and each worker $i$ responds with ${\bf S}_{i}{\bf A}{\bf w}+{\bf e}_{i}$, where ${\bf e}_{i}={\bf 0}$ if worker $i$ is honest, otherwise can be arbitrary. Once master receives $\{{\bf S}_{i}{\bf A}{\bf w}+{\bf e}_{i}\}_{i=1}^{m}$, it can run the error correcting procedure to recover ${\bf A}{\bf w}$. To apply this in our setting, we take ${\bf A}$ to be the identity matrix, such that ${\bf S}_{i}{\bf A}={\bf S}_{i}$, and the master can recover ${\bf w}$ from $\{{\bf r}_{i}={\bf S}_{i}{\bf w}+{\bf e}_{i}\}_{i=1}^{m}$, if at most $t$ of the ${\bf e}_{i}$’s are non-zero. For convenience, we rewrite the encoding matrix ${\bf S}_{i}$ for the $i$’th worker node from Section 4.2 below:

[TABLE]

Here $q=(m-2t)$ and $l=d-(p-1)q$, where $p=\lceil\frac{d}{q}\rceil$. Note that $1\leq l<q$, and if $q$ divides $d$, then $l=q$. All the unspecified entries of ${\bf S}_{i}$ are zero. By stacking up the ${\bf S}_{i}$’s gives us our desired encoding matrix ${\bf S}=[{\bf S}_{1}^{T}\ {\bf S}_{2}^{T}\ \ldots\ {\bf S}_{m}^{T}]^{T}$. Note that $b_{1i},b_{2i},\ldots,b_{qi}$ are such that if we let ${\bf b}_{i}=[b_{i1}\ b_{i2}\ldots b_{im}]^{T}$ for every $i\in[q]$, then $\{{\bf b}_{1},{\bf b}_{2},\ldots,{\bf b}_{q}\}$ is a set of orthonormal vectors. This implies that ${\bf S}$ is orthonormal, and, therefore, ${\bf S}^{+}={\bf S}^{T}$. By taking ${\bf R}={\bf S}^{T}$, we have ${\bf R}^{+}={\bf S}$. Now we show that ${\bf S}$ satisfies P.1-P.2.

Our Encoding Satisfies P.1.

We need to show a map $f:[p]\to\mathcal{P}([d])$ that satisfies P.1. Let us define the function $f$ as follows, where $(q=m-2t)$ and $p=\lceil\frac{d}{q}\rceil$:

[TABLE]

and for any $\mathcal{U}\subseteq[p]$, we define $f(\mathcal{U}):=\cup_{i\in\mathcal{U}}f(i)$. It is clear from the definition of $f$ that (i) $\{f(i):i\in[p]\}$ partitions $[d]$; (ii) for every $i\in[p-1]$ we have $|f(i)|=q$, and that $|f(p)|\leq q$. Recall that $q=m-2t$. For the 3rd property, note that, for any $\mathcal{U}\subseteq[p]$, all the columns of ${\bf S}_{i\mathcal{U}}$ whose indices belong to $[d]\setminus f(\mathcal{U})$ are identically zero, which implies that we have

[TABLE]

which in turn implies that

[TABLE]

Since ${\bf S}^{+}={\bf S}^{T}$, we have ${\bf S}_{i\mathcal{U}}^{+}={\bf S}_{i\mathcal{U}}^{T}$ for every $i\in[m]$ and every $\mathcal{U}\subseteq[p]$. With these, our update rule ${\bf v}_{i\mathcal{U}}^{t+1}={\bf S}_{i\mathcal{U}}{\bf w}^{t}-\alpha_{t}{\bf S}_{i\mathcal{U}}{\bf X}^{T}\phi^{\prime}({\bf X}{\bf w}^{t})$191919We emphasize that we used ${\bf S}^{+}={\bf S}^{T}$ crucially to equivalently write our update rule ${\bf v}_{i\mathcal{U}}^{t+1}={\bf R}_{i\mathcal{U}}^{+}{\bf w}^{t}-\alpha{\bf R}_{i\mathcal{U}}^{T}{\bf X}^{T}\phi^{\prime}({\bf X}{\bf w}^{t})$ from (18) as ${\bf v}_{i\mathcal{U}}^{t+1}={\bf S}_{i\mathcal{U}}{\bf w}^{t}-\alpha_{t}{\bf S}_{i\mathcal{U}}{\bf X}^{T}\phi^{\prime}({\bf X}{\bf w}^{t})$. This follows because ${\bf S}^{+}={\bf S}^{T}$ and we take ${\bf R}^{+}={\bf S}$, which together imply that ${\bf R}_{i\mathcal{U}}^{+}={\bf R}_{i\mathcal{U}}^{T}={\bf S}_{i\mathcal{U}}$. can equivalently be written as

[TABLE]

where

[TABLE]

Observe that (25) is the same update rule as (19), which implies that if each worker $i$ updates ${\bf v}_{i\mathcal{U}}$ according to the CD update rule, then the collective update at all the worker nodes automatically updates ${\bf w}_{f(\mathcal{U})}$ according the CD update rule. Now we show that our invariant ${\bf v}^{t+1}={\bf S}{\bf w}^{t+1}$ is maintained. We show this by induction. Base case ${\bf v}^{0}={\bf S}{\bf w}^{0}$ holds by construction. For the inductive case, assume that ${\bf v}^{t}={\bf S}{\bf w}^{t}$ holds at time $t$ and we show ${\bf v}^{t+1}={\bf S}{\bf w}^{t+1}$ holds at time $t+1$.

Define $\overline{\mathcal{U}}:=[p]\setminus\mathcal{U}$ and $\overline{f(\mathcal{U})}:=[d]\setminus f(\mathcal{U})$. Since we did not update ${\bf v}_{i\overline{\mathcal{U}}}^{t}$’s, we have ${\bf v}_{i\overline{\mathcal{U}}}^{t+1}={\bf v}_{i\overline{\mathcal{U}}}^{t}$ for every $i\in[m]$. This, together with the inductive hypothesis (i.e., ${\bf v}^{t}={\bf S}{\bf w}^{t}$), implies that

[TABLE]

Since $f(\overline{\mathcal{U}})=\overline{f(\mathcal{U})}$, we have from (22) that

[TABLE]

It is clear from (25) that ${\bf w}_{\overline{f(\mathcal{U})}}^{t}$ did not get an update when we updated ${\bf v}_{i\mathcal{U}}^{t}$’s, which implies that ${\bf w}_{\overline{f(\mathcal{U})}}^{t+1}={\bf w}_{\overline{f(\mathcal{U})}}^{t}$. Substituting this in (27) gives ${\bf S}_{i\overline{\mathcal{U}}}{\bf w}^{t}={\bf S}_{i\overline{\mathcal{U}},\overline{f(\mathcal{U})}}{\bf w}_{\overline{f(\mathcal{U})}}^{t+1}$, which, by (22), yields ${\bf S}_{i\overline{\mathcal{U}}}{\bf w}^{t}={\bf S}_{i\overline{\mathcal{U}}}{\bf w}^{t+1}$. This, together with (26), implies

[TABLE]

We already have from (22) and (24) that

[TABLE]

Since (28) and (29) hold for every $i\in[m]$, we have ${\bf v}^{t+1}={\bf S}{\bf w}^{t+1}$. Hence, the invariant is maintained.

Our Encoding Satisfies P.2.

If we let

[TABLE]

then the collective update (24) from all the workers can be written as

[TABLE]

It is easy to verify that for every choice of $\mathcal{U}\subseteq[p]$, ${\bf S}_{[m]\mathcal{U},f(\mathcal{U})}$ is a full column-rank matrix, which implies that we can in principle recover the updated ${\bf w}_{f(\mathcal{U})}^{t+1}$ from ${\bf v}_{[m]\mathcal{U}}^{t+1}={\bf S}_{[m]\mathcal{U},f(\mathcal{U})}{\bf w}_{f(\mathcal{U})}^{t+1}$. Now we show that not only can we recover ${\bf w}_{f(\mathcal{U})}^{t+1}$ from $\{{\bf S}_{i\mathcal{U},f(\mathcal{U})}{\bf w}_{f(\mathcal{U})}^{t+1}\}_{i=1}^{m}$, but also efficiently recover ${\bf w}_{f(\mathcal{U})}^{t+1}$ from the erroneous vectors $\{{\bf S}_{i\mathcal{U},f(\mathcal{U})}{\bf w}_{f(\mathcal{U})}^{t+1}+{\bf e}_{i\mathcal{U}}\}_{i=1}^{m}$, where at most $t$ out of $m$ error vectors $\{{\bf e}_{i\mathcal{U}}\}_{i=1}^{m}$ are non-zero and may have arbitrary entries. Let $\mathcal{U}=\{j_{1},j_{2},\ldots,j_{|\mathcal{U}|}\}$, and for every $i\in[m]$, let ${\bf e}_{i\mathcal{U}}=[e_{ij_{1}}e_{ij_{2}}\ldots e_{ij_{|\mathcal{U}|}}]^{T}$. Master equivalently writes $\{{\bf S}_{i\mathcal{U},f(\mathcal{U})}{\bf w}_{f(\mathcal{U})}^{t+1}+{\bf e}_{i\mathcal{U}}\}_{i=1}^{m}$ as $|\mathcal{U}|$ systems of linear equations.

[TABLE]

where, for every $i\in\mathcal{U}$, $\tilde{\bf e}_{i}=[e_{1i},e_{2i},\ldots,e_{mi}]^{T}$ and $\tilde{\bf S}_{i,f(\mathcal{U})}$ is an $m\times|f(\mathcal{U})|$ matrix whose $j$’th row is equal to the $i$’th row of ${\bf S}_{j\mathcal{U}}$, for every $j\in[m]$. Note that at most $t$ entries in each $\tilde{\bf e}_{i}$ are non-zero. Observe that $\{{\bf S}_{i\mathcal{U},f(\mathcal{U})}{\bf w}_{f(\mathcal{U})}^{t+1}+{\bf e}_{i\mathcal{U}}\}_{i=1}^{m}$ and $\{\tilde{\bf S}_{i,f(\mathcal{U})}{\bf w}_{f(\mathcal{U})}^{t+1}+\tilde{\bf e}_{i}\}_{i\in\mathcal{U}}$ are equivalent systems of linear equations, and we can get one from the other. Observe that (31) is similar to (8): $\tilde{\bf S}_{i,f(\mathcal{U})}$ is equal to $\tilde{\bf S}_{i}$ (for the same $i$) with some of its zero columns removed; and adding zero columns to $\tilde{\bf S}_{i,f(\mathcal{U})}$ will not change the value of $\tilde{h}_{i}({\bf w}_{f(\mathcal{U})}^{t+1})$. Now, using the machinery developed in Section 4 we can recover ${\bf w}_{f(\mathcal{U})}^{t+1}$ from (31) in $O(|\mathcal{U}|m^{2})$ time.

5.2 Resource Requirement Analysis

In this section, first we give our algorithm developed for distributed coordinate descent in the presence of $t$ (out of $m$) adversarial worker nodes, whose pictorial description is given in Figure 2.

We use two encoding matrices ${\bf L}\in\mathbb{R}^{(p_{1}m)\times n}$ and ${\bf R}\in\mathbb{R}^{d\times(p_{2}m)}$. Let ${\bf L}=[{\bf L}_{1}^{T}\ {\bf L}_{2}^{T}\ \ldots\ {\bf L}_{m}^{T}]^{T}$ and ${\bf R}=[{\bf R}_{1}\ {\bf R}_{2}\ \ldots\ {\bf R}_{m}]$, where each ${\bf L}_{i}$ is a $p_{1}\times n$ matrix with $p_{1}=\lceil\frac{n}{m-2t}\rceil$ and each ${\bf R}_{i}$ is a $d\times p_{2}$ matrix with $p_{2}=\lceil\frac{d}{m-2t}\rceil$. Worker $i$ stores both $\widetilde{\bf X}_{i}^{L}={\bf L}_{i}{\bf X}$ and $\widetilde{\bf X}_{i}^{R}={\bf X}{\bf R}_{i}$. Roughly, ${\bf L}$ is used to recover ${\bf X}{\bf w}$ from the erroneous $\{{\bf L}_{i}{\bf X}{\bf w}+{\bf e}_{i}\}_{i=1}^{m}$, and ${\bf R}$ is used to update the parameter vector reliably despite errors. Here ${\bf L}$ is a full column-rank matrix and ${\bf R}$ is a full row-rank matrix. Initialize with an arbitrary ${\bf w}^{0}$ and let ${\bf v}^{0}={\bf R}^{+}{\bf w}^{0}$. Repeat the following until convergence:

At iteration $t$, master sends $({\bf w}_{f(\mathcal{U})}^{t-1}-{\bf w}_{f(\mathcal{U})}^{t})$202020Observe that master need not send the locations $f(\mathcal{U})$, because workers can compute those by themselves, as they know both $\mathcal{U}$ and the function $f$. to all the workers (at $t=0$, master sends ${\bf w}_{0}$), where $\mathcal{U}\subseteq[p_{2}]$ is the set of indices used for updating ${\bf v}_{i\mathcal{U}}^{t-1}$’s in the previous iteration, which in turn updated ${\bf w}_{f(\mathcal{U})}^{t-1}$; see (24) and (25) in Section 5.1. 2. 2.

Worker $i$ computes $\widetilde{\bf X}_{i}^{L}({\bf w}_{f(\mathcal{U})}^{t-1}-{\bf w}_{f(\mathcal{U})}^{t})={\bf L}_{i}{\bf X}({\bf w}_{f(\mathcal{U})}^{t-1}-{\bf w}_{f(\mathcal{U})}^{t})$ and sends it to the master.212121With some abuse of notation, when we write ${\bf X}{\bf w}_{f(\mathcal{U})}$, we implicitly assume that ${\bf w}_{f(\mathcal{U})}$ is a length $d$ vector, which has [math]’s in the indices that lie in $\overline{f(\mathcal{U})}$. Upon receiving $\{\widetilde{\bf X}_{i}^{L}({\bf w}_{f(\mathcal{U})}^{t-1}-{\bf w}_{f(\mathcal{U})}^{t})+{\bf e}_{i}\}_{i=1}^{m}$, where at most $t$ of the ${\bf e}_{i}$’s are non-zero and may have arbitrary entries, the master applies the decoding procedure of Section 4 and recovers ${\bf X}({\bf w}_{f(\mathcal{U})}^{t-1}-{\bf w}_{f(\mathcal{U})}^{t})$. We assume that master keeps ${\bf X}{\bf w}^{t-1}$ from the previous iteration (which is equal to ${\bf 0}$ if $t=0$), it can compute ${\bf X}{\bf w}^{t}={\bf X}{\bf w}^{t-1}-{\bf X}({\bf w}_{f(\mathcal{U})}^{t-1}-{\bf w}_{f(\mathcal{U})}^{t})$. Note that if $t=0$, this is equal to ${\bf X}{\bf w}^{0}$. 3. 3.

After obtaining ${\bf X}{\bf w}^{t}$, master computes $\phi^{\prime}({\bf X}{\bf w}^{t})$, picks a subset $\mathcal{U}\subseteq[p_{2}]$ (randomly or in a round robin fashion to cover $[p_{2}]$ in a few iterations), and sends $(\phi^{\prime}({\bf X}{\bf w}^{t}),\mathcal{U})$ to all the workers. 4. 4.

Each worker node $i\in[m]$ updates ${\bf v}_{i\mathcal{U}}^{t+1}\leftarrow{\bf v}_{i\mathcal{U}}^{t}-\alpha_{t}\nabla_{i\mathcal{U}}\phi(\widetilde{\bf X}{\bf v}^{t})={\bf v}_{i\mathcal{U}}^{t}-\alpha_{t}(\widetilde{\bf X}_{i\mathcal{U}}^{R})^{T}\phi^{\prime}({\bf X}{\bf w}^{t})$, while keeping the other coordinates of ${\bf v}_{i}^{t}$ unchanged. Worker $i$ sends ${\bf v}_{i\mathcal{U}}^{t+1}$ to the master. Note that ${\bf v}_{i\mathcal{U}}^{t+1}={\bf R}_{i\mathcal{U},f(\mathcal{U})}^{+}{\bf w}_{f(\mathcal{U})}^{t+1}$, where ${\bf w}_{f(\mathcal{U})}^{t+1}=[{\bf w}_{f(\mathcal{U})}^{t}-\alpha{\bf X}_{f(\mathcal{U})}^{T}\phi^{\prime}({\bf X}{\bf w}^{t})]$; see (24) and (25) in Section 5.1. 5. 5.

Upon receiving $\{{\bf v}_{i\mathcal{U}}^{t+1}+{\bf e}_{i\mathcal{U}}\}_{i=1}^{m}$, where at most $t$ of the $\{{\bf e}_{i\mathcal{U}}\}_{i=1}^{m}$’s are non-zero and may have arbitrary entries, master applies the decoding procedure (since our encoding satisfies P.2) and recovers ${\bf w}_{f(\mathcal{U})}^{t+1}$.

Now we analyze the total amount of resources (storage, computation, and communication) required by the above algorithm and prove Theorem 2. Fix an $\epsilon>0$. Let the corruption threshold $t$ satisfy $t\leq\lfloor(\epsilon/(1+\epsilon))\cdot(m/2)\rfloor$.

5.2.1 Storage Requirement:

By a similar analysis done in Section 4.5, we can show that the total storage at all worker nodes is roughly equal to $2(1+\epsilon)|{\bf X}|$.

5.2.2 Computational Complexity:

We can divide the computational complexity of our scheme as follows:

• •

Encoding the data matrix. By a similar analysis done in Section 4.5, we can show that the total encoding time is $O\left((\frac{\epsilon}{1+\epsilon}m+1)nd\right)$. Note that this encoding is to be done only once.

• •

Computation at each worker node. Suppose that in each iteration of our algorithm, all the workers update $\tau$ coordinates of ${\bf v}_{i}$’s. Fix an iteration $t$ and assume that at iteration $(t-1)$, workers updated the coordinates in the set $\mathcal{U}\subseteq[p_{2}]$, where $|\mathcal{U}|=\tau$. Recall from P.1 that updating $\tau=|\mathcal{U}|$ coordinates of each ${\bf v}_{i}^{t-1}$ automatically updates ${\bf w}_{f(\mathcal{U})}^{t-1}$. Upon receiving $({\bf w}_{f(\mathcal{U})}^{t-1}-{\bf w}_{f(\mathcal{U})}^{t})$ from the master node, each worker $i$ computes $\widetilde{\bf X}_{i}^{\bf L}({\bf w}_{f(\mathcal{U})}^{t-1}-{\bf w}_{f(\mathcal{U})}^{t})$, and reports back the resulting vector. Note that $({\bf w}_{f(\mathcal{U})}^{t-1}-{\bf w}_{f(\mathcal{U})}^{t})$ has at most $|f(\mathcal{U})|=\frac{\tau m}{1+\epsilon}$ non-zero elements, which together with that $\widetilde{\bf X}_{i}^{\bf L}$ is a $p_{1}\times d$ matrix, implies that computing $\widetilde{\bf X}_{i}^{\bf L}({\bf w}_{f(\mathcal{U})}^{t-1}-{\bf w}_{f(\mathcal{U})}^{t})$ takes $O(p_{1}\cdot|f(\mathcal{U})|)=O(n\tau)$ time.222222Note that in the very first iteration, master sends ${\bf w}^{0}$, which may be a dense length $d$ vector, and computing $\widetilde{\bf X}_{i}{\bf L}{\bf w}^{0}$ at the $i$’th worker can take $O(p_{1}d)=O((1+\epsilon)\frac{nd}{m})$ time. This is only for the first iteration. In the second round, given $\phi^{\prime}({\bf X}{\bf w}^{t})$, since $(\widetilde{\bf X}_{i\mathcal{U}}^{R})^{T}$ is of dimension $n\times\tau$, updating ${\bf v}_{i\mathcal{U}}^{t}$ requires $O(n\tau)$ time, where $\tau=|\mathcal{U}|$. So, the total time taken by each worker is $O(n\tau)$.

• •

Computation at the master node. Once master receives $\{{\bf L}_{i}{\bf X}({\bf w}_{f(\mathcal{U})}^{t-1}-{\bf w}_{f(\mathcal{U})}^{t})+{\bf e}_{i}\}_{i=1}^{m}$, applying the decoding procedure of Section 4 to obtain ${\bf X}({\bf w}_{f(\mathcal{U})}^{t-1}-{\bf w}_{f(\mathcal{U})}^{t})$ from these erroneous vectors requires $O(p_{1}m^{2})=O((1+\epsilon)nm)$ time. After that obtaining ${\bf X}{\bf w}^{t}$ takes another $O(n)$ time. Given ${\bf X}{\bf w}^{t}$, computing $\phi^{\prime}({\bf X}{\bf w}^{t})$ takes $O(n)$ time, assuming that computing $\ell^{\prime}(\langle{\bf x}_{i},{\bf w}^{t}\rangle;y_{i})$ requires unit time, where $\langle{\bf x}_{i},{\bf w}^{t}\rangle$ is equal to the $i$’th entry of ${\bf X}{\bf w}^{t}$. Upon receiving $\{{\bf v}_{i\mathcal{U}}^{t+1}+{\bf e}_{i\mathcal{U}}\}_{i=1}^{m}$, where ${\bf v}_{i\mathcal{U}}^{t+1}={\bf R}_{i\mathcal{U},f(\mathcal{U})}^{+}{\bf w}_{f(\mathcal{U})}^{t+1}$, for all $i\in[m]$, recovering ${\bf w}_{f(\mathcal{U})}^{t+1}$ requires $O(\tau m^{2})$ time. So, the total time taken by the master node is $O((1+\epsilon)nm+\tau m^{2})$.

5.2.3 Communication Complexity:

Suppose workers update $\tau$ coordinates of ${\bf v}_{i}$’s in each iteration. Then (i) master broadcasts $\left(\frac{\tau m}{1+\epsilon}+n\right)$ real numbers, $\frac{\tau m}{1+\epsilon}$ in the first round to represent ${\bf w}_{f(\mathcal{U})}^{t}$ and $n$ in the second round to represent $\phi^{\prime}({\bf X}{\bf w}^{t})$; and (ii) each worker sends $\left(\tau+(1+\epsilon)\frac{n}{m}\right)$ real numbers, $(1+\epsilon)\frac{n}{m}$ in the first round for computing ${\bf X}{\bf w}^{t}$ at the master node and $\tau$ in the second iteration to represent ${\bf v}_{i\mathcal{U}}^{t}$.

6 Extensions

In this section, we give a few important extensions of our coding scheme developed earlier in Section 4. First we give a Byzantine-resilient and communication-efficient method for stochastic gradient descent (SGD). Second we show how to exploit the specific structure of our encoding matrix to efficiently extend our coding technique to the streaming data model. In the end, we give a few more important applications, where our method can be applied constructively.

6.1 Stochastic Gradient Descent

Stochastic gradient descent (SGD) [HM51] is another alternative if full gradients are too costly to compute. In each iteration of SGD, we sample a data point uniformly at random, compute a gradient on that sample, and update the parameter vector based on that. We start with an arbitrary/random parameter vector ${\bf w}_{0}\in\mathbb{R}^{d}$ and update it according the following update rule:

[TABLE]

where $r_{t}$ is sampled uniformly at random from $\{1,2,\ldots,n\}$. This ensures that the expected value of the gradient is equal to the true gradient. Due to its simplicity and remarkable empirical performance, SGD has become arguably the most widely-used optimization algorithm in many large-scale applications, especially in deep learning [Bot10, RSS12, DCM*+*12]. We want to run SGD in a distributed setup, where data is distributed among $m$ worker nodes and at most $t$ of them can be corrupt; see Section 2.3 for details on our adversary model.

Our solution.

In the plain SGD, we sample a data point randomly and compute its gradient. So, we give a method in which, at any iteration $t$, master picks a random number $r_{t}$ in $\{1,2,\ldots,n\}$, broadcasts it, and recovers the $r_{t}$’th data point ${\bf x}_{r_{t}}$. Once the master has obtained ${\bf x}_{r_{t}}$, it can compute a gradient on it and updates the parameter vector. Since master recovers the data points, we can optimize for non-convex problems also; essentially, we could optimize anything that the plain SGD can. Our method is described below.

We encode ${\bf X}^{T}$ using the $\left\lceil d/(m-2t)\right\rceil\times d$ encoding matrix ${\bf S}^{(2)}$, which has been defined in Section 4.5. For simplicity, we denote ${\bf S}^{(2)}$ by ${\bf S}$. Let ${\bf S}=[{\bf S}_{1}^{T}\ {\bf S}_{2}^{T}\ \ldots\ {\bf S}_{m}^{T}]^{T}$. Note that the $j$’th worker stores ${\bf S}_{j}{\bf X}^{T}$. Let $\widetilde{\bf X}:={\bf S}{\bf X}^{T}$, which is a $\left\lceil d/(m-2t)\right\rceil\times n$ matrix, whose $i$’th column is the encoding $\widetilde{\bf x}_{i}:={\bf S}{\bf x}_{i}$ of the $i$’th data point ${\bf x}_{i}$. Using the method developed in Section 4, given $\{{\bf S}_{j}{\bf x}_{i}+{\bf e}_{j}\}_{j=1}^{m}$, where ${\bf e}_{j}={\bf 0}$ if the $j$’th worker is honest, otherwise can be arbitrary, master can recover ${\bf x}_{i}$ exactly in $O((1+\epsilon)md)$ time. Our main theorem is stated below, a proof of which trivially follows from Section 4.

Theorem 3** (Stochastic Gradient Descent).**

Let ${\bf X}\in\mathbb{R}^{n\times d}$ denote the data matrix. Let $m$ denote the total number of worker nodes. We can compute a stochastic gradient in a distributed manner in the presence of $t$ corrupt worker nodes and $s$ stragglers, with the following guarantees, where $\epsilon>0$ is a free parameter.

• •

$(s+t)\leq\left\lfloor\frac{\epsilon}{1+\epsilon}\cdot\frac{m}{2}\right\rfloor$.

• •

Total storage requirement is roughly $(1+\epsilon)|{\bf X}|$.

• •

Computational complexity for each stochastic gradient computation:

• –

at each worker node is $O((1+\epsilon)\frac{d}{m})$.

• –

at the master node is $O((1+\epsilon)dm)$.

• •

Communication complexity for each stochastic gradient computation:

• –

each worker sends $\left((1+\epsilon)\frac{d}{m}\right)$ real numbers.

• –

master broadcasts $\lceil\log n\rceil$ bits.

• •

Total encoding time is $O\left(nd\left(\frac{\epsilon}{1+\epsilon}m+1\right)\right)$.

Observe the distributed gain of our method in the communication exchanged between the workers and the master: (i) master only broadcasts an index in $\{1,2,\ldots,n\}$, which only takes $\lceil\log n\rceil$ bits; and (ii) each worker sends roughly $\frac{1+\epsilon}{m}$ fraction of the total dimension $d$. Hence, this method is particularly useful in distributed settings with communication-constrained and band-limited links. The Remarks 2, 4, 5 are also applicable for Theorem 3.

Remark 10** (One-round vs. two-round approach).**

Unlike the two-round approach taken for gradient computation in PGD and also for CD, we give a one-round approach for each iteration of SGD. This is because in each SGD iteration we need to compute the gradient on only one data point (not the entire dataset, as in the case for each PGD iteration). Because of this, recovering a (random) data-point itself at the master and then computing a gradient on it locally (which is what we do) would be far more efficient than computing gradient on a single data point in a distributed manner. This is in contrast to each gradient computation for PGD, which requires computation of the full gradient (which is the summation of gradients on all $n$ data points). In principle, we can use the one-round for each PGD iteration also in which first we recover all the $n$ data points at master and then compute the full gradient locally, but this approach would defeat the purpose of distributed computation both in terms of storage and computational complexity. Note that our two-round approach for PGD is significantly more efficient than this.

The reason behind taking the two-round approach for CD is because in order to update the local parameter vectors in the $t$’th iteration, workers need access to the MV multiplication $\widetilde{\bf X}^{R}{\bf v}^{t}={\bf X}{\bf w}^{t}$ (see the paragraph after (18) for more details), and in order to provide that we use an extra round – the first round is used for computing ${\bf X}{\bf w}$ and the second round is used for updating the local parameter vectors. Again, for CD also, we could adopt a one-round approach where master recovers all the $n$ data points and then do the parameter update, but that would be highly inefficient and defeat the purpose of distributed computation.

One of the main advantages of the one-round approach for SGD is that since we are recovering the data point itself at the master, we can use it to optimize any function, both convex and non-convex. This is in contrast to the two-round approach, which can only be used for generalized linear models only.

6.2 Encoding in The Streaming Data Model

An attractive property of our encoding scheme is that it is very easy to update with new data points. More specifically, our encoding requires the same amount of time, irrespective of whether we get all the data at once, or we get each sample point one by one, as in the online/streaming model. This setting encompasses a more realistic scenario, in which we design our coding scheme with the initial set of data points and distribute the encoded data among the workers. Later on, when we get some more samples, we can easily incorporate them into our existing encoded data. We show that updating $(m-2t)$ new data points in $\mathbb{R}^{d}$ requires $O\left((m-2t)\left((2t+1)d\right)\right)$ time in total, i.e., $O\left((2t+1)d\right)$ amortized-time per data point. This is the best one can hope for, since the offline encoding of $n$ data points requires $O\left((2t+1)nd\right)$ total time. At the end of the update, the final encoded matrix that we get is the same as the one we would have got had we had all the $n+1$ data points in the beginning. Therefore, the decoding is not affected by this method at all. Note that we use the same encoding matrices both for gradient computation as well as for coordinate descent. So, it suffices to prove our result in the streaming model for any one of them, and we show it for gradient computation below.

Theorem 4**.**

The total time complexity in encoding all the data points at once, i.e., when encoding is done offline, is the same as the total time complexity in encoding the data points one by one as they come in the streaming model, i.e., when encoding is done online.

Proof.

Let ${\bf S}^{(1)}$ and ${\bf S}^{(2)}$ denote the encoding matrices for encoding ${\bf X}$ and ${\bf X}^{T}$, respectively; see Section 4.2. For convenience, we copy over the corresponding encoding matrices ${\bf S}_{i}^{(1)}$ and ${\bf S}_{i}^{(2)}$ from (11) for the $i$’th worker node in Figure 3.

Suppose at some point of time we have encoded $n$ data points each lying in $\mathbb{R}^{d}$ and distributed the encoded data among the $m$ worker nodes. Now a new data sample ${\bf x}\in\mathbb{R}^{d}$ comes in. We will show how to incorporate it in the existing scheme in $O\left((2t+1)d\right)$ time on average.

Updating the encoding matrices.

Fix an arbitrary worker $i\in[m]$. Note that the new data matrix ${\bf X}$ has dimension $(n+1)\times d$. So, the new encoding matrix ${\bf S}_{i}^{(1)}$ should have $(n+1)$ columns, and we have to add one more column to ${\bf S}_{i}^{(1)}$. By examining the repetitive structure of ${\bf S}_{i}^{(1)}$, it is obvious which column to add: if $l_{1}<q$, then we add the $p_{1}$-dimensional vector $[0,0,\ldots,0,b_{(l_{1}+1)i}]^{T}$ as the last column; otherwise, if $l_{1}=q$, then we add the $(p_{1}+1)$-dimensional vector $[0,0,\ldots,0,b_{1i}]^{T}$ as the last column. In the second case, the number of rows of ${\bf S}_{i}^{(1)}$ increases by one – the last row has all zeros, except for the last element, which is equal to $b_{1i}$. Note that ${\bf S}_{i}^{(2)}$ does not change at all. Observe that if the $i$’th worker performs this update, then it does not have to store its entire encoding matrix ${\bf S}_{i}^{(1)}$, it only needs to store $n$, $q=(m-k)$, and the $q$ real numbers $b_{1i},b_{2i},\ldots,b_{qi}$, where $q=m-k$, which could be much smaller as compared to $n$ and $d$, and are enough to define ${\bf S}_{i}^{(1)}$ and ${\bf S}_{i}^{(2)}$.

Updating the encoded data.

Now we show how to update the encoded data with the new sample ${\bf x}$. We need to update both ${\bf S}_{i}^{(1)}{\bf X}$ as well as ${\bf S}_{i}^{(2)}{\bf X}^{T}$ for every worker $i\in[m]$.

• •

Updating ${\bf S}_{i}^{(1)}{\bf X}$. If $l_{1}<q$, then we add $b_{(l_{1}+1)i}{\bf x}^{T}$ to the last row of ${\bf S}_{i}^{(1)}{\bf X}$; otherwise, if $l_{1}=q$, then we add $b_{1i}{\bf x}$ as a new row in ${\bf S}_{i}^{(1)}{\bf X}$. In the first case, the resulting matrix still has $p_{1}$ rows, whose first $p_{1}-1$ rows are same as before, and the last row is the sum of the previous row and $b_{(l_{1}+1)i}{\bf x}^{T}$. In the second case, the resulting matrix has $(p_{1}+1)$ rows, whose first $p_{1}$ rows are the same as before and the last row is equal to $b_{1i}{\bf x}^{T}$. Note that each row of ${\bf S}_{i}^{(1)}$ for $i\leq 2t$, has at most $(m-2t)$ non-zero elements; whereas, for $i>2t$, each row of ${\bf S}_{i}^{(1)}$ has exactly one non-zero entry. Since there are $p_{1}=\left\lceil n/(m-2t)\right\rceil$ rows in each ${\bf S}_{i}^{(1)}$, updating ${\bf S}_{i}^{(1)}{\bf X}$ for every $i\leq 2t$ takes $O(d)$ time; and for $i>2t$, update in ${\bf S}_{i}^{(1)}{\bf X}$ happens only once in $(m-2t)$ new data points (whenever the second case occurs and the resulting ${\bf S}_{i}^{(1)}$ has $(p_{1}+1)$ rows). So, updating $(m-2t)$ data points at all $m$ worker nodes require $O\left(2t*(m-2t)d+(m-2t)*d\right)=O((m-2t)(2t+1)d)$ time, i.e., $O\left((2t+1)d\right)$ time per data point.

• •

Updating ${\bf S}_{i}^{(2)}{\bf X}^{T}$. Note that ${\bf X}^{T}$ is a $d\times(n+1)$ matrix whose last column is equal to the new data sample ${\bf x}$. Now, to update ${\bf S}_{i}^{(2)}{\bf X}^{T}$, we add ${\bf S}_{i}^{(2)}{\bf x}$ as an extra column. The resulting matrix is of size $p_{2}\times(n+1)$, whose first $n$ columns are the same as before and the last column is equal to ${\bf S}_{i}^{(2)}{\bf x}$. Since total number of non-zero entries in ${\bf S}_{i}^{(2)}$ is equal to $d$ if $i\leq 2t$ and equal to $p_{2}=\left\lceil d/(m-2t)\right\rceil$ if $i>2t$, the total time required to update a new data point is $O(2t*d+(m-2t)*p_{2})=O\left((2t+1)d\right)$.

Observe that at the end of this local update at each worker node, the final encoded matrix that we get is the same as the one we would have got had we had all the $n+1$ data points in the beginning. The decoding is not affected by this method at all. This completes the proof of Theorem 4. ∎

Remark 11** (Updating the encoded data efficiently with new features).**

Observe that since we encode both $X$ and $X^{T}$ in an analogous fashion, it follows by symmetry that we can not only update efficiently upon receiving a new data sample, but can also update efficiently if we decide to enlarge the dimension $d$ of each of the $n$ data samples at some point of time – maybe we figure out some new features of the data to get a more accurate model to overcome under-fitting. In these situations, we don’t need to encode the entire dataset all over again, just a simple update is enough to incorporate the changes.

Remark 12** (What allows our encoding to be efficient for streaming data?).**

The efficient update property of our coding scheme is made possible by the repetitive structure of our encoding matrix (see Figure 3), together with the fact that this structure is independent of the number of data points $n$ and the dimension $d$ – it only depends on the number of worker nodes $m$ and the corruption threshold $t$. We remark that other data encoding methods in literature, even for weaker models, do not support efficient update. For example, the encoding of [KSDY17], which was designed for mitigating stragglers, depends on the dimensions $n$ and $d$ of the data matrix. So, it may not efficiently update if a new data point comes in.

6.3 More Applications.

There are many iterative algorithms, other than the gradient descent for learning GLMs, which use repeated MV multiplication. Some of them include (i) the power method for computing the largest eigenvalue of a diagonalizable matrix, which is used in Google’s PageRank algorithm [ISW06], Twitter’s recommendation system [GGL*+*13], etc.; (ii) iterative methods for solving sparse linear systems [Saa03]; (iii) many graph algorithms, where the graph is represented by a fixed adjacency matrix, [KG11]. In large-scale implementation of these systems, where Byzantine faults are inevitable, the method described in this paper can be of interest.

In most of these applications, the underlying matrix ${\bf A}$ is generally sparse, which is exploited to gain computational efficiency. So, it is desired not to lose sparsity even if we want resiliency against Byzantine attacks. Fortunately, our encoding matrix ${\bf S}$ is sparse (see (11)), which ensures that the encoded matrix ${\bf S}{\bf A}$ will not lose the sparsity of ${\bf A}$: Each of the first $pk$ rows of ${\bf S}$ has at most $(m-k)$ (where $k=2t$) non-zero elements, and each of the remaining rows has exactly one 1. Since $m$ is the number of worker nodes, which may be small, and we can take $t$ to be up to $\lfloor\frac{m-1}{2}\rfloor$, we may have a few non-zero entries in each row of ${\bf S}$ (in the extreme case when $2t=m-1$, each row of ${\bf S}$ has only one non-zero entry). In a sense, we are getting Byzantine-resiliency almost for free without compromising the computational efficiency that is made possible by the sparsity of the matrix.

7 Numerical Experiments

In this section, we validate the efficacy of our proposed methods by numerical experiments. We run distributed gradient descent (GD) and coordinate descent (CD) for linear regression $\arg\min_{{\bf w}\in\mathbb{R}^{d}}\|{\bf X}{\bf w}-{\bf y}\|_{2}^{2}$. As mentioned in Section 2.1, for linear regression (which is equal to ridge regression when $h=0$), the projected gradient descent (PGD) reduces to gradient descent (GD). Since we are doing exact computation (computing the gradients exactly in the case of GD and updating the coordinates exactly in the case of CD), (i) there is no need to check the convergence, and (ii) our algorithm will perform exactly the same whether we are working with synthetic datasets or real datasets, hence, we will work with a synthetic dataset. We run our algorithms232323We implement our algorithm in Python, and run it on an iMac machine with 3.8 GHz Quad-Core Intel Core i5 processor and 16 GB 2400 MHz DDR4 memory. with $m=15$ worker nodes on two datasets: $(n=10,000,d=250)$ and $(n=20,000,d=22,000)$. For both the datasets, we generate $({\bf X},{\bf y})$ by sampling ${\bf X}\leftarrow\mathcal{N}(0,I)$ and ${\bf y}={\bf X}\theta+{\bf z}$, where $\theta\in\mathbb{R}^{d}$ has $d/3$ non-zero entries, all of them are i.i.d. according to $\mathcal{N}(0,4)$, and each entry of ${\bf z}\in\mathbb{R}^{n}$ is sampled from $\mathcal{N}(0,1)$ i.i.d. In each round of the gradient computation, the adversary picks $t$ worker nodes uniformly at random, and adds independent random vectors of appropriate length as errors, whose entries are sampled from $\mathcal{N}(0,\sigma^{2})$ i.i.d. with $\sigma=100$, to the true vectors.

7.1 $n=10,000,d=250,m=15$

In Figure 4, we plot the total time taken (which is the sum of the maximum time taken by any single worker node and the time taken by the master node in both rounds) for updating different number of coordinates in one CD iteration, with varying number of corrupt worker nodes from $t=1$ to $t=7$. We plot the time needed for updating $\gamma$-fraction of $d$ coordinates for four different values of $\gamma$ (i.e., $\gamma=0.1,0.25,0.5,1$) and we denote it by CD$(\gamma d)$ for $\gamma=0.1,0.25,0.5,1$. Recall that CD$(d)$ is equivalent to full gradient computation as in the case of GD. Note that, when $t=7$, we have $\epsilon=m-1$, which is the main cause behind the significant increment in time for $t=7$.

7.2 $n=20,000,d=22,000,m=15$

In Figure 5, we report separately, the maximum time taken by any single worker node and the time taken by the master node (together in both the rounds) in one CD iteration for updating different number of coordinates and also for GD, with varying number of corrupt worker nodes from $t=1$ to $t=6$. As in the above case, we report the time needed for updating $\gamma$-fraction of $d$ coordinates for four different values of $\gamma$. Observe that the time taken by the master node is orders of magnitude less than the time taken by the worker nodes. We can also observe that with the running time in a worker node per iteration for CD($0.1d$) is 95% less than that for GD, while this time saving in the master node is more than 40%.

Acknowledgements

The work of Deepesh Data and Suhas Diggavi was partially supported by the Army Research Laboratory under Cooperative Agreement W911NF-17-2-0196, by the UC-NL grant LFR-18-548554, and by the NSF award 1740047. The work of Linqi Song was partially supported by the NSF awards 1527550, 1514531, by the City University of Hong Kong grant 7200594, and by the Hong Kong RGC ECS 21212419. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Laboratory or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation here on.