{'fraction_non_alphanumeric': 0.04199060408265081, 'fraction_numerical': 0.013740489751798113, 'mean_word_length': 4.24384129060388, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Experimental Evaluation of Multi-Round Matrix Multiplication on MapReduce

A common approach in the design of MapReduce algorithms is to minimize the number of rounds. Indeed, there are many examples in the literature of monolithic MapReduce algorithms, which are algorithms requiring just one or two rounds. However, we claim that the design of monolithic algorithms may not be the best approach in cloud systems. Indeed, multi-round algorithms may exploit some features of cloud platforms by suitably setting the round number according to the execution context. In this paper we carry out an experimental study of multi-round MapReduce algorithms aiming at investigating the performance of the multi-round approach. We use matrix multiplication as a case study. We first propose a scalable Hadoop library, named M$_3$, for matrix multiplication in the dense and sparse cases which allows to tradeoff round number with the amount of data shuffled in each round and the amount of memory required by reduce functions. Then, we present an extensive study of this library on an in-house cluster and on Amazon Web Services aiming at showing its performance and at comparing monolithic and multi-round approaches. The experiments show that, even without a low level optimization, it is possible to design multi-round algorithms with a small running time overhead.


I. INTRODUCTION
MapReduce is a computational paradigm 1 for processing large-scale data sets in a sequence of rounds executed on conglomerates of commodity servers [6]. This paradigm, and in particular its open source version Hadoop, has emerged as a de facto standard and has been widely adopted by a number of large Web companies (e.g., Google, Yahoo!, Amazon, Microsoft) and universities (e.g., CMU, Cornell, and the University of Maryland) [11]. MapReduce was initially introduced for log processing and web indexing but it has also been successfully used for other applications, including machine learning [14], data mining [2], scientific computing [21], and bioinformatics [22]. Informally, a MapReduce algorithm transforms an input multiset of key-value pairs into an output multiset of key-value pairs in a sequence of rounds. Each round consists of three steps: each input pair is individually transformed into a multiset of new pairs by a map function (map step); then the new pairs are grouped by key (shuffle step); finally, each group of pairs with the same key is processed, separately for each key, by a reduce function that produces the next new set of key-value pairs (reduce step).
The MapReduce paradigm has a functional flavor, in that it merely requires that the algorithm designer specifies the computation in terms of map and reduce functions. This enables parallelism without forcing an algorithm to cater for the explicit allocation of processing resources. Nevertheless, the paradigm implicitly posits the existence of an underlying unstructured and possibly heterogeneous parallel infrastructure, where the computation is eventually run. For these reasons, the MapReduce paradigm has gained popularity in recent years in cloud services for the development of large scale computations. Indeed, MapReduce is currently offered as a service by prominent cloud providers, including Amazon Web Services and Microsoft Azure [10].
As MapReduce is increasingly used for solving computational hungry problems on large datasets, it is crucial to design efficient and scalable algorithms. Many research efforts have been dedicated to capture efficiency in MapReduce algorithms (e.g., [2], [8], [11], [13], [17]). The major source of inefficiency is the communication requirements for transporting data from mappers to reducers in the shuffle step: since communication is a major factor determining the performance of algorithms on current computing systems, the amount of shuffled data should be minimized. Another issue that is usually taken into account for improving performance is the number of rounds of a MapReduce algorithm since the initial setups of a round and of the shuffle step are very costly operations. Thus, several MapReduce algorithms have been designed that require a very small number, usually one or two, of rounds (e.g., [2], [11]). We denote with monolithic algorithm a MapReduce algorithm requiring one or two rounds. A common approach for obtaining a monolithic algorithm is to decompose the problem into small subproblems which are executed concurrently in a single round, with each subproblem solved by a single application of the reduce function.
However, we claim that the design of monolithic algorithms may not be the best approach for long running MapReduce algorithms in cloud systems. Although it is true that the best performance is usually reached by reducing the round number, monolithic MapReduce algorithms may not exploit some features of cloud computing. Some examples follow. a) Service market: Some cloud providers offer a market where it is possible to bid on the cost of a service and to use it only when the actual price is below the bid. For reducing computing costs of long but low-priority computations, it would be desirable to develop MapReduce algorithms that can be stopped and restarted according to the price of the service. Unfortunately, current implementations of the MapReduce paradigm do not allow to stop a computation at an arbitrary point and then restart it. However, it is possible, like in Hadoop, to restart a computation from the beginning of the round that has been interrupted, losing the work that was already executed in that round. This clearly penalizes monolithic algorithms that consist of just a few long rounds: stopping a round and waiting for a better price might not be convenient since a significant part of previous work has to be executed again. It would be then interesting to develop MapReduce algorithms featuring a larger number of short rounds in order to reduce the amount of discarded work. b) Resource requirements: Some MapReduce algorithms need strong resource requirements in order to be monolithic. For instance, the multiplication of two dense matrices of size √ n× √ n in two rounds must [17] exchange approximately n n/m pairs during each shuffle, where m denotes an upper bound to the memory that can be used by a map/reduce function. This amount of data is linear in the input size only if m ∼ n, that is, when the multiplication can be almost solved sequentially. Similar requirements are also needed for joining relations [2] and enumerating triangles in a graph [16]. In a big-data era, the required local memory may exceed system resources and the huge amount of data in the shuffle step could arise issues related to network performance and system failure. Indeed, the network may be subject to congestion due to the large amount of data created within a small time interval, penalizing the scalability and fault tolerance of the network. It is then desirable to design algorithms that can tradeoff round number and resource requirements. We also observe that distributing a large computation among different rounds may help to checkpoint the computation and thus to restore it if the system completely fails or is heavily penalized by a fault.
It is then interesting to investigate multi-round algorithms where the round number can be set according to the execution context. However it is challenging to guarantee performance similar to the ones provided by monolithic algorithms. The first issue is that, by distributing the computation among different rounds, the total amount of communication in the shuffle steps should not increase with respect to the monolithic version since, as we have already mentioned, communication is the main bottleneck. Moreover, the latency required by each new round should be amortized by a sufficient amount of computation and communication performed within the round.
In this paper we carry out an experimentally study of multi-round MapReduce algorithms for matrix multiplication aiming at investigating the performance of the multi-round approach. Matrix multiplication is an important building block for many problems arising in different contexts, in particular scientific computing and graph processing [7]. Since an increasing number of computational problems is solved in the cloud and MapReduce is a common tool, it is then of interest to investigate efficient MapReduce solutions to matrix multiplication. Moreover, the communication and computation intensive nature of this problem strongly suggests a multiround approach.
The objective of this paper is therefore twofold. Firstly, it proposes a scalable Hadoop library, named M 3 (Matrix Multiplication in MapReduce), for matrix multiplication in the dense and sparse cases which allows to tradeoff round number with resource requirements. Secondly, it presents an extensive study of this library on an in-house cluster and on a cloud provider aiming at showing its performance but also at comparing monolithic and multi-round approaches. The experiments show that, even without a low level optimization, it is possible to exploit the benefit of multi-round algorithms with a small decrease of the performance.

