Oguz Kaya
email: oguz.kaya@ens-lyon.fr

Bora Uçar
email: bora.ucar@ens-lyon.fr

Scalable Sparse Tensor Decompositions in Distributed Memory Systems [Technical Paper]

Keywords: G.1.0 [Numerical Analysis]: General-Numerical algorithms, Parallel algorithms, G.2.2 [Discrete Mathematics]: Graph Theory-Hypergraphs, G.4 [Mathematical Software]: Algorithm design and analysis, Parallel and vector implementations

We investigate an efficient parallelization of the most common iterative sparse tensor decomposition algorithms on distributed memory systems. A key operation in each iteration of these algorithms is the matricized tensor times Khatri-Rao product (MTTKRP). This operation amounts to element-wise vector multiplication and reduction depending on the sparsity of the tensor. We investigate a fine and a coarse-grain task definition for this operation, and propose hypergraph partitioning-based methods for these task definitions to achieve the load balance as well as reduce the communication requirements. We also design a distributed memory sparse tensor library, HyperTensor, which implements a well-known algorithm for the CANDECOMP-/PARAFAC (CP) tensor decomposition using the task definitions and the associated partitioning methods. We use this library to test the proposed implementation of MTTKRP in CP decomposition context, and report scalability results up to 1024 MPI ranks. We observed up to 194 fold speedups using 512 MPI processes on a well-known real world data, and significantly better performance results with respect to a state of the art implementation.

INTRODUCTION

Tensors or multi-dimensional arrays arise in many fields, including analysis of Web graphs [START_REF] Kolda | The TOPHITS model for higher-order web link analysis[END_REF], knowledge bases [START_REF] Carlson | Toward an architecture for never-ending language learning[END_REF], product reviews at online stores [START_REF] Bennett | The netflix prize[END_REF], chemometrics [START_REF] Appellof | Strategies for analyzing data from video fluorometric monitoring of liquid chromatographic effluents[END_REF], signal processing [START_REF] Lathauwer | From matrix to tensor: Multilinear algebra and signal processing[END_REF], computer vision [START_REF] Vasilescu | Multilinear analysis of image ensembles: Tensorfaces[END_REF], and more. Tensor decomposition algorithms are used as an important tool to understand the tensors and glean hidden or latent information SC '15, November 15-20, 2015, Austin, TX, USA ACM ISBN 978-1-4503-2138-9. DOI: http://dx.doi.org/10.1145/2807591.2807624 from them. Considerable efforts are being put in designing numerical algorithms for different tensor decomposition problems (see a short [START_REF] Comon | Tensors: A brief introduction[END_REF] and a long survey [START_REF] Kolda | Tensor decompositions and applications[END_REF]), and algorithmic and software contributions go hand in hand with these efforts [START_REF] Andersson | The N-way toolbox for MATLAB[END_REF][START_REF] Bader | Matlab tensor toolbox version 2.6[END_REF][START_REF] Choi | DFacTo: Distributed factorization of tensors[END_REF][START_REF] Kang | GigaTensor: Scaling tensor analysis up by 100 times -Algorithms and discoveries[END_REF][START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF].

