Data Science Capstone Project Ideas With Dataset

Data Encoding for Byzantine-Resilient Distributed Optimization∗

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.

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.

⌋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-Offbetween

the corruption threshold and the resource requirements (storage, computational, and communication

Complexity). As An Example, For T ≤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.

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].

Arxiv:1907.02664V2 [Cs.Dc] 4 Nov 2020

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 computa- tion. 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 adver- saries 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 mini- mization, 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 com- pute 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 parti- tioned 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,1 which also has been shown to be very effective for solving generalized lin- ear 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 knowl- edge, ours is the first paper that studies CD against Byzantine attacks and provides an efficient solution for that.

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 op- erating 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, 1Alternatively, 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.

2

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 9) for PGD and Theorem 2 (on page 11) for CD, and demonstrate a trade-offbetween the Byzantine resilience (in terms of the number of adversarial nodes) and the resource requirement (storage, computational, and communication complexity).

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

2M

m−2t (which is a constant, even if t is a constant (< 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.

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).

2Storage redundancy is defined as the ratio of the size of the encoded matrix and the size of the raw data matrix.

Notation

We denote vectors by bold small letters (e.g., x, y, z, etc.)

And Matrices By Bold Capital Letters (E.G.,

A, F, S, X, etc.). We denote the amount of storage required by a matrix X by |X|.

For Any Positive

integer n ∈N, we denote the set {1, 2, . . . , n} by [n]. For n1, n2 ∈N, where n1 ≤n2, we write [n1 : n2] to denote the set {n1, n1 + 1, . . . , n2}. For any vector u ∈Rn and any set S ⊂[n], we write uS to denote the |S|-length vector, which is the restriction of u to the coordinates in the set S. The support of a vector u ∈Rn is defined by supp(u) := {i ∈[n] : ui̸ = 0}. We say that a vector u ∈Rn is t-sparse if |supp(u)| ≤t.

While stating our results, we assume that performing the basic arithmetic operations (addition, subtraction, multiplication, and division) on real numbers takes unit time.

Problem Setting And Our Results

Given a dataset consisting of n labelled data points (xi, yi) ∈Rd × R, i ∈[n], we want to learn a model/parameter vector w ∈Rd, which is a minimizer of the following empirical risk minimization problem:

(1)

where fi(w), i = 1, 2, . . . , n, denotes the risk associated with the i’th data point with respect to w and

Pn

i=1 fi(w) the average empirical risk associated with the n data points with respect to w. Our main focus in this paper is on generalized linear models (GLM), where fi(w) = ℓ(⟨xi, w⟩; yi) for some differentiable loss function ℓ. Here, each fi : Rd →R is differentiable, h : Rd →R is convex but not necessarily differentiable, and ⟨xi, w⟩is the dot product of xi and w. We do not necessarily need each fi to be convex, but we require f(w) to be a convex function. Note that f(w) + h(w) is a convex function. In the following we study different algorithms for solving (1) to learn a GLM.

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 w0 ∈Rd, and then update the parameter vector according to the following

Wt+1 = Proxh,Αt(Wt −Αt∇F(Wt)),

t = 1, 2, 3, . . .

(2)

where α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 α, the proximal

2Α∥Z −W∥2

2 + h(z).

(3)

Observe that if h = 0, then proxh,α(w) = w for every w ∈Rd, 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 Fi(W) = 1

2(⟨xi, w⟩−yi)2 and h(w) = λ∥w∥1. It turns out that proxh,α(z) for Lasso is equal to the soft-thresholding operator Sλα(z) [Tib15], which, for j ∈[d], is defined as

Zj −Λα

if zj > λα.

4

• 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 minw∈C f(w), where C ⊆Rd is a closed, convex set. Define an indicator function IC for C as follows: IC(w) := 0, if w ∈C; and IC(w) := ∞, otherwise. Now, observe the following equivalence

W∈C F(W) ⇐⇒Min

w∈Rd f(w) + IC(w). 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 C [Tib15]. So, the proximal gradient update step is to compute the usual gradient and then project it back onto the set C.

• Logistic regression. Here fi is the logistic function, defined as

,

where ui = ⟨xi, w⟩, and h = 0. As noted earlier, since h = 0, PGD reduces to GD for logistic regression.

2. Since Fi’S And H Are Differentiable,

we can alternatively solve this simply using GD. Let X ∈Rn×d denote the data matrix, whose i’th row is equal to the i’th data point xi. For simplicity,

M +J. In

a distributed setup, all the data is distributed among m worker machines (worker i has Xi) and master updates the parameter vector using the update rule (2). At the t’th iteration, master sends wt to all the workers; worker i computes the gradient (denoted by ∇if(wt)) on its local data and sends it to the master; master aggregates all the received m local gradients to obtain the global gradient

I=1

∇if(wt).

(4)

Now, master updates the parameter vector according to (2) and obtains wt+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.

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 X ∈Rn×d denote the data matrix and y ∈Rn the corresponding label vector. To make it distinct from the last section, we denote the objective function by φ and write it as φ(Xw; y) to emphasize that we want to learn a GLM, where the objective function depends

5

on the data points only through their inner products with the parameter vector. Formally, we want to

!

.

(5)

For U ⊆[d], we write ∇Uφ(Xw; y) to denote the gradient of φ(Xw; y) with respect to wU, where wU denotes the |U|-length vector obtained by restricting w to the coordinates in U. To make the notation less cluttered, let φ′(Xw; y) denote the n-length vector, whose i’th entry is equal to ℓ′(⟨xi, w⟩; yi) :=

∂

∂uℓ(u; yi)|u=⟨xi,w⟩. Note that ∇φ(Xw; y) = XT φ′(Xw; y) and that ∇Uφ(Xw; y) = XT

Uφ′(Xw; Y), Where Xu Denotes The N×|U|

matrix obtained by restricting the column indices of X to the elements in 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 X = [X1 X2 . . . Xm] and w = [wT

M Vector. Each Worker I Stores Xi And Is

responsible for updating (a few coordinates of) wi – hence the terminology, model-parallelism. We store the label vector 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 S ⊂[m] of workers and asks them to update their wi’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 wi’s. Master can (deterministically or randomly) pick a subset U (which may or may not be different for all workers) of f ≤d/m coordinates and asks each worker to updates only those coordinates. If master picks U deterministically, it can cycle through and update all coordinates of the parameter vector in ⌈d/mf⌉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 U. We will adopt this approach in our method to make the distributed CD Byzantine-resilient. Let r =

M. For Any I ∈[M], Let

wi = [wi1 wi2 . . . wir]T and Xi = [Xi1 Xi2 . . . Xir], where Xij is the j’th column of Xi. For any i ∈[m] and U ⊆[r], let wiU denote the |U|-length vector that is obtained from wi by restricting its entries to the coordinates in U; similarly, let XiU denote the n × |U| matrix obtained by restricting the column indices of Xi to the elements in U.

In Algorithm 1, for each worker i to update wi according to (6), where the partial gradient of φ with

Iuφ′(Pm

j=1 Xjwj; y) and worker i has only (Xi, wi), every other worker j sends Xjwj to the master, who computes φ′(Pm

J=1 Xjwj; Y)5 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.6 Similarly, corrupt workers can send adversarially chosen information to affect the stopping criterion.

3Here we are not optimizing the average of loss functions – since n is a fixed number, this does not affect the solution space. 4After the 1st iteration, worker i need not multiply Xi with wi to obtain Xiwi in every iteration; as only a few coordinates of wi are updated, it only needs to multiply those columns of Xi that corresponds to the updated coordinates of wi.

5Note that even after computing Xw, master needs access to the labels yi, i = 1, 2, . . . , n to compute φ′(Xw; y). Since y ∈Rn is just a vector, we can either store that at master, or, alternatively, we can encode 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 y.

6Specifically, suppose the i’th worker is corrupt and the adversary wants master to compute φ′(Xw + e; y) for any arbitrary vector e ∈Rn of its choice, then the i’th worker can send Xiwi + e to the master.

Algorithm 1 Distributed Coordinate Descent

1: Initialize. Each worker i ∈[m] starts with an arbitrary/random wi ∈Rr, where r =

M And, For

simplicity, we assume that m divides d. 2: while (until the stopping criteria at master is not satisfied) do

4:

Worker i computes Xiwi and sends it to the master node.4

5:

Worker i receives (U ⊆[r], φ′(Xw; y)) from the master node.

6:

Worker i updates its local parameter vector as (where ∇iUφ(Xw; y) = XT

(6)

while keeping the other coordinates of wi unchanged, and sends the updated wi to the master.

8:

Master receives {Xiwi}i∈[m] from the m workers.

Master First Computes Xw = Pm

i=1 Xiwi and then computes φ′(Xw; y).

10:

Master picks U ⊆[r] (where U can be picked either randomly or in a round-robin fashion) and sends (U ⊆[r], φ′(Xw; y)) to all workers.

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 < m

2 Worker Nodes7, 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 ui+ei, where ui is the true vector and ei is the error vector, where ei = 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 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].

7Our 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.

Our Approach To Gradient Computation

Recall that fi(w) = ℓ(⟨xi, w⟩; yi) for some differentiable loss function ℓ, and the gradient of fi at w is equal to ∇fi(w) = (xi)T ℓ′(⟨xi, w⟩; yi), where ℓ′(⟨xi, w⟩; yi) :=

∂Uℓ(U; Yi)|U=⟨Xi,W⟩. Note That ∇Fi(W) ∈Rd Is A

column vector. Let f ′(w) denote the n-length vector whose i’th entry is equal to ℓ′(⟨xi, w⟩; yi). With this

I=1 Fi(W), We Have ∇F(W) = 1

nXT f ′(w). Since n is a constant, it is enough to compute XT f ′(w). So, for simplicity, in the rest of the paper we write

∇F(W) = Xt F ′(W),

∀w ∈Rd.

(7)

A natural approach to computing the gradient ∇f(w) is to compute it in two rounds: (i) compute f ′(w) in the 1st round by first multiplying X with w and then master locally computes f ′(w) from Xw (master can do this locally, because Xw is an n-dimensional vector whose i’th entry is equal to ⟨xi, w⟩and (f ′(w))i = ℓ′(⟨xi, w⟩; yi));8 and then (ii) compute ∇f(w) = XT f ′(w) in the 2nd round by multiplying XT with f ′(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 17 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 infor- mation 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 X using a sparse encoding matrix S(1) = [(S(1)

I X

at the i’th worker node; and for the 2nd round, we encode XT using another sparse encoding matrix

M )T ]T , And Store S(2)

i XT at the i’th worker node. Now, in the 1st round of the gradient computation at w, the master node broadcasts w and the i’th worker node replies with S(1)

I Xw

(a corrupt worker may report an arbitrary vector); upon receiving all the vectors, the master node applies error-correction procedure to recover Xw and then locally computes f ′(w) as described above. In the 2nd round, the master node broadcasts f ′(w) and similarly can recover XT f ′(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 v to a fixed matrix 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.

8Note that even after computing Xw, master needs access to the labels yi, i = 1, 2, . . . , n to compute f′(w). See Footnote 5 for a discussion on how master can get access to the labels.

8

A trivial approach. Take a generator matrix G of any real-error correcting linear code. Encode A as AT G =: B. Divide the columns of B into m groups as B = [B1 B2 . . . Bm], where worker i stores Bi. Master broadcasts v and each worker i responds with vT Bi + eT

I , Where Ei = 0 If The I’Th Worker Is Honest,

otherwise can be arbitrary. Note that at most t of the ei’s can be non-zero. Responses from the workers can be combined as vT B + eT . Since every row of B is a codeword, vT B = vT AT G is also a codeword. Therefore, one can take any off-the-shelf decoding algorithm for the code whose generator matrix is G and obtain vT AT . 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 < 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 Ω(n2 + d2) 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 ∼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 X ∈Rn×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 ϵ > 0 is a free parameter.

K

. • Total storage requirement is roughly 2(1 + ϵ)|X|. • Computational complexity for each gradient computation:

– At Each Worker Node Is O((1 + Ε) Nd

m ). – at the master node is O((1 + ϵ)(n + d)m). • Communication complexity for each gradient computation:

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



. Remark 2. The statement of Theorem 1 allows for any s and t as long as (s + t) ≤

. As We Are

handling both erasures and errors in the same way9 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 I ⊂[m] with |I| ≤t, with the understanding that this can also work with stragglers.

In Theorem 1, ϵ is a design choice and a free parameter that can take any value in the interval [0, m−1], where ϵ = 0 implies no corruption and ϵ = m −1 implies that corruption threshold t can be anything up to

M−1

2 . If we want to tolerate t corrupt workers, then ϵ must satisfy ϵ ≥

M−2T.10

9When there are only stragglers, one can design an encoding scheme where both the master and the worker nodes oper- ate oblivious to encoding, while solving a slightly altered optimization problem [KSDY17], in which gradients are computed approximately, leading to more efficient straggler-tolerant GD.

10We could have written everything in terms of t, m, n, d, but we chose to introduce another variable ϵ which, in our opinion, clearly brings out the tradeoffbetween the corruption threshold and the resource requirements without cluttering the expressions.

9

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 w to all the workers; upon receiving w, all workers compute the gradients on their local data in O( 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 + ϵ) factor more;11 (see also Remark 4) (ii) the amount of computation at each worker node is O(1 + ϵ) factor more; (iii) the amount of computation at the

Master Node Is O((1 + Ε)(1 + N

d )) factor more, which is comparable in cases where n is not much bigger than

D; (Iv) Master Broadcasts (1 + N

d ) factor more data, which is comparable if n is not much bigger than d; and



factor more data, which is O(1 + ϵ) – 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 δ > 0, if we want to get δ close to m/2, i.e., t = m/2 −δ, then we must have (1 + ϵ) ≥m/2δ. In particular, at ϵ = 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 δ is a constant, i.e., t is close to

If T = M−1

2 , then ϵ = m −1. In this case, our storage redundancy factor is O(m). In contrast, the trivial scheme (see “trivial approach” on page 9) 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 ≤C · M−1

2 , where c < 1 is a constant, and use the trivial scheme if t is close to m−1 2 .

Time. Note That O(Nd) Is Equal To The Time

required for distributing the data matrix 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 re-

Quirement), But It Can Also Tolerate Up To ⌊M−1

2 ⌋corrupt worker nodes (by taking ϵ = m −1 in Theorem 1). It is not hard to prove that this bound is information-theoretically optimal, i.e., no algorithm can tolerate

⌈M

2 ⌉corrupt worker nodes, and at the same time correctly computes the gradient.

Our Approach To Coordinate Descent

We use data encoding and add redundancy to enlarge the parameter space. Specifically, we encode the data matrix X using an encoding matrix R = [R1 R2 . . . Rm], where each Ri is a d×p matrix (with pm ≥d), and store XRi at the i’th worker. Define eXR := XR. Now, instead of solving (5), we solve the encoded problem arg minv∈Rpm φ( eXRv; y) using Algorithm 1 (together with decoding at the master); see Figure 2 on page 25 for a pictorial description of our algorithm. We design the encoding matrix R such that at every iteration of our algorithm, updating any (small) subset of coordinates of vi’s (let v = [vT

M]) Automatically

updates some (small) subset of coordinates of w; and, furthermore, by updating those coordinates of vi’s, we can efficiently recover the correspondingly updated coordinates of w, despite the errors injected by the adversary. In fact, at any iteration t, the encoded parameter vector vt and the original parameter vector wt satisfies vt = R+wt, where R+ := RT (RRT )−1 is the Moore-Penrose pseudo-inverse of R, and wt evolves in the same way as if we are running Algorithm 1 on the original problem.

11For example, by taking ϵ = 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.

10

We will be effectively updating the coordinates of the parameter vector 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 vi, then k(m −2t) coordinates of w will get updated. For comparison, Algorithm 1 updates km coordinates of the parameter vector 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 wi according to (6), worker i needs access to φ′(Xw; y), which master computes after receiving {Xjwj}j∈[m] from the workers. In our Byzantine-resilient algorithm for CD also master will need to compute Xw 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 Xw, we encode X using an encoding

2 . . . Lt

m]T , where each Li is a p′×n matrix (with p′m ≥n) and worker i stores eXL i = LiX. Note that in order to compute Xw, in the first round of gradient computation as described in Section 2.4, master broadcasts w to all the workers and each worker i computes eXL

I W And Sends It The The Master (Corrupt

workers may report arbitrary vectors), who then decodes and obtains Xw. However, in coordinate descent, though master wants to compute Xw in each CD iteration, we can significantly improve the computation required at each worker: since only a few coordinates of the original parameter vector 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 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 ϵ > 0 is a free parameter.

K

. • Total storage requirement is roughly 2(1 + ϵ)|X|. • If each worker i updates τ coordinates of vi, then

Τm

1+ϵ coordinates of the corresponding w gets updated.

– The Computational Complexity In Each Iteration

∗at each worker node is O(nτ). ∗at the master node is O((1 + ϵ)nm + τm2).

real numbers.



real numbers.



. 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 protec- tion. Let ϵ be any number in the interval [0, m −1] – for illustration, we can take ϵ = 2, which means t ≤m

3

workers are corrupt. In Algorithm 1, if each worker i updates

1+Ε Coordinates

of w) in each iteration, then (i) each worker requires O( nτ

Wi; (Ii) Master Requires O(Nm) Time To Compute Pm

i=1 Xiwi from {Xiwi}i∈[m]; (iii) each worker sends n real numbers (required for Xiwi) to master; and (iv) master broadcasts n real numbers (required for φ′(Xw; y)).

11

In our scheme (i) the total storage requirement is O(1+ϵ) factor more; (ii) the amount of computation at each worker node is O(1+ϵ) factor more; (iii) the amount of computation at the master node is O((1+ϵ)+ τm

N )

factor more – typically, since τ is a constant and number of workers is much less than n, this again could be



factor more data, which could be a constant if τm is smaller



factor more data, where the 1st term is much smaller than 1 as τ is typically a constant, and the 2nd term is close to zero as (1 + ϵ) 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

2T+1 Groups Of Size (2T + 1)

each and also divide the data matrix as X = [X1 X2 . . . X

M

2t+1 ] (assume, for simplicity, that (2t+1) divides m). Now, store the i’th block Xi at the (2t + 1) workers in the i’th group of workers. Let the parameter

. . . Wtm

2t+1 ]T . In each CD iteration, the local parameter updates in any wi 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 = m

3 , With A

constant overhead in storage and computational complexity. The Remarks 2, 4, 5 are also applicable for Theorem 2.

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 compu- tation [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 up- dates 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 distri- bution) 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 gra- dient 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 mas- ter 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 Ω(

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

12

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(m3d) 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.

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.

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 + ϵ) =

Constant (< 1

2) fraction of m.12 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

M

m−2t) more computation at worker, which is a constant even if t is a constant (< 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(

M−2T(1 + N

d )) more than Draco. In high-dimensional settings, where n is not much bigger than d, and

T Is A Constant (< 1

2) fraction of m, this overhead is constant. Overall, for a constant fraction of corruption,

Say, T = M

3 , Draco requires Ω(t) times more storage and computation at workers than our scheme (which could be significant in large-scale settings), and requires Ω(1 + n

D ) Times Less Computation At Master. Note

that the computation time at workers scales at least as Ω( 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.

12To highlight the storage redundancy gain of our method over that of Draco, consider the following two concrete scenarios, where the data matrix X ∈Rn×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 × nd real numbers, whereas Draco stores 201 × 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 ≈17 more than ours.

13

Yu et al. [YLR+19] (which is a concurrent work13) 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 m−1

M

1+2δ −1 in our setting; here δ can take any value in [0 : m−1

M

δ+1, which is roughly the same as ours. For example,

3 Corrupt Workers (I.E., Δ = M−3

6 ), the storage overhead of our scheme and of [YLR+19] is a

6

1+3/m ≈6, respectively. (ii) The encoding time complexity of our scheme is

O(Nd(M−2Δ)), Whereas, It Is O(M Log2(M) Nd

δ+1) in [YLR+19]. Note that for constant δ (i.e., corruption close to 1/2), the encoding time of our scheme is much less (by a factor of O(m log2(m))) than that of [YLR+19],

Log2(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(

1+Δ) Time. (Iv) The

decoding time complexity per gradient computation in [YLR+19] is O(m log2(m)d), whereas, ours requires O((1 + ϵ)(n + d)m) time. Note that when n is not much bigger than d and we want a constant fraction of

Corruption, Say, 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,

N+D

1+2δ and d real numbers in ours and the scheme in [YLR+19]. Note that if n ≤dm and to tolerate a constant fraction of corruption, say, 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 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

M

log2(m)).

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 compu- tation 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 13Yu 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.

14

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.

Matrix-Vector Multiplication

For the task of a more fundamental problem of matrix-vector (MV) multiplication in the presence of Byzan- tine 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]14 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 require- ment, 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, M

3 corrupt workers, our method outperforms the one in [YLR+19] in terms of the decoding time complexity by a factor of O(log2(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 14The 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.

15

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.

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 A ∈Rnr×nc and a vector v ∈Rnc, we want to compute Av 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 A is encoded and distributed among all the worker nodes, and the master node recovers the MV product Av using real-error correction; see Figure 1. We will think of our encoding matrix as S = [ST

M], Where Each Si Is

a p × nr matrix and pm ≥nr. We will derive the matrix S in Section 4.2. For the value of p, looking

3 , We Would Have P = 3N

m ). For i ∈[m], we store the matrix SiA at the i’th worker node. As described in Section 2, the computation proceeds as follows: The master sends v to all the worker nodes and

Receives {Siav + Ei}M

i=1 back from them. Let ei = [ei1, ei2, . . . , eip]T for every i ∈[p]. Note that ei = 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 {Siav + Ei}M

i=1 as p systems of linear equations.

(8)

where, for every i ∈[p], ˜ei = [e1i, e2i, . . . , emi]T , and ˜Si is an m × nr matrix whose j’th row is equal to the i’th row of Sj, for every j ∈[m]. Note that at most t entries in each ˜ei are non-zero. Observe that

I=1 And {˜Siav + ˜Ei}P

i=1 are equivalent systems of linear equations, and we can get one from the other. Note that ˜Si’s constitute the encoding matrix S, which we have to design. In the following, we will design these matrices ˜Si’s (which in turn will determine the encoding matrix S), with the help of another matrix F, which will be used to find the error locations, i.e., identities of the compromised worker nodes. We will design the matrix F (of dimension k × 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 ˜Si’s such that C.1 F˜Si = 0 for every i ∈[p].

C.2 For any t-sparse u ∈Rm, we can efficiently find all the non-zero locations of u from Fu. C.3 For any T ⊂[m] such that |T | ≥(m −t), let ST denote the |T |p × nr matrix obtained from S by restricting it to all the Si’s for which i ∈T . We want ST to be of full column rank.

If we can find such matrices, then we can recover the desired MV multiplication Av 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 ST Av, where ST is the |T |p × nr matrix obtained from S by restricting it to Si’s for all i ∈T , where T is the set of all honest worker nodes.

Now, by C.3, since ST is of full column rank, we can recover Av from ST Av exactly. Details follow. Suppose we have matrices F and ˜Si’s such that C.1 holds. Now, multiplying (8) by F yields

(9)

for every i ∈[p], where ∥˜ei∥0 ≤t. In Section 4.1, we give our approach for finding all the corrupt worker nodes with the help of any error locator matrix F. Then, in Section 4.2, we give a generic construction for

W ←−Proxh,Α(W −Α∇F(W))

Figure 1 This figure shows our 2-round approach to the Byzantine-resilient distributed gradient descent to optimize (1) for learning a generalized linear model. Since the gradient at w is equal to ∇f(w) = XT f′(w) (see (7)), we compute it in 2 rounds, using a matrix-vector (MV) multiplication as a subroutine in each round. In the 1st round, first we compute Xw, and then compute f′(w) from Xw – since the j’th entry of Xw is equal to ⟨xj, w⟩, we can compute f′(w) from Xw (see Section 2.4).

In the 2nd round we compute XT f′(w) – which is equal to ∇f(w) – using another application of MV multiplication. For a matrix A and a vector v, to make our distributed MV multiplication Av Byzantine-resilient, we encode A using a sparse matrix

M . . . St

m]T and distribute SiA to worker i (denoted by Wi). Note that in the first round, we have A = X, v = w, and we encode X using S(1), and in the second round, we have A = XT , v = f′(w), and encode XT using S(2). The adversary can corrupt at most t workers (the compromised ones are denoted in red color), potentially different sets of t workers in different rounds. The master node (denoted by M) broadcasts v to all the workers. Each worker performs the local MV product and sends it back to M. If Wi is corrupt, then it can send an arbitrary vector. Once the master has received all the vectors (out of which t may be erroneous), it sends them to the decoder (denoted by Dec), which outputs the correct MV product Av.

designing ˜Si’s (and, in turn, our encoding matrix S) such that C.1 and C.3 hold. In Section 4.3, we show how to compute the desired matrix-vector product Av 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 F that we use in our construction.

Remark 8. As we will see in Section 4.2, the structure of our encoding matrix S is independent of our error locator matrix F. Specifically, the repetitive structure of the non-zero entries of S as well as their locations will not change irrespective of what the F matrix is. This makes our construction very generic, as we can choose whichever 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!

Finding The Corrupt Worker Nodes

Observe that supp(˜ei) may not be the same for all i ∈[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 I = Sp

I=1 Supp(˜Ei), Which Is The Set

of all corrupt worker nodes. Note that |I| ≤t. We want to find this set I efficiently, and for that we note the following crucial observation. Since the non-zero entries of all the error vectors ˜ei’s occur in the same set I, a random linear combination of ˜ei’s has support equal to 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

17

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, dis- crete 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.

I=1 Αi˜Ei, Where Αi’S Are Sampled I.I.D. From An

absolutely continuous distribution. Then with probability 1, we have supp(ˆe) = I. From (9) we have fi = F˜ei for every i ∈[p]. Take a random linear combination of fi’s with coefficients αi’s chosen i.i.d. from an absolutely continuous distribution, for example, the Gaussian distribution. Let

I=1 Αi˜Ei. Note That, With Probability

1, supp(˜e) is equal to the set of all corrupt worker nodes, and we want to find this set efficiently. In other words, given F˜e, we want to find supp(˜e) efficiently. For this, we need to design a k × m matrix F (where k < m) such that for any sparse error vector e ∈Rm, we can efficiently find supp(e) from f = Fe. 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 S. In particular, we will use a k × m Vandermonde matrix along with the Reed-Solomon type decoding, which can correct up to k/2 errors and has decoding complexity of O(m2); 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 F˜ei for every i ∈[p]: Note that we can get F˜ei by multiplying (8) with F. Since F is a k × m matrix, and we compute F˜hi(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(m2) time. Since p is much bigger than m, the total time complexity is O(pkm).

Designing The Encoding Matrix S

Now we give a generic construction for designing ˜Si’s such that C.1 and C.3 hold. Fix any k × m matrix F such that we can efficiently find e from Fe, provided e is sufficiently sparse. We can assume, without loss of generality, that F has full row-rank; otherwise, there will be redundant observations in Fe that we can discard and make F smaller by discarding the redundant rows. Let N(F) ⊂Rm denote the null-space of F. Since rank(F) = k, dimension of N(F) is q = (m −k). Let {b1, b2, . . . , bq} be a basis of N(F), and let bi = [bi1 bi2 . . . bim]T , for every i ∈[q]. We set bi’s the columns of the following matrix F⊥:

B21

. . .

B22

. . .

Bq2

... ... ...

...

B2M

. . .

(10)

The following property of F⊥will be used for recovering the MV product in Section 4.3. Claim 1. For any subset T ⊂[m], such that |T | ≥(m −t), let F⊥

The Restriction Of F⊥To The Rows In T . Then F⊥

T is of full column rank.

18

Proof. Note that q = m −k, where k = 2t. So, if we show that any q rows of F⊥are linearly independent, then, this in turn will imply that for every T ⊂[m] with |T | ≥(m −t), the sub-matrix F⊥

T Will Have Full

column rank. In the following we show that any q rows of F⊥are linearly independent. To the contrary, suppose not; and let T ′ ⊂[m] with |T ′| = q be such that the q ×q matrix F⊥

T ′ Is Not A Full Rank Matrix. This

implies that there exists a non-zero c′ ∈Rq such that F⊥ T ′c′ = 0. Let b = F⊥c′. Note that b̸ = 0 (because columns of F⊥are linearly independent) and also that ∥b∥0 ≤m −q = k. Now, since FF⊥= 0, we have Fb = 0, which contradicts the fact that any k columns of F are linearly independent.

Now we design ˜Si’s. For i ∈[p], we set ˜Si as follows:

0

. . .

B21

. . .

0

. . .

0

0

. . .

B22

. . .

0

. . .

0

... ... ...

... ... ...

... ... ...

...

0

. . .

B2M

. . .

0

. . .



where l = q if i < p; otherwise l = nr −(p −1)q. The first (i −1)q and the last nr −[(i −1)q + l] columns of ˜Si are zero. This also implies that the number of rows in each Si is p = ⌈nr/q⌉. Claim 2. For every i ∈[p], we have F˜Si = 0.

Proof. By construction, the null-space of F is N(F) = span{b1, b2, . . . , bq}, which implies that Fbi = 0, for every i ∈[q]. Since all the columns of ˜Si’s are either 0 or bj for some j ∈[q], the claim follows. The above constructed matrices ˜Si’s give the following encoding matrix Si for the i’th worker node:

B1I . . . Bqi

...

(11)

All the unspecified entries of Si are zero. The matrix Si is for encoding the data for worker i. By stacking up the Si’s on top of each other gives us our desired encoding matrix S. To get efficient encoding, we want S to be as sparse as possible. Since S is completely determined by F⊥, whose columns are the basis vectors of N(F), it suffices to find a sparse basis for N(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 N(F) by reducing F to its row-reduced-echelon-form (RREF) using the Gaussian elimination [HK71]. This will result in F⊥whose last q rows forms a q × 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 F⊥that we obtain by the RREF will be very sparse – only the first 2t rows may be dense. Since computing S is equivalent to computing F⊥, and we can compute F⊥in O(k2m) time using the Gaussian elimination, the time complexity of computing S is also O(k2m).

Now we prove an important property of the encoding matrix S that will be crucial for recovery of the desired matrix-vector product. Claim 3. For any T ⊂[m] such that |T | ≥(m −t), let ST denote the |T |p × nr matrix obtained from S by restricting it to all the blocks Si’s for which i ∈T . Then ST is of full column rank.

Proof. For i ∈[p −1], let Bi = [(i −1)q + 1 : iq] and Bp = [(p −1)q + 1 : nr −(p −1)q], where we see Bi’s as a collection of some column indices. Consider any two distinct i, j ∈[p]. It is clear that for any two vectors u1 ∈Bi, u2 ∈Bj, we have supp(u1) ∩supp(u2) = φ, which means that all the columns in distinct Bi’s are linearly independent. So, to prove the claim, we only need to show that the columns within the same Bi’s are linearly independent. Fix any i ∈[p], and consider the |T |p × q sub-matrix S(i)