A. Our results
The results provided in the paper are the following. 1) We propose the M 3 library for performing dense and sparse matrix multiplication in Hadoop. The library is based on the multi-round algorithms proposed in [17] which exploit a 3D decomposition of the problem. In this paper, we integrate the high level explanation in [17] by providing a detailed description, including the pseudocode of the algorithms, filling the gap between the theoretical results and the actual implementation. For the sake of completeness, we also present a MapReduce algorithm for dense matrix multiplication based on a 2D approach. The algorithms exhibit a tradeoff among round number, the amount of data in the shuffle steps and the amount of memory required by each reduce function. The library is publicly available at http://www.dei.unipd.it/m3.

2)
We carry out an extensive experimental evaluation of the M 3 library on an in-house cluster and on Amazon Elastic MapReduce (EMR), and show that the algorithms efficiently scale with increasing input size and processor number. We compare the performance of the multi-round approach with the monolithic one by measuring the performance of our algorithms when the round number increases. The experiments on the in-house cluster show that there is an average performance loss of 7% for each new round added by a multi-round approach with respect to the monolithic two-round approach (the loss is 18% on EMR). We give evidence that this factor is mainly due to the distributed file system HDFS used by Hadoop for storing pairs.
Matrix multiplication is an important building block for many problems and an experimental study of MapReduce algorithms is of independent interests, in particular considering recent efforts for perform scientific computing in the cloud. Moreover, we consider matrix multiplication as an interesting problem for studying multi-round algorithms since its high, but still tractable, communication and computation requirements allow to significantly load the system and to assess the performance under stress. We remark that the performance of our algorithms should not be compared with state-of-theart research on parallel sparse/dense matrix multiplication in High Performance Computing (e.g., [4]). Indeed, we are interested in investigating algorithm performance in cloud systems with MapReduce. In these settings, abstraction layers add a significant burden to algorithm performance. Research on the reduction of this load is crucial (see, e.g. [18], [23]), but it is out of the scope of this paper.

B. Previous Work
The MapReduce [6] paradigm has been widely studied and we refer to the book by Lin and Dyer [12] for a survey. The design of efficient MapReduce algorithms has been investigated from practical and theoretical perspectives. For instance, best practices in designing large-scale algorithms in MapReduce are proposed in [1], [13], while theoretical models for analyzing MapReduce are studied in [2], [8], [11], [17]. The majority of the algorithms requires just a few rounds, although monolithic solutions are not known for some problems (e.g., computing the diameter of a graph [5] or matrix inversion [17]).
Multi-round algorithms that trade round number with resource requirements have been proposed, with a theoretical approach, in [17] for some linear algebra problems, including matrix multiplication, in [8] for sorting and searching, and in [16] for triangle enumeration. To the best of our knowledge, the only experimental analysis of multi-round algorithms has been provided in [16] for triangle enumeration and shows that the performance of a multi-round approach is equivalent to a monolithic one.
Matrix multiplication is one of the most studied problems in the literature and we refer to the book by Golub and Van Loan [7] for further insights. However, to the best of our knowledge, this problem has been scarcely studied in MapReduce: ignoring naive implementations available on the web, the only scientific works are the aforementioned [17], and [19] which introduces HAMA, a library including MapReduce algorithms for matrix multiplication and conjugate gradient for matrices stored in the HBase database.
Comparison with previous work: This work completes the theoretical investigations in [17] by integrating the original description with more details that allow to fill the gap between theory and practice, and by proposing an Hadoop implementation. The M 3 library includes as special case the algorithms in HAMA [19]: the iterative approach proposed in the framework requires √ n rounds for multiplying a matrix of size √ n × √ n and is a special case of the algorithm based on the 2D approach described in Section III-A (i.e., it suffices to set m = √ n and ρ = 1). Since the HAMA library assumes a very different input representation and our experiments show that the aforementioned 2D approach is slower than the 3D approach, we do not take into account HAMA in this paper (however, further investigations will be carried out in the full version of this paper). Moreover, it should be said that the algorithms in [19] break the functional approach of MapReduce by allowing concurrent accesses to the distributed file system within each map and reduce function, and that the experiments have been carried out on small inputs (input side at most 5000, while we run our algorithms with input side 16000 and 32000).