One of the most common tensor decomposition approaches is the CANDECOMP/PARAFAC (CP) formulation, which approximates a given tensor as a sum of rank-one tensors. Two most common methods for computing a CP decomposition are (i) CP-ALS [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via an N-way generalization of "Eckart-Young" decomposition[END_REF][START_REF] Harshman | Foundations of the PARAFAC procedure: Models and conditions for an "explanatory" multi-modal factor analysis[END_REF], which is based on the alternating least squares method; and (ii) CP-OPT [START_REF] Acar | A scalable optimization approach for fitting canonical tensor decompositions[END_REF], which is based on the gradient descent method. Both of these methods are iterative, where the computational core of each iteration is a special operation called the matricized tensor times Khatri-Rao product (MTTKRP). When the input tensor is sparse and N dimensional, the MTTKRP operation amounts to element-wise multiplication of N -1 vectors and scaled reduction of those products according to the sparsity structure of the tensor. This computationally involved operation has received recent interest for efficient execution in different settings such as Matlab [START_REF] Andersson | The N-way toolbox for MATLAB[END_REF][START_REF] Bader | Matlab tensor toolbox version 2.6[END_REF], MapReduce [START_REF] Kang | GigaTensor: Scaling tensor analysis up by 100 times -Algorithms and discoveries[END_REF], shared memory [START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF], and distributed memory [START_REF] Choi | DFacTo: Distributed factorization of tensors[END_REF].

We investigate an efficient parallelization of the MTTKRP operation in distributed memory environments for sparse tensors in the context of the CP decomposition methods CP-ALS and CP-OPT. For this purpose, we formulate two task definitions, a coarse-grain and a fine-grain one. These definitions are given by applying the owner-computes rule to a coarse and a fine-grain partition of the tensor nonzeros. We define the coarse-grain partition of a tensor as a partition of one of its dimensions. In matrix terms, a coarse-grain partition corresponds to a row-wise or a column-wise partition. Two very recent parallel algorithms DFacTo [START_REF] Choi | DFacTo: Distributed factorization of tensors[END_REF] and SPLATT [START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF] have coarse-grain tensor partitions and hence have coarse-grain tasks. We define the fine-grain partition of a tensor as a partition of its nonzeros. This has the same significance in matrix terms. Based on these two task granularities, we present two parallel algorithms for the MTTKRP operation. We address the computational load balance and communication cost reduction problems for the two algorithms, and present hypergraph partitioning-based models to tackle these problems with off-the-shelf partitioning tools.

The MTTKRP operation also arises in the close variants of CP-ALS and CP-OPT for some other tensor decomposition methods [START_REF] Liavas | Parallel algorithms for constrained tensor factorization via the alternating direction method of multipliers[END_REF]; hence, it has been implemented as a standalone routine [START_REF] Bader | Efficient MATLAB computations with sparse and factored tensors[END_REF] to enable algorithm development. Once this operation is efficiently done, the other parts of the de-composition algorithms are usually straightforward. For this reason, most of the related work on high performance tensor decomposition algorithms focus on this particular operation. As hinted above, the majority of our contribution is also on the efficiency of the MTTKRP operation. Nonetheless, we design a library for the parallel CP-ALS algorithm to test the proposed MTTKRP algorithms in a suitable context and give experimental results using this library.

The organization of the rest of the paper is as follows. In the next section, we give the notation, describe the MT-TKRP operation and CP-ALS method, and review the related work that influenced our efforts. Hypergraph theoretical definitions are also given in this section. In Section 3, we describe the coarse and fine-grain MTTKRP algorithms, analyze their efficient parallelization requirements, and present hypergraph models for reducing the parallelization overhead. Section 4 contains experimental results, where we report speedups and perform comparisons with a state of the art distributed memory CP-ALS implementation.

BACKGROUND AND NOTATION

Bold, upper case Roman letters are used for matrices as in A. Matrix elements are shown with the corresponding lowercase letters as in ai,j. Matrix sizes are sometimes shown in the lower right corner, e.g., AI×J . Matlab notation is used to refer to the entire rows and columns of a matrix, e.g., A(i, :) and A(:, j) refer to the ith row and jth column of A respectively.

Tensors

We use calligraphic font to refer to tensors, e.g., X. The order of a tensor is the number of its dimensions, which we denote with N . For the sake of simplicity of the notation and the discussion, we describe all the notation and the algorithms for N = 3, even though our algorithms and implementations have no such restriction. We explicitly generalize the discussion to general order-N tensors whenever we find necessary. As in matrices, an element of a tensor is denoted by a lowercase letter and subscripts corresponding to the indices of the element, e.g., the element (i, j, k) of a third-order tensor is x i,j,k . A fiber in a tensor is defined by fixing every index but one, e.g., if X is a third-order tensor, X :,j,k is a mode-1 fiber and Xi,j,: is a mode-3 fiber. A slice in a tensor is defined by fixing only one index, e.g., Xi,:,:, refers to the ith slice of X in mode 1. We use |Xi,:,:| to denote the number of nonzeros in Xi,:,:.

Tensors can be matricized in any mode. This is achieved by identifying a subset of the modes of a given tensor X as the rows and the other modes of X as the columns of a matrix and appropriately mapping the elements of X to those of the resulting matrix. We will be exclusively dealing with the matricizations of tensors along a single mode. For example, take X ∈ R I 1 ו••×I N . Then X (1) denotes the mode-1 matricization of X in such a way that the rows of X (1) corresponds to the first mode of X and the columns corresponds to the remaining modes. The tensor element xi 1 ,...,i N corresponds to the element i1, 1 + N j=2 (ij -1) j-1 k=1 I k of X [START_REF] Acar | A scalable optimization approach for fitting canonical tensor decompositions[END_REF] . Specifically, each column of the matrix X (1) becomes a mode-1 fiber of the tensor X. Matricizations in the other modes are defined similarly.

Given two matrices AI 1 ×J 1 and BI 2 ×J 2 , the Kronecker product is defined as

A ⊗ B =    a1,1B • • • a1,J 1 B . . . . . . . . . aI 1 ,1B • • • aI 1 ,J 1 B    .
For AI 1 ×J and BI 2 ×J , the Khatri-Rao product is defined as

A B = A(:, 1) ⊗ B(:, 1) • • • A(:, J) ⊗ B(:, J) , which is of size I1I2 × J.
For AI×J and BI×J , the Hadamard product is defined as

A * B =    a1,1b1,1 • • • a1,J b1,J . . . . . . . . . aI,1bI,1 • • • aI,J bI,J    .
The CP-decomposition of rank R (or with R components) of a given tensor X factorizes X into a sum of R rank-one tensors. For X ∈ R I×J×K , it yields to x i,j,k ≈ R r=1 airbjrc kr and X ≈ R r=1 ar •br •cr, for ar ∈ R I , br ∈ R J and cr ∈ R K , where • is the outer product of the vectors. Here the matrices A = [a1, . . . , aR], B = [b1, . . . , bR], and C = [c1, . . . , cR] are called the factor matrices, or factors. For N -mode tensors, we use U1, . . . , UN to refer to the factor matrices.

We are now equipped with the notation to present the Alternating Least Squares (ALS) method for obtaining a rank-R approximation of a tensor X with the CP-decomposition. A common formulation of CP-ALS is shown in Algorithm 1 for the 3rd order tensors. At each iteration, each factor matrix is recomputed while fixing the other two; e.g., A ← X

(1) (C B)(B T B * C T C) † .
This operation is performed in the following order: MA = X (1) (C B), V = (B T B * C T C) † , and then A ← MAV. Here V is a dense matrix of size R×R and is easy to compute. The important issue is the efficient computation of the MTTKRP operations yielding MA, similarly MB = X (2) (A C) and MC = X (3) (B A).

The sheer size of the Khatri-Rao products makes them impossible to compute explicitly; hence, efficient MTTKRP algorithms find other means to carry out the MTTKRP operation (see the next subsection).

Algorithm 1: CP-ALS for the 3rd order tensors Input : X: A 3rd order tensor R: The rank of approximation

Output: CP decomposition [[λ; A, B, C]] repeat A ← X (1) (C B)(B T B * C T C) † Normalize columns of A B ← X (2) (C A)(A T A * C T C) † Normalize columns of B C ← X (3) (B A)(A T A * B T B) †
Normalize columns of C and store the norms as λ until no improvement or maximum iterations reached

Related work

SPLATT [START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF] is an efficient implementation of the MT-TKRP operation for sparse tensors on shared memory systems. It is our understanding that the codes are implemented for the 3-mode tensors-there is no experimental results with the higher order tensors. Their discussion includes the generalization of the techniques to higher order tensors whenever relevant. SPLATT implements the MT-TKRP operation based on the slices of the dimension in which the factor is updated, e.g., on the mode-1 slices when computing A ← X (1) (C B)(B T B * C T C) † . Nonzeros of the fibers in a slice are multiplied with the corresponding rows of B and the results are accumulated to be later scaled with the corresponding row of C to compute the row of A corresponding to the slice. Parallelization is done using OpenMP directives, and the load balance (in terms of the number of nonzeros in the slices of the mode for which MTTKRP is computed) is achieved by using the dynamic scheduling policy. Hypergraph models are used to optimize cache performance by reducing the number of times a row of A, B, and C are accessed. Smith et al. [START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF] also use N -partite graphs (where N is the order of the tensor) to reorder the tensors for all dimensions. Experiments are conducted on an HP ProLiant BL280c G6 server with dual 8core E5-2670 Xeon processors running at 2.6GHz. Smith et al. implement a sparse-tensor vector product algorithm called TVec for carrying out the MTTKRP operation and report speedups with respect to this algorithm. They report 3.7x speedup on the serial execution and 29.8x speedup on 16-way parallel execution of the SPLATT with respect to the TVec.

GigaTensor [START_REF] Kang | GigaTensor: Scaling tensor analysis up by 100 times -Algorithms and discoveries[END_REF] is an implementation of CP-ALS which follows the Map-Reduce paradigm. All important steps (the MTTKRPs and the computations of the B T B * C T C) are performed using this paradigm. A distinct advantage of Gi-gaTensor is that thanks to Map-Reduce, the issues of faulttolerance, load balance, and out of core tensor data are automatically handled. On a real world data, speedup studies with upto 100 machines (each machine has 2 quad-core Intel 2.83 GHz CPUs) are presented, where the speedup with 100 machines is 1.4 times the speedup with 25 machines. One iteration of CP-ALS as implemented in GigaTensor takes more than 10 3 seconds for a random tensor of size 10 5 ×10 5 ×10 5 with 10 5 /50 nonzeros on 35 machines, eventually reaching between 10 4 and 10 5 seconds for 10 9 ×10 9 ×10 9 with 10 9 /50 nonzeros. The presentation [START_REF] Kang | GigaTensor: Scaling tensor analysis up by 100 times -Algorithms and discoveries[END_REF] of GigaTensor focuses on three-mode tensors and expresses the map and the reduce functions for this case. To the best of our understanding, additional map and reduce functions are needed for higher order tensors, which would incur overheads.

DFacTo [START_REF] Choi | DFacTo: Distributed factorization of tensors[END_REF] is a distributed memory implementation of the MTTKRP operation. It performs two successive sparse matrix-vector multiplies (SpMV) to compute a column of the product X (1) (C B). A crucial observation (made also elsewhere [START_REF] Kolda | Tensor decompositions and applications[END_REF]) is that this operation can be implemented as X T

(2) B(:, r), which can be reshaped into a matrix to be multiplied with C(:, r) to form the rth column of the MTTKRP. Although SpMV is a well-investigated operation, there is a peculiarity here: the result of the first SpMV forms the values of the sparse matrix used in the second one. Therefore, there are sophisticated data dependencies between the two SpMVs. Notice that DFacTo is rich in SpMV operations: there are two SpMVs per factorization rank per dimension of the input tensor. DFacTo needs to store the tensor matricized in all dimensions, i.e., X (1) , . . . , X (N ) . In low dimensions, this can be a slight memory overhead; yet in higher dimensions the overhead could be non-negligible. DFacTo uses MPI for parallelization yet fully store the fac-tor matrices in all MPI ranks. The rows of X (1) , . . . , X (N ) are blockwise distributed (statically). With this partition, each process computes the corresponding rows of the factor matrices. Finally, DFacTo performs an MPI Allgatherv operation to communicate the new results to all processes, which results in (In/P ) log 2 P communication volume per process (assuming a hypercube algorithm) when computing the n th factor matrix having In rows using P processes. Experiments are presented on machines equipped with two 2.1 GHz 12-core AMD 6172 processors where up to 32 machines are used. In sequential runs, DFacTo is shown to be 5 times faster than GigaTensor and 10 times faster than a MATLAB implementation [START_REF] Bader | Efficient MATLAB computations with sparse and factored tensors[END_REF]. On a real world data set, DFacTo obtains about 3.5 speedup on 32 machines with respect to an execution on four machines.

Tensor Toolbox [START_REF] Bader | Matlab tensor toolbox version 2.6[END_REF] is a MATLAB toolbox for handling tensors. It provides many essential operations and enables fast and efficient realizations of complex algorithms in MAT-LAB for sparse tensors [START_REF] Bader | Efficient MATLAB computations with sparse and factored tensors[END_REF]. Among those operations, MT-TKRP implementations are provided and used in CP-ALS method. Here, each column of the output is computed by performing N -1 sparse tensor vector multiplications. Another well-known MATLAB toolbox is the N -way toolbox [START_REF] Andersson | The N-way toolbox for MATLAB[END_REF] which is essentially for dense tensors and incorporates now support for sparse tensors [START_REF] Acar | A scalable optimization approach for fitting canonical tensor decompositions[END_REF] through Tensor Toolbox. Tensor Toolbox and the related software provide excellent means for rapid prototyping of algorithms and also efficient programs for tensor operations that can be handled within MATLAB.

Hypergraphs and hypergraph partitioning

A hypergraph H = (V, E) is defined as a set of vertices V and a set of hyperedges E. Each hyperedge is a set of vertices. The vertices of a hypergraph can be associated with weights denoted by w[•], and the hyperedges can be associated with costs denoted by c[•]. For a given integer K ≥ 2, a K-way vertex partition of a hypergraph H = (V, E) is denoted as Π = {V1, . . . , VK }, where the parts are non-empty, mutually exclusive,

V k ∩ V = ∅ for k = ; and collectively exhaustive, V = V k . Let W k = v∈V k be the total weight in V k and Wavg = v∈V w[v]/K be the average part weight. If each part V k ∈ Π satisfies the balance criterion W k ≤ Wavg(1 + ε), for k = 1, 2, . . . , K (1)
we say that Π is balanced where ε represents the maximum allowed imbalance ratio. In a partition Π, a hyperedge that has at least one vertex in a part is said to connect that part. The number of parts connected by a hyperedge h, i.e., connectivity, is denoted as λ h . Given a vertex partition Π of a hypergraph H = (V, E), one can measure the size of the cut induced by Π as

χ(Π) = h∈E c[h](λ h -1) . (2)
This cut measure is called the connectivity-1 cutsize metric. Given ε > 0 and an integer K > 1, the standard hypergraph partitioning problem is defined as the task of finding a balanced partition Π with K parts such that χ(Π) is minimized. The hypergraph partitioning problem is NPhard [START_REF] Lengauer | Combinatorial Algorithms for Integrated Circuit Layout[END_REF].

A recent variant of the above problem is the multi-constraint hypergraph partitioning [START_REF] Çatalyürek | Hypergraph Models for Sparse Matrix Partitioning and Reordering[END_REF][START_REF] Karypis | Multilevel algorithms for multi-constraint hypergraph partitioning[END_REF]. In this variant, each vertex has an associated vector of weights. The partitioning objective is the same as above, and the partitioning constraint is to satisfy a balancing constraint for each weight. Let w[v, i] denote the C weights of a vertex v for i = 1, . . . , C. In this variant, the balance criterion (1) is rewritten as

W k,i ≤ Wavg,i (1+ε) for k = 1, . . . , K and i = 1, . . . , C (3)
where the ith weight W k,i of a part V k is defined as the sum of the ith weights of the vertices in that part (i.e., W k,i = v∈V k w[v, i]), Wavg,i is the average part weight for the ith weight of all vertices (i.e., Wavg,i = v∈V w[v, i]/K), and ε again represents the allowed imbalance ratio.

PARALLELIZATION

A common approach in implementing the MTTKRP is to explicitly matricize a tensor across all modes, and then perform the Khatri-Rao product using the matricized tensors [START_REF] Choi | DFacTo: Distributed factorization of tensors[END_REF][START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF]. Matricizing a tensor in a mode i requires column index values up to N k =i I k , which can exceed the integer value limits supported by modern architectures when using tensors of higher order and very large dimensions. Also, matricizing across all modes results in N replications of a tensor; which can exceed the memory limitations. Hence, in order to be able to handle large tensors we store them in coordinate format for the MTTKRP operation, which is also the method of choice in Tensor Toolbox [START_REF] Bader | Matlab tensor toolbox version 2.6[END_REF].

With a tensor stored in the coordinate format, MTTKRP operation can be performed as shown in Algorithm 2. As seen on Line 1 of this algorithm, a row of B and a row of C are retrieved, and their Hadamard product is computed and scaled with a tensor entry to update a row of MA. In general, for an N -mode tensor

MU 1 (i1, :) ← MU 1 (i1, :)+xi 1 ,i 2 ...,i N [U2(i2, :) * • • • * UN (iN , :)]
: M A = X (1) (B C) Initialize M A to zeros of size I A × R foreach x i,j,k ∈ X do 1 M A (i, :) ← M A (i, :) + x i,j,k [B(j, :) * C(k, :)]
As factor matrices are accessed row-wise, we define computational units in terms of the rows of factor matrices. It follows naturally to partition all factor matrices row-wise and use the same partition for the MTTKRP operation for each mode of an input tensor across all CP-ALS iterations to prevent extra communication. A crucial issue is the task definitions, as this pertains to the issues of load balancing and communication. We identify a coarse-grain and a finegrain task definition for this computational kernel.

In the coarse-grain task definition, ith atomic task consists of computing the row MA(i, :) using the nonzeros in the tensor slice Xi,:,: and the rows of B and C corresponding to the nonzeros in that slice. The input tensor X does not change throughout the iterations of tensor decomposition algorithms; hence it is viable to make the whole slice Xi,:,: available to the process holding MA(i, :) so that the MTTKRP operation can be performed by only communicating the rows of B and C. Yet, as CP-ALS requires the MTTKRP in all modes, and each nonzero x i,j,k belongs to slices X(i, :, :), X(:, j, :), and X(:, :, k), we need to replicate tensor entries in the owner processes of these slices. This may require up to N times replication of the tensor, depending on its partitioning. Note that an explicit matricization always requires exactly N replications of tensor entries.

In the fine-grain task definition, an atomic task corresponds to the multiplication of a tensor entry with the Hadamard product of the corresponding rows of B and C. Here, tensor nonzeros are partitioned among processes with no replication to induce a task partition by following the owner-computes rule. This necessitates communicating the rows of B and C that are needed by these atomic tasks. Furthermore, partial results on the rows of MA need to be communicated, as without duplicating tensor entries, we cannot in general compute all contributions to a row of MA. Here, the partition of X should be useful in all modes, as the CP-ALS method requires the MTTKRP in all modes.

The coarse-grain task definition resembles to the one-dimensional (1D) row-wise (or column-wise) partitioning of sparse matrices, whereas the fine-grain one resembles to the twodimensional (nonzero-based) partitioning of sparse matrices for parallel sparse matrix-vector multiply (SpMV) operations. As is confirmed for SpMV in modern applications, 1D partitioning usually leads to harder problems of load balancing and communication cost reduction. The same phenomenon is likely to be observed in tensors as well. Nonetheless, we cover the coarse-grain task definition, as it is used in the state of the art parallel MTTKRP [START_REF] Choi | DFacTo: Distributed factorization of tensors[END_REF][START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF] methods, which partition the matricized tensor row-wise (or equivalently, partition the input tensor by slices).

Coarse-grain task model

In the coarse-grain task model, computing the rows of MA, MB, and MC are defined as the atomic tasks which are partitioned across all processes. Let µA denote the partition of the first mode's indices among the processes, i.e., µA(i) = p, if the process p is responsible for computing MA(i, :). Similarly, let µB and µC define the partition of the second and the third mode indices. The process owning MA(i, :) needs the entire tensor slice Xi,:,:; similarly, the process owning MB(j) needs X:,j,: and the owner of MC (k, :) needs X :,:,k . This necessitates duplication of some tensor nonzeros to prevent unnecessary communication.

One needs to take the context of CP-ALS into account when parallelizing the MTTKRP method. First, the output MA of MTTKRP is transformed into A. Since A(i, :) is computed simply by multiplying MA(i, :) with the matrix (B T B * C T C) † , we make the process which owns MA(i, :) responsible for computing A(i, :). Second, N MTTKRP operations follow one another in an iteration. Assuming that every process has the required rows of the factor matrices while executing MTTKRP for the first mode, it is advisable to implement the MTTKRP in such a way that its output MA, after transformed into A, is communicated. This way, all processes would have the necessary data for executing the MTTKRP for the next mode. With these in mind, the coarse-grain parallel MTTKRP method executes Algorithm 3 at each process p. Jp and Kp correspond to the unique second and third mode indices in X Ip,:,: On exit : A(Ip, :) is computed, its rows are sent to processes needing them;

A T A is available Initialize M A (Ip, :) to all zeros of size |Ip| × R. foreach i ∈ Ip do foreach x i,j,k ∈ X i,:,: do 1 M A (i, :) ← M A (i, :) + x i,j,k [B(j, :) * C(k, :)] 2 A(Ip, :) ← M A (Ip, :)(B T B * C T C) † foreach i ∈ Ip do 3
Send A(i, :) to all processes having nonzeros in X i,:,: .

Receive A(i, :) from µ A (i) for each owned nonzero x i,j,k .

As seen in Algorithm 3, the process p computes MA(i, :) for all i with µA(i) = p on Line 1. This is possible without any communication, if we assume the preconditions that B T B, C T C, and the required rows of B and C are available at each process. We need to satisfy this precondition for the MTTKRP operations in the remaining modes. Once MA(Ip, :) is computed, it is multiplied on Line 2 with (B T B * C T C) † to obtain A(Ip, :). Then, the process p sends the rows of A(Ip, :) to other processes who will need them (messages destined to the same process are combined into a single message), then receive the ones that it will need. More precisely, the process p sends A(i, :) to the process r, where µB(j) = r or µC (k) = r for a nonzero x i,j,k ∈ X, thereby satisfying the stated precondition for the remaining MTTKRP operations. Since computing A(i, :) necessitates that B T B and C T C are available at each process, we also need to compute A T A and make the result available to all processes. We perform this by computing local multiplications A(Ip, :) T A(Ip, :) and all-reducing the partial results on Line 4.

We now investigate the problems of obtaining load balance and reducing the communication cost in an iteration of CP-ALS given in Algorithm 3; that is, when Algorithm 3 is applied N times, once for each mode. Line 1 is the computational core; the process p performs i∈Ip |Xi,:,:| many Hadamard products (of vectors of size R) without any communication. In order to achieve load balance, one should partition the slices of the first mode equitably by taking the number of nonzeros into account. Line 2 does not involve communication, where load balance can be achieved by assigning almost equal number of slices to the processes. Line 4 requires almost equal number of slices to the processes for load balance. The operations at Lines 2 and 4 can be performed very efficiently using BLAS3 routines; therefore, their costs are generally negligible in compare to the cost of the sparse irregular computations at Line 1. There are two lines, 3 and 4, requiring communication. Line 4 requires the collective all-reduce communication, and one can rely on efficient MPI implementations. The communication at Line 3 is irregular and would be costly. Therefore, the problems of computational load balance and communication cost need to be addressed with regards to Lines 1 and 3, respectively.

Let ai be the task of computing MA(i, :) at Line 1. Let also TA be the set of tasks ai for i = 1, . . . , IA, and bj, c k , TB, and TC be defined similarly. In the CP-ALS algorithm, there are dependencies among the triplets ai, bj, and c k of the tasks for each nonzero x i,j,k . Specifically, the task ai depends on the tasks bj and c k . Similarly, bj depends on the tasks ai and c k , and c k depends on the tasks ai and bj. These dependencies are the source of the communication at Line 3; e.g., A(i, :) should be sent to the processes that own B(j, :) and C(k, :) for the each nonzero x i,j,k in the slice Xi,:,:. These dependencies can be modeled using graphs or hypergraphs by identifying the tasks with vertices, and their dependencies with edges and hyperedges, respectively. As is well-known in similar parallelization contexts [START_REF] Çatalyürek | Hypergraph-partitioning-based decomposition for parallel sparse-matrix vector multiplication[END_REF][START_REF] Hendrickson | Load balancing fictions, falsehoods and fallacies[END_REF][START_REF] Hendrickson | Graph partitioning models for parallel computing[END_REF], the standard graph partition related metric of "edge cut" would loosely relate to the communication cost. We therefore propose a hypergraph model for modeling the communication and computational requirements of the MTTKRP operation in the context of CP-ALS.

For a given tensor X, we build a hypergraph H = (V, E) with the vertex set V and hyperedge set E as follows. The vertex set V corresponds to the set of tasks. We abuse the notation and use ai to refer to a task and its corresponding vertex in V ; the meaning should be clear from the context. We can then set V = TA ∪ TB ∪ TC . Since one needs load balance on Line 1 for each mode, we use vectors of size N for vertex weights. We set w[ai, 1] = |Xi,:,:|, w[bj, 2] = |X:,j,:|, and w[c k , 3] = |X :,:,k | as the vertex weights in the indicated modes, and assign weight 0 in all other modes. The dependencies causing communication at Line 3 are oneto-many: for a nonzero x i,j,k , any one of ai, bj, and c k is computed by using the data associated with the other two. Following the guidelines in building the hypergraph models [START_REF] Uçar | Revisiting hypergraph models for sparse matrix partitioning[END_REF], we add a hyperedge corresponding to each task (or each row of the factor matrices) in TA ∪ TB ∪ TC , in order to model the data dependencies to that specific task. We use n a i to denote the hyperedge corresponding to the task ai, and define n b i and n c k similarly. For each nonzero x i,j,k , we add the vertices ai, bj, and c k to the hyperedges n a i , n b j , and n c k . Let Ji and Ki be all unique 2nd and 3rd mode indices, respectively, of the nonzeros in Xi,:,:. Then, the hyperedge n a i contains the vertex ai, all vertices bj for j ∈ Ji, and all vertices c k for k ∈ Ki.

Figure 1a demonstrates the hypergraph model for the coarse-grain task decomposition of a sample 3 × 3 × 3 tensor with nonzeros (1, 2, 3), (2, 3, 1), and (3, 1, 2). Labels for the hyperedges, shown as small blue circles, are not provided; yet it should be clear from the context, e.g., ai ∈ n a i is shown with a non-curved black line between the vertex and the hyperedge. For the first, second, and third nonzeros, we add the corresponding vertices to related hyperedges, which are shown with green, magenta, and orange curves, respectively. It is interesting to note that the N -partite graph of Smith et al. [START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF] and this hypergraph model are related; their adjacency structure when laid as a matrix give the same matrix.

Consider now a P -way partition of the vertices of H = (V, E) and the association of each part with a unique process for obtaining a P -way parallel CP-ALS. The task ai incurs |Xi,:,:| Hadamard products in the first mode in the process µA(i), and has no other work when computing the MTTKRPs MB = X (2) (A C) and MC = X (3) (A B) in the other modes. This observation (combined with the corresponding one for bjs and c k s) shows that any partition satisfying the balance condition for each weight component (3) corresponds to a balanced partition of Hadamard products (Line 1) in each mode. Now consider the messages, which are of size R, sent by the process p regarding A(i, :). These messages are needed by the processes that have slices in modes 2 and 3 that have nonzeros whose first mode index is i. These processes correspond to the parts connected by the hyperedge n a i . For example, in Figure 1a due to the nonzero (3, 1, 2), a3 has a dependency to b1 and c2, which are modeled by the curves from a3 to n b 1 and n c 2 . Since a3, b1, and c2 reside in different parts, a3 contributes one to the connectivity of n b 1 and n c 2 , which accurately captures the total communication volume by increasing the total connectivity-1 metric by two.

Notice that the part µA(i) is always connected by n a iassuming that the slice Xi,:,: is not empty. Therefore, minimizing the connectivity-1 metric (2) exactly amounts to minimizing the total communication volume required for MT-TKRP in all computational phases of a CP-ALS iteration, if we set c[•] = R for all hyperedges.

In the gradient-descent based algorithm CP-OPT [START_REF] Acar | A scalable optimization approach for fitting canonical tensor decompositions[END_REF], all factors are updated simultaneously based on their values in the previous iteration; hence, one can use the same hypergraph with a slight modification in vertex weights. In this case, we need to use a single-constraint hypergraph partitioning by setting w

Fine-grain task model

In the fine-grain task model we assume a partition of the nonzeros of the tensor; hence there is no duplication. Let µ X denote the partition of the tensor nonzeros, e.g., µ X (i, j, k) = p if the process p holds x i,j,k . Whenever µ X (i, j, k) = p, the process p performs the Hadamard product of B(j, :) and C(k, :) in the first mode, A(i, :) and C(k, :) in the second mode, and A(i, :) and B(j, :) in the third mode; then scales the result with x i,j,k . Let again µA, µB, and µC denote the partition on the indices of the corresponding modes; that is, µA(i) = p denotes that the process p is responsible for computing MA(i, :). As in the coarse-grain model, we take the CP-ALS context into account while designing the MT-TKRP algorithm: (i) the process which is responsible for MA(i, :) is also held responsible for A(i, :); (ii) at the beginning, all rows of the factor matrices B(j, :) and C(k, :) where µ X (i, j, k) = p are available to the process p; (iii) at the end, all processes get the R × R matrix A T A; and (iv) each process has all the rows of A that are required for the upcoming MTTKRP operations in the second or the third modes. With these in mind, the fine-grain parallel MTTKRP method executes Algorithm 4 at each process p. Jp and Kp correspond to the unique second and third mode indices in Xp On exit : A(Ip, :) is computed, its rows are sent to processes needing them; A(Fp, :) is updated; and

A T A is available Initialize M A (Fp, :) to all zeros of size |Fp| × R foreach x i,j,k ∈ Xp do 1 M A (i, :) ← M A (i, :) + x i,j,k [B(j, :) * C(k, :)] foreach i ∈ Fp \ Ip do 2 Send M A (i, :) to µ A (i)
3 Receive contributions to M A (i, :) for i ∈ Ip and add them up

4 A(Ip, :) ← M A (Ip, :)(B T B * C T C) † foreach i ∈ Ip do 5
Send A(i, :) to all processes having nonzeros in X i,:,: .

foreach i ∈ Fp \ Ip do 6
Receive A(i, :) from µ A (i)

7 Locally compute A(Ip, :) T A(Ip, :) and all-reduce the results to form A T A

In Algorithm 4, Xp and Ip denote the set of nonzeros and the set of first mode indices assigned to the process p. As seen in the algorithm, the process p has the required data to carry out all Hadamard products corresponding to all x i,j,k ∈ Xp on Line 1. After scaling by x i,j,k , the Hadamard products are locally accumulated in MA, which contains a row for all i where Xp has a nonzero on the slice Xi,:,:. Notice that a process generates a partial result for each slice of the first mode in which it has a nonzero. The relevant row indices (the set of unique first mode indices in Xp) are denoted by Fp. Some indices in Fp are not owned by the process p. For such an i ∈ Fp \ Ip, the process p sends its contribution to MA(i, ; ) to the owner process µA(i) at Line 2. Again, messages of the process p destined to the same process are combined. Then, the process p receives contributions to MA(i, :) where µA(i) = p at Line 3. Each such contribution is accumulated, and the new A(i, :) is computed by a multiplication with (B T B * C T C) † at Line 4. Finally, the updated A is communicated to satisfy the precondition of the upcoming MTTKRP operations. Note that the communications at Lines 5 and 6 are the duals of those at Lines 3 and 2, respectively, in the sense that the directions of the messages are reversed and the messages always pertain to the same row indices. For example, MA(i, :) is sent to the process µA(i) at Line 2 andA(i, :) is received from µA(i) at Line 6. Hence, the process p generates partial results for each slice Xi,:,: on which it has a nonzero; either keeps it or sends it to the owner of µA(i) (at Line 2), and if so, receives the updated row A(i, :) later (at Line 6). Finally, the algorithm finishes by an all-reduce operation to make A T A available to all processes for the upcoming MTTKRP operations.

We now investigate the computational and communication requirements of Algorithm 4. As before, the computational core is at Line 1, where the Hadamard products of the row vectors of B(j, :) and C(k, :) are computed for each nonzero x i,j,k that a process owns, and the result is added to MA(i, :). Therefore, each nonzero x i,j,k incurs three Hadamard products (one per mode) in its owner process µ X (i, j, k) in an iteration of the CP-ALS algorithm. Other computations are either efficiently handled with BLAS3 routines (Lines 4 and 7), or they (Line 3) are not as many as the number of Hadamard products. The significant communication operations are the sends at Lines 2 and 5, and the corresponding receives at Lines 3 and 6. Again, the allreduce operation at Line 7 would be handled efficiently by the MPI implementation. Therefore, the problems of computational load balance and reduced communication cost need to be addressed by partitioning the tensor nonzeros equitably among processes (for load balance at Line 1) while trying to confine all slices of all modes to a small set of processes to reduce the communication (to reduce communication at Lines 2 and 6).

We propose a hypergraph model for the computational load and communication volume of the fine-grain parallel MTTKRP operation in the context of CP-ALS. For a given tensor, the hypergraph H = (V, E) with the vertex set V and hyperedge set E is defined as follows. Since we need to partition the nonzeros and the indices along each mode, V has four types of vertices: ai, for i = 1, . . . , IA; bj, for j = 1, . . . , IB; c k for k = 1, . . . , IC ; and the vertex x i,j,k we define for each nonzero in X by abusing the notation. The vertex x i,j,k represents the Hadamard products B(j, :) * C(k, :), A(i, :) * C(k, :), and A(i, :) * B(j, :) to be performed during the MTTKRP operations at different modes. There are three types of hyperedges in E, and they represent the dependencies of the Hadamard products to the rows of the factor matrices. The hyperedges are defined as follows: E = EA ∪ EB ∪ EC where EA contains a hyperedge n a i for each A(i, :); EB contains a hyperedge n b j for each B(j, :); and EC contains a hyperedge n c k for each C(k, :). Initially, n a i , n b j and n c k contain only the corresponding vertices ai, bj, and c k . Then, for each nonzero x i,j,k ∈ X, we add the x i,j,k to the hyperedges n a i , n b j and n c k to model the data dependency to the corresponding rows of A, B, and C. As in the previous subsection, one needs a multi-constraint formulation for three reasons. First, since x i,j,k vertices represent the Hadamard products, they should be evenly distributed among the processes. Second, as the owners of the rows of A, B, and C are possible destinations of partial results (at Line 2), it is advisable to spread them evenly. Third, almost equal number of rows per process implies balanced memory usage and BLAS3 operations at Lines 4 and 7; even though this operation is not our main concern for the computational load. With these constraints, we associate a weight vector of size N + 1 with each vertex. For a vertex x i,j,k we set w[x i,j,k ] = [0, 0, 0, 3], for vertices ai, bj and c k , we set

w[ai] = [1, 0, 0, 0], w[bj] = [0, 1, 0, 0], w[c k ] = [0, 0, 1, 0].
In Figure 1b, we demonstrate the fine-grain hypergraph model for the sample tensor of the previous section. We exclude the vertex weights from the figure. Similar to the coarse-grain model, dependency to the rows of A, B and C are modeled by connecting each vertex x i,j,k to the hyperedges n a i , n b j , and n c k to capture the communication cost of the parallel fine-grain MTTKRP.

Consider now a P -way partition of the vertices of H = (V, E) and associate each part with a unique process to obtain P -way parallel CP-ALS with the fine-grain MTTKRP operation. The partition of the x i,j,k vertices uniquely partition the Hadamard products. Satisfying the fourth balance constraint, as in (3), balances the computational loads of the processes. Similarly, satisfying the first three constraints leads to an equitable partition of the rows of factor matrices among processes. Let µ X (i, j, k) = p and µA(i) = ; that is, the vertex x i,j,k and ai are in different parts. Then, the process p sends a partial result on MA(i, :) to the process at Line 2 of Algorithm 4. By looking at these messages carefully, we see that all processes whose corresponding parts are in the connectivity set of n i A will send a message to µA(i). Since µA(i) is also in this set, we see that the total volume of send operations concerning MA(i, :) at Line 2 is equal to λ n a i -1. Therefore, the connectivity-1 cut size metric (2) over the hyperedges in EA encodes the total volume of messages sent at Line 2, if we set c[•] = R. Since the send operations at Line 5 are duals of the send operations at Line 2, total volume of messages sent at Line 5 for the first mode is also equal to this number. By extending this reasoning to all hyperedges, we see that the cumulative (over all modes) volume of communication is equal to the connectivity-1 cut-size metric [START_REF] Andersson | The N-way toolbox for MATLAB[END_REF].

The presence of multiple constraints usually causes difficulties for the current state of the art partitioning tools (Aykanat et al. [START_REF] Aykanat | Multi-level direct K-way hypergraph partitioning with multiple constraints and fixed vertices[END_REF] discuss the issue for PaToH [START_REF] Çatalyürek | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF]). In preliminary experiments, we observed that this was the case for us as well. In an attempt to circumvent this issue, we designed the following alternative. We start with the same hypergraph, and remove the vertices corresponding to the rows of factor matrices. Then use a single constraint partitioning on the resulting hypergraph to partition the vertices of nonzeros. We are then faced with the problem of assigning the rows of factor matrices, given the partition on the tensor nonzeros. Any objective function (relevant to our parallelization context) in trying to solve this problem would likely be a variant of the NP-complete Bin-Packing Problem [17, p.124], as in the SpMV case [START_REF] Bisseling | Communication balancing in parallel sparse matrix-vector multiplication[END_REF][START_REF] Uçar | Encapsulating multiple communication-cost metrics in partitioning sparse rectangular matrices for parallel matrix-vector multiplies[END_REF]. Therefore, we heuristically handle this problem, by making the following observations: (i) the modes can be partitioned one at a time; (ii) each row corresponds to a removed vertex for which there is a corresponding hyperedge; (iii) the partition of x i,j,k vertices define a connectivity set for each hyperedge. This is essentially multi-constraint partitioning which takes advantage of the listed observations. We use connectivity of corresponding hyperedges as key value of rows to impose an order. Then, rows are visited in decreasing order of key values. Each visited row is assigned to the least loaded part in the connectivity set of its corresponding hyperedge, if that part is not overloaded. Otherwise, the visited vertex is assigned to the least loaded part at that moment. Note that in the second case we add one to the total communication volume; otherwise the total communication volume obtained by partitioning x i,j,k vertices remains the same. Note also that key values correspond to the volume of receives at Line 3 of Algorithm 4; hence, this method also aims to balance the volume of receives at Line 3, and the volume of dual send operations at Line 5.

