We consider a matrix multiplication problem where the goal is to compute C = A T B given two input matrices A ∈ R s×r and B ∈ R s×t . Many applications in optimization and machine learning require multiplications of large matrices of dimension of the order of 10 5 × 10 5 . These large-scale matrix multiplications cannot be carried out in a single machine mainly due to low-latency requirement in many applications. This requirement can be met by dividing input matrices A and B into m and n sub-matrices, respectively, and distributing the tasks of computing the product of these sub-matrices among P worker nodes [1]- [12]. A master node then collects the results of the worker nodes, and aggregates them appropriately to obtain C. Since the master node requires the results of all worker nodes so as to successfully compute C, a few slow worker nodes, referred to as stragglers, can significantly increase the computational delay.
In [6], Lee et al. proposed a scheme, referred to as 1-D maximum distance separable (MDS) codes, to mitigate the effect of stragglers for the case of n = 1, in which the sub-matrices of the matrix A are encoded using a (P, m) MDS code, and the task of computing A T i B, for i ∈ [P ], is distributed among P worker nodes. This scheme can be extended to matrix-matrix multiplication in a natural way by considering each column of B as a vector; however, as shown in [2], the recovery threshold of this scheme -defined as the minimum number of worker nodes that the master node needs to wait for to compute C, will not be optimal. In [2], Yu et al. proposed a coding scheme, referred to as Polynomial codes, which divides both A and B into submatrices. As was shown in [2], the Polynomial codes achieve the minimum possible recovery threshold, i.e., K = mn. The main drawback of Polynomial codes is that the decoding process requires interpolating polynomials of degree K, which is equivalent to inverting a K × K Vandermonde matrix. However, inverting real-valued Vandermonde matrices is highly numerically unstable even for moderate K [13].
In [9], Lee et al. proposed a coding scheme using Product codes, where both the input matrices A and B are split (column-wise) into m = n sub-matrices, and encoded using an ( √ P , m) MDS code, unlike the scheme in [6]. Product codes can be implemented using several component codes such as Polynomial codes [1], OrthoPoly codes [14], Random Khatri-Rao-Product (RKRP) codes [15], etc. [4], [16]. Product codes are generally better than Polynomial codes in terms of numerical stability. For example, decoding of a Product code of dimension
whereas decoding of a Polynomial code of dimension K would require inversion of a K × K Vandermonde matrix. The main drawback of Product codes is that their recovery threshold is not optimal unlike Polynomial codes. In [17], Baharav et al. proposed a scheme, referred to as d-dimensional Product codes for matrix multiplication, by spreading component matrices of A and B over d/2 dimensions each, and encoding them using a d-dimensional Product code. While d-dimensional Product codes can perform better than 2-dimensional Product codes for certain regimes of P and K, the recovery threshold of d-dimensional Product codes is still not optimal when the number of stragglers is linear in P . Moreover, (P, K) d-dimensional Product codes of rate R are built upon (P 1/d , K 1/d ) MDS component codes of rate R 1/d . Unless P is very large, it is impractical to use large values of d since the set of possible MDS component codes becomes trivial for large d.
In [18], Mallick et al. proposed a scheme for matrix-vector multiplication using Luby Transform (LT) codes [19], a type of rateless codes, which has several desired properties: (i) high numerical stability; (ii) linear decoding complexity, and (iii) recovery threshold of P (1 -α), where α ∈ [0, 1) is the fraction of worker nodes that are stragglers. However, this scheme is not directly extendable to the matrix-matrix multiplication problem. In [20], Wang et al. proposed the use of LT codes for distributed matrix multiplication. The scheme of [20] has all the properties (i)-(iii); however, for this scheme, both the communication and the computation at each worker node are substantially more expensive than those of Polynomial codes and Product codes.
In this work, we propose two novel encoding schemes. The first scheme is based on LT codes, referred to as Factored LT (FLT) codes, which is better in terms of numerical stability as well as decoding complexity, when compared to Polynomial codes. In particular, the decoding complexity of FLT codes is O(rt log K), whereas the decoding complexity of Polynomial codes is O(rt log 2 K log log K), where r and t are the number of columns in A and B, respectively, and K = mn where m and n are the number of sub-matrices into which A and B are split, respectively. As in the case of LT codes for erasure channels, the performance of FLT codes can be improved for finite lengths, particularly for moderate values of K. For this regime, we propose a Raptor code based scheme, termed Factored Raptor (FR) codes, which performs well even for moderately large K. Our simulation results show that the recovery threshold of FR codes is better than that of Product codes. For example, a (10000, 6400) FR code has recovery threshold 7060, whereas the recovery threshold of a (21,18)
We use boldface capital letters for matrices and underlined variables to represent vectors. We denote the (i, j)th element of the matrix A by A ij , and represent the ith row and jth column of the matrix A by A i,: and A :,j , respectively. We assume that vectors without transposes are column vectors unless stated otherwise. We denote {1, . . . , i} by [i], and denote the cardinality of a set S by |S|.
We consider a system where there is one master node and P worker nodes. The goal of the master node is to compute C = A T B in a distributed fashion using P worker nodes. Worker nodes can only communicate with the master node and cannot communicate among themselves. To distribute the computation among workers, A and B are divided along columns to m and n sub-matrices of size s × r m and s × t n , respectively, as shown below.
Alternatively, the output matrix C can be written in terms of the components of A and B as follows:
Then, computing C is equivalent to computing mn blocks
where a p i for i ∈ I p A and b p j for j ∈ I p B are randomly sampled from a Gaussian distribution with zero mean and unit variance; and we choose a p i = b p j = 0 for i / ∈ I p A and j / ∈ I p B . The master node computes ( A p ) T and B p and sends them to the worker node p. The worker node p computes C p = ( A p ) T B p and returns C p to the master node. We refer to
as the encoding function. The master node collects C p 's from a subset of worker nodes, referred to as non-stragglers, and attempts to recover the matrix C from the results of the non-straggling worker nodes using a decoding function f . Given an encoding function C and a decoding function f , the recovery threshold is defined as the minimum integer N such that the master node can recover the matrix C from the results of any N (non-straggling) worker nodes.
The goal is to design an encoding function and a decoding function so as to minimize the recovery threshold.
In this section, we propose Factored LT (FLT) codes, which, in this paper, will be used as a component of the distributed matrix multiplication scheme proposed in Section III-B. The application of FLT codes as a standalone scheme for distributed matrix multiplication is left for future work.
In the following, we briefly describe the encoding and decoding of FLT codes.
1) Encoding: LT codes, introduced by Luby in [19], are a class of rateless erasure codes that can be used to generate a (potentially infinite) sequence of output symbols from K source symbols. The number of source symbols involved in generating an output symbol is referred to as the degree of the output symbol. The degree of output symbols follow a degree distribution Ω(x). In [19], it was shown that in the case of a single erasure channel, the source symbols can be recovered from any N = K(1 + ) output symbols with high probability, where is the overhead. In particular, as shown in [19], for some proper choice of degree distribution Ω(x), the overhead vanishes as K grows unbounded. To apply LT codes to the task of matrix-matrix multiplication, we treat A T i B j 's for i ∈ [m] and j ∈ [n] as source symbols, where A i 's and B j 's are component matrices of A and B, respectively. As described in Section II, the matrix-matrix multiplication problem requires each output symbol to be the product of A and B, where A and B are sum of randomly chosen chunks of A and B, respectively. Unlike the case of a single erasure channel, for the matrix-matrix multiplication problem a degree-d output symbol cannot be generated from any arbitrary d source symbols. Thus, the encoding of LT codes is not directly applicable to our setting.
Let Ω(x) = K i=1 Ω i x i be a degree distribution, where Ω i = 0 for all prime i > max(m, n). The encoding process performed by the master node for each worker node p ∈ [P ] is described as follows:
1. Randomly choose a degree d by sampling from the degree distribution Ω(x).
The master node computes ( A p ) T = A T 2 + A T 3 and B p = B 1 + B 3 , and sends them to the worker node p. The worker node p is tasked to compute
) and send it back to the master node. Notice that
3 is a linear combination of d = 4 source symbols, as desired. The encoding process is illustrated in Fig. 1.
2) Decoding: Without loss of generality, suppose that the master node collects results from the first N workers with N ≤ P . Given the above encoding scheme, we have
We use M ∈ R N ×mn to represent the above coefficient matrix. To guarantee decodability, the master node should collect results from enough number of workers such that the coefficient matrix M has full column rank. Given such M , the master node then uses a peeling decoding algorithm, which can be described using a bipartite graph as follows.
The bipartite graph is constructed with one part of the nodes (referred to as the left nodes) being the set of source symbols
and the other part (referred to as the right nodes) being the set of output symbols C p for p ∈ [N ]. There is an edge between the vertices A T i B j and C p if A T i B j is involved in the computation of C p . The peeling decoding algorithm works in iterations. In each iteration, the master node searches for a right node of degree 1. If such a node cannot be found, the decoding algorithm terminates.
Otherwise, the master node recovers the (unique) left node connected to the found degree-1 right node, and removes all edges adjacent to the recovered left node. The master node runs this process iteratively until all left nodes are recovered or no degree-1 right node can be found. An early termination of the decoding process (before all left nodes are recovered) yields a decoding failure; otherwise, the decoding process is dubbed successful.
In the FLT coding scheme presented in section III-A.1, the output symbols are generated by taking a linear combination of randomly chosen source symbols. Therefore, when the number of encoded symbols are close to the number of source symbols, a small fraction of source symbols remain uncovered by any output symbols, which means there exists a fraction of all zero columns in the coefficient matrix M. These uncovered source symbols cannot be recovered at the end of the peeling decoding process. This implies that FLT codes do not perform well for moderate values of N . This is a well known issue with LT codes [19]. To address this issue, we propose a Raptor code based scheme, termed Factored Raptor (FR) codes, which is described below.
1) Encoding: An FR code is an FLT code concatenated with an outer code. In the case of a single erasure channel, the source symbols are encoded using the outer code, and the input symbols of the LT code are formed by a codeword of the outer code. Unlike the case of a single erasure channel, the outer code of FR codes cannot be chosen arbitrarily, due to the structure imposed by the matrix-matrix multiplication problem. In distributed matrix multiplication problem, the worker nodes are asked to compute product of a linear combination of chunks of one matrix with that of the other matrix. Collection of these products must form a codeword of the FLT code, and the input symbols corresponding to the codeword of the FLT code must form a codeword of the outer code. This requirement is met by encoding A and B with an MDS code and sending the linear combination of chunks of the respective encoded matrix to the worker nodes. We briefly describe the encoding process below.
1. For m ≥ m and ñ ≥ n such that P = mñ, the master node encodes the matrix
respectively. 2. The master node then encodes A and B according to the degree distribution of an FLT code as described in Section III-A.1. By the above construction, C is a codeword of the FLT code, and A i B j 's for i ∈ [ m] and j ∈ [ñ] form a codeword of an ( m, m) × (ñ, n) outer Product code.
2) Decoding: Decoding of FR codes consists of alternating decoding iterations on the FLT code and the outer Product code. At the master node, encoding of FR codes induces a Tanner graph which has A i B j 's for i ∈ [ m] and j ∈ [ñ] as the input symbols (i.e., left nodes) and C p 's as the output symbols (i.e., right nodes). Decoding iterations of
Worker node p
Fig. 1: Encoding of a degree-4 node the FLT code proceed on this Tanner graph as described in Section III-A.
respectively. Using the same procedure as in Example 1 for generating an output symbol of an FLT code, for this example suppose each worker node p for p ∈ [10] is requested to compute C p as follows:
Example 3: In this example, we describe the decoding algorithm at the master node using the code described in Example 2. For this example, consider that the master node collects results from the worker nodes {1, 3, 5, 7}. The peeling decoding is run on the subgraph, denoted by G 0 , induced by the worker nodes {1, 3, 5, 7}. In iteration 0, the source symbol A T 1 B 2 is recovered from worker 1 since C 1 is a degree-1 node in G 0 . Peel the edges connected to the source symbol A T 1 B 2 in G 0 and denote the residual graph by G 1 . In iteration 1, the source symbol
Similarly, the source symbol A T 2 B 1 is recovered from worker 3 in iteration 3. Since there are no degree-1 symbols left after iteration 3, the decoder of the FLT code cannot proceed further. Decoding of the FR code can proceed further by decoding the outer Product code. Arrange the input symbols to the FLT code in a 3 × 2 matrix as shown in Fig. 2e. Both first row and first column have two computation results as shown in Fig. 2e. In iteration 4, the missing computation results in both first row and first column are recovered by decoding the (3, 2) MDS code. The first few iterations of the decoding process for this example are illustrated in Fig. 2.
3) Optimal Decoding by Inactivation Decoding: Maximum likelihood decoding of FLT and FR codes can be implemented using a low-complexity algorithm known as inactivation decoding [21], which is a combination of peeling decoding and matrix inversion. The inactivation decoding starts with peeling decoding, which continues until it encounters a stopping set. At this point, the decoder assigns a variable to the value of an unrecovered symbol, referred to as an inactivated symbol. Then, the peeling decoder starts again and computes the unrecovered symbols in terms of the variable representing the value of the inactivated symbol. This procedure is repeated again, when the peeling decoder is stuck, by choosing another inactivated symbol from the unrecovered input symbols. The decoding stops when all the input symbols are either recovered or inactivated. Finally, optimal decoding of the inactivated symbols is performed via Gaussian elimination, and the decoded values are back-substituted into the decoded input symbols which depend on them.
In our encoding scheme, the master node sends two matrices A p and B p to each worker node p, and each worker node computes only one product. Therefore, the computation and communication costs of our proposed scheme are the same as those of Polynomial codes and Product codes. Unlike the FLT coding scheme, the LT coding scheme proposed in [20] requires the master node to send (on average) log K chunks of the input matrices to each worker node, and each worker node must compute (on average) log K products.
In the following, we briefly describe the decoding complexity of FR codes. The peeling of every edge in the Tanner graph corresponds to performing rt K operations. Let d avg be the average degree of output degree distribution. Then, each block A T i B j will be involved in O( davgN K ) such operations on average. We are interested in the case when m -m and ñ -n are very small and hence, the decoding complexity of the outer code is negligible when compared to that of the LT part of the FR code. In addition, we are interested in a regime where r and t are very large and hence, the approximate complexity of the decoding algorithm is O(rt davgN K ). 242
A T 2 B 1 ? ? ? ? ?
(e) Iteration 4
First, we present simulation results to compare the recovery threshold of FR codes and 3-D Product codes. For the FR code, we set m = n = 80 and choose the entries of A and B to be realizations of an i.i.d. Gaussian random variable with zero mean and unit variance. We encode the source symbols A T i , for i ∈ [m] and j ∈ [n], using an (82, 80)×(82, 80) Product code, where an (82, 80) systematic RKRP code [15] is used as the row and column code. The output symbols from the encoder of the Product code are then encoded using a (10000, 6724) FLT code with right degree distribution Ω(x) = 0.013x + 0.5x 2 + 0.1661x 3 + 0.0726x 4 + 0.0826x 5 + 0.0581x 8 + 0.0340x 9 + 0.0576x 18 + 0.0160x 66 . As discussed in Section III-A.1, the number of source symbols is K = mn = 6400. For the chosen values of K = 6400 and P = 10000, the maximum fraction of straggling workers that any scheme can tolerate is given by P -K P = 0.36 [1]. Note that Polynomial codes achieve this threshold, and have zero probability of decoding failure for any α <= 0.36 where α = P -N P denotes the fraction of straggling workers. In our simulations, for each value of α being considered, the αP straggling workers were chosen uniformly at random.
To illustrate the superiority of the encoding scheme in Section III-A. the probability that the master node is unable to recover the output matrix C after receiving computations from αP fraction of workers. We have also plotted the decoding failure of a (21, 18) × (22, 19) × (22, 19) Product code in Fig. 3 for comparison. We observe that the encoding scheme described in Section III-A.1 (Scheme III) has better recovery threshold when compared to the other two encoding schemes (Schemes I and II). In addition, the recovery threshold of FR codes (Scheme III) is higher than that of 3-D Product codes. We have also plotted the error performance of the FR codes under inactivation decoding as described in Section III-B.3. One can see that the recovery threshold under inactivation decoding, instead of using only peeling decoding, is higher.
Next, for the chosen parameters above, we demonstrate the numerical stability of the proposed FR codes by simulating the average relative error defined as η ave := E ||C -C|| 2 /||C|| 2 , where C is the estimate of the output matrix C. In Fig. 4, we plot the average relative error versus fraction of straggling workers. While computing average relative error, we do not account for instances of decoding failure due to stopping sets. It can be seen that the average relative error of the simulated FR code is around 10 -9 and does not increase as the fraction of straggling workers increases. It is also observed that 3-D product codes and FR codes have comparable average relative error, see Fig. 4; whereas the FR codes have a better recovery threshold than 3-D product codes, see Fig. 3. In Fig. 4, it can also be seen that the average relative error of Polynomial codes is substantially higher than that of product codes and FR codes.