C. Paper organization
The remaining sections of the paper are organized as follow. Section II introduces MapReduce, the matrix multiplication problem and the experimental settings. Section III describes and analyzes the proposed MapReduce algorithms. Section IV shows some additional technicalities required by the Hadoop implementation. Section V provides and explains the experimental results. Finally, Section VI gives some final remarks.

A. MapReduce and Hadoop
A MapReduce algorithm consists of a sequence of rounds. The input/output of a round is a multiset of key-value pairs k; v , where k and v denote the key and the value respectively. The input of a round may contain the output of the previous rounds but also other pairs. A round is organized in three steps: 1) Map step: each input pair is individually given in input to a map function, which returns a new multiset of pairs. The pairs returned by all applications of the map function are called intermediate pairs. 2) Shuffle step: the intermediate pairs are grouped by key. 3) Reduce step: each group of pairs associated with the same key is individually given to a reduce function, which returns a new multiset of pairs. The pairs returned by all applications of the reduce function are the output pairs of the round.
For simplicity we refer to a single application of the map (resp., reduce) function with mapper (resp., reducer). For a MapReduce algorithm, we denote with round number the number of rounds, with shuffle size the maximum number of intermediate pairs in each round, and with reducer size the maximum amount of memory words required by each reducer. 2 Hadoop 3 is the most used open source implementation of the MapReduce paradigm. A MapReduce algorithm is executed by Hadoop on p processors roughly as follows. In each round, each processor is associated with a fixed number of map tasks and reduce tasks. In the map step, the runtime system evenly distributes the input pairs to the map tasks and then each map task applies the map function to each single pair. In the reduce step, all groups of same-key pairs are randomly assigned to each reduce task and then each reduce task applies the reduce function to each single group. The input/output pairs are stored in the distributed file system HDFS [20], while intermediate pairs are temporary stored in the disk of each machine. Hadoop allows to personalize different settings of the framework, in particular it allows to specify how groups are distributed among reduce tasks by defining a partitioner. A partitioner receives as input a key k and the total number T of reduce tasks, and it returns a value in [0, T ) univocally denoting the reduce task that will apply the reduce function to the group associated with key k.
We recall that there exist other efficient implementations of the MapReduce paradigm, including Spark [23], MapReduce-MPI [18], Sphere [9]. Our algorithms can be executed on these frameworks with minor programming changes.

B. Matrix multiplication
The focus of the paper is matrix multiplication in a general semiring, that is we are ruling out Strassen-like algorithms. For the sake of simplicity, we focus on square matrices of size √ n × √ n. We denote with A and B the two input matrices and with C = A · B the output matrix.
In the sparse case, we denote with 0 ≤ δ ≤ 1 the density of non-zero entries in a given matrix of size √ n × √ n, that is the number of non-zero entries is δ n. A random Erdös-Rényi matrix of size √ n × √ n and parameter 0 ≤ δ ≤ 1 is a matrix where each entry is non-zero with probability δ ; the expected density is clearly δ . When δ << 1/n 1/4 , the expected density of the product of two Erdös-Rényi matrices is δ 2 √ n [3].

C. Experimental Framework
The experiments have been carried out on an in-house cluster and on Amazon Web Services. The in-house cluster consists of 16 nodes, each equipped with a 4-core Intel i7 processor Nehalem @ 3.07GHz, 12 GB of memory, 6 disks of 1TB and 7200 RPM in RAID0. The interconnection system is a 10 gigabit Ethernet network. The operating system is Debian/Squeeze with kernel 2.6.32. All experiments have been executed with Hadoop 2.4.0 on 16 nodes (1 master + 16 slaves; a node was both master and slave).
The experiments in the cloud have been executed with Amazon Elastic MapReduce (EMR) on instances of type c3.8xlarge and i2.xlarge in the US East Region. 4 The c3.8xlarge is a compute optimized instance featuring 32 virtual cores on a physical processor Intel Xeon E5-2680 @2.8GHz, 64 GB of memory, 2 solid state disks of 320 GB; instances are connected by a 10 gigabit network. The i2.xlarge is a storage optimized instance featuring 4 virtual cores on a physical processor Intel Xeon E5-2670 @2.5GHz, 32 GB of memory, 1 solid state disk of 800 GB optimized for very high random I/O performance; instances are connected by a moderate-performance network. All experiments have been run on clusters with 9 (1 master + 8 slave) instances of the same type, with a Debian/Squeeze operating system and Hadoop 2.2.0.