Discussion

The pros and cons of the coarse and fine-grain models resemble to those of the 1D and 2D partitionings in the SpMV context. There are two advantages of the coarse-grain model over the fine-grain one. First, there is only one communication/computation step in the coarse-grain model; whereas fine-grain model requires two computation and communication phases. Second, the coarse-grain hypergraph model is smaller than the fine-grain hypergraph model (one vertex per index of each dimension versus additionally having one vertex per nonzero). However, one major limitation with the coarse-grain parallel decomposition is that restricting all nonzeros in an (N -1) dimensional slice of an N dimensional tensor to reside in the same process poses a significant constraint towards communication minimization; an (N -1)dimensional data typically has very large "surface" which necessitates gathering numerous rows of the factor matrices. As N increases, one might expect this phenomenon to increase the communication volume even further. Also, if the tensor X is not large enough in one its modes, this poses a granularity problem which potentially causes load imbalance and limits the parallelism. None of these issues exists in the fine-grain task decomposition.

EXPERIMENTS

We conducted our experiments on the IDRIS Ada cluster, which consists of 304 nodes each having 128 GBs of memory and four 8-core Intel Sandy Bridge E5-4650 processors running at 2.7 GHz. We ran our experiments up to 1024 cores, and assigned one MPI process per core. The cluster allowed runs up to 2048 cores; but none of the methods we experimented with scaled to that point. All codes we used in our benchmarks were compiled using the Intel C++ compiler (version 14.0.1.106) with -O3 option for compiler optimizations, -axAVX,SSE4.2 to enable SSE optimizations whenever possible, and -mkl option to use the Intel MKL library (version 11.1.1) for LAPACK, BLAS and VML (Vector Math Library) routines. We also compiled DFacTo code using the Eigen (version 3.2.4) library [START_REF] Guennebaud | [END_REF] with EIGEN USE MKL ALL flag set to ensure that it utilizes the routines provided in the Intel MKL library for the best performance. We used PaToH [START_REF] Çatalyürek | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF] (version 3.2) with default options to partition the hypergraphs that we formed. We created all partitions offline, and ran our experiments on these partitioned tensors on the cluster. We do not report timings for partitioning hypergraphs, which is quite costly, for two reasons. First, in most applications (see the provenance of our data set below), the tensors from the real-world data are built incrementally; hence, a partition for the updated tensor can be formed by refining the partition of the previous tensor. Also, one can decompose a tensor multiple times with different ranks of approximation. Thereby, the cost of an initial partitioning of a tensor can be amortized across multiple runs. Second, we had to partition the data on a system different from the one used for CP-ALS runs. It would not be very useful to compare the runtimes from different systems and repeating the same runs on different system would not add much to the presented results. We also note that the proposed coarse and fine-grain MTTKRP formulations are independent from the partitioning method. Our dataset for experiments consists of four sparse tensors, whose sizes are given in Table 1. The first tensor has user × movie × time dimensions, and is formed from the Netflix Prize competition [START_REF] Bennett | The netflix prize[END_REF]. In this tensor, the nonzeros correspond to the user reviews of movies, and review date extends the data to the third dimension. The values of the nonzeros are determined by the corresponding review scores given by the users. The tensor NELL-B comes from the Never Ending Language Learning (NELL) knowledge database of the "Read the Web" project [START_REF] Carlson | Toward an architecture for never-ending language learning[END_REF], which consists of tuples of the form (entity, relation, entity) such as ('Chopin','plays musical instrument','piano'). The nonzeros of this tensor correspond to these entries discovered by NELL from the web, and the values are set to be the "belief" scores given by the learning algorithms used in NELL. Delicious and Flickr are the datasets for the web-crawl of Delicious.com and Flickr.com during 2006 and 2007, which is formed by Görlitz et al. [START_REF] Görlitz | Pints: Peer-to-peer infrastructure for tagging systems[END_REF]. These datasets consist of tuples of the form (users × resources × tags × time); hence we naturally form 4-mode tensors out of these tuples.