III. ALGORITHMS
The M 3 library currently contains two MapReduce algorithms for performing matrix multiplication which are based on the 3D approach proposed in [17]. The first algorithm multiplies two square dense matrices and is described in Section III-B. The second algorithm multiplies two sparse matrices and is given in Section III-C; although it applies to any sparse input, we analyze it in the case of random Erdös-Rényi matrices. The algorithms exhibit tradeoff among round number, shuffle size and reducer size, and rely on the wellknow 3D decomposition of the lattice representing the n 3/2 elementary products (some of which are zeros in the sparse case) into cubes of a given size. In other words, the problem is decomposed into multiplications of smaller square matrices that are solved sequentially by the reducers. The subproblems are evenly distributed among the round in order to guarantee the specified tradeoff.
For the sake of comparison, we also provide in Section III-A an algorithm for dense matrix multiplication which is based on the 2D decomposition of the lattice: the problem is decomposed into the multiplication of smaller rectangular matrices where the longest side is √ n. This approach is quite common in naive implementations of matrix multiplication in MapReduce available on the web: however, as we will see, this approach is inefficient. The algorithm exhibits a tradeoff among round number, shuffle size, and reducer size. This result does not appear in [17].
The tradeoffs in our algorithms are highlighted by expressing the round number as a function of two parameters: the replication factor ρ and the subproblem size m. The replication factor gives an estimate of the desired volume of intermediate 4 http://aws.amazon.com/ec2/instance-types. Algorithm 1: Map and reduce functions of the r-th round of the 2D algorithm, with 0 ≤ r < n/(ρm).
data in each round, while the subproblem size bounds the memory requirements of each reducer: indeed, our algorithms guarantee that the shuffle size is Θ (ρ) times the input size and the reducer size is Θ (m).

A. 2D Algorithm for Dense Matrix Multiplication
The 2D MapReduce algorithm divides the input matrix A into n/m submatrices of size m/ √ n × √ n, the input matrix B into n/m submatrices of size √ n × m/ √ n, and the output matrix into (n/m) 2 submatrices of size m/ √ n × m/ √ n. We label each submatrix of A and B with A i and B i and each submatrix of C with C i, j , for any 0 ≤ i, j < n/m. We clearly have C i, j = A i · B j . The input matrix A is stored as a collection of pairs (i, −1); A i , for any 0 ≤ i < n/m, where −1 is used as dummy value; similarly for matrix B. The output matrix is stored as a collection of pairs (i, j);C i, j , for any 0 ≤ i, j < n/m.

For any
√ n ≤ m ≤ n and 1 ≤ ρ ≤ n/m, the algorithm requires R = n/(ρm) rounds, the shuffle size is 2ρn, and the reducer size is 3m. For the sake of simplicity, we assume that √ n is an integer and that √ n divides m. In the r-th round, with 0 ≤ r < R, the algorithm computes the products C i, j = A i · B j with j = i + + rρ mod (n/m), for any 0 ≤ i < n/m and 0 ≤ < ρ. Specifically, at the beginning of the r-th round, the input pairs are (i, −1); A i , (−1, j); B j . Then, the map functions emits for each pair (i, −1); A i the pairs (i, j); A i with j = (i + + rρ) mod (n/m) for any 0 ≤ < ρ; similarly, for each pair (−1, j); B j , the mapper emits the pairs (i, j); B j with i = ( j − − rρ) mod (n/m) for any 0 ≤ < ρ. Then, each reducer associated with key (i, j) receives A i and B j and emits the pair (i, j); A i · B j . The pseudocode of the map and reduce functions is in Algorithm 1.
Theorem 1. The above algorithm requires R = n/(ρm) rounds, the shuffle size is 2ρn, and the reducer size is 3m.
Proof: Subproblem (i, j) is executed with = ( j − i) mod ρ in round r = ( j − i − )/ρ. The reducer associated with key (i, j) correctly receives matrices A i and B j in round r since the mapper with input (i, j); A i emits the pair (i, i + + rρ); A i where i + + rρ = j for the above values of and r (a similar argument applies for B j ). The shuffle size is 2ρn since there are at most ρ copies of each submatrix A i and B j . Each reducer just requires at most 3m memory words for computing C i, j = A i · B j .
Our algorithm guarantees that each mapper performs the same amount of work: indeed, each submatrix of A and B is replicated exactly ρ times. We observe that a naive way to distribute subproblems among rounds may significantly unbalance the work performed by mappers: in the worst case, there can be one submatrix of A replicated n/m times, while the remaining submatrices of A are replicated once (similarly for B); then one mapper requires O (n/m) work, while all the remaining mappers perform O (1) work.

B. 3D Algorithm for Dense Matrix Multiplication
The 3D MapReduce algorithm for dense matrix multiplication was initially proposed in [17]. The paper provides an high level description of the algorithm in the MR computational model, however it does not provide important details such as a description of the map and reduce functions that are required for actually implementing the algorithm in MapReduce. In this section, we close this gap by providing a detailed description of the algorithm, including a pseudocode for the map and reduce functions.
For any 1 ≤ m ≤ n and 1 ≤ ρ ≤ n/m, the algorithm requires R = √ n/(ρ √ m) + 1 rounds, the shuffle size is 3ρn, and the reducer size is 3m. For the sake of simplicity, we assume that √ n and √ m are integers and that √ m divides √ n. The input matrices A, B and the output matrix C are divided into submatrices of size √ m × √ m, and we denote these submatrices with A i, j , B i, j and C i, j , respectively, for 0 ≤ i, j < n/m. Clearly, we have C i, j = ∑ √ n/m−1 h=0 A i,h · B h, j . The input matrix A is stored as a collection of pairs (i, −1, j); A i, j , where −1 is used as dummy value. Matrices B and C are stored similarly. The above division of the input matrices implies (n/m) 3/2 products between submatrices, which are partitioned into n/m groups as follows: group G , with 0 ≤ < n/m, contains the products A i,h · B h, j for h = (i + j + ) mod n/m and for every 0 ≤ i, j < n/m. We note that each submatrix of A and B appears exactly once in each group.
The algorithm works in R rounds. In the r-th round, with 0 ≤ r < R − 1, the algorithm computes all products in G and the sum C i, j = ∑ r−1 k=0 A i,i+ j+ +kρ B i+ j+ +kρ, j for every rρ ≤ < (r + 1)ρ and 0 ≤ i, j < n/m; each product A i,h · B h, j is computed within the reducer associated with key (i, h, j). Eventually, in the last round r = R − 1, the matrix C is created by computing C i, j = ∑ ρ−1 =0 C i, j for every 0 ≤ i, j < n/m. A more detailed explanation follows. At the beginning of the rth round, the input pairs are (i, −1, j); A i, j , (i, −1, j); B i, j . Starting from the second round (i.e., r ≥ 1) the algorithm also receives as input the pairs (i, , j);C i, j , for every 0 ≤ i, j < n/m and 0 ≤ < ρ, that have been outputted by the previous round. In the map step, each input pair is replicated ρ times since it is required in ρ reducers. Specifically, for each pair (i, −1, j); A i, j , the map function emits the pairs (i, j, h); A i, j with h = j − i − mod n/m for every 0 ≤ < ρ. Similarly, for each pair (i, −1, j); B i, j , the map function emits the pairs (h, i, j); B i, j with h = i − j − mod n/m for every 0 ≤ < ρ. On the other hand, the map function Algorithm 2: Map and reduce functions of the r-th round of the 3D dense algorithm, with 0 ≤ r < √ n/(ρ √ m) + 1.
emit (h, i, j); D 10: case D is C i, j : 11: if r is the last round then emits for each pair (i, , j);C i, j the pair (i, h, j);C i, j with h = i+ j + +rρ mod n/m. Finally, each reducer associated with key (i, h, j) with h = i + j + , for any 0 ≤ i, j < n/m and 0 ≤ < ρ, receives as input A i,h , B h, j and C i, j , computes C i, j = C i, j + A i,h B h, j and emits (i, , j),C i, j . The pseudocode of the map and reduce functions is described in Algorithm 2.
Theorem 2. The above algorithm requires R = √ n/(ρ √ m)+1 rounds, the shuffle size is 3ρn, and the reducer size is 3m.
Proof: The correctness of the algorithm is proved in [17]. However, we have to show that each reducer receives the correct set of submatrices since pair distribution is not taken into account in the original paper. Each product A i,k · B k, j is computed when = (k − i − j) mod ρ and in round r = (k − i − j − )/ρ. In this round, the mapper associated with pair (i, k); A i,k emits the pair (i, k, k − i − − rρ); A i,k where k − i − − rρ = j for the above values of and r. A similarly argument applies to B h, j and C i, j .
Since each submatrix of A and B is replicated ρ times and there are ρn/m matrices C i, j , the shuffle size is 3ρn. Each reducer requires at most 3m memory words for computing C i, j = C i, j + A i,i+ j+ +rρ · B i+ j+ +rρ, j .
We observe that there are 3n n/m pairs in each shuffle in a two-round algorithm for matrix multiplication (i.e, ρ = n/m), and it is possible to get a one-round only if the problem is solved sequentially by a single reducer (i.e., m = n). As shown in [17], this is the best we can get if only semiring operations are allowed. Finally, we note that the algorithm guarantees that almost all mappers perform the same amount of work because each submatrix of A and B is replicated ρ times.