We now briefly describe the methods compared in the experiments. DFacTo is the CP-ALS routine we compare our methods against, which uses block partitioning of the rows of the matricized tensors to distribute matrix nonzeros equally among processes. In the method ht-coarsegrain-block, we operate the coarse-grain CP-ALS implementation in Hy-perTensor on a tensor whose slices in each mode are partitioned consecutively to ensure equal number of nonzeros among processes, similarly to the DFacTo's approach. The method ht-coarsegrain-hp runs the same computational kernel, but it uses the hypergraph partitioning of the tensor for the coarse-grain task definition given in Section 3.1. The method ht-finegrain-random partitions the tensor nonzeros as well as the rows of the factor matrices randomly to establish load balance. It uses the fine-grain CP-ALS implementation in HyperTensor. Finally, the method ht-finegrain-hp corresponds to the fine-grain CP-ALS implementation in HyperTensor operating on the tensor partitioned according to the hypergraph model proposed for the finegrain task decomposition in Section 3.2. We benchmark our methods against DFacTo on Netflix and NELL-B datasets only; as DFacTo can only process 3-mode tensors. Hence, for Flickr and Delicious tensors we only compare our methods using different partitioning schemes. We let the CP-ALS implementations run for 20 iterations on each data with R = 10, and record the average time spent per iteration.