C. 3D Algorithm for Sparse Matrix Multiplication
An algorithm for sparse matrix multiplication easily follows from the previous 3D algorithm by exploiting the density of the input and output matrices. We assume that the input matrices are random Erdös-Rényi matrices with size √ n × √ n and expected density δ << n 1/4 ; we will later see how to extend the result to the general case. We recall that the output density of the product of two random Erdös-Rényi matrices is The input matrices A and B and the output matrix C are partitioned into blocks of size . Each submatrix is represented as a list of non-zero entries, although other formats can be used like the Compressed Row Storage. The sparse algorithm follows by applying the 3D dense algorithm in Section III-B to the submatrices of size √ m × √ m . The sparse algorithm exploits the fact that each submatrix A i, j or B i, j (resp., C i, j ) is a random Erdös-Rényi matrix with expected density δ (resp., δ O ). Then the expected space required for computing A i,h · B h, j is (2δ + δ 2 √ n)m ∼ 3m, and thus the reducer size is 3m even for this algorithm. The pseudocode can be easily derived from the pseudocode for the dense case in Algorithm 2 by replacing m with m . Theorem 3. The above algorithm requires R = δ n 3/4 /(ρ √ m) + 1 rounds, the expected shuffle size is 3ρδ 2 n 3/2 , and the expected reducer size is 3m.
Proof: The claim easily follows from Theorem 2 and by the sparsity of the input/output matrices.
Suppose that A and B are general sparse matrices with density δ << 1 and that an approximation of the density size of the output matrix is knownδ O . Let δ M = max{δ ,δ 0 }, then it is sufficient to split each input/output matrix into blocks of size √ m × √ m with m = m/δ M and apply the previous sparse algorithm. For improving the load balancing among reducers columns and rows of the input matrices should be randomly permuted [17]. In this case it is easy to show that the algorithm requires R = √ δ M n/(ρ √ m) + 1 rounds, the expected shuffle size is 3ρδ M n, and the expected reducer size is 3m. A good approximation of the output matrices can be computed with a scan of the input matrices (see, e.g., [15]). Finally, we observe that if the output size is not known, then it suffices to set m = m/δ but the bounds on shuffle and reducer sizes do not apply anymore.

IV. IMPLEMENTATION
The algorithms described in previous sections are implemented as Hadoop jobs, one for each round. Each job comprises a single map and a single reduce function. Matrices are represented as SequenceFiles where keys are triplets or pairs (for 3D and 2D multiplication, respectively) and values are serialized objects representing blocks. The representation used for blocks will vary, depending on the matrix being dense or sparse. For dense matrices, each serialized block is the sequence of elements of the matrix in row-major order. As for sparse matrices, only non-zero elements are serialized along with their indexes. The entries of the matrices are doubles.

A. Map and reduce implementation
Map functions are straightforward implementations of the ones described in Algorithms 1 and 2. Reduce functions, while being adherent to the specification given in pseudocode, present a few caveats from both a correctness and a performance perspective.
Reduce functions in Hadoop take input values as instances of the Iterable interface, giving Hadoop developers the freedom of changing the actual underlying implementation without breaking client code. If one needs all the input values at the same time to perform the reduce computation (as is our case), then a deep copy of each value must be saved in local variables or collections to ensure correctness. Simply storing the reference does not work, since the implementation of the Iterable interface used by Hadoop returns a mutable object for each next call 5 . Failing to perform a deep copy of the input will make the reducer use the same block three times instead of the proper blocks from matrices A, B, and C. Note that by performing this copy we incur an unavoidable performance penalty.
As for the performance in general, it is crucial to use a fast linear algebra library for local multiplications. For the dense case, we use the JBLAS library 6 , that provides Java bindings to the fast BLAS/LAPACK routines. Unfortunately, JBLAS does not support sparse matrices, for which we use MTJ 7 . Even though this is the fastest library for sparse matrix multiplication among the tested ones, it is orders of magnitude slower than JBLAS on inputs of the same size thus preventing us to perform complete experiments on sparse matrices.

B. Hadoop configuration
As mentioned in Section II, we used both a in-house cluster and Amazon Web Services (AWS) to carry out our experiments. Our cluster has been configured as follows. HDFS replication has been set to one, that is, redundancy has been disabled. We found that the default replication, 3, degraded performance of ≈ 5% while not providing benefits in our experimental setting. Each machine of the cluster runs two mappers and two reducers, each with 3GB of memory, in order to allow us to use matrix blocks with many elements.
As for the configuration Hadoop instances provided by AWS, we used the default configuration, customized by Amazon for each machine instance type. The rationale is that AWS provides Hadoop as-a-service, enabling users to concentrate on their application rather than the configuration and management of an Hadoop cluster.

C. Partitioner
The keys of the dense 3D algorithms proposed in the Section III-B are triplets (i, j, k). A partitioner receives as input a key and the number T of reduce tasks in the system and 5 At least in the most recent version, 2.4.0 6 http://jblas.org/ 7 https://github.com/fommil/matrix-toolkits-java/ returns the identifier t ∈ [0, T ) of the reduce tasks that will be responsible of the key. When the key is a triplet (i, j, k), a simple partitioner for these keys is t = (31 2 i+31 j +k) mod T .
However, as Figure 1 shows, this function is not able to evenly distribute reduces among the T reduce tasks.
We propose an alternative partitioner. In the r-th round there are ρ(n/m) reduces denoted by keys (i, h, j) with 0 ≤ i, j < n/m and (i + j + rρ) mod n/m ≤ h < (i + j + (r + 1)ρ) mod n/m. The partitioner univocally maps each key (i, h, j) in the range [0, ρn/m) by setting z = iρn/m + jρ + h with h = h mod ρ. We note that the mapping consists of a row-major ordering on the first, third and second coordinates, where the second coordinate has been adjusted to be in the range [0, ρ). Keys mapped in [0, T ρn/(mT ) ) are evenly distributed among the reduce tasks, while the remaining at most T − 1 keys are randomly distributed among the tasks. The pseudocode is given in Algorithm 3. This partitioner is also used in the 3D sparse algorithm and a slightly different approach is used for the 2D dense algorithm.