In Figures 2a and2b, we show the time spent per CP-ALS iteration on Netflix and NELL-B tensors for all algorithms. In Figures 2c and2d, we show the speedups per CP-ALS iteration on Flickr and Delicious tensors for methods excluding DFacTo. While performing the comparison with DFacTo, we preferred to show the time spent per iteration instead of the speedup, as the sequential running time of our methods differs from that of DFacTo. Also in Table 2, we show the statistics regarding the communication and computation requirements of the 512-way partitions of the Netflix tensor for all proposed methods.

We first observe in Figure 2a that on the Netflix tensor ht-finegrain-hp clearly outperforms all other methods by achieving a speedup of 194x with 512 cores over a sequential execution, whereas ht-coarsegrain-hp, ht-coarsegrain-block, DFacTo, and ht-finegrain-random could only yield to 69x, 63x, 49x, and 40x speedups, respectively. Table 2 shows that on 512-way partitioning of the same tensor, ht-finegrainrandom, ht-coarsegrain-block, and ht-coarsegrain-hp result in total send communication volumes of 142M, 80M, and 77M units respectively; whereas ht-finegrain-hp partitioning requires only 7.6M units; which explains the superior parallel performance of ht-finegrain-hp. Similarly in Figures 2b, 2c, and 2d ht-finegrain-hp achieves 81x, 129x, and 123x speedups on NELL-B, Flickr, and Delicious tensors, while the best of all other methods could only obtain 39x, 55x, and 65x, respectively, upto 1024 cores. As seen from these figures, ht-finegrain-hp is the fastest of all proposed methods. It experiences slow-downs at 1024 cores for all instances. Even though other methods could scale to 1024 cores in some instances, they still remained to be significantly slower than ht-finegrain-hp.

One point to note in Figure 2b is that our MTTKRP kernels get notably more efficient than DFacTo on NELL-B, despite that DFacTo uses optimized kernels from the Eigen and MKL libraries. We believe that the main reason for this difference is due to DFacTo's column-by-column way of computing the factor matrices in order to be able to use these column vectors in the SpMV kernels. We instead operate our sparse computations on the row vectors of the factor matrices, and compute all columns of the factor matrices simultaneously. This vector-wise mode of operation on the rows of the factors can significantly improve the data locality and cache efficiency which leads to this performance difference. Another point in the same figure is that Hyper-Tensor running with ht-coarsegrain-block partitioning scales better than DFacTo, even though they use similar consecutive row or slice partitioning of the matricized tensor or the tensor. This is mainly due to the fact that DFacTo uses MPI Allgatherv to communicate local rows of factor matrices to all processes, which incurs significant increase in the communication volume, whereas our coarse-grain computation kernel performs point-to-point communications only with processes that need these rows. Finally, in Figure 2a, DFacTo has a slight edge over our kernels in sequential execution. This could be due to the fact that it requires 3(|X| + F ) multiply/add operations across all modes, where F is the average number of non-empty fibers and is upperbounded by |X|, whereas our implementation always taks 6|X| operations.

Partitioning metrics in the Table 2 show that ht-finegrainhp is able to successfully reduce the total communication volume, which brings about its improved scalability. However, we do not see the same effect when using hypergraph partitioning for the coarse-grain task model, due to inherent limitations of 1D partitionings. Another point to note is that despite reducing the total communication volume explicitly, imbalance in the communication volume still exists among the processes, as this objective is not directly captured by the hypergraph models. However, the most important limitation on further scalability is the communication latency. The results in Table 2 show that each process communicates with almost all others especially on small dimensions of the tensors, and the number of messages is doubled for the fine-grain computation kernel-recall that it has two communication phases. To alleviate this issue, one can try to limit the latency following previous work [START_REF] Selvitopi | A novel method for scaling iterative solvers: Avoiding latency overhead of parallel sparse-matrix vector multiplies[END_REF][START_REF] Uçar | Encapsulating multiple communication-cost metrics in partitioning sparse rectangular matrices for parallel matrix-vector multiplies[END_REF], and we plan to exploit this possibility in the future work.