A. Dense matrices
We start our experimental analysis by evaluating the dense algorithm described in Section III-B.
1) Impact of subproblem size m: The total amount of shuffled data in all rounds is O n n/m . This quantity is independent of the replication factor ρ, but increases as m decreases. Since the performance of an algorithm strongly relies on the total amount of shuffled data, it is reasonable to select a large m and then decompose the problem into larger subproblems. On the other hand, a large value of m reduces parallelism since there are (n/m) 3/2 independent subproblems. Moreover, it increases the load on each machine since a matrix multiplication of size √ m × √ m is computed by each reducer with Θ m 3/2 work and Θ (m) memory, which may exceed hardware limits. This fact is highlighted in Figure 2(a) which shows the running time with input matrix side √ n ∈ {16000, 32000} and different values of the subproblem size √ m ∈ {1000, 2000, 4000}; the experiments have been carried out with ρ = n/m and with ρ = 1, that is with a monolithic approach and an extreme case of multi-round. All experiments show that the performance improves with larger values of m, but the gain decreases for larger values: with input side 32000 and ρ = n/m, the performance gain when √ m moves from 1000 to 2000 is 1.99, while it is 1.12 when it moves from 2000 to 4000. We observe that the monolithic approach is slightly less sensible to variations of m since it can exploit much more parallelism in each round than a multi-round approach. The value √ m = 8000 is missing in the experiments since all executions failed due to an out-of-memory error. Although each submatrix requires only 488MB, the actual memory requirements are much more due to the internal structure of Hadoop, which may keep in the memory of each processor the submatrices outputted by map tasks as well as the submatrices required in input by the reduce tasks. Advanced optimizations can increase the maximum value of m, but do not significantly change the proposed results.
Similar experiments have been done on EMR that show a similar behavior on c3.8xlarge and i2.xlarge instances, with √ m = 4000 being the optimal choice. Interestingly, smaller instance types require √ m = 2000 for avoiding memory errors. All the following experiments assume √ m = 4000.
2) Impact of replication factor ρ: We are ready to study the performance of a multi-round approach. We recall that, for any replication factor 1 ≤ ρ ≤ n/m, the number of rounds is √ n/(ρ √ m) + 1. Figures 2(b) and 2(e) investigate the running times with different values of the replication factor and input size √ n = 16000 and √ n = 32000, respectively. For ρ = n/m, we get a monolithic approach (i.e., two rounds). Without any surprise, the best running time is reached by the monolithic approach, therefore supporting the common practice to minimize round number. However, we observe that a multi-round approach increases the running time with respect to the monolithic two-round approach by an average factor 7% for each additional round. We believe that the main reason of this increase is due to the distributed file system HDFS, which is optimized for writing and reading large files. When ρ = n/m, each reduce task writes two large chunks of output pairs (one per round) on the HDFS containing all outputs of the associated reducers in a given round. On the other hand, for values of ρ smaller than n/m, each reduce task writes a larger number of chunks of smaller size. Although the total size of the written files is the same in each approach, the system is not able to exploit the features of HDFS with a multiround approach. We conjecture that this gap can be closed by improving the access to the distributed file system. (a) Time vs subproblem size with √ n ∈ {16000, 32000}. The label max (resp., min) means that it has been used maximum (resp., minimum) replication, that is ρ = n/m (resp., ρ = 1).   √ n ∈ {16000, 32000}, respectively. Running times significantly increase with respect to experiments on the inhouse cluster: the running times with √ n = 16000 (resp., √ n = 32000) on EMR are on the average 4.7 (resp., 1.4) times larger than the ones on the cluster, even if the computational resources of the two systems are somewhat similar. It is interesting to note that the gap decreases with larger input sizes. The average performance loss with respect to the monolithic two-round approach is 17% for each additional round.
Finally, we observe that the algorithm efficiently scales with the input size. The running time, for any replication factor, increases on the in-house cluster by a factor ∼ 8 when the input side doubles, which agrees with the cubic (in the matrix side) complexity of the algorithm. However, the factor is ∼ 5 on EMR: the smaller factor is due to high fixed costs which are not efficiently amortized with small inputs.
3) Cost analysis: We now analyze how communication, computation and infrastructure affect the running time. We define the following three costs: a) Infrastructure cost T in f r : it is the time for setting up the Hadoop system and each round. It is given by the time of the algorithm when no computation is done in the reduce functions and the value of each pair is replaced by an empty matrix (i.e., the amount of shuffle data is negligible). b) Computation cost T comp : it is the time for performing the local computation. LetT comp be the time of the algorithm when any local multiplication is performed on two locally created matrices (the cost of the creation is negligible) and the value of each pair is replaced by an empty matrix. Then, the computation cost isT comp − T in f r . c) Communication cost T comm : it is the time for exchanging pairs between mappers and reducers, and includes the costs of reading/writing data on HDFS and of shuffle steps. LetT comm be the time of the algorithm when no computation is done in the reduce function (the output pair is a copy of an input pair). Then, the communication cost isT comm − T in f r .
We think that the infrastructure cost is a good approximation of the true time for setting up the system and each round. However, this is not true for the communication and computation costs since the procedure does not take into account the concurrency between the three components. Nevertheless, these costs sill provide an interesting overview of the main factors determining the running time of an algorithm.     (f) Comparison between the 2D and 3D approaches with √ n = 16000. The replication factor is ρ ∈ {1, 2, 4} for the 3D approach, and ρ ∈ {1, 2, 4, 8, 16} for the 2D approach. nents for √ n = 16000 and ρ ∈ {1, 2, 4}, and for √ n = 32000 and ρ ∈ {1, 2, 4, 8}, respectively. The infrastructure cost increases linearly with the round number, and the average fixed cost of a round is 17 seconds. The computation cost for a given input size is independent of the replication value, which confirms that the work is evenly distributed among cluster nodes. The communication cost increases when ρ increases as already noted for Figures 2(b) and 2(e).
The cost analysis on EMR with c3.8xlarge instances is reported in Figures 3(b) and 3(e) for √ n = 16000 and ρ ∈ {1, 2, 4}, and for √ n = 32000 and ρ ∈ {1, 2, 4, 8}, respectively. It should be noted that the average computation cost is similar to the respective component in the in-house cluster, although there is a larger variance due to the unpredictable load of the physical machines of EMR. The average infrastructure cost is 30 seconds.
Finally, we report the cost analysis with √ n = 16000 and ρ ∈ {1, 2, 4} on i2.xlarge instances. A i2.xlarge instance has faster disks optimized for random I/Os but slower network than a c3.8xlarge instance. However, we observe that the communication costs are smaller than the respective costs in Figure 3(b) for c3.8xlarge instances. This fact supports the claim in Section V-A2 that the main bottleneck is the inability of HDFS to efficiently read/write small chunks of data.