Since DFacTo implements the MTTKRP operation in a series of parallel SpMV operations, one may consider using graph/hypergraph partitioning techniques to speed up the parallel execution. We did not investigate these alternatives for two reasons. First, our results on the parallel performance along with the statistical description of the htcoarsegrain-block and ht-coarsegrain-hp partitions in Table 2 show that there is not much hope to obtain decent parallel performance from 1D partitioning of tensors (or row-wise partitioning of matricized tensors). Second, as stated in Section 2.2, there are sophisticated data dependencies between the successive SpMV operations in DFacTo, which requires further effort to apply the graph/hypergraph models.

CONCLUSION

We have investigated the efficient parallelization of the matricized tensor times Khatri-Rao product (MTTKRP) operation in the context of the alternating least squares (ALS) method for the CANDECOMP/PARAFAC (CP) decomposition. This operation is the computational core of the CP-ALS algorithm as well as many of its variants. We have taken a close look at the computational tasks and formulated a coarse and a fine-grain task definition. The coarse-grain task definition is based on a one-dimensional partitioning of the tensor. Such partitioning are inherently limited in scalability-the same limits as in the one-dimensional partitioning based SpMV operations. The fine-grain task definition is based on a nonzero partitioning of the tensor, and overcomes these limitations. We have identified the parallelization requirements of the coarse-grain and fine-grain MTTKRP methods, and presented hypergraph models to meet those needs. We have designed a library for implementing the CP-ALS algorithm and presented scalability results

  is computed. Here, indices of the corresponding rows of the factor matrices and MU 1 coincide with indices of the unique tensor entry of the operation. Algorithm 2: MTTKRP for the 3rd order tensors Input : X: tensor B, C: Factor matrices in all modes except the first I A : Number of rows of the factor A R: Rank of the factors Output