4) Scalability:
We now study the scalability of the algorithm by running the dense algorithm on a smaller number of nodes. For this experiment, we run the dense algorithm with input size 16000 on 4, 8 and 16 nodes of the in-house cluster. The results are given in Figure 2(d) with different replication factor ρ ∈ {1, 2, 4}. The plot shows that, for any replication factor, the algorithm scales efficiently although there is a small reduction in the speed up with 16 nodes. This may be due to a loss in data locality and to a larger cost of the shuffle step when the algorithm is run on 16 nodes.

5) 2D vs 3D algorithms:
We now move to analyze the performance of the 3D multiplication strategy compared with the 2D strategy. Figure 3(f) clearly shows that the 3D approach has a significant performance advantage. This is due to the fact that the total shuffle size is, for the 3D approach, O n n/m , whereas for the 2D approach it is O n 2 /m . Since the major bottleneck in MapReduce is communication, the 2D approach incurs a significant penalty.

B. Sparse matrices
In this section, we investigate the performance of the 3D sparse algorithm and show that it improves performance by exploiting input sparsity. Figure V with different replication factors of the 3D sparse algorithm with two input Erdös-Rényi matrices with √ n ∈ {2 20 , 2 22 , 2 24 }, and an average of 8 non-zero entries per row and per column (i.e., δ ∈ {1/2 17 , 1/2 19 , 1/2 21 }, respectively). Since the output matrices have expected density δ O ∈ {1/2 14 , 1/2 16 , 1/2 18 } (see Section II-B), the input matrices are partitioned into submatrices of size 2 18 , 2 19 and 2 20 , respectively, so that the expected number of non-zero elements in each submatrix of the output is comparable with the subproblem size of the dense case. This shows that, by exploiting the sparseness of the input matrix, we can tackle much bigger problems, under the same memory constraints. As mentioned in Section IV, we avoided the computation of local products, because of the lack of an efficient Java implementation of sequential sparse matrix multiplication, that will be the subject of future work. However, as we have experimentally shown, the time spent computing local products is a fixed additive factor, whereas the communication is the dominant component, thus these experiments clearly show the tradeoffs between number of rounds and time for sparse matrix multiplication.

VI. CONCLUSION
In this paper we have proposed the Hadoop library M 3 for performing dense and sparse matrix multiplication in MapReduce by exploiting the theoretical results in [17]. We have carried out an extensive experimental study on an inhouse cluster and on a cloud provider which offers MapReduce as a service (Amazon Elastic MapReduce). The results show that M 3 efficiently scales with processor number and input size. The proposed algorithms exhibit a tradeoff among round number, the amount of shuffled data and the amount of space required by each reducer. We exploit this tradeoff for assessing the overhead of a multi-round approach with respect to a monolithic approach. We have shown that the running time increases by an average 7% in an in-house cluster (18% on Amazon) for each additional round with respect to the monolithic two-round approach. As we have observed, the main source of inefficiency as the round number increases is the distributed file system HDFS where input/output pairs are stored. As future work, we plan to further investigate sparse matrix multiplication on MapReduce and to test our algorithms on other implementations of the MapReduce paradigm, in particular Spark, where the management of the input/output pairs of each round is more efficient.

ACKNOWLEDGMENTS
This paper was supported in part by: MIUR of Italy under project AMANDA; University of Padova under projects CPDA121378 and AACSE; Amazon in Education Grant; equipment donation by Samsung. Silvestri is also grateful to Paolo Rodeghiero for initial discussions on Hadoop.