Algorithm 3 :

 3 Coarse-grain MTTKRP for the first mode of 3rd order tensors at process p within CP-ALS Input : Ip, indices where µ A (i) = p X Ip,:,: , tensor slices B T B and C T C, R × R matrices B(Jp, :), C(Kp, :), Rows of the factor matrices, where

Figure 1 :

 1 Figure 1: Coarse and fine-grain hypergraph models for the 3 × 3 × 3 tensor X = {(1, 2, 3), (2, 3, 1), (3, 1, 2)}.

  [ai] = |Xi,:,:|, w[bj] = |X:,j,:|, and w[c k ] = |X :,:,k | to induce load balance by partitioning the hypergraph.

Algorithm 4 :

 4 Fine-grain MTTKRP for the first mode of 3rd order tensors at process p within CP-ALS Input : Xp, tensor nonzeros where µ X (i, j, k) = p Ip, the first mode indices where µ A (Ip) = p Fp, the set of unique first mode indices in Xp B T B and C T C: R × R matrices B(Jp, :), C(Kp, :): Rows of the factor matrices, where

Table 1 :

 1 Size of tensors used in the experiments

	Tensor	I 1	I 2	I 3	I 4	#nonzeros
	Netflix	480K 17K	2K	-	100M
	NELL-B	3.2M	301	638K	-	78M
	Delicious 530K 17M 2.4M	1K	140M
	Flickr	319K 28M 1.6M 713	112M

Locally compute A(Ip, :) T A(Ip, :) and all-reduce the results to form A T A.

ACKNOWLEDGMENTS

This work was performed using HPC resources from GENCI-[TGCC/CINES/IDRIS] (Grant 2015-100570). Additional computational resources were used from the PSMN com-puting center at ENS Lyon.

on up to 1024 cores. The experiments showed that the proposed fine-grain MTTKRP can achieve the best performance with respect to other alternatives with a good partitioning, reaching up to 194x speedups on 512 cores. In our analysis and experiments, we identified the communication latency as the dominant hindrance for further scalability of the fastest proposed method. We will investigate this in the future. We also note that the size of the hypergraphs that we build can cause discomfort to all existing partitioning tools. Methods that partition huge hypergraphs efficiently and effectively are needed for handling larger tensors than those treated in this work.

During the revision process, we have come across to a new study called DMS on distributed memory tensor factorization [START_REF] Smith | DMS: Distributed sparse tensor factorization with alternating least squares[END_REF] which uses SPLATT formulation. The communication is based on all-to-all primitives. The current versions of the two codes (ours using MPI, DMS using OpenMP+MPI) do not allow a useful comparison. We plan to update our codes and do a comparison in the near future.