Keywords:

computational scheme and data structure; and an effective distributed memory parallelization together with scalable data partitioning schemes, obtained through novel hypergraph models, for the parallel algorithm. This chapter lays the foundations of this thesis, as these techniques are properly adopted by methods in the subsequent chapters. In Chapter 4, we discuss the computation and parallelization of a "modified" CP decomposition, in which we put the constraint of having only nonzero entries in CP decomposition. This enables a very fast parallel algorithm whose complexity does not increase with tensor dimensionality, and provides comparable results to the standard CP decomposition in approximation quality. Next, we investigate the parallelization of Tucker decomposition algorithms in shared and distributed memory environments for sparse tensors in Chapter 5, also adopt the new computational scheme to reduce the cost of the expensive TTMc step. Afterwards, Chapter 6 discusses the distributed memory parallelization of a related problem, non-negative matrix factorization, to which most of our parallelization and partitioning techniques in tensor decompositions gracefully apply. Chapter 7 investigates an optimal use of our new computational scheme in computing dense CP and Tucker decompositions with an in-depth theoretical analysis. Finally, Chapter 8 summarizes the available high performance sparse tensor software, and exemplifies the software developed during the thesis by highlighting its capabilities. All these chapters are summarized more in detail in what follows.

Introduction

In the age of big data, the volume of collected data grows rapidly not only due to sheer increase in the number of data sources, such as the number of users of social media, but also as a natural result of the insatiate need of obtaining a more precise and in-depth understanding of the big datasets arising in real-world applications. In parallel with the growth in data size, the number of data features also escalates, which in turn increases the dimensionality of data. In response to this, tensors have been increasingly employed in the recent past owing to their ability to naturally model such high dimensional datasets, and tensor decomposition methods have proven to be an effective tool for gaining insights from the analysis of such data.

Tensor decompositions have been used in a vast range of application domains, including the 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], recommender systems [START_REF] Lars | Pairwise interaction tensor factorization for personalized tag recommendation[END_REF][START_REF] Rendle | Learning optimal ranking with tensor factorization for tag recommendation[END_REF][START_REF] Ymeonidis | Tag recommendations based on tensor dimensionality reduction[END_REF], signal processing [START_REF] Lathauwer A N D B | From matrix to tensor: Multilinear algebra and signal processing[END_REF], computer vision [START_REF] Vasilescu | Multilinear analysis of image ensembles: Ten-sorFaces[END_REF], health care [START_REF] Perros | Sparse hierarchical Tucker factorization and its application to healthcare[END_REF], chemometrics [START_REF] Dav I D S O N | Strategies for analyzing data from video fluorometric monitoring of liquid chromatographic effluents[END_REF], forensic data analysis [START_REF] Maruhashi | MultiAspectForensics: Pattern mining on large-scale heterogeneous networks with tensor analysis[END_REF], and many others [START_REF] Cichocki | Tensor decompositions for signal processing applications: From two-way to multiway component analysis[END_REF][START_REF]Tensor decompositions and applications[END_REF][START_REF] Sidiropoulos | Tensor decomposition for signal processing and machine learning[END_REF]. In these applications, tensor decomposition are used to find latent relations or predict missing elements in data using its low rank structure. For this aptitude of tensor decompositions, there have been considerable efforts in designing numerical algorithms for different tensor decomposition problems [START_REF]Tensor decompositions and applications[END_REF], and algorithmic and software contributions go hand in hand with these efforts [3,[START_REF] Bader | Matlab tensor toolbox version 2[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] Kaya A N D B | Scalable sparse tensor decompositions in distributed memory systems[END_REF][START_REF]High performance parallel algorithms for the Tucker decomposition of sparse tensors[END_REF][START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF][START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF]. These decompositions can be computed to express both dense and sparse tensors, while an efficient computation of the latter has only very recently drawn significant interest with the emerging applications using high dimensional big sparse datasets.

The two prominent tensor decomposition methods that are widely used in the literature are CAN-DECOMP/PARAFAC (CP) decomposition (or equivalently, canonical polyadic decomposition (CPD)) and Tucker decomposition. Both these formulations have uses in various applications; in particular, CP formulation is deemed useful for understanding latent components [START_REF] Kolda | The TOPHITS model for higher-order web link analysis[END_REF], whereas Tucker formulation is considered to be more appropriate for compression [5], identifying relations among the factors [START_REF] Jeon | Haten2: Billion-scale tensor decompositions[END_REF], and predicting missing data entries [START_REF] Rendle | Learning optimal ranking with tensor factorization for tag recommendation[END_REF].

CP formulation approximates a given tensor as a sum of rank-one tensors. The standard algorithm for computing CP decomposition is 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, whereas other variants also exist [1]. Most of these algorithms are iterative, in which the computational core of each iteration involves 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 the element-wise multiplication of N -1 matrix row vectors and their scaled sum reduction according to the nonzero structure of the tensor. When the dimensionality of the ten-2 CHAPTER 1. INTRODUCTION sor increases, this operation gets computationally more expensive; hence, efficiently carrying out MTTKRP for higher dimensional tensors is of our particular interest due to the needs of emerging applications [START_REF] Perros | Sparse hierarchical Tucker factorization and its application to healthcare[END_REF]. This operation has received recent attraction for efficient execution in different settings such as MATLAB [3,[START_REF] Bader | Matlab tensor toolbox version 2[END_REF], MapReduce [START_REF] Kang | GigaTensor: Scaling tensor analysis up by 100 times -Algorithms and discoveries[END_REF], shared memory [START_REF] Kaya | A parallel nonzero CP decomposition algorithm for higher order sparse data analysis[END_REF][START_REF]Parallel CP decomposition of sparse tensors using dimension trees[END_REF][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][START_REF] Kaya A N D B | Scalable sparse tensor decompositions in distributed memory systems[END_REF][START_REF]Parallel CP decomposition of sparse tensors using dimension trees[END_REF][START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF].

Tucker formulation gives another decomposition that has been employed in many data analysis problems [START_REF] Rendle | Learning optimal ranking with tensor factorization for tag recommendation[END_REF][START_REF] Ymeonidis | Tag recommendations based on tensor dimensionality reduction[END_REF][START_REF] Zhang | Cubic analysis of social bookmarking for personalized recommendation[END_REF]. One of the most popular algorithms for computing Tucker decomposition an alternating least squares (ALS) based method called Higher Order Orthogonal Iteration (HOOI) [START_REF]On the best rank-1 and rank-(R 1 ,R 2 ,...,R N ) approximation of higher-order tensors[END_REF].

A key operation in HOOI algorithm involves the multiplication of an N -dimensional tensor with N -1 matrices, which we call tensor-times-matrix chain product (TTMc). In many applications involving Tucker decomposition, the tensor formed from the data is sparse and big, both in terms of dimension sizes and the number of nonzero entries. The sparsity needs to be adequately exploited so as to compute Tucker decomposition efficiently.

The main goal of this thesis is to address such computational challenges in computing these two tensor decompositions from a multitude of perspectives, in an effort to render these tensor decomposition methods affordable to use in analyzing massive scale datasets. A particular focus is given on effectively carrying out MTTKRP and TTMc operations, which constitute the most costly step in many tensor decomposition algorithms. These approaches can be summarized as follows:

• Shared and distributed memory parallelizations with novel partitioning routines for load balancing and communication reduction.

• Tensor data structures for effectively carrying out tensor operations in parallel.

• New computational schemes and algorithms providing asymptotic computational gains in computing standard tensor decompositions.

• Faster numerical methods for handling high dimensional tensors.

• Complexity analysis pertaining to finding optimal computational schemes for tensor decompositions.

Other directions are also taken in the literature in accelerating the computation of tensor decompositions. The first approach is using randomization techniques that operate on a smaller sample of the tensor [START_REF] Krishnamurthy | Low-rank matrix and tensor completion via adaptive sampling[END_REF][START_REF] Vervliet | A randomized block sampling approach to canonical polyadic decomposition of large-scale tensors[END_REF]. The second direction involves new generalized tensor decomposition schemes and algorithms yielding faster convergence to the solution [START_REF] Sorber | Optimization-based algorithms for tensor decompositions: Canonical polyadic decomposition, decomposition in rank-$(l_r,l_r,1)$ terms, and a new generalization[END_REF]. Additionally, there are also techniques that exploit an underlying structure on the factor matrices implied by the problem or application at hand [START_REF] Sidiropoulos | Tensor decomposition for signal processing and machine learning[END_REF][START_REF] Vervliet | Tensorlab 3.0 -numerical optimization strategies for large-scale constrained and coupled matrix/tensor factorization[END_REF][START_REF] Vervliet | Breaking the curse of dimensionality using decompositions of incomplete tensors: Tensor-based scientific computing in big data analysis[END_REF]. Indeed, the focus of this thesis, using parallelization and HPC techniques for tensor decompositions, is an orthogonal research direction that can be coupled with any of these approaches, hence remains pertinent in all cases.

The approaches employed in the thesis for computing CP/Tucker decomposition of big sparse and dense tensors are organized in the following chapters as follows. In Chapter 3, we discuss novel techniques for a fast computation of the CP decomposition of big sparse tensors. This efficient computation is achieved through a multitude of approaches: a new computational scheme for carrying out the most expensive step (MTTKRP) of CP decomposition algorithms; a novel sparse tensor data structure that enables a scalable representation of high dimensional sparse tensors with significant computational and memory gains; a shared memory parallelization of CP decomposition algorithms using this performs O(N ) and O(log N ) less work in computing a CP decomposition over two efficient stateof-the-art methods, the latter being the algorithm introduced in Chapter 3. With these algorithmic contributions and a highly tuned parallel implementation, we achieve up to 16.7x speedup in sequential, and up to 10.5x speedup in parallel executions over these libraries on a 28-core workstation. In doing so, we incur with up to 24x less preprocessing time, and use up to O(log N ) less memory for storing intermediate computations. We show using a real-world tensor that the accuracy of our method is comparable to the standard CP decomposition.

Chapter 5: Parallel Tucker Decomposition of Sparse Tensors. In this chapter, we investigate an efficient parallelization of an algorithm for computing the well-known Tucker decomposition of general N -dimensional sparse tensors. The algorithm is iterative and uses the alternating least squares method. At each iteration, for each dimension of an N -dimensional input tensor, the following operations are performed: (i) the tensor is multiplied with (N -1) matrices (TTMc step); (ii) the product is then converted to a matrix; and (iii) a few leading left singular vectors of the resulting matrix are computed (truncated SVD step (TRSVD)) to update one of the matrices for the next TTMc step. We propose an efficient parallelization of these algorithms for the current parallel platforms with multicore nodes. We discuss a set of preprocessing steps which takes all computational decisions out of the main iteration of the algorithm and provides an intuitive shared-memory parallelism for the TTMc and TRSVD steps. We show how the dimension tree-based framework can be employed to perform the TTMC step faster. We propose a coarse-and a fine-grain distributed memory parallel algorithms, investigate their data dependencies, and identify efficient communication schemes. We demonstrate the way the computation of singular vectors in the TRSVD step can be carried out efficiently following the TTMc step. Finally, we develop a hybrid MPI-OpenMP implementation of the overall algorithm, and report scalability results up to 4096 cores of an IBM BlueGene/Q supercomputer.

Chapter 6: Parallel Non-negative Sparse Matrix Factorization. Non-negative matrix factorization (NMF) is the problem of finding two non-negative low rank factors W and H for a given input matrix A such that A ⇡ WH. NMF is a useful tool for many internet applications such as topic modeling in text mining and community detection in social networks. The algorithms for computing NMF have a notable similarity to tensor decomposition methods, which opens up the possibility to apply all our parallelization techniques to the computation of NMF for sparse matrices. To this end, the focus of this chapter is on scaling distributed NMF to very large sparse datasets by employing effective algorithms, communications, and partitioning schemes that leverage the sparsity of the input matrix A. Our method employs an efficient point-to-point communication scheme, and removes any redundant communications that exist in state-of-the-art algorithms. Unlike the state-of-the-art, our method also permits arbitrary partitions of the sparse input matrix as well as factor matrices, which enables us to employ effective partitioning schemes for better scalability. We experiment on a cluster using up to 3072 cores, and report significant better scalability with respect to existing methods owing to our efficient communication and partitioning strategies. We also run our algorithms on an IBM BlueGene/Q supercomputer, and achieve scalability up to 32768 processors with a 32x speedup over the execution using 256 processors.

Chapter 7: Computing Dense Tensor Decompositions using Dimension Trees. The focus of this chapter is on computing CP and Tucker decompositions of dense tensors using the dimension treebased computational scheme. We extend our earlier algorithms using dimension trees to the dense case, and propose efficient algorithms for computing both CP and Tucker decompositions. These algorithms are oblivious to the structure of the tree, which determines the computational cost in these algorithms. This raises the question of finding an optimal tree structure. We prove that finding an optimal dimension tree for computing both CP and Tucker decompositions is NP-hard. For CP decomposition, we show that an optimal tree has to be binary, whereas in the Tucker case we provide a counter-example for which the optimal dimension tree is not binary. Finally, we introduce an optimal greedy algorithm that finds the ordering minimizing the cost of a series of TTMs in computing Tucker decomposition.

Chapter 8: High Performance Tensor Software. In this chapter, we give an overview of the existing high performance software for sparse tensor decompositions, and briefly discuss the capabilities of two software libraries developed during this thesis. The first software package is HYPERTENSOR, which involves parallel algorithms for computing sparse CP and Tucker decompositions, and employs an novel sparse tensor data structure and dimension tree-based computational scheme in these parallelizations. The second software package is PACOS (Partitioning and communication framework for sparse irregular applications) that drives all partitioning and communication routines of HYPER-TENSOR for distributed memory parallelism. Even though PACOS is employed in HYPERTENSOR for parallel tensor factorization, the provided routines are generic, which enables using PACOS in the parallelization of many sparse (multi)-linear algebra routines including sparse matrix vector/matrix multiplication and various graph algorithms. We have recently integrated PACOS into NMFLIBRARY, a high performance parallel non-negative matrix factorization library, which provided much better scalability owing to effective partitioning strategies and efficient communication routines provided in PACOS. À l'âge du "big data", le volume de données recueillies augmente rapidement non seulement en raison de l'augmentation du nombre de sources de données, comme le nombre d'utilisateurs des réseaux sociaux, mais aussi comme un résultat naturel du besoin insatiable d'obtenir une compréhension plus précise et plus complète des grands ensembles de données découlant d'applications du monde réel. Parallèlement à la croissance de la taille des données, le nombre de fonctions de données augmente également, ce qui augmente la dimensionnalité des données. En réponse à cela, les tenseurs ont été de plus en plus employés dans le passé récent en raison de leur capacité à modeler naturellement de tels ensembles de données de grande dimension, et les méthodes de décomposition de tenseur se sont révélées être un outil efficace pour tirer des idées de l'analyse de ces données.

Les décompositions de tenseur ont été utilisées dans une vaste gamme de domaines d'application, y compris l'analyse des graphiques Web [START_REF] Kolda | The TOPHITS model for higher-order web link analysis[END_REF], des bases de connaissances [START_REF] Carlson | Toward an architecture for never-ending language learning[END_REF], les systèmes de recommandation [START_REF] Lars | Pairwise interaction tensor factorization for personalized tag recommendation[END_REF][START_REF] Rendle | Learning optimal ranking with tensor factorization for tag recommendation[END_REF][START_REF] Ymeonidis | Tag recommendations based on tensor dimensionality reduction[END_REF], le traitement des signaux [START_REF] Lathauwer A N D B | From matrix to tensor: Multilinear algebra and signal processing[END_REF], la vision par ordinateur [START_REF] Vasilescu | Multilinear analysis of image ensembles: Ten-sorFaces[END_REF], les soins de santé [START_REF] Perros | Sparse hierarchical Tucker factorization and its application to healthcare[END_REF], la chimiométri [START_REF] Dav I D S O N | Strategies for analyzing data from video fluorometric monitoring of liquid chromatographic effluents[END_REF], l'analyse de données médico-légales [START_REF] Maruhashi | MultiAspectForensics: Pattern mining on large-scale heterogeneous networks with tensor analysis[END_REF] et bien d'autres [START_REF] Cichocki | Tensor decompositions for signal processing applications: From two-way to multiway component analysis[END_REF][START_REF]Tensor decompositions and applications[END_REF][START_REF] Sidiropoulos | Tensor decomposition for signal processing and machine learning[END_REF]. Dans ces applications, la décomposition de tenseur sert à trouver des relations latentes ou à prédire les éléments manquants dans les données en utilisant sa structure de rang réduit. Pour cette aptitude des décompositions tensorielles, des efforts considérables ont été déployés dans la conception d'algorithmes numériques pour différents problèmes de décomposition de tenseur [START_REF]Tensor decompositions and applications[END_REF], et les contributions algorithmiques et logicielles vont de pair avec ces efforts [3,[START_REF] Bader | Matlab tensor toolbox version 2[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] Kaya A N D B | Scalable sparse tensor decompositions in distributed memory systems[END_REF][START_REF]High performance parallel algorithms for the Tucker decomposition of sparse tensors[END_REF][START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF][START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF]. Ces décompositions peuvent être calculées pour exprimer des tenseurs denses et creux, alors qu'un calcul efficace de ce dernier n'a que très récemment suscité un intérêt significatif pour les applications émergentes utilisant des ensembles de données à grande dimension et à grande échelle.

Les deux méthodes de décomposition de tenseur proéminentes qui sont largement utilisées dans la littérature sont la décomposition CANDECOMP/PARAFAC (CP) (ou équivalent, la décomposition polyadique canonique (CPD)) et la décomposition Tucker. Les deux formulations ont des utilisations dans diverses applications; en particulier, la formulation de CP est jugée utile pour comprendre les composants latents, alors que la formulation Tucker est considérée plus appropriée pour la compression, l'identification des relations entre les facteurs et la prédiction des entrées de données manquantes.

La formulation CP se rapproche d'un tenseur donné en tant que somme des tenseurs de premier rang. L'algorithme standard pour le calcul de la décomposition CP est 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], qui est basé sur la méthode des moindres carrés alternés, tandis que d'autres variantes existent également [1]. La plupart de ces algorithmes sont itératifs dans lesquels le noyau de calcul de chaque itération implique une opération spéciale appelée le produit de tenseur matricizé et le produit Khatri-Rao (MTTKRP). Lorsque le tenseur d'entrée est creux et N -dimensionnel, l'opération MTTKRP équivaut à la multiplication élémentaire des vecteurs de rangée de matrice N -1 et à leur réduction de somme à l'échelle selon la structure nonzéro du tenseur. Lorsque la dimensionnalité du tenseur augmente, cette opération devient plus coûteuse. Par conséquent, la réalisation efficace de MTTKRP pour des tenseurs dimensionnels plus élevés est de notre intérêt particulier en raison des besoins des applications émergentes [START_REF] Perros | Sparse hierarchical Tucker factorization and its application to healthcare[END_REF]. Cette opération a reçu une attraction récente pour une exécution efficace dans différents domaines tels que MATLAB [3,[START_REF] Bader | Matlab tensor toolbox version 2[END_REF], MapReduce [START_REF] Kang | GigaTensor: Scaling tensor analysis up by 100 times -Algorithms and discoveries[END_REF], mémoire partagée [START_REF] Kaya | A parallel nonzero CP decomposition algorithm for higher order sparse data analysis[END_REF][START_REF]Parallel CP decomposition of sparse tensors using dimension trees[END_REF][START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF] et mémoire distribuée [START_REF] Choi | DFacTo: Distributed factorization of tensors[END_REF][START_REF] Kaya A N D B | Scalable sparse tensor decompositions in distributed memory systems[END_REF][START_REF]Parallel CP decomposition of sparse tensors using dimension trees[END_REF][START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF].

La formulation Tucker donne une autre décomposition qui a été utilisée dans de nombreux problèmes d'analyse de données [START_REF] Rendle | Learning optimal ranking with tensor factorization for tag recommendation[END_REF][START_REF] Ymeonidis | Tag recommendations based on tensor dimensionality reduction[END_REF][START_REF] Zhang | Cubic analysis of social bookmarking for personalized recommendation[END_REF]. L'un des algorithmes les plus populaires pour calculer la décomposition Tucker est une méthode basée sur les moindres carrés alternatifs (ALS) appelée iteration orthogonale à ordre supérieur (HOOI) [START_REF]On the best rank-1 and rank-(R 1 ,R 2 ,...,R N ) approximation of higher-order tensors[END_REF]. Une opération clé dans l'algorithme HOOI implique la multiplication d'un tenseur N -dimensionnel avec des matrices N -1, que nous appelons un produit à chaîne matricielle tensorielle (TTMc). Dans de nombreuses applications impliquant la décomposition Tucker, le tenseur formé à partir des données est creux et grand, tant en termes de dimensions que de nombre d'entrées nonzéro. La sparticité doit être exploitée de manière adéquate afin de calculer efficacement la décomposition Tucker. L'objectif principal de cette thèse est d'aborder ces défis informatiques dans le calcul de ces deux décompositions tensorielles à partir d'une multitude de points de vue, dans le but de rendre ces méthodes de décomposition de tenseur abordables à utiliser dans l'analyse d'ensembles de données à grande échelle. Une attention particulière est accordée à la réalisation efficace des opérations MTTKRP et TTMc, qui constituent l'étape la plus coûteuse dans de nombreux algorithmes de décomposition de tenseur. Ces approches peuvent être résumées comme suit:

• Parallélisations en mémoire partagées et distribuées avec de nouvelles routines de partitionnement pour l'équilibrage de charge et la réduction de la communication.

• Structures de données de tenseur pour effectuer efficacement des opérations tensorielles en parallèle.

• Nouveaux schémas et algorithmes de calcul fournissant des gains de calcul asymptotiques dans le calcul des décompositions de tenseur standard.

• Méthodes numériques plus rapides pour la manipulation de tenseurs dimensionnels élevés.

• Analyse de complexité concernant la recherche de schémas informatiques optimaux pour les décompositions des tenseurs.

D'autres directions sont également prises dans la littérature pour accélérer le calcul des décompositions des tenseurs. La première approche consiste à utiliser des techniques de randomisation qui fonctionnent sur un échantillon plus petit du tenseur [START_REF] Krishnamurthy | Low-rank matrix and tensor completion via adaptive sampling[END_REF][START_REF] Vervliet | A randomized block sampling approach to canonical polyadic decomposition of large-scale tensors[END_REF]. La deuxième direction implique de nouveaux schémas généraux de décomposition de tenseur et des algorithmes donnant une convergence plus rapide à la solution [START_REF] Sorber | Optimization-based algorithms for tensor decompositions: Canonical polyadic decomposition, decomposition in rank-$(l_r,l_r,1)$ terms, and a new generalization[END_REF]. En outre, il existe également des techniques qui exploitent une structure sous-jacente sur les matrices factorisées impliquées par le problème ou la demande en cours [START_REF] Sidiropoulos | Tensor decomposition for signal processing and machine learning[END_REF][START_REF] Vervliet | Tensorlab 3.0 -numerical optimization strategies for large-scale constrained and coupled matrix/tensor factorization[END_REF][START_REF] Vervliet | Breaking the curse of dimensionality using decompositions of incomplete tensors: Tensor-based scientific computing in big data analysis[END_REF]. En effet, l'objectif de cette thèse, en utilisant la parallélisation et les techniques HPC pour les décompositions des tenseurs, est une direction de recherche orthogonale qui peut être couplée à l'une de ces approches, donc reste pertinente dans tous les cas.

Les approches employées dans la thèse pour le calcul de la décomposition CP/Tucker de grands tenseurs creux et denses sont organisées dans les chapitres suivants comme suit. Dans le Chapitre 3, nous discutons de nouvelles techniques pour un calcul rapide de la décomposition du CP de grands tensorants creux. Ce calcul efficace est réalisé grâce à une multitude d'approches: un nouveau schéma de calcul pour effectuer l'étape la plus coûteuse (MTTKRP) des algorithmes de décomposition du CP; une nouvelle structure de données de tenseur creux qui permet une représentation évolutive de tenseurs creux à haute dimension avec des gains significatifs de calcul et de mémoire; une parallélisation en mémoire partagée des algorithmes de décomposition du CP à l'aide de ce schéma de calcul et de la structure de données; et une parallélisation efficace en mémoire distribuée avec des schémas de partitionnement de données évolutifs, obtenus grâce à de nouveaux modèles d'hypergraphe, pour l'algorithme parallèle. Ce chapitre pose les bases pour les chapitre suivants car les techniques utilise dans ce chapitre sont adequatement adopte aux problèmes dans les chapitres suivants. Dans le CHAPTER 1. INTRODUCTION Chapitre 4, nous discutons du calcul et de la parallélisation d'une décomposition CP "modifiée", dans laquelle nous mettons la contrainte d'avoir seulement des entrées nonzéro dans la décomposition du CP. Ceci permet un algorithme parallèle très rapide dont la complexité n'augmente pas avec la dimensionnalité du tenseur et fournit des résultats comparables à la décomposition standard de CP en qualité d'approximation. Ensuite, nous étudions la parallélisation des algorithmes de décomposition Tucker dans les environnements de mémoire partagés et distribués pour les tensorants creux au Chapitre 5, adoptons également le nouveau schéma de calcul pour réduire le coût de l'étape TTMc coûteuse. Par la suite, le Chapitre 6 analyse la parallélisation en mémoire distribuée d'un problème associé, la factorisation des matrices non-négatives, à laquelle la plupart de nos techniques de parallélisation et de partitionnement des décompositions des tenseurs sont gracieusement appliquées. Le Chapitre 7 étudie une utilisation optimale de notre nouveau schéma informatique dans le calcul des décompositions denses de CP et de Tucker avec une analyse théorique approfondie. Enfin, le Chapitre 8 résume les logiciels haute performance disponibles de tenseurs creux et illustre le logiciel développé au cours de la thèse en mettant en évidence ses capacités. Tous ces chapitres sont résumés plus en détail dans ce qui suit.

Vue d'Ensemble de la Thèse

Chapitre 2: Contexte. Ce chapitre fournit la plupart du contexte requis par les chapitres suivants. Il introduit la notation tensorielle, décrit les algorithmes de décomposition de tenseur, fournit une description formelle d'une structure de données arborescentes et le problème de partitionnement des hypergraphes et explique les ensembles de données et les environnements de calcul parallèles utilisés dans les expériences de la plupart des chapitres.

Chapitre 3: Décomposition CP Parallèle des Tenseurs Creux Utilisant des Arbres de Dimension. Ce chapitre se concentre sur le calcul de la décomposition CP des tenseurs creux, qui a été appliquée avec succès à de nombreux problèmes bien connus dans la recherche sur le Web, l'analyse graphique, les systèmes de recommandation, les analyses de données sur les soins de santé et bien d'autres domaines. Dans ces applications, un calcul efficace de la décomposition CP des tenseurs creux est essentiel pour pouvoir traiter et analyser des données à grande échelle. Pour ce faire, nous étudions un calcul efficace et une parallélisation de la décomposition CP des tenseurs creux dans ce chapitre. Nous proposons un nouveau schéma de calcul basé sur un arbre en utilisant une structure de données appelée arbre de dimension pour réduire considérablement le coût du MTTKRP dans CP-ALS. Ensuite, nous parallélisons efficacement ce schéma de calcul dans le contexte CP-ALS dans des environnements de mémoire partagés et distribués, et proposons des modèles de distribution de données et de tâches pour une meilleure scalabilité. Nous comparons nos implémentations parallèles avec une bibliothèque de factorisation de tenseur parallèle courant en utilisant des tenseurs formés à partir d'ensembles de données du monde réel et synthétiques. Grâce à nos contributions et implémentations algorithmiques, nous rapportons jusqu'à 5,96x, 5,65x et 3,9x d'accélération dans l'exécution séquentiel et les exécutions parallèles en mémoire partagées en distribuée et obtenons une grande scalabilité allant jusqu'à 4096 coeurs sur un superordinateur IBM BlueGene/Q. Chapitre 5: Décomposition Tucker Parallèle des Tenseurs Creux. Dans ce chapitre, nous étudions une parallélisation efficace d'un algorithme pour calculer la décomposition bien connue de Tucker des tenseurs creux N -dimensionnels généraux. L'algorithme est itératif et utilise la méthode des moindres carrés alternés. À chaque itération, pour chaque dimension d'un tenseur d'entrée N -dimensionnel, les opérations suivantes sont effectuées: (i) le tenseur est multiplié par des matrices (N -1) (étape TTMc); (ii) le produit est ensuite converti en une matrice; et (iii) quelques vecteurs primaires de gauche de la matrice résultante sont calculés (étape SVD tronquée (TRSVD)) pour mettre à jour une des matrices pour la prochaine étape TTMc. Nous proposons une parallélisation efficace de ces algorithmes pour les plates-formes parallèles actuelles avec des noeuds multi-coeurs. Nous discutons d'un ensemble d'étapes de prétraitement qui supprime toutes les décisions de calcul de l'itération principale de l'algorithme et fournit un parallélisme intuitif en mémoire partagée pour les étapes TTMc et TRSVD. Nous montrons la manière dont le cadre arborescent des dimensions peut être utilisé pour effectuer l'étape TTMc plus rapidement. Nous proposons des algorithmes parallèles en mémoire distribuée, étudions leurs dépendances de données et identifions des systèmes de communication efficaces. Nous démontrons la manière dont le calcul des vecteurs singuliers dans l'étape TRSVD peut être effectué efficacement à la suite de l'étape TTMc. Enfin, nous développons une implémentation hybride MPI-OpenMP de l'algorithme et rapportons des résultats de scalabilité jusqu'à 4096 coeurs d'un superordinateur IBM BlueGene/Q. Chapitre 6: Factorisation Non-Négative Parallèle des Matrices Creuses. La factorisation matricielle non négative (NMF) est le problème consistant à trouver deux facteurs négatifs non négatifs W et H pour une matrice d'entrée donnée A telle que A = WH. NMF est un outil utile pour de nombreuses applications Internet telles que la modélisation de sujet dans l'exploration de texte et la détection communautaire dans les réseaux sociaux. Les algorithmes de calcul de NMF ont une similitude notable avec les méthodes de décomposition de tenseur, ce qui ouvre la possibilité d'appliquer toutes nos techniques de parallélisation au calcul de NMF pour les matrices dispersées. À cette fin, l'accent est mis sur la mise à l'échelle de NMF distribué sur de très grands ensembles de données dispersés en utilisant des algorithmes efficaces, des communications et des schémas de partitionnement qui tirent parti de la matrice d'entrée A. Notre méthode utilise une communication point-à-point efficace Schéma, et supprime toutes les communications redondantes qui existent dans les algorithmes courant. Contrairement à l'état de l'art, notre méthode permet également des partitions arbitraires de la matrice d'entrée sparse ainsi que des matrices de facteurs, ce qui nous permet d'utiliser des schémas de partitionnement efficaces pour une meilleure évolutivité. Nous expérimentons sur un cluster en utilisant jusqu'à 3072 noyaux et signalons une meilleure évolutivité significative par rapport aux méthodes existantes grâce à nos stratégies efficaces de communication et de partitionnement. Nous exécutons également nos algorithmes sur un superordinateur IBM BlueGene/Q et réalisons une évolutivité jusqu'à 32768 processeurs avec une accélération de 32x sur l'exécution à l'aide de 256 processeurs.

CHAPTER 1. INTRODUCTION

Chapitre 7: Calcule de Décompositions des Tenseurs Denses Utilisant des Arbres de Dimension. L'objectif de ce chapitre est de calculer les décompositions CP et Tucker des tenseurs denses en utilisant le schéma de calcul basé sur les arborescences. Nous étendons nos algorithmes antérieurs à l'aide d'arbres de dimension dans le cas dense et proposons des algorithmes efficaces pour calculer les décompositions de CP et de Tucker. Ces algorithmes sont indifférents á la structure de l'arbre, ce qui détermine le coût de calcul de ces algorithmes. Cela soulève la question de trouver une structure arborescente optimale. Nous prouvons que trouver un arbre de dimension optimale pour calculer les décompositions CP et Tucker est NP-difficile. Pour la décomposition du CP, nous montrons qu'un arbre optimal doit être binaire, alors que dans le cas de Tucker, nous fournissons un contre-exemple pour lequel l'arborescence des dimensions optimales n'est pas binaire. Enfin, nous introduisons un algorithme glouton optimal qui trouve l'ordre minimisant le coût d'une série de TTM dans le calcul de la décomposition Tucker.

Chapitre 8: Logiciels Haute Performance de Tenseur. Dans ce chapitre, nous donnons un aperçu du logiciel de haute performance existant pour les décompositions des tenseurs creux et discutons brièvement les capacités de deux bibliothèques de logiciels développées au cours de cette thèse. Le premier paquet de logiciels est HYPERTENSOR, qui implique des algorithmes parallèles pour calculer des décompositions CP et Tucker des tenseurs creux, et utilise une nouvelle structure de données tensorielle et un schéma de calcul basé sur les arborescences dans ces parallélisations. Le deuxième logiciel est PACOS (partitionnement et cadre de communication pour les applications irrégulières creuses) qui déploie toutes les routines de partitionnement et de communication de HYPERTENSOR pour le parallélisme en mémoire distribuée. Même si PACOS est utilisé dans HYPERTENSOR pour la factorisation du tenseur en parallèle, les routines fournies sont génériques, ce qui permet l'utilisation de PACOS dans la parallélisation de nombreuses routines d'algèbre linéaire ou multilinéaire, y compris le produit des matrices creuses et vecteurs et divers algorithmes des graphes. Nous avons récemment intégré PACOS dans NMFLIBRARY, une bibliothèque parallèle haute performance de factorisation des matrices non-négative, qui a permis une meilleure évolutivité grâce à des stratégies de partitionnement efficaces et à des routines de communication efficaces fournies dans PACOS.

Chapter 2 Background

Tensor notation and operations

We denote the set {1,...,N} of integers as N N for N 2 Z + . For vectors, we use bold lowercase Roman letters, as in x. For matrices, we use bold uppercase Roman letters, e.g., X. For tensors, we generally follow the notation in Kolda and Bader's survey [START_REF]Tensor decompositions and applications[END_REF]. We represent tensors using bold calligraphic fonts, e.g., X . The order of a tensor is defined as the number of its dimensions, or equivalently, modes, which we denote by N . We use italic lowercase letters with corresponding indices to represent vector, matrix, and tensor elements, e.g., x i for a vector x, x i,j for a matrix X , and x i,j,k for a 3-dimensional tensor X . For the column vectors of a matrix, we use the same letter in lowercase and with a subscript corresponding to the column index, e.g., x i to denote X (:,i).Aslice of a tensor in the nth mode is a set of tensor elements obtained by fixing the index only along the nth mode. A fiber of a tensor is defined by fixing every index but one, and a slice of a tensor is obtained by fixing only one index. We use the MATLAB notation to refer to matrix rows and columns as well as tensors fibers and slices, e.g., X (i, :) and X (:,j) are the ith row and the jth column of X , whereas X (i, j, :) and X (:, :,k) represent a fiber and a slice of X in the third dimension, respectively.

The multiplication of an

N -dimensional tensor X 2 R I 1 ⇥•••⇥I N along a dimension n 2 N N with a vector v 2 R In is a tensor Y 2 R I 1 ⇥•••⇥I n-1 ⇥1⇥I n+1 ⇥•••⇥I N with elements y i 1 ,...,i n-1 ,1,i n+1 ,...,i N = In X j=1 v j x i 1 ,...,i n-1 ,j,i n+1 ,...,i N . (2.1)
This operation is called tensor-times-vector multiply (TTV) and is denoted by

Y = X ⇥ n v. The cost of this operation is O( Q i2N N I i ). X can also be multiplied with a matrix U 2 R K⇥In in a mode n, which results in a tensor Y = X ⇥ n U, Y 2 R I 1 ⇥•••⇥I n-1 ⇥K⇥I n+1 ⇥•••⇥I N with elements y i 1 ,...,i n-1 ,k,i n+1 ,...,i N = In X j=1 u k,j x i 1 ,...,i n-1 ,j,i n+1 ,...,i N .
(2.2)

In other words, TTV of X with each row vector of U forms the kth slice of Y in the nth dimension. This operation is called tensor-times-matrix multiply (TTM) or n-mode product and has the cost O(K

Q i2N N I i ).
The order of TTVs or TTMs in a set of distinct modes is irrelevant, i.e., X ⇥ i u ⇥ j w = X ⇥ j w ⇥ i u for u 2 R I i , w 2 R I j , i 6 = j, and i, j 2 N N .

CHAPTER 2. BACKGROUND

A tensor X can be matricized in some modes, meaning that a matrix X can be associated with X by identifying a subset of its modes to correspond to the rows of X, and the rest of the modes to correspond to the columns of X . This involves a mapping of the elements of X to those of the matricization X of the tensor. We will be exclusively dealing with the matricizations of tensors along a single mode, meaning that a single mode is mapped to the rows of the resulting matrix, and the rest of the modes correspond to its columns. We use X (n) to denote matricization along a mode n, e.g., for

X 2 R I 1 ⇥•••⇥I N , the matrix X (1) 2 R I 1 ⇥I 2 I 3 .
..I N denotes the mode-1 matricization of X . Specifically, the tensor element x i 1 ,...,i N corresponds to the element

⇣ i 1 ,i 2 + P N j=3 h (i j -1) Q j-1 k=2 I k i⌘ of X (1)
in this matricization. Matricizations in other modes are defined likewise.

The Hadamard product of two vectors u, v 2 R I is a vector w = u ⇤ v, w 2 R I with elements

w i = u i • v i .
The outer product of N>1 vectors u (1) ,...,u (N ) of respective sizes I 1 ,...,I N is denoted by

X = u (1) •••••u (N ) where X 2 R I 1 ⇥•••⇥I N is a N -dimensional tensor with elements x i 1 ,...,i N = Q t2N N u (t)
it . The Kronecker product of vectors u 2 R I and v 2 R J results in a vector w = u ⌦ v, w 2 R IJ , defined as

w = u ⌦ v = 2 6 6 6 4 u 1 v u 2 v . . . u I v 3 7 7 7 5
.

For matrices U 2 R I⇥K and V 2 R J⇥K , their Khatri-Rao product corresponds to the Kronecker product of their corresponding columns, i.e.,

W = U V = [u 1 ⌦ v 1 ,...,u K ⌦ v K ] , (2.3)
where W 2 R IJ⇥K .

We use the shorthand notation • i6 =n u (i) to denote operation • on a set {u

1 ,...,u (N ) } of operands, e.g., X ⇥ i6 =n u (N ) , and ⇤ n2N N u (n) gives u (1) ⇤•••⇤u (N ) .

(i) denotes X ⇥ 1 u (1) ⇥ 2 •••⇥ n-1 u (n-1) ⇥ n+1 u (n+1) ⇥ n+2 •••⇥ N u
We next describe the algorithms for computing Tucker and CP decompositions. Throughout the rest of the manuscript, we assume that the input tensor

X 2 R I 1 ⇥ 1 •••⇥ N I N has dimensions I n ≥ 2
for n 2 N N , as we can otherwise remove all dimensions of size 1, and execute the algorithms on the lower dimensional tensor.

CP decomposition

The rank-R CP decomposition of a tensor X expresses or approximates X as a sum of R rank-1 tensors. For instance, a rank-R CP decomposition yields

X ⇡ P R r=1 a r • b r • c r for X 2 R I⇥J⇥K with a r 2 R I , b r 2 R J ,
and c r 2 R K . This decomposition provides an element-wise approximation (or equality) x i,j,k ⇡ P R r=1 a ir b jr c kr , as shown in Figure 2.1. The minimum R value rendering this approximation an equality is called the rank (or CP-rank) of the tensor X , and computing this rank is NP-hard [START_REF] Håstad | Tensor rank is np-complete[END_REF]. Here, the matrices

A =[a 1 ,...,a R ], B =[b 1 ,...,b R ], and C =[c 1 ,...,c R ]
are called factor matrices, or factors. For N -mode tensors, we use U (1) ,...,U (N ) to refer to factor matrices having I 1 ,...,I N rows and R columns, u (n) j to refer to the jth column of U (n) , and u (n) i,j to denote the element of U (n) with the index (i, j).

≈ + = + ... I ⊗ ⊗ a 1 b 1 c 1 ⊗ a R b R c R C x i,j,k X B A i J K I J K I J K I J K R R R j k Figure 2.1: CP decomposition of a 3rd order tensor X 2 R I⇥J⇥K .
The traditional algorithm for computing CP decomposition is CP-ALS; an iterative algorithm, shown in Algorithm 1, that progressively updates factors U (n) in an alternating fashion starting from an initial guess. CP-ALS continues until it can no longer improve the solution, or it reaches the allowed maximum number of iterations. Factor matrices can be initialized randomly, or using the truncated SVD of the matricizations of X [START_REF]Tensor decompositions and applications[END_REF]. Each iteration of CP-ALS consists of N subiterations, where in the nth subiteration U (n) is updated using X as well as the current values of all other factor matrices.

Algorithm 1 CP-ALS: ALS algorithm for computing CP decomposition

Input: X : An N -dimensional tensor, X 2 R I1⇥•••⇥I N R:
The rank of CP decomposition U (1) ,...,U (N ) : Initial factor matrices Output: [[ λ; U (1) ,...,U (N ) ]] : A rank-R CP decomposition of X 1: for n =1,...,N do I Initialization 2:

W (n) U (n) T U (n) 3: repeat I Main iteration 4:
for n =1,...,N do I Subiterations 5:

M (n) X (n) ( i6 =n U (i) ) I MTTKRP 6: H (n) ⇤ i6 =n W (i)
7:

U (n) M (n) H (n) † I H (n) † is the pseudoinverse of H (n) . 8: λ COLUMN-NORMALIZE(U (n) )
I Normalize columns and store the norms in λ.

9:

W (n) U (n) T U (n) 10:
until convergence or the maximum number of iterations 11: return [[ λ; U (1) ,...,U (N ) ]]

Computing the matrix M (n) 2 R In⇥R at Line 5 of Algorithm 1 is the sole part involving the tensor X , and it is the most expensive computational step of CP-ALS. The operation

X (n) ( i6 =n U (i) )
is called the matricized tensor-times Khatri-Rao product (MTTKRP). The Khatri-Rao product of the involved U (n) s defines a matrix of size (

Q i6 =n I i ) ⇥ R according to (2.
3), and can get very costly in terms of computational and memory requirements when I i or N is large, which is the case for many real-world sparse tensors. To alleviate this, various methods are proposed in the literature that enable performing MTTKRP without forming the Khatri-Rao product. One such formulation [6] expresses MTTKRP in terms of a series of TTVs, and computes the resulting matrix M (n) column by column. With this formulation, the rth column of M (n) can be computed using N -1 TTVs as in

M (n) (:,r) X ⇥ i6 =n u (i) r .
(2.4)

CHAPTER 2. BACKGROUND

Another equivalent formulation computes MTTKRP as a summed reduction of Hadamard products of factor matrix rows scaled with tensor entries, which is given in Algorithm 2. For a sparse tensor Algorithm 2 Performing MTTKRP for a tensor X in a mode n.

Input: X : An N -dimensional tensor, X 2 R I1,...,I N U (1) ...U (N ) : Factor matrices n: The dimension of matricization for MTTKRP Output:

M (n) : MTTKRP result matrix, M (n) 2 R In⇥R 1: M (n) 0 2: for x i1,...,in,...,i N 2 X do 3: M (n) (i n , :) += x i1,...,in,...,i N ⇥ ⇤ j6 =n U (j) (i j , :) ⇤
having nnz(X ) nonzero entries, the overall cost of this algorithm is O(nnz(X )NR) as each nonzero induces an Hadamard product of N -1 row vectors followed by an addition of a row vector of size R. This algorithm is used is commonly used by many existing implementations [6,[START_REF] Kang | GigaTensor: Scaling tensor analysis up by 100 times -Algorithms and discoveries[END_REF].

Once M (n) is obtained, the Hadamard product of the matrices U (i) T U (n) of size R ⇥ R is computed at Line 6 to form the matrix H (n) 2 R R⇥R . Note that within the subiteration n, only U (n) is updated among all factor matrices. Therefore, for efficiency, one can precompute all matrices

W (i) = U (i) T U (i) of size R ⇥ R for 1  i  N , performed at Line 2, then update U (n) T U (n) once U (n) changes as in Line 9.
As the rank R of approximation is typically much smaller than tensor dimensions I n in practice, performing these Hadamard products to compute H (n) and the matrixmatrix multiplication to compute U (n) T U (n) is cheaper compared with the TTV step. Once both M (n) and H (n) are computed, another matrix-matrix multiplication is performed using M (n) and the pseudoinverse of H (n) in order to update the matrix U (n) , which is not very expensive for the same reason. Finally, U (n) is normalized column-wise, and the norms of column vectors are stored in a vector λ 2 R R . Convergence is achieved when the relative improvement in the error norm, i.e.,

|X -

P R r=1 λ r (u

(1) r •••••u (N )
r )|, is deemed small. The cost of this computation is insignificant; hence, we skip its details.

Tucker decomposition

Tucker decomposition expresses a given tensor X 2 R I 1 ⇥•••⇥I N with a smaller core tensor G multiplied by a factor matrix U (n) of size

I n ⇥ R n in each mode n 2 N N .
Here, the tuple (R 1 ,...,R N ) forms the requested rank of the decomposition across different modes. This decomposition is denoted by [[ G; U (1) ,...,U (N ) ]] , which expresses (or approximates)

X as G ⇥ 1 U (1) ⇥ 2 •••⇥ N U (N ) . For example, an (R 1 ,R 2 ,R 3 )-rank Tucker decomposition [[ G; U (1) , U (2) , U (3) ]] of X 2 R I 1 ⇥I 2 ⇥I 3 pro- vides x i,j,k ⇡ P R 1 p=1 P R 2 q=1 P R 3 r=1 g p,q,r u (1)
i,p u

(2) A well-known algorithm for computing Tucker decomposition is Higher Order Orthogonal Iteration (HOOI) [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF], which is shown in Algorithm 3. In this algorithm, factor matrices are initialized first. This initialization can be done randomly or using the higher-order SVD of X [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF]. Then, the "repeat-until" loop applies the ALS method. Here, in each mode n, X ⇥ i6 =n U (i) T is computed at Line 4, and we call this operation tensor-times-matrix chain product (TTMc). This produces a tensor of size

j,q u (3) k,r .
R 1 ⇥ R 2 ⇥•••⇥R n-1 ⇥ I n ⇥ R n+1 ⇥•••⇥R N , which is then matricized along the nth mode U (1) ⊗ ≈ X G I 1 U (2) U ( 3 ) I 2 I 3 I 1 I 2 I 3 R 1 R 2 R 3 Figure 2.2: Tucker decomposition of a 3rd order tensor X 2 R I 1 ⇥I 2 ⇥I 3 as a core tensor G 2 R R 1 ⇥R 2 ⇥R 3 multiplied by matrices U (1) 2 R I 1 ⇥R 1 , U (2) 2 R I 2 ⇥R 2 and U (3) 2 R I 3 ⇥R 3 in different modes. In the CP-decomposition, G is a diagonal tensor having the same size R 1 = R 2 = R 3 = R
along each dimension, and U (1) , U (2) and U (3) have the same number of columns.

into the matrix

Y (n) 2 R In⇥Π i6 =n R i .
Next, the leading R n left singular vectors of Y (n) are computed to form the new U (n) at Line 5. After all matrices U (n) are updated, the core tensor G is formed at Line 6, and the fit (|X |-|G|)/|X | is measured to check convergence at the end of each iteration.

Algorithm 3 HOOI: ALS algorithm for computing Tucker decomposition

Input:

X : An N -dimensional tensor, X 2 R I1⇥•••⇥I N R 1 ,...,R N :
The rank of Tucker decomposition in each dimension Output: [[ G; U (1) ,...,U (N ) ]] : A rank-(R 1 ,...,R N ) Tucker decomposition of X 1: Initialize the matrix

U (n) 2 R In⇥Rn for n 2 N N . 2: repeat 3:
for n =1,...,N do 4:

Y X ⇥ i6 =n U (i) T 5: U (n) R n leading left singular vectors of Y (n) 6: G X ⇥ 1 U (1) T ⇥ 2 •••⇥ N U (N ) T 7:
until convergence or the maximum number of iterations 8: return [[ G; U (1) ,...,U (N ) ]]

Dimension tree

A dimension tree is a data structure that partitions the mode indices of an N -dimensional tensor in a hierarchical manner for computing tensor decompositions efficiently. It was first used in the hierarchical Tucker format representing the hierarchical Tucker decomposition of a tensor [START_REF] Grasedyck | Hierarchical singular value decomposition of tensors[END_REF], which was introduced as a computationally feasible alternative to the original Tucker decomposition for higher order tensors. We provide the formal definition of a dimension tree along with some basic properties as follows.

Definition 1. A dimension tree T for N dimensions is a tree with a root, denoted by ROOT(T ), and N leaf nodes, denoted by the set LEAVES(T ). In a dimension tree T , each non-leaf node has at least CHAPTER 2. BACKGROUND two children, and each node t 2 T is associated with a mode set µ(t) ✓ N N satisfying the following properties:

1. µ(ROOT(T )) = N N .
2. For each non-leaf node t 2 T , the mode sets of its children partition µ(t).

3. The nth leaf node, denoted by L n 2 LEAVES(T ), has µ(L n )={n}.

For the simplicity of the presentation, we assume w.l.g that the sequence L 1 ,...,L N corresponds to a relative ordering of the leaf nodes in a post-order traversal of the dimension tree. If this is not the case, we can relabel tensor modes accordingly. We define the inverse mode set of a node t as µ 0 (t)=N N \ µ(t). For each node t with a parent P (t), µ(t) ⇢ µ(P (t)) holds due to the second property, which in turn yields µ 0 (t) ⊃ µ 0 (P (t)).

Hypergraph partitioning

A hypergraph H =(V, E) consists of a set V of vertices and a set E of hyperedges. Each hyperedge is a subset of V . 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,aK-way vertex partition of a hypergraph H =( V, E) is denoted by Π={V 1 ,...,V K }, where the parts are non-empty, i.e., V k 6 = ; for k 2 N K ; mutually exclusive, i.e., V k \ V `= ; for k 6 = `; and collectively exhaustive, i.e., V = S V k .

Let W k = P v2V k w[v] be the total vertex weight in V k , and W avg = P v2V w[v]/K denote the average part weight. If each part V k 2 Π satisfies the balance criterion

W k  W avg (1 + "), for k 2 N K , (2.5)
we say that Π is balanced where " represents the allowed maximum 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 is called its connectivity, and is denoted by λ h . Given a vertex partition Π of a hypergraph H =(V, E), one can measure the cutsize metric induced by Π as

χ(Π) = X h2E c[h](λ h -1) . (2.6)
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 NP-hard [START_REF] Lengauer | Combinatorial Algorithms for Integrated Circuit Layout[END_REF].

A common 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 2 N C . In this variant, the balance criterion (2.5) is rewritten as

W k,i  W avg,i (1 + ") for k 2 N K and i 2 N C , (2.7)
2.6. TENSOR DATASETS 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 = P v2V k w[v, i], and W avg,i represents the average part weight for the ith weight of all vertices, i.e., W avg,i = P v2V w[v, i]/K.

Tensor datasets

Most experiments throughout the manuscript involve six real-world sparse tensors whose sizes are shown in Table 2.1. Netflix tensor has user ⇥ movie ⇥ time dimensions, which we formed from the data of the Netflix Prize competition [START_REF] Bennett | The netflix prize[END_REF]. In this tensor, nonzeros correspond to the user reviews for movies, and the 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. We obtained the NELL tensor 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 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. There is also a smaller version of this dataset from which we construct another tensor named NELL2. Delicious and Flickr are the datasets for the web-crawl of Delicious.com and Flickr.com during 2006 and 2007, which are 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 form (time ⇥ users ⇥ resources ⇥ tags); hence, we form 4-mode tensors out of these tuples. We obtained the Amazon review dataset from SNAP [START_REF] Leskovec | SNAP Datasets: Stanford large network dataset collection[END_REF], which contains product review texts by users. We first processed this dataset with the standard text processing routines. We used the nltk package [START_REF] Bird | Natural Language Processing with Python[END_REF] in Python to tokenize the review text, to discard the stop words, to apply Porter stemmer, and to keep the words that are in the US, GB, or CA dictionaries. Afterwards, we retained only the words with at least five occurrences in the whole review set. Then, we created a three dimensional tensor whose dimensions correspond to the users, products, and retained words. Numerical values are set to the frequency of a word in a review.

We also used synthetic tensors created randomly having 4, 8, 16, and 32 dimensions. In these random tensors, each dimension is of size 10M, and there are 100M nonzeros with a uniform random distribution of indices. Using these tensors, we measure the effect of tensor dimensionality on the performance.

CHAPTER 2. BACKGROUND

Parallel compute environments

Most of the experiments in the manuscript were conducted on a shared memory and a separate distributed memory system, which are described in what follows. The shared memory system has two CPU sockets (Intel(R) Xeon(R) E5-2695 v3) each having 14 cores at a clock speed of 2.30GHz with Turbo Boost disabled. The system has a total memory of size 768GB and L1, L2, L3 caches of sizes of 32KB, 256KB, and 35MB, respectively. We call this system Grunch throughout the manuscript. The distributed memory system is an IBM Blue Gene/Q cluster. This system consists of 6 racks of 1024 nodes with each node having 16 GBs of memory and a 16-core IBM PowerPC A2 processor running at 1.6 GHz. This system is referred to as Blue Gene/Q in the rest of the document. Other compute environments used specific to a chapter are described within that chapter.

Part I Parallel Sparse Tensor and Matrix Decompositions

Chapter 3

Parallel CP Decomposition of Sparse Tensors Using Dimension Trees

In this chapter, we investigate the parallelization of CP-ALS algorithm in shared and distributed memory environments for sparse tensors. For the shared memory computations, we propose a novel computational scheme that significantly reduces the computational cost of MTTKRP while providing an effective parallelization. We then perform theoretical analyses corresponding to the computational gains and the memory utilization of this scheme, which are also validated with experiments. Next, we propose a fine-grain distributed memory parallel CP-ALS algorithm, and compare it against a medium-grain variant [START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF]. Finally, we discuss effective partitioning routines for these algorithms.

The organization of the rest of the chapter is as follows. We explain the use of dimension trees to carry out MTTKRPs within CP-ALS in the next section. In Section 3.2, we discuss an effective shared and distributed memory parallelization of CP-ALS iterations. Afterwards, we introduce partitioning problems pertaining to distributed memory parallel performance, and use these partitioning methods in our experiments using real-world tensors. Section 3.3 provides an overview of the existing literature. Finally, we present experimental results to demonstrate the performance gains of our algorithms for shared and distributed memory parallelism over an efficient state of the art implementation in Section 3.4.

Computing CP decomposition using dimension trees

As introduced in Section 2.4, a dimension tree partitions the mode indices of an N -dimensional tensor in a hierarchical manner for computing tensor decompositions efficiently. It was first used in the hierarchical Tucker format representing the hierarchical Tucker decomposition and SVD of a tensor [START_REF] Grasedyck | Hierarchical singular value decomposition of tensors[END_REF].

Here, we propose a novel way of using dimension trees for computing the standard CP decomposition of tensors with a formulation that asymptotically reduces the computational cost. In doing so, we do not alter the original CP decomposition in any way. The reduction in the computational cost is made possible by storing partial TTV results, and hence by trading off more memory. A similar idea of reusing partial results without the use of a tree framework was moderately explored by Baskaran et al. [START_REF] Baskaran | Efficient and scalable computations with sparse tensors[END_REF] for computing the Tucker decomposition of sparse tensors and by Phan et al. for computing the CP decomposition of dense tensors [START_REF] Phan | Fast alternating LS algorithms for high order CANDECOMP/PARAFAC tensor factorizations[END_REF]. We adopt the dimension tree data structure to fully utilise the potential of reducing the cost of MTTKRP operations in computing the CP decomposition. If a dimension tree has the height dlog(N )e with its first blog(N )c levels forming a complete binary tree, we call it a balanced binary dimension tree (BBDT). In Figure 3.1, we show a BBDT for 4 dimensions (associated with a sparse tensor described later).

Using dimension trees to perform successive TTVs

At each subiteration of the CP-ALS algorithm, using (2.4), X is multiplied with the column vectors of matrices in N -1 modes in performing the MTTKRP. Some of these TTVs involve the same matrices as the preceding subiterations. As a series of TTVs can be done in any order, this opens up the possibility to factor out and reuse TTV results that are common in consecutive subiterations for reducing the computational cost. For instance, in the first subiteration of CP-ALS using a 4-mode tensor X , we compute X ⇥ 2 u

(2)

r ⇥ 3 u (3) r ⇥ 4 u (4)
r and eventually update u

(1)
r for each r 2 N R , whereas in the second subiteration we compute X ⇥ 1 u

(1)

r ⇥ 3 u (3) r ⇥ 4 u (4)
r and update u

(2) r . In these two subiterations, the matrices U (3) and U (4) remain unchanged, and both TTV steps involve the TTV of X with u

r and Y r ⇥ 1 u (1)
r to obtain the final TTV results for these modes. We use dimension trees to systematically detect and reuse such partial results by associating a tree T with an N -dimensional tensor X as follows. With each node t 2 T , we associate R tensors

T 1 (t),...,T R (t). T r (t) corresponds to the TTV result T r (t)=X ⇥ m 1 u (m 1 ) r ⇥ m 2 ••• ⇥ m |µ 0 (t)| u (m |µ 0 (t)| ) r
, where m i 2 µ 0 (t) are distinct elements of the inverse mode set µ 0 (t). Therefore, the inverse mode set µ 0 (t) corresponds to the set of modes in which TTV is performed on X to form T r (t). For the root of the tree, µ 0 (ROOT(T )) = ;; thus, all tensors of the root node correspond to the original tensor X . Since µ 0 (t) ⊃ µ 0 (P (t)) for a node t and its parent P (t), tensors of P (t) can be used as partial results to update the tensors of t. Let δ(t)=µ 0 (t)\µ 0 (P (t)). We can then compute each tensor of t from its parent's as T r (t)=T r (P (t))⇥ d 1 u

(d 1 ) r ⇥ d 2 •••⇥ d |δ(t)| u (d |δ(t)| ) r for δ(t)={d 1 ,...,d |δ(t)| }.
This procedure is called DTREE-TTV and is shown in Algorithm 4.D TREE-TTV first checks if the tensors of t are already computed, and immediately returns if so. This happens, for example, when two children of t call DTREE-TTV on t consecutively, in which case for the second DTREE-TTV call, the tensors of t would be already computed. If the tensors of t are not already computed, DTREE-TTV on P (t) is called first to make sure that P (t)'s tensors are up-to-date. Then, each T r (t) is computed by performing a TTV on the corresponding tensor T r (P (t)) of the parent. We use the notation T : (t) to denote all R tensors of a node t. T r (t) T r (P (t)) ⇥ d1 u

(d1) r ⇥ d2 •••⇥ d |δ(t)| u (d |δ(t)| ) r CHAPTER 3.

PARALLEL CP DECOMPOSITION OF SPARSE TENSORS USING DIMENSION TREES

Algorithm 5 DTREE-CP-ALS: Dimension tree-based CP-ALS algorithm Input: X : An N -mode tensor R: The rank of CP decomposition Output: [λ; U (1) ,...,U (N ) ]: Rank-R CP decomposition of X 1: T CONSTRUCT-DIMENSION-TREE(X ) I The tree has leaves {L 1 ,...,L N }. 2: for n =2...N do 3:

W (n) U (n) T U (n) 4: repeat 5:
for n =1,...,N do

6:

for all t 2 T do 7:

if n 2 µ 0 (t) then 8:

DESTROY(T : (t))

I Destroy all tensors that are multiplied by U (n) . 9:

DTREE-TTV(L n ) I Perform the TTVs for the leaf node tensors.

10:

for r =1,...,Rdo I Form M (n) column-by-column (done implicitly).

11:

M (n) (:,r) T r (L n ) I T r (L n
) is a vector of size I n for the leaf node L n .

12:

H (n) ⇤ i6 =n W (i)
13:

U (n) M (n) H (n) † 14: λ COLUMN-NORMALIZE(U (n) )
15:

W (n) U (n) T U (n) 16:
until converge is achieved or the maximum number of iterations is reached. 17: return [[ λ; U (1) ,...,U (N ) ]]

Dimension tree-based CP-ALS algorithm

At the nth subiteration of Algorithm 1, we need to compute X ⇥ i6 =n u (i) r for all r 2 N R in order to form M (n) . Using a dimension tree T with leaves L 1 ,...,L N , we can perform this simply by executing DTREE-TTV(L n ), after which the rth tensor of L n provides the rth column of M (n) . Once M (n) is formed, the remaining steps follow as before. We show the whole CP-ALS using a dimension tree in Algorithm 5. At Line 1, we construct a dimension tree T with the leaf order L 1 ,...,L N obtained from a post-order traversal of T . This tree can be constructed in any way that respects the properties of a dimension tree described in Section 3.1; but for our purposes we assume that it is formed as a BBDT. At Line 8 within the nth subiteration, we destroy all tensors of a node t if its set of multiplied modes µ 0 (t) involve n, as in this case its tensors involve multiplication using the old value of U (n) which is about to change. Note that this step destroys the tensors of all nodes not lying on a path from L n to the root. Afterwards, DTREE-TTV is called at Line 9 for the leaf node L n to compute its tensors. This step computes (or reuses) the tensors of all nodes from the path from L n to the root. Next, the rth column of M (n) is formed using T r (L n ) for r 2 N R . Once M (n) is ready, H (n) and U (n) are computed as before, after which U (n) is normalized.

Performing TTVs in CP-ALS using a BBDT in this manner provides significant computational gains with a moderate increase in the memory cost. We now state two theorems pertaining to the computational and memory efficiency of DTREE-CP-ALS.

Theorem 1. Let X be an N -mode tensor. The total number of TTVs at each iteration of Algorithm 5 using a BBDT is at most RN dlog N e.

Proof. As we assume that the sequence L 1 ,...,L N is obtained from a post-order traversal of T , for each internal node t, the subtree rooted at t has the leaves L i , L i+1 ,...,L i+k-1 corresponding to k consecutive mode indices for some positive integers i and k. As we have µ(L i )=i for the ith leaf node, we obtain µ(t)={i, i +1,...,i+ k -1} due to the second property of dimension trees. As a result, for each leaf node L i+k 0 for 0  k 0 <k , we have i + k 0 2 µ(t); hence i + k 0 / 2 µ 0 (t) as µ 0 (t)=N N \ µ(t). Therefore, within an iteration of Algorithm 5 the tensors of t get computed at the ith subiteration, stay valid (not destroyed) and get reused until the i + k -1th subiteration, and finally get destroyed in the following subiteration. Once destroyed, the tensors of t are never recomputed in the same iteration, as all the modes associated with the leaf tensors in its subtree are processed (which are the only nodes that can reuse the tensors of t). As a result, in one CP-ALS iteration, tensors of every tree node (except the root) get to be computed and destroyed exactly once. Therefore, the total number of TTVs in one iteration becomes the sum of the number of TTVs performed to compute the tensors of each node in the tree once. In computing its tensors, every node t has R tensors, and for each tensor it performs TTVs for each dimension in the set δ(t) in Algorithm 4, except the root node, whose tensors all are equal to X and never change. Therefore, we can express the total number of TTVs performed within a CP-ALS iteration due to one of these R tensors as

X t2T \{ROOT(T )} |δ(t)| = X t2T \{ROOT(T )} |µ(P (t)) \ µ(t)|. ( 3
X t2T \{ROOT(T )} |µ(P (t)) \ µ(t)| = X t2T \LEAVES(T ) |µ(t)| . (3.2)
Note that in constructing a BBDT, at the root node we start with the mode set µ(ROOT(T )) = N N .

Then, at each level l>0, we form the mode sets of the nodes at level l by partitioning the mode sets of their parents at level l -1. As a result, at each level l, each dimension n 2 N N can appear in only one set µ(t) for a node t belonging to the level l of the BBDT. With this observation in mind, rewriting (3.2) by iterating over nodes by levels of the BBDT yields

X t2T \LEAVES(T ) |µ(t)| = dlog N e X l=1 X t2T \LEAVES(T ),LEVEL(t)=l |µ(t)|  dlog N e X l=1 N = N dlog N e .
As there are R tensors in each BBDT tree node, the overall cost becomes RN dlog N e TTVs for a CP-ALS iteration.

Note in comparison that the traditional scheme [START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF] incurs R(N -1) TTVs in each mode, and RN (N -1) TTVs in total in a CP-ALS iteration. This yields a factor of (N -1)/ log N reduction in the number performed of TTVs using dimension trees. We would like to stress that in terms of the actual computational cost, this corresponds to a lower bound on the expected speedup for the following reason. As the tensor is multiplied in different dimensions, resulting tensors are expected to have many index overlaps, effectively reducing their number of nonzeros and rendering subsequent

CHAPTER 3. PARALLEL CP DECOMPOSITION OF SPARSE TENSORS USING DIMENSION TREES

TTVs significantly cheaper. This renders using a BBDT much more effective as it avoids repeating such expensive TTVs at the higher levels of the dimension tree by reusing partial results. In the light of this observation, we argue the actual obtained speedup the following way. On one extreme, multiplying tensor in dimensions might create no or few index overlaps, which makes the cost of each TTV approximately equal, yielding the speedup factor of approximately N/ log N . On the other extreme, the first TTVs performed the original tensor may drastically reduce the number of tensor elements so that the cost of the subsequent TTVs becomes negligible. The traditional scheme multiplies the original tensor N times, once per dimension, wheres a BBDT suffices with 2 such TTVs as the root node has only two children. Therefore, we expect in practice the actual speedup to be between these two extremes depending on the sparsity of the tensor, and having more speedup with higher index overlap after multiplications.

For sparse tensors, one key idea we use for obtaining high performance is performing TTVs for all R tensors T : (t) of a node t 2 T in a vectorized manner. We illustrate this on a 4-dimensional tensor X and a BBDT, and for clarity, we put the mode set µ(t) of each tree node t in the subscript, as in t 1234 . Let t 1234 represent the root of the BBDT with T r (t 1234 )=X for all r 2 N R . The two children of t 1234 are t 12 and t 34 with the corresponding tensors T r (t 12 )=T r (t 1234 ) ⇥ 3 u

(3) r ⇥ 4 u (4) r and T r (t 34 )=T r (t 1234 )⇥ 1 u (1) r ⇥ 2 u (2)
r , respectively. Since T r (t 1234 )s are identical for all r 2 N R , the nonzero pattern of all tensors T r (t 12 ) are identical for all r 2 N R . The same is true for T r (t 34 )s. The same argument applies to the children t 1 and t 2 of t 12 , as well as t 3 and t 4 of t 34 . As a result, each node in the tree involves R tensors with identical nonzero patterns. This opens up two possibilities in terms of efficiency. First, it is sufficient to compute only one set of nonzero indices for each node t 2 T to represent the nonzero structure of all of its tensors, which reduces the computational and memory cost by a factor of R. Second, we can perform the TTVs for all tensors at once in a vectorized manner by modifying (2.1) to perform R TTVs of the form Y (r)

X (r) ⇥ d v r for V =[ v 1 | ...|v R ] 2 R I d ⇥R as follows: y (:) i 1 ,...,i d-1 ,1,i d+1 ,...,i N = I d X j=1 V(j, :) ⇤ x (:) i 1 ,...,i d-1 ,j,i d+1 ,...,i N , (3.3)
where y (:) i 1 ,...,i d-1 ,1,i d+1 ,...,i N and x (:) i 1 ,...,i d-1 ,j,i d+1 ,...,i N are vectors of size R with elements y (r) i 1 ,...,i d-1 ,1,i d+1 ,...,i N and x (r) i 1 ,...,i d-1 ,j,i d+1 ,...,i N in Y (r) and X (r) , for all r 2 N R . We call this operation tensor-timesmultiple-vector multiplication (TTMV) as R column vectors of V are multiplied simultaneously with R tensors of identical nonzero patterns. We can similarly extend this formula to the multiplication

Z (r) Y (r) ⇥ d 2 W =( X (r) ⇥ d 1 v r ) ⇥ d 2 w r in two modes d 1 and d 2 , d 1 <d 2 , with matrices V 2 R I d 1 ⇥R and W 2 R I d 2 ⇥R as z (:) i 1 ,...,i d 1 -1 ,1,i d 1 +1 ,...,i d 2 -1 ,1,i d 2 +1 ,...,i N = I d 2 X j 2 =1 W(j, :) ⇤ y (:) i 1 ,...,i d 1 -1 ,1,i d 1 +1 ,...,i d 2 -1 ,j 2 ,i d 2 +1 ,...,i N = (I d 1 ,I d 2 ) X (j 1 ,j 2 )=(1,1) V(j 1 , :) ⇤ W(j 2 , :) ⇤ x (:) i 1 ,...,i d 1 -1 ,j 1 ,i d 1 +1 ,...,i d 2 -1 ,j 2 ,i d 2 +1 ,...,i N . (3.4)
The formula similarly generalizes to d dimensions for d>2 where each dimension adds another Hadamard product with a corresponding matrix row. This "thick" mode of operation provides a significant performance gain thanks to the increase in locality. Also, performing TTVs in this manner in CP-ALS effectively reduces the RN dlog N e TTVs required in Theorem 1 to N dlog N e TTMV calls within an iteration. In our approach, for each node t 2 T , we store a single set I t containing index tuples of the form (i 1 ,...,i N ) representing the nonzeros x (r) i 1 ,...,i N 2 T r (t) for all r 2 N R . Also, for each such index tuple we hold a vector of size R corresponding to the values of this nonzero in each one of R tensors, and we denote this value vector as V t (i 1 ,...,i N ).

The sparsity of input tensor X determines the sparsity of the tensors of the nodes of the dimension tree. Computing this sparsity structure for each TTV can get very expensive, and is redundant. As X stays fixed, we can compute the sparsity of each tensor once, and reuse it in all CP-ALS iterations. To this end, we need a data structure that can express this sparsity, while exposing parallelism to update the tensor elements in numerical TTV computations. We now describe computing this data structure in what we call the symbolic TTV step.

Symbolic TTV

For simplicity, we proceed with describing how to perform the symbolic TTV using the same 4dimensional tensor X and BBDT. However, the approach naturally generalizes to any N -dimensional tensor and dimension tree.

The first information we need is the set of nonzero indices I t for each node t in a dimension tree, which we determine as follows. As mentioned, the two children t 12 and t 34 of t 1234 have the corresponding tensors T r (t 12 )=T r (t 1234 ) ⇥ 3 u

(3) r ⇥ 4 u (4) r and T r (t 34 )=T r (t 1234 ) ⇥ 1 u (1) r ⇥ 2 u (2)
r , respectively. Using (3.4), the nonzero indices of T r (t 12 ) and T r (t 34 ) take the form (i, j, 1, 1) and (1, 1,k,l), respectively (in practice, tensor indices in the multiplied dimensions are always 1; hence we omit storing them), and such a nonzero index exist in these tensors only if there exists a nonzero x i,j,k,l 2 T r (t 1234 ). Determining the set I t 12 (or I t 34 ) can simply be done by starting with a list I t 1234 of tuples, then replacing each tuple (i, j, k, l) in the list with the tuple (i, j) (or with (k, l) for I t 34 ). This list may contain duplicates, which can efficiently be eliminated by sorting. Once the set of nonzeros for t 12 (or t 34 ) is determined, we proceed to detect the nonzero patterns of its children t 1 and t 2 (or, t 3 and t 4 ).

To be able to carry out the numerical calculations (3.3) and (3.4) for each nonzero at a node t, we need to identify the set of nonzeros of the parent node P (t)'s tensors that contribute to this nonzero. Specifically, at t 12 , for each nonzero index (i, j) 2I t 12 we need to bookmark all nonzero index tuples of t 1234 of the form (i, j, k, l) 2I t 1234 , as it is this set of nonzeros of T r (t 1234 ) that contribute to the nonzero of T r (t 12 ) with index (i, j, 1, 1) in (3.4). Therefore, for each such index tuple of t 12 we need a reduction set R t 12 (i, j) which keeps track of all such index tuples of the parent. We determine these sets simultaneously with I t . In Figure 3.1, we illustrate a sample BBDT for a 4-dimensional sparse tensor with I t , shown with arrays, and R t , shown using arrows.

Next, we provide the following theorem to help us analyze the computational and the memory cost of symbolic TTV using a BBDT for sparse tensors.

Theorem 2. Let X be an N -mode sparse tensor. The total number of index arrays in a BBDT of X is at most N (dlog N e + 1).

Proof. Each node t in the dimension tree holds an index array for each mode in its mode set µ(t). As stated in the proof of Theorem 1, the total size of mode sets at each tree level is at most N . Therefore, the total number of index arrays cannot exceed (dlog N e + 1)N in a BBDT. Theorem 2 shows that the storage requirement for the tensor indices of a BBDT cannot exceed (dlog N e + 1)-times the size of the original tensor in the coordinate format (which has N index arrays of size nnz(X )), yielding the overall worst-case memory cost O(nnz(X )N (dlog N e + 1)). Indeed, this is a pessimistic estimate for real-world tensors, as with significant index overlap in non-root tree tensors, the total memory cost could reduce to as low as O(nnz(X )N ). This renders the approach very suitable for higher dimensional tensors. In computing the symbolic TTV, we sort |µ(t)| index arrays for each node t 2 T . In addition, for each node t of a BBDT, we sort an extra array to determine the reduction set R t . In the worst case, each array can have up to nnz(X ) elements. Therefore, combining with the number of index arrays as given in Theorem 2, the overall worst case cost of sorting becomes O(N (dlog N e + 1) + 2N -1)nnz(X ) log(nnz(X )) = O(N log N nnz(X ) log(nnz(X )). We note, however, that both the total index array size and sorting cost are pessimistic overestimates, since the nonzero structure of real-world tensors exhibits significant locality in indices. For example, on two tensors from our experiments (Delicious and Flickr), we observed a reduction factor of 2.57 and 5.5 in the number of nonzeros of the children of the root of the BBDT. Consequently, the number of nonzeros in a node's tensors reduces dramatically as we approach to the leaves. In comparison, existing approaches [START_REF]Tensor-matrix products with a compressed sparse tensor[END_REF] sort the original tensor once with a cost of O(N nnz(X ) log(nnz(X ))) at the expense of computing TTVs from scratch in each CP-ALS iteration.

Symbolic TTV is a one-time computation whose cost is amortized. Normally, choosing an appropriate rank R for a sparse tensor X requires several executions of CP-ALS. Also, CP-ALS is known to be sensitive to the initialization of factor matrices; therefore, it is often executed with multiple initializations [START_REF]Tensor decompositions and applications[END_REF]. In all of these use cases, the tensor X is fixed; therefore, the symbolic TTV is required only once. Moreover, CP-ALS usually have a number of iterations which involve many costly numeric TTV calls. As a result, the cost of the subsequent numeric TTV calls over many iterations and many CP-ALS executions easily amortizes that of this symbolic preprocessing. Nevertheless, in case of need, this step can efficiently be parallelized in multiple ways. First, symbolic TTV is essentially a sorting of multiple index arrays; hence, one can use parallel sorting methods. Second, the BBDT structure naturally exposes a coarser level of parallelism; once a node's symbolic TTV is computed, one can proceed with those of its children in parallel, and process the whole tree in this way. Finally, in a distributed memory setting where we partition the tensor to multiple processes, each process can perform the symbolic TTV on its local tensor in parallel. We benefit only from this parallelism in our implementation.

After symbolic TTV is performed, index arrays of all nodes in the tree stay fixed and are kept throughout CP-ALS iterations. However, at each subiteration n, only the values of tensors which are necessary to compute DTREE-TTV(L n ) are kept. The following theorem provides an upper bound on the number of such value arrays, which gives an upper bound on the memory usage for tensor values.

Theorem 3. For an N -mode tensor X , the total number of tree nodes whose values are not destroyed is at most dlog N e at any instant of Algorithm 5 using a BBDT.

Proof. Note that at the beginning of the nth subiteration of Algorithm 5, the tensors of each node t 2 T involving n in µ 0 (t) are destroyed at Line 8. These are exactly the tensors that do not lie in the path from the leaf L n to the root, as they do not involve n in their mode set µ. The TTMV result for L n depends only on the nodes on this path from L n to the root; therefore, at the end of the nth subiteration, only the tensors of the nodes on this path will be computed using Algorithm 4. As this path length cannot exceed dlog N e in a BBDT, the number of nodes whose values are not destroyed cannot exceed Rdlog N e at any instant of Algorithm 5.

CHAPTER 3. PARALLEL CP DECOMPOSITION OF SPARSE TENSORS USING DIMENSION TREES

Theorem 3 puts a nice upper bound of dlog N e on the maximum number of allocated value arrays, thus on the maximum memory utilization due to tensor values, in DTREE-CP-ALS. Note that in the worst case, each tensor value array may have up to nnz(X ) elements, requiring nnz(X )Rdlog N e memory to store its values in total. Indeed, having significant nonzero index overlaps after TTVs may significantly reduce this cost, as the size of a value array is proportional to the number of nonzero elements of the associated node's tensors.

Parallel CP-ALS for sparse tensors using dimension trees

We first present shared memory parallel algorithms involving efficient parallelization of the dimension tree-based TTMVs in Section 3.2.1. Later, in Section 3.2.2, we present distributed memory parallel algorithms that use this shared memory parallelization.

Shared memory parallelism

Algorithm 6 SMP-DTREE-TTV

Input: t: A dimension tree node/tensor Output: Numerical values V t are computed. 1: if EXISTS(V t ) then 2:
return I Numerical values V t are already computed. 3: SMP-DTREE-TTV(P (t))

I Compute the parent's values V P (t) first. 4: parallel for all (i 1 ,...,i N ) 2I t do I Compute each V t (i 1 ,...,i N ) in parallel.

5:

V t (i 1 ,...,i N ) zeros(1,R) I Initialize with a zero vector of size 1 ⇥ R.

6:

for all (j 1 ,...,j N ) 2R(i 1 ,...,i N ) do I Perform updates using elements in R.

7:

r V P (t) (j 1 ,...,j N ) 8:
for all d 2 δ(t) do I Multiply the vector r with corresponding matrix rows. 9:

r r ⇤ U (d) (j d , :) 10: V t (i 1 ,...,i N ) V t (i 1 ,...,i N )+r I Do the update due to element (j 1 ,...,j N ).
After forming the dimension tree T for a tensor X with symbolic structures I t and R t for all tree nodes, we can perform numeric TTMV computations in parallel. In Algorithm 6, we provide the shared memory parallel TTMV algorithm, called SMP-DTREE-TTV, for a node t of a dimension tree. The goal of SMP-DTREE-TTV is to compute the tensor values V t for a given node t. Similar to Algorithm 4, it starts by checking if V t is already computed, and returns immediately in that case. Otherwise, it calls SMP-DTREE-TTV on the parent node P (t) to make sure that parent's tensor values V P (t) are available. Once V P (t) is ready, the algorithm proceeds with computing V t for each nonzero index (i 1 ,...,i N ) 2I t . As for each index (i 1 ,...,i N ) 2I t , the reduction set is defined during the symbolic TTV, the V t can independently be updated in parallel. In performing this update, for each element V P (t) (j 1 ,...,j N ) of the parent, the algorithm multiplies this vector with the rows of the corresponding matrices of the TTMV in δ(t) of t, then adds it to V t (i 1 ,...,i N ).

For shared memory parallel CP-ALS, we replace Line 9 with the call SMP-DTREE-TTV(L n ) in Algorithm 5. The parallelization of the rest of the computations is trivial. In computing the matrices W (n) and U (n) at Lines 3, 13 and 15, we use parallel dense BLAS kernels. Computing the matrix H (n) at Line 12 and normalizing the columns of U (n) are embarrassingly parallel element-wise matrix operations. We skip the details of the parallel convergence check whose cost is negligible.

Distributed memory parallelism

Parallelizing CP-ALS in a distributed memory setting involves defining unit parallel tasks, data elements, and their interdependencies. Following to this definition, we partition and distribute tensor elements as well as factor matrices to all available processes. We provide a fine-grain and a mediumgrain parallel task model together with the associated distributed memory parallel algorithms. We start the discussion with the following straightforward lemma that enables us to distribute tensor nonzeros for parallelization.

Lemma 1 (Distributive property of tensor-times-vector multiplies). Let X , Y, and Z be tensors in

R I 1 ⇥•••⇥I N with X = Y + Z. Then, for any n 2 N N and u 2 R In X ⇥ n u = Y ⇥ n u + Z ⇥ n u holds.
Proof. Using (2.1) we express the element-wise result of X ⇥ n u as

(X ⇥ n u) i 1 ,...,1,...,i N = In X j=1 u j (x i 1 ,...,j,...,i N -z i 1 ,...,j,...,i N + z i 1 ,...,j,...,i N ) = In X j=1 u j (x i 1 ,...,j,...,i N -z i 1 ,...,j,...,i N )+ In X j=1 u j z i 1 ,...,j,...,i N =(Y ⇥ n u) i 1 ,...,1,...,i N +(Z ⇥ n u) i 1 ,...,1,...,i N .
By extending the previous lemma to P summands and all but one mode TTV, we obtain the next corollary.

Corollary 1. Let X and X 1 ,...,X P be tensors in

R I 1 ⇥•••⇥I N with P P i=1 X i = X . Then, for any n 2 N N and u (i) 2 R I i for i 2 N N \{n}, we obtain X ⇥ i6 =n u (i) = X 1 ⇥ i6 =n u (i) +•••+X P ⇥ i6 =n u (i) .
Proof. Multiplying the tensors X and X 1 ,...,X P in any mode n 0 6 = n in the equation gives tensors

X 0 = X ⇥ n 0 u (n 0 ) and X 0 i = X i ⇥ n 0 u (n 0 )
, and X 0 = P P i=1 X 0 i holds due to the distributive property. The same process is repeated in the remaining modes to obtain the desired result.

Fine-grain parallelism

Corollary 1 allows us to partition the tensor X in the sum form X 1 + •••+ X P for any P>1, then perform TTMVs on each tensor partition X p independently, finally sum up these results to obtain the TTMV result for X multiplied in N -1 modes. As X is sparse, an intuitive way to achieve this decomposition is by partitioning its nonzeros to P tensors where P is the number of available distributed processes. This way, for any dimension n, we can perform the TTMV of X p with the columns of the set of factor matrices {U (1) ,...,U (N ) } in all modes except n. This yields a "local" matrix M (n) p at each process p, and all these local matrices must subsequently be "assembled" by summing up their rows corresponding to the same row indices. In order to perform this assembling of rows, we also partition the rows of matrices M (n) so that each row is "owned" by a process that is responsible for holding the final row sum. We represent this partition with a vector σ (n) 2 R In where σ (n) (i)=p implies that the final value of M (n) (i, :) resides at the process p. We assume the same partition σ (n) on the corresponding factor matrices U (n) , as this enables each process to compute the rows of U (n) that it owns using the rows of M (n) belonging to that process without incurring any communication.

Algorithm 7 DMP-DTREE-CP-ALS: Dimension tree-based CP-ALS algorithm

Input: X p : An N -mode tensor

I (n) p : The index set with elements i 2 I (n) p where σ (n) (i)=p F (n) p :
The set of all unique indices of X p in mode n U (1) (F

p , :),...,U (N ) (F (N ) p

, :): Distributed initial matrices (rows needed by process p) R: The rank of CP decomposition Output: [λ; U (1) ,...,U (N ) ]:

Rank-R CP decomposition of X p with distributed U (n) 1: T = {L 1 ,...,L N , I 1 ,...} CONSTRUCT-DIMENSION-TREE(X p ) 2: for n =2...N do 3: W (n) ALL-REDUCE([U (n) (I (n) p , :)] T U (n) (I (n) p , :)) 4: repeat 5:
for n =1...N do

6:

for all t 2 T do I Invalidate tree tensors that are multiplied in mode n.

7:

if n 2 µ 0 (t) then DESTROY(V t ) I Destroy all tensors that are multiplied by U (n) . 9:

SMP-DTREE-TTV(L n ) I Perform the TTVs for the leaf node tensors.

10:

M (n) (F (n) p , :) V Ln (:) I Form M (n) using leaf tensors (done implicitly). 11: COMM-FACTOR-MATRIX(M (n) , "fold",F (d) p ,I (d)
p ,σ (d) ) I Assemble M (n) (I (n)
p , :).

12:

H (n) ⇤ i6 =n W (i)
13:

U (n) (I (n) p , :) M (n) (I (n) p , :)H (n) † 14: λ COLUMN-NORMALIZE(U (n) ) 15:
COMM-FACTOR-MATRIX(U (n) , "expand",F

p ,I

p ,σ (d) ) I Send/Receive U (n) .

16:

W (n) ALL-REDUCE([U (n) (I (n) p , :)] T U (n) (I (n)
p , :)) 17: until converge is achieved or the maximum number of iterations is reached.

This approach amounts to a fine-grain parallelism where each fine-grain computational task corresponds to performing TTMV operations due to a nonzero element x i 1 ,...,i N 2 X . Specifically, according to (3.4), the process p needs the matrix rows U (1) (i 1 , :),...,U (N ) (i N , :) for each nonzero x i 1 ,...,i N in its local tensor X p in order to perform its local TTMVs. For each dimension n, we represent the union of all these "required" row indices for the process p by F (n) p . Similarly, we represent the set of "owned" rows by the process p by I (n) p . In this situation, the set

F (n) p \I (n)
p correspond to the rows of M (n) for which the process p generates a partial TTMV result, which need to be sent to their owner processes. Equally, it represents the set of rows of U (n) that are not owned by the process p and are needed in its local TTMVs according to (3.4). These rows of U (n) are similarly to be received from their owners in order to carry out the TTMVs at process p. Hence, a "good" partition in general involves a significant overlap of F p of rows of each factor matrix U (n) that it needs for its local computations. Similar to Algorithm 5, at Line 1 we start by forming the dimension tree for the local tensor X p . We then compute the matrices W (n) corresponding to U (n) T U (n) using the initial factor matrices. We do this step in parallel in which each process computes the local contribution

[U (n) (I (n) p , :)] T U (n) (I (n)
p , :) due to its owned rows. Afterwards, we perform an ALL-REDUCE communication to sum up these local results to obtain a copy of W (n) at each process. The cost of this communication is typically negligible as W (n) is a small matrix of size R ⇥ R. The main CP-ALS subiteration for mode n begins with destroying tensors in the tree that will become invalid after updating U (n) . Next, we perform SMP-DTREE-TTV on the leaf node L n , and obtain the "local" matrix for all i 2 F p \ I p do I Send non-owned rows to their owners.

3:

Send M(i, :) to the process σ(i).

4:

for all i 2 I p do I Gather all partial results of owned rows together.

5:

Receive and sum up all partial results for M(i, :). 6: else if comm = "expand" then for all i 2 I p do I Send owned rows to all processes in need.

8:

Send M(i, :) to the all processes p 0 with i 2 F p 0 .

9:

for all i 2 F p \ I p do I Receive rows that are needed locally. Receive M(i, :) from the process σ(i).

M (n) . Then, the partial results for the rows of M (n) are communicated to be assembled at their owner processes. We name this as the fold communication step following the convention from the fine-grain parallel sparse matrix computations. Afterwards, we form the matrix H (n) locally at each process p in order to compute the owned part

U (n) (I (n) p , :) using the recently assembled M (n) (I (n)
p , :). Once the new distributed U (n) is computed, we normalize it column-wise and obtain the vector λ of norms. The computational and the communication costs of this step are negligible. The new U (n) is finalized after the normalization, and we then perform an expand communication step in which we send the rows of U (n) from the owner processes to all others in need. This is essentially the inverse of the fold communication step in the sense that each process p that sends a partial row result of M (n) (i, :) to another process q in the fold step receives the final result for the corresponding row U (n) (i :) from the process q in the expand communication. Finally, we update the matrix W (n) using the new U (n) in parallel.

The expand and the fold communications at Lines 11 and 15 constitute the most expensive communication steps. We outline these communications in Algorithm 8. In the expand communication, the process p sends the partial results for the set F p \ I p of rows to their owner processes, while similarly receiving all partial results for its set I p of owned rows and summing them up. Symmetrically, in the fold communication, the process p sends the rows with indices I p , and receives the rows with indices F p \ I p . The exact set of row indices that needs to be communicated in fold and expand steps depends on the partitioning of X and the factor matrices. As this partition does not change once determined, the communicated rows between p and q stays the same in CP-ALS iterations. Therefore, in our implementation we determine this row set once outside the main CP-ALS iteration, and reuse it at each iteration.

CHAPTER 3. PARALLEL CP DECOMPOSITION OF SPARSE TENSORS USING DIMENSION TREES

Medium-grain parallelism

For an N -mode tensor X 2 R I 1 ⇥•••⇥I N and using P = Q N i=1 P i processes, the medium-grain decomposition imposes a partition with P 1 ⇥•••⇥P N Cartesian topology on the dimensions of X . Specifically, for each dimension n, the index set N In is partitioned into P n sets

I (n) 1 ,...,I (n)
Pn . With this partition, the process with the index (p 1 ,...,p N ) 2

P 1 ⇥•••⇥P N gets X (I (1) p 1 ,...,I (n)
p N ) as its local tensor. Each factor matrix U (n) is also partitioned following this topology where the set of rows

U (n) (I (n) j , :
) is owned by the processes with index (p 1 ,...,p N ) where j 2 N P N and p n = j, even though these rows are to be further partitioned among the processes having p n = j. As a result, one advantage of the medium-grain partition is that only the processes with p n = j need to communicate with each other in mode n. This does not necessarily reduce the volume of communication, but it can reduce the number of communicated messages by a factor of P n in the nth dimension.

One can design an algorithm specifically for the medium-grain decomposition [START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF]. However, using the fine-grain algorithm on a medium-grain partition effectively provides a medium-grain algorithm. For this reason, we do not need nor provide a separate algorithm for the medium-grain task model, and use the fine-grain algorithm with a proper medium-grain partition instead, which equally benefits from the topology.

Partitioning

The distributed memory algorithms that we described require partitioning the data and the computations, as in any distributed memory algorithm. In order to reason about their computational load balance and communication cost, we use hypergraph models. Once the models are built, different hypergraph partitioning methods can be used to partition the data and the computations. We discuss a few partitioning alternatives.

Partitioning for the fine-grain parallelism

We propose a hypergraph model to capture the computational load and the communication volume of the fine-grain parallelization given in Algorithm 7. For the simplicity of the discussion, we present the model for a 3rd order tensor X 2 R I 1 ⇥I 2 ⇥I 3 and factor matrices U (1) 2 R I 1 ⇥R , U (2) 2 R I 2 ⇥R , and U (3) 2 R I 3 ⇥R . For these inputs, we construct a hypergraph H =( V, E) with the vertex set V and the hyperedge set E. The generalization of the model to higher order tensors should be clear from this construction.

The vertex set V = V (1) [V (2) [V (3) [V (X ) of the hypergraph involves four types of vertices. The first three types correspond to the rows of the matrices U (1) , U (2) , and U (3) . In particular, we have vertices v

(1) 3) for k 2 N I 3 . These vertices represent the "ownership" of the corresponding matrix rows, and we assign unit weight to each such vertex. The fourth type of vertices are denoted by v (X ) i,j,k , which we define for each nonzero x i,j,k 2 X . This vertex type relates to the number of operations performed in TTMV due to the nonzero element x i,j,k 2 X in using (3.4) in all modes. In the N -dimensional case, this includes up to N vector Hadamard products involving the value of the nonzero x i,j,k , and the corresponding matrix rows. The exact number of performed Hadamard products depends on how nonzero indices coincide as TTVs are carried out, and cannot be determined before a partitioning takes place. In our earlier work [START_REF] Kaya A N D B | Scalable sparse tensor decompositions in distributed memory systems[END_REF], this cost was exactly N Hadamard products per nonzero, as the MTTKRPs were computed without reusing partial results and without index compression after each TTMV. In the current case, we assign a cost of N to each vertex v (X ) i,j,k to represent an upper bound on the computational cost, and expect this to lead to a good load balance in practice. With these vertex definitions, one can use multi-constraint partitioning (2.7) with one constraint per vertex type. In this case, the first, the second, and the third types have unit weights in the first, second, and third constraints, respectively, and zero weight in all other constraints. The fourth vertex type also gets a unit weight (N , or equivalently, 1) in the fourth constraint, and zero weight for others. Here, balancing the first three constraints corresponds to balancing the number of matrix rows at each process (which provides the memory balance as well as the computational balance in dense matrix operations), whereas balancing the fourth type corresponds to balancing the computational load due to TTMVs.

i 2 V (1) for i 2 N I 1 , v (2) j 2 V (2) for j 2 N I 2 , and v (3) k 2 V (
v 1 (1) v 1 (2) v 1 (3) v 1,2,3 (Χ) v 2 (1) v 2 (2) v 2 (3) v 2,3,1 (Χ) v v 3 (2) v 3 (1) v (Χ) n 2 (1) n 2 (2) n 2 (3) n 1 (1) n 1 (2) n 1 (3) n 3 (3) n 3 (2) n 3 (1) 3,1,2 3 (3)
As TTMVs are carried out using (3.4), data dependencies to the rows of U (1) (i, :), U (2) (j, :), and U (3) (k, :) take place when performing Hadamard products due to each nonzero x i,j,k . We introduce three types of hyperedges in E = E (1) [ E (2) [ E (3) to represent these dependencies as follows:

E (1) contains a hyperedge n (1) i for each matrix row U (1) (i, :), E (2) contains a hyperedge n (2)
j for each row U (2) (j, :), and E (3) contains a hyperedge n

(3) k for each row U (3) (k, :). Initially, n (1) i , n (2) j and n (3) k contain the corresponding vertices v (1) i , v (2)
j , and v

(3) k , as the owner of a matrix row has a dependency to it by default. In computing the MTTKRP using (3.4), each nonzero x i,j,k requires access to U (1) (i, :), U (2) (j, :), and U (3) (k, :). Therefore, we add the vertex v (X ) i,j,k to the hyperedges n Consider now a P -way partition of the vertices of H =( V, E) where each part is associated with a unique process to obtain a P -way parallel execution of Algorithm 7. We consider the first subiteration of Algorithm 7 that updates U (1) , and assume that each process already has all data

CHAPTER 3. PARALLEL CP DECOMPOSITION OF SPARSE TENSORS USING DIMENSION TREES

elements to carry out the local TTMVs at Line 9. Now suppose that the nonzero x i,j,k is owned by the process p and the matrix row U (1) (i, :) is owned by the process q. Then, the process p computes a partial result for M (1) (i, :) which needs to be sent to the process q at Line 3 of Algorithm 8. By construction of the hypergraph, we have v

(1) i 2 n (1)
i which resides at the process q, and due to the nonzero x i,j,k we have v

(X ) i,j,k 2 n (1)
i which resides at the process p; therefore, this communication is accurately represented in the connectivity λ n (1) i of the hyperedge n (1) i . In general, the hyperedge n

(1) i incurs λ n (1)
i -1 messages to transfer the partial results for the matrix row M (1) (i, :) to the process q at Line 3. Therefore, the connectivity-1 cutsize metric (2.6) over the hyperedges exactly encodes the total volume of messages sent at Line 3, if we set c[•]=R. Since the send operations at Line 8 are duals of the send operations at Line 3, the total volume of messages sent at Line 8 for the first mode is also equal to this number. By extending this reasoning to all other modes, we obtain that the cumulative (over all modes) volume of communication in one iteration of Algorithm 7 equals to the connectivity-1 cut-size metric. As the communication due to each mode take place in different stages, one might alternatively use a multi-objective hypergraph model to minimize the communication volume due to each mode (or equivalently, hyperedge type) independently.

As discussed above, the proper model for partitioning the data and the computations for the finegrain parallelism calls for a multi-constraint and a multi-objective partitioning formulation to achieve the load balance and minimize the communication cost with a single call to a hypergraph partitioning routine. Since these formulations are expensive, we follow a two-step approach. In the first step, we partition only the nonzeros of the tensor on the hypergraph H =( V (X ) ,E) using just one load constraint due to the vertices in V (X ) , and we thereby avoid multi-constraint partitioning. We also avoid multi-objective partitioning by treating all hyperedge types as the same, and thereby aim to minimize the total communication volume across all dimensions, which works well in practice. Once the nonzero partitioning is settled, we partition the rows of the factor matrices in a way to balance the communication, which is not achievable using standard partitioning tools.

We now discuss three methods for partitioning the described hypergraph.

Random: This approach visits the vertices of the hypergraph and assigns each visited vertex to a part chosen uniformly at random. It is expected to balance the TTMV work assigned to each process while ignoring the cost of communication. We use random partitioning only as a "worst case" point of reference for other methods.

Standard: In this standard approach, we feed the hypergraph to a standard hypergraph partitioning tool to obtain balance on the number of tensor nonzeros and the amount of TTMV work assigned to a process, while minimizing the communication volume. This approach promises significant reductions in communication cost with respect to the others, yet imposes high computational and memory requirements.

Label propagation-like: Given that the standard partitioning approach is too costly in practice, we developed a fast hypergraph partitioning heuristic which has reasonable memory and computational costs. The method is based on the balanced label-propagation algorithm [START_REF] Slota | PuLP: Scalable multi-objective multi-constraint partitioning for small-world networks[END_REF][START_REF] Backstrom | Balanced label propagation for partitioning massive graphs[END_REF], and includes some additional adaptations to handle hypergraphs [START_REF] Buurlage | Self-improving sparse matrix partitioning and bulk-synchronous pseudostreaming[END_REF][START_REF] Enne | Label propagation for hypergraph partitioning[END_REF]. The heuristic starts with an initial assignment of vertices to parts, and then proceeds with multiple passes over the hypergraph. At each pass, the vertices are visited in an order, and are possibly moved to other parts in order to reduce the cutsize while respecting the balance constraints.

For the heuristic to be efficient on hypergraphs, some adaptations are needed. Each pass involves two types of updates. In the first step, each hyperedge chooses a "preferred part" by considering the current part of its vertices. Next, each vertex updates its part according to the preferred parts of the hyperedges that include the vertex. In both steps, the most dominant part index is chosen for the update. The heuristic runs in linear time on the size of the hypergraph per iteration, and requires a memory of 2|V | + |E| +4P . Running the algorithm for a few iterations provides reasonably good partitions. This basic algorithm can have many variants. In one variant, we visit the vertices in an order imposed by an increasing ordering by size of the hyperedges. This variant has an overhead of sorting the hyperedges. In another variant, we reweigh the preference of a hyperedge of size s by the multiplier (1 -1 P ) s-1 . This last variant has a memory overhead for storing the weights for efficiency purposes; for each size s, the value (1 -1 P ) s-1 is needed.

Partitioning for the medium-grain parallelism

We propose a hypergraph model to capture the communication requirements of the medium-grain algorithm. It is an adaptation of the checkerboard partitioning designed for matrices [START_REF] Çatalyürek A N D C | A hypergraph-partitioning approach for coarse-grain decomposition[END_REF][START_REF] Çatalyürek | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF] to tensors.

The partitioning heuristic for the medium-grain parallelism proceeds in N steps for an N -dimensional tensor, and for simplicity, we discuss the three dimensional case in which we partition the tensor for a P 1 ⇥ P 2 ⇥ P 3 process topology where P = P 1 P 2 P 3 . Each step partitions the tensor slices in the corresponding mode. In the first step, a hypergraph H 1 =(V (1) ,E (1) ) is built containing a vertex v

(1) i

for each tensor slice X (i, :, :) in the first mode, and hyperedges n

(2) j and n

k for each index j 2 N I 2 and k 2 N I 3 in the second and the third modes. Next, each vertex v

(1) i is included in the hyperedge n (2) j if X (i, j, :) contains a nonzero; similarly, v (1)
i is included in the hyperedge n (3) k if X (i, :,k) is not empty. The resulting hypergraph is shown in Figure 3.3a for the previous sample tensor. The vertices have weights corresponding to the number of nonzeros in the associated slices. This hypergraph is then partitioned into P 1 parts. A similar hypergraph is built in the second step for the second mode, as shown in Figure 3.3b. In this hypergraph, each vertex v (2) j has P 1 weights designating the number of nonzeros of the jth slice of the second mode of X that are contained in the parts identified in the first step. After having the second hypergraph partitioned into P 2 parts, in the third step, the hypergraph associated with the third mode is built similarly, as shown in Figure 3.3c, and partitioned into P 3 parts. Here, each vertex v (3) k has P 1 P 2 weights representing the number of nonzeros of the kth slice of the third mode of X that are contained in the Cartesian product of the parts identified in the first two steps. The resulting partitions are used to assign tensor nonzeros to the P processes. The underlying P 1 ⇥ P 2 ⇥ P 3 topology confines the communication in each mode to a 2D grid of processors, e.g., in the first mode, each processor communicates with P 2 P 3 processors. This bounds the maximum number of messages sent per process by P 2 P 3 , but does not does not necessarily decrease the actual volume of communication.

Nevertheless, as P can be large, the number of constraints P 1 and P 1 P 2 can similarly get large, in which case the state of the art partitioners are known to not perform well both in terms of partition quality and speed. For higher dimensional tensors, this situation only gets worse. That is why explicit communication reduction using hypergraph partitioning for the medium-grain algorithm is not feasible in practice. Hence, we use the partitioning heuristic by Smith and Karypis [START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF] to partition medium-grain hypergraphs for load balance, and to expect a communication reduction due to partition topology indirectly. We also determine the partition topology by choosing P 1 , P 2 , and P 3 proportional to the tensor dimensions I 1 , I 2 , and I 3 .

CHAPTER 3. PARALLEL CP DECOMPOSITION OF SPARSE TENSORS USING DIMENSION TREES

v 1 (1) v 2 (1) v 3 (1)
n 2

(2)

n 3 (3) n 3 (2) n 1 (3) n 1 (2) n 2 (3) (a) v 1 (2) v 2 (2) v 3 (2) n 3 (1) n 2 (3) n 1 (1) n 3 (3)
n 2

(1)

n 1 (3) (b) v 1 (3) v 2 (3) v 3 (3) n 2 (1) n 3 (2) n 3 (1) n 1 (2) n 1 (1) n 2 (2)
(c)

Mode partitioning

Once the nonzero partitioning is obtained for the given fine-or medium-grain parallelism, we proceed with partitioning the mode indices (or, equivalently, the rows of the factor matrices) using a similar heuristic common in similar work [START_REF] Kaya A N D B | Scalable sparse tensor decompositions in distributed memory systems[END_REF][START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF]. For each matrix row i in dimension n, we identify the processes that have a data dependency to that row. These are exactly the processes which have at least one nonzero with index i in the nth dimension. Next, all row indices are sorted in increasing order of the number of dependent processes. Finally, each row is greedily assigned to the process having the minimum total communication volume among all processes dependent to that row.

Related work

There has been many recent advances in the efficient computations of tensor factorizations in general, and CP decomposition in particular. We briefly mention these here and refer the reader to the original sources for details. In [6], Bader and Kolda show how to efficiently carry out MTTKRP as well as other fundamental tensor operations on sparse tensors in MATLAB. GigaTensor [START_REF] Kang | GigaTensor: Scaling tensor analysis up by 100 times -Algorithms and discoveries[END_REF] is a parallel implementation of CP-ALS using the Map-Reduce framework. DFacTo [START_REF] Choi | DFacTo: Distributed factorization of tensors[END_REF] is a C++ implementation with distributed memory parallelism using MPI, and it uses a particular formulation of MTTKRP using sparse matrix-vector multiplication. SPLATT [START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF] is an efficient parallelization of MTTKRP and CP-ALS both in shared [START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF] and distributed memory environments [START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF] using OpenMP and MPI, and is implemented in C. It uses a medium-grain distributed parallelism with a Cartesian partitioning of the tensor, and generalizes this technique to the tensor completion problem [START_REF] Smith | An exploration of optimization algorithms for high performance tensor completion[END_REF]. It is the fastest publicly available CP-ALS implementation in the existing literature, and their approach translates to performing N (N -1) TTMVs in performing the MTTKRP in the main CP-ALS iteration. Karlsson et al. similarly discuss a parallel computation of the tensor completion problem using CP formulation [START_REF] Karlsson | Parallel algorithms for tensor completion in the CP format[END_REF] in which they replicate the entire factor matrix U (n) among MPI processes unlike our approach, and report scalability results only up to 100 cores. For computing the CP decomposition of dense tensors, Phan et al. [START_REF] Phan | Fast alternating LS algorithms for high order CANDECOMP/PARAFAC tensor factorizations[END_REF] propose a scheme that divides the tensor modes into two sets, precomputes the TTMVs for each mode set, and finally reuses these partial results to obtain the final MTTKRP result in each mode. This provides a factor of 2 improvement in the number of TTMVs over the traditional approach, and our dimension tree-based framework can be considered as the generalization of this ap-proach that provides a factor of N/ log N improvement. Finally, Chi and Kolda present an alternative Alternating Poisson Regression (CP-APR) algorithm for computing the CP decomposition of large scale sparse datasets in [START_REF] Chi | On tensors, sparsity, and nonnegative factorizations[END_REF].

While this work was under evaluation, another paper appeared [START_REF] Li | Model-driven sparse CP decomposition for higher-order tensors[END_REF]. Li et al. use the same idea of storing intermediate tensors but use a different formulation based on tensor times tensor multiplication and a tensor times matrix through Hadamard products for shared memory systems. The overall approach is similar to that by Phan et al. [START_REF] Phan | Fast alternating LS algorithms for high order CANDECOMP/PARAFAC tensor factorizations[END_REF], where the difference lies in the application of the method to sparse tensors and auto-tuning to better control memory use and gains in the operation counts.

Aside from CP decomposition, Baskaran et al. [START_REF] Baskaran | Efficient and scalable computations with sparse tensors[END_REF] provide a shared memory parallel implementation for the Tucker decomposition of sparse tensors. Kaya et al. [START_REF]High performance parallel algorithms for the Tucker decomposition of sparse tensors[END_REF] provide efficient shared and distributed memory parallelization of the Tucker decomposition for sparse tensors using OpenMP and MPI. Austin et al. [5] discuss a high performance distributed memory parallelization of dense Tucker factorization in the context of data compression. Finally, Perros et al. [START_REF] Perros | Sparse hierarchical Tucker factorization and its application to healthcare[END_REF] investigate an efficient computation of hierarchical Tucker decomposition for sparse tensors.

Experiments

We first investigate the way CP-ALS implementations compare using a single thread to assess the algorithmic impact of using a BBDT in the same implementation. Then, we compare these implementations using multiple threads to evaluate their shared memory parallel performance. Finally, we investigate the medium-and the fine-grain distributed memory parallel algorithms.

Dataset and environment

We experimented with five real-world tensors, namely Netflix, NELL, Flickr, Delicious, and Amazon, which are described more in detail in Section 2.6 and whose sizes are given in Table 2.1. We also used synthetic tensors having 4, 8, 16, and 32 dimensions, which are described in the same section, to measure the effect of increasing dimensionality on performance.

We provide the dimension-tree based CP-ALS implementation in our tensor factorization library called HYPERTENSOR. It is a C++11 implementation providing shared and distributed memory parallelism through OpenMP and MPI libraries. We compared our code against SPLATT v1.1.1 [START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF], a C code with OpenMP and MPI parallelizations.

For shared memory experiments, we used Grunch environment described in Section 2.7 using 1 and 14 threads (single socket). All codes were compiled with gcc/g++-5.3.0 using OpenMP directives and compiler options -O3, -ffast-math, -funroll-loops, -ftree-vectorize, -fstrict-aliasing on this shared memory system. Distributed memory experiments were done on Blue Gene/Q using up to 256 nodes (4096 cores). Each core of Blue Gene/Q can handle an arithmetic and a memory operation simultaneously; therefore, we assigned 32 threads per node (2 threads per core) for better performance. On this system, all codes were compiled using Clang C++ compiler (version 3.5.2) with IBM MPI wrapper using the same optimization flags, and linked against IBM ESSL library for LAPACK and BLAS routines.

Shared memory experiments

We compare the shared memory performance of the dimension tree-based CP-ALS algorithm with the state of the art. We experimented with four methods called ht-tree2, ht-tree3, ht-tree, and splatt.

The ht-tree method, implemented in HYPERTENSOR, uses a full BBDT to carry out TTMVs. The ht-tree2 method is the same implementation as ht-tree except that it uses a 2-level flat dimension tree. In this tree, N leaf nodes are directly connected to the root, hence no intermediate results are generated. However, TTMVs are performed one mode at a time to benefit from the index compression to reduce the operation count. As a result, this method performs N -1 TTMVs for each mode in an iteration just as SPLATT; we thus expect comparable performance. ht-tree3 uses a three level BBDT, whose second level has two nodes holding the partial TTMV results corresponding to the first and the second half of the set of dimensions. This is analogous to the approach by [START_REF] Phan | Fast alternating LS algorithms for high order CANDECOMP/PARAFAC tensor factorizations[END_REF] for computing dense CP decompositions, but it also employs our data structure and shared memory parallelization for sparse tensors. Note that for 3-and 4-dimensional tensors, ht-tree3 and ht-tree use identical trees, thus give the same results. These results for ht-tree3 are indicated with an asterisk in the tables. Finally, splatt corresponds to the parallel CP-ALS implementation in SPLATT. We ran all algorithms for 20 iterations with the rank of approximation R = 20 (except for the sequential execution of 16and 32-dimensional random tensors, which are run for 2 iterations due to their cost), and recorded the average time spent per CP-ALS iteration. Test instances in which a method gets out of memory are indicated with a dash symbol.

Sequential execution

In Table 3.1, we give the sequential per-iteration run time of all three methods. We report the run time in seconds for splatt, and the relative speedup with respect to splatt for the other three methods.

We first note that ht-tree2 runs slightly slower than splatt on three dimensional Amazon and NELL tensors (0.99x and 0.87x), and notably slower on Netflix tensor (0.60x). This is so because SPLATT has a specially tuned implementation for 3-dimensional tensors, whereas we use a single code for all dimensions. On all higher dimensional tensors, ht-tree2 performs significantly better than splatt, up to 2.08x on Random8D, which shows the efficiency of our implementation for N -dimensional tensors even before using a BBDT. The gap between splatt and ht-tree2 narrows using Random16D, as ht-tree2 depletes the memory in one NUMA node, which is discussed more in Section 3.4.2, and starts accessing the distant memory in the other NUMA node. ht-tree2 gets out of memory using Random32D, as it stores O(N 2 ) index arrays.

We now measure the effect of dimension trees by comparing ht-tree with and ht-tree2 in Table 3.1. These two methods use the same TTMV implementation, whereas ht-tree uses a full BBDT. On Delicious, Flickr, Netflix, and NELL, ht-tree obtains 1.78x, 1.61x, 1.63x, and 1.47x speedup over ht-tree2 thanks to the BBDT. Likewise, on random tensors, we observed 1.43x, 1.91x, 3.01x speedup on tensors Random4D, Random8D, and Random16D, respectively, using ht-tree. This validates our performance expectation (Theorem 1) that as the dimensionality of the tensor increases, a BBDT results in significantly fewer TTMVs and better performance.

Comparing ht-tree with splatt similarly yields a speedup of 1.99x, 1.97x, 1.26x, 2.05x, and 3.95x on tensors Delicious, Flickr, NELL, Random4D, Random8D, and Random16D, respectively, which similarly meets our expectation of performance gain from Theorem 1. On Amazon, ht-tree was only 2% faster, whereas on Netflix, splatt was only 2% faster, which was the only instance in which splatt had a slight edge over ht-tree.

Finally, we note that the performance gap between ht-tree and ht-tree3 widens significantly as the tensor gets higher dimensional. Using Random8D, Random16D, and Random32D, ht-tree is 1.47x, 1.95x, and 2.50x faster than ht-tree3 as it incurs significantly fewer TTMVs. These results suggest that using a full BBDT is indeed the ideal method of choice.

Shared memory parallel execution

In Table 3.2, we give the run time results of all methods with shared memory parallelism using 14 threads. We first note that in HYPERTENSOR, using dimension trees consistently yields better execution times. Using ht-tree, we obtain 2.15x, 1.66x, 1.58x, 1.32x, and 1.10x speedup over ht-tree2 on Delicious, Flickr, Netflix, NELL, and Amazon tensors, respectively. For random tensors, we get 1.62x, 2.69x, and 4.23x speedup on Random4D, Random8D, and Random16D. Comparing ht-tree with splatt, we observe a speedup of 2.00x, 1.81x, 1.46x, 1.47x, 2.18x, and 3.47x on tensors Delicious, Flickr, NELL, Random4D, Random8D, and Random16D, respectively. This demonstrates that the use of a BBDT in CP-ALS computations can be effectively parallelized in a shared memory Table 3.3: Symbolic precomputation timings. We report the exact timing for splatt in seconds, and relative timing with respect to splatt for other methods, indicating the ratio at which splatt is faster than these methods this precomputation. setting, on top of significantly reducing the amount of TTMV work. On three dimensional Amazon and Netflix tensors, splatt has a slight edge over ht-tree, which runs 5% and 14% slower, respectively. Another point to note is that splatt has somewhat better parallel speedup in general (over its own sequential run time) than ht-tree2 and ht-tree. This is mostly due to the fact that TTMV is a memory-bound computation; hence, once the memory bandwidth is fully utilized, one cannot expect further speedup through multi-threading. When performing TTMVs, our implementation makes slightly more memory accesses due to extra pointer arrays involved in the dimension tree nodes, which saturates the bandwidth earlier, and in turn somewhat affects the parallel speedup. Nevertheless, using ht-tree we achieve up to 5.65x faster runs over splatt in a shared memory parallel execution.

Conformally with the sequential case, ht-tree gets significantly faster than ht-tree3 as the tensor dimensionality increases, up to 2.59x using Random32D, which demonstrates the effectiveness of using a full BBDT.

Preprocessing cost

In Table 3.3, we provide symbolic TTV costs of our methods as well as the precomputation cost of splatt for setting up its data structures. All runtimes are for a sequential execution; neither SPLATT nor HYPERTENSOR parallelizes this step in their current version. We first note that ht-tree incurs significantly less cost than ht-tree2 in all instances. This is expected, as ht-tree2 has O(N 2 ) index arrays to be sorted, whereas ht-tree has only O(N log N ) of them. For the same reason, ht-tree gets notably faster than ht-tree3 as the tensor dimensionality increases to 32. The cost is comparable between splatt and ht-tree for Delicious, Flickr, Netflix, NELL, and Amazon tensors, but splatt takes significantly less time as the dimensionality increases, up to 3.97x on Random16D, as it sorts only O(N ) arrays.

Comparing these timings with the iteration times in Table 3.1, we see that this precomputation is amortized in a few iterations, except for Netflix and NELL. In practice, CP-ALS is typically executed multiple times with different initial matrices and ranks of approximation using the same symbolic dimension tree construct, which should render this preprocessing cost less important even for these two tensors.

Memory usage

We provide the memory consumption of all methods in Table 3.4. The first column corresponds to the amount of memory used to store factor matrices, which is common to all methods. We provide the memory usage for storing index arrays for all methods, and also give the buffer memory consumption for storing the value arrays of intermediate tensors for ht-tree3 and ht-tree. We first note that ht-tree2 uses the highest amount of memory to store index arrays as expected, up to 8.31 times more than splatt using Random16D. In all tensors, the amount of index memory used by ht-tree is only slightly higher than splatt, the worst cases being Flickr and Random32D tensors for which ht-tree consumes 1.35 and 1.31 times more memory in comparison. Even though splatt uses only O(N ) index arrays, we realized upon inspecting the implementation that it uses two different representations of a tensor for faster execution, effectively doubling the memory requirements. Aside from this, we note that the amount of buffer memory used to store the values of intermediate tensors is reasonable, Random8D being the only exception in which the buffer size exceeds the index size. Finally, ht-tree3 uses more memory for index storage than ht-tree for high dimensional tensors, as it uses more index arrays. It uses less buffer, however, as it has only one level that stores intermediate tensors. All in all, we conclude upon considering these results that ht-tree provides remarkable performance improvements with a reasonable increase in the memory usage.

Distributed memory experiments

We compare the performance and the scalability of the fine-and the medium-grain parallel CP-ALS algorithms. In these experiments, we do not use SPLATT software to benchmark medium-grain parallelization for two reasons. First, we would like to compare the effect of load balance and communication cost in different algorithms using different partitionings, while isolating the effects of the efficiency of local CP-ALS computations. Since SPLATT's medium-grain implementation does not use BBDTs for local TTMVs, and is notably slower, comparing it against HYPERTENSOR's finegrain implementation which has faster local TTMVs would not be fair, nor would correctly reveal the effect of different partitioning strategies. Second, we were not able to get SPLATT to work on our distributed system despite our full efforts. Therefore, we instead performed medium-grain partitioning of tensors following the description of SPLATT's heuristic [START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF], and ran HYPERTENSOR on these partitions which incurs the same cost in terms of the communication volume and the number of messages as SPLATT, while using more efficient TTMV kernels. For local CP-ALS computations, we use the BBDT-based method ht-tree for shared memory parallelism, as it gives the best performance. This way, the experiments become more precise in terms of measuring the influence of medium-and fine-grain algorithms and associated partitionings on parallel scalability.

We investigate the performance in two tables. In Table 3.5, we give the run time results of the medium-and the fine-grain algorithms up to 256 MPI ranks using 4096 cores. Since we achieved the maximum scalability in most tensors with 256 MPI ranks and 4096 cores, the discussion is mostly confined to this case. Especially, in Table 3.6, we give the detailed load balance and communication cost metrics just for this case.

In Table 3.5, we compare the execution of Algorithm 7 with three partitioning methods. The fine-hp and fine-lb correspond, respectively, to the standard hypergraph partitioning and the labelpropagation-like heuristic of Section 3.2.2. For fine-hp, we used PaToH [START_REF] Çatalyürek A N D C | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF] with the default settings. On Amazon tensor, we could not obtain results for fine-hp as the tensor was too big for PaToH to handle. For fine-lb, we ran the three alternatives, each for three passes, and chose the partitioning with the smallest cut. The med-gd method corresponds to the medium-grain partitioning heuristic [START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF]. The fine-rd method refers to the random partitioning of the fine-grain hypergraph, which is given as a reference to illustrate the impact of a good partitioning. Due to memory constraints, we were not able to execute Algorithm 7 on a single node, as the original tensors are large. Therefore, for each tensor, we give the results starting from the minimum number of nodes needed, and for the same instance we also give the single threaded results. In passing from one core per node to 16 core per node, we see some speedup over 16 in Table 3.5; this is because we use two threads per core (see Section 3.4.1).

In order to be able to discuss the speedup results, we give the number of tensor nonzeros per part, computational load (the number of Hadamard products), communication volume, and the number of messages incurred by these three partitionings, using 256 MPI ranks in Table 3.6. For the four performance metrics, we give the maximum and the average value observed across all processes. We see in all instances except fine-hp on NELL that balancing the number of nonzeros per part gracefully translates into balancing the actual computational load. For each dataset, we use the run time of the single threaded execution of the fastest method as our baseline in computing the parallel speedup.

As seen in Table 3.5, using 256 MPI ranks on Delicious, fine-hp and fine-lb yield 2.5x and 1.9x speedup, respectively, over med-gd. We observe in Table 3.6 that this is due to better minimization of the total and the maximum communication volume. On Flickr, fine-hp is 3.9x faster than med-gd at 256 MPI ranks with 14.7x and 10.6x less total and maximum communication volume, while fine-lb shows a speedup of 3.2x over med-gd with 6.4x and 5.4x less total and maximum communication cost. In both tensors, med-gd results in about the half the communication volume of fine-rd. In overall, on Delicious fine-hp and fine-lb obtain 186x and 142x speedup using 4096 cores, whereas med-gd and fine-rd give 75x and 40x speedup for the same tensor. On Flickr, fine-hp and fine-lb similarly yield 178x and 148x speedup using 4096 cores, while med-gd and fine-rd could achieve 46x and 23x speedup. For the Delicious and Flickr tensors, while passing from 8 nodes (with 8 ⇥ 16 cores) to 256 nodes (with 256 ⇥ 16 cores), med-gd results in 5.29x and 4.32x speedup. The finehp and fine-lb result in 11.33x and 10.11x speedup for Delicious, and 11.62x and 10.50x speedup for Flickr in the same scenario. fine-rd is significantly slower than the other methods, incurring the highest communication as shown in Table 3.6.

On Netflix and NELL, med-gd yields 321x and 197x speedup using 4096 cores. fine-hp shows a comparable performance with 260x and 164x speedup, whereas fine-lb is slightly slower than fine-hp with 236x and 157x speedup. We first note that in passing from 4 nodes (with 4 ⇥ 16 cores) to 256 nodes (with 256 ⇥ 16 cores), med-gd results in 12.73x and 12.03x speedup for Netflix and NELL, respectively. The fine-hp and fine-lb partitioning result in 10.18x and 11.18x speedup for Netflix, and 9.72x and 9.32x speedup for NELL in the same scenario. fine-rd is similarly the slowest of all methods giving 67x and 62x speedup for these two tensors. Using Netflix and NELL, med-gd gets 23% and 20% faster than fine-hp, and 36% and 25% faster than fine-lb. In Table 3.6, we see that this is due to med-gd incurring smaller maximum communication volume and having less number of messages. We investigated this outcome and observed that when a tensor is long in one mode and short in all others, the communication due to the long mode dominates the overall cost, and the communication for the small modes remains negligible in comparison. On Netflix with P = 256 MPI ranks, using med-gd with P = p ⇥ q ⇥ r topology, the worst case communication volume for the first mode is upper-bounded by 480K(qr -1), as I 1 = 480K indices are distributed to p process "slices" each with qr processes. Similarly, for the second and the third modes, the worst case communication volumes are 17K(pr -1) and 2K(pq -1). In such cases, choosing a large p and smaller q and r significantly reduces the worst-case communication cost in the first mode, while the cost in other modes stays low. The medium-grain heuristic achieves this. Specifically, on Netflix, the mediumgrain heuristic chooses a grid size of 64 ⇥ 4 ⇥ 1. This advantage is lost when there are at least two Figure 3.4: Running time dissection of a parallel CP-ALS iteration using fine-hp and hp-tree scheme.

The legends "ttv", "dmat", "comm", and "fit" correspond to the time spent for TTMVs, dense matrix operations following the TTMVs, communication, and fit computation. The fit computation's time is not discernible in the plot.

long dimensions (see Delicious and Flickr), as using more processes in one long mode can increase the communication significantly in the other long modes.

On Amazon tensor, med-gd starts to lose scalability at 256 MPI ranks. In Table 3.6 we observe that this is due to the load imbalance. Amazon tensor has some relatively "dense" slices that make load balancing difficult for the medium-grain heuristic. This problem never arises in the fine-grain partitioning due to finer granularity of tasks; as a result, fine-lb runs 1.85x faster than med-gd using 256 MPI ranks. In this tensor, from 128 nodes to 256 nodes, med-gd displays a speedup of 1.08. With fine-lb, the parallel algorithm enjoys 1.86x speedup in passing from 64 nodes to 128 nodes, and 3.2x speedup in passing from 64 to 256 nodes.

In Figure 3.4, we present the dissection of the parallel run time for a CP-ALS iteration on Flickr tensor using 256 MPI ranks. We choose Flickr as representative, as it includes the highest proportion of dense matrix operations in comparison to all other tensors. Despite this fact and using a BBDT for faster TTMVs, the TTMV step still remains to be the dominant computational cost. In this figure, we first observe that the workload due to TTMV and dense matrix computations decrease with the increasing number of processes. Second, we expect in general that having more processes increases the total communication volume; yet we observe in the plot that the communication cost declines until 128 MPI ranks. This is because a good partitioning can reduce the communication volume per process (while increasing the total communication volume). At 256 MPI ranks, however, communication cost starts to increase and become the bottleneck. The fit computation takes a negligible amount of time hence is not discernible in the plot.

On Flickr tensor, the three variants of the fine-lb took 58.38, 89.66, and 65.04 seconds to partition the hypergraph, med-gd took 190 seconds, and fine-hp took 207 minutes. In all data instances, fine-lb gives good results while being a fast partitioning heuristic. fine-hp consistently provides better partitions than fine-lb in all instances, yet the partitioning cost might render it impractical to use in real-world scenarios. med-gd heuristic is only effective when the tensor nonzeros are homogeneously distributed, and the tensor has only one large dimension. One might consider reducing the communication volume on a medium-grain topology using hypergraph partitioning, yet the high number of constraints prevents this approach from being amenable. Therefore, we believe that fine-lb serves well in most practical situations.

Conclusion

We investigated an efficient computation of successive tensor-times-vector multiplication in the context of the well-known CP-ALS algorithm for sparse tensor factorization. We introduced a computational scheme using dimension trees that asymptotically reduces the computational cost of the TTMV operations for higher order tensors while using a reasonable amount of memory. Our technique provides performance benefits for lower order tensors, and gets progressively better as the dimensionality of the tensor increases in comparison to the state of the art. We proposed an effective shared memory parallelization of this method with a precomputation step in order to efficiently carry out numerical computations within the CP-ALS iterations. We introduced a fine-grain parallelization approach in the distributed memory setting, compared it against a recently proposed medium-grain variant, discussed good partitionings for both approaches, and validated these findings with experiments on real-world tensors. The proposed computational scheme can be applied to both dense and sparse tensors as well as other tensor decomposition algorithms involving successive tensor-times-vector and -matrix multiplications. We are planning to investigate this potential in our future work.

Chapter 4 Parallel Nonzero CP Decomposition of Sparse Tensors

Similar to the previous chapter, the focus of this chapter is on investigating an efficient parallel computation of CP-ALS for high dimensional big sparse tensors, which is particularly motivated by emerging big data applications [START_REF] Perros | Sparse hierarchical Tucker factorization and its application to healthcare[END_REF]. The techniques introduced in the previous chapter provide substantial benefits for high dimensional tensors, reducing the cost of a CP-ALS iteration to O(N log N ) TTMVs from O(N 2 ). We investigate a new parallel algorithm for performing MTTKRPs, which further improves computational gains for high dimensional tensors, and its effective parallelization in shared memory.

Instead of the TTV-based formulation for MTTKRPs used in the previous chapter, this chapter employs the equivalent nonzero-based formulation for MTTKRPs provided in Algorithm 2, and discusses the cost analysis accordingly. Even though the MTTKRP implementation in DFACTO and SPLATT translates to Algorithm 2, they aim to reduce the cost of this algorithm with the help index overlaps after multiplication in a dimension. To illustrate, for a 3-dimensional tensor X having nonzeros x i,j,k 1 and x i,j,k 2 , one can first compute x i,j,k 1 U (3) (k 1 , :) + x i,j,k 2 U (3) (k 2 , :), then multiply this result with U (2) (j, :) to obtain the final contribution to M (1) (i, :). Even though this approach can substantially reduce the number of Hadamard multiplications, the worst-case complexity of this approach stays the same, i.e., O(nnz(X )NR). Throughout the chapter, we let HYPERTENSOR represent our dimension tree-based algorithm described in the previous chapter, in which MTTKRP is expressed as a series of R tensor-times-vector multiply operations whose worst-case amortized cost translates to O(nnz(X ) log NR) for each tensor dimension, which provides significant performance gains as the tensor dimensionality increases in compare to SPLATT. The goal is chapter is to further reduce this cost down to O(nnz(X )R) while providing effective shared memory parallelizations.

We summarize the contributions in this chapter as following:

• We introduce a method for efficiently computing MTTKRP for high dimensional sparse tensors.

This scheme asymptotically reduces the computational cost of MTTKRP as well as a common precomputation step performed in SPLATT and the dimension tree-based method discussed in Chapter 3.

• We provide an efficient parallelization of this computational scheme that runs up to 10.5 times faster than SPLATT and 3.28 times faster than the dimension tree-based CP-ALS on a 28-core 48 CHAPTER 4. PARALLEL NONZERO CP DECOMPOSITION workstation.

The rest of the chapter is organized as follows. In Section 4.1, we describe our method for computing MTTKRP and CP-ALS, which imposes a nonzero constraint on factor matrices to reduce the computational costs. We compare the complexity of our approach with the state of the art, and discuss a carefully tuned parallelization of this approach for a shared memory NUMA architecture. Finally, we provide in Section 4.2 a comparison of sequential and parallel executions of our method with two state-of-the-art implementations.

Parallel CP decomposition using nonzero factors

Here, we first introduce our approach for efficiently performing MTTKRP with the assumption that factor matrices do not involve zeros or very small entries. When a such entry with |U (n) (i, j)| <✏ is encountered in the course of CP-ALS, we slightly "perturb" the factor matrix by replacing it with sign(U (n) (i, j))✏ where sign(x) equals to -1 if x is negative, and 1 otherwise. In practice such a perturbation is expected to have a negligible impact on the quality of solution for a sufficiently small ✏. Next, we introduce a shared memory parallelization of this scheme, and argue how to establish load balance among processes. Finally, we discuss optimization strategies for better parallel performance on a NUMA architecture.

Computing CP decomposition with nonzero factors

The cost of the algorithm in Algorithm 1 is dominated by the MTTKRP step at Line 5 that involves the multiplication of the elements of the sparse tensor X with the rows of N -1 factor matrices at each subiteration. This amounts to performing N -1 vector Hadamard products and a vector addition for each nonzero element of the tensor, as pointed out at Line 5 of Algorithm 2. Here, we present a new technique for efficiently performing this costly step with the nonzero factor matrix assumption. In this case, for each nonzero x i 1 ,...,i N 2 X , instead of performing the Hadamard product of N -1 row vectors, one can precompute a vector z i 1 ,...,i N 2 R R as z i 1 ,...,i N = ⇤ n2N N U (n) (i n , :), then perform the MTTKRP update due to this nonzero as M (n) (i n , :) += z i 1 ,...,i N ↵ U (n) (i n , :). A similar idea is also employed in the CP-APR algorithm for handling sparse tensors [START_REF] Chi | On tensors, sparsity, and nonnegative factorizations[END_REF]. Here, the cost per nonzero reduces to a single Hadamard division, which can always be performed since U (n) (i n ,j) 6 =0 . Once new U (n) is computed using M (n) at Line 7, z i 1 ,...,i N needs to be updated accordingly with the new U (n) (i n , :). This can be done by dividing it with the old value of U (n) (i n , :), then multiplying by its new value, which amounts to a single Hadamard multiplication and division. This way, instead of N -1 vector Hadamard products, we perform a Hadamard multiplication and two Hadamard divisions for each tensor element, which effectively reduces the cost of MTTKRP to O(nnz(X )R). In contrast, SPLATT and DFACTO require up to N -1 vector Hadamard products per tensor nonzero, yielding the worst-case complexity O(nnz(X )NR), and HYPERTENSOR takes O(nnz(X ) log NR) time. Therefore, our approach provides significant computational gains over all these methods particularly as X gets higher dimensional. In doing so, we use only U (n) for executing CP-ALS in dimension n, whereas SPLATT and DFACTO access all N -1 factor matrices except U (n) ; hence, our method also yields a better memory footprint.

In our method, we need to store the matrix Z which takes O(nnz(X )R) space. In contrast, HY-PERTENSOR uses O(log N ) buffers each taking up to O(nnz(X )R) space. SPLATT uses only O(PR)

Parallelization

Performing MTTKRP in a mode n amounts to performing a divide-add operation for each vector z i 1 ,...,i N to eventually form the matrix row M (n) (i n , :). Similarly, after the new U (n) (i n , :) is computed, one needs to update the vector z i 1 ,...,i N with a Hadamard multiplication and a division. Therefore, all nonzero elements of X whose nth index equals to i n contributes a summand to M (n) (i n , :), and the corresponding vectors in Z needs to be updated subsequently using the old and new values of U (n) (i n , :). To perform this, for each dimension n 2{1,...,N} and for each matrix row i n 2 I n we compute a reduction list of tensor nonzeros whose nth index is i n , which we denote as rl (n) (i n ). This way, each row M (n) (i n , :) can be computed in parallel by performing |rl (n) (i n )| vector operations. Similarly, once U (n) is computed, one can process each row i n in parallel to update the corresponding vectors z i 1 ,...,i N for each contributing nonzero x i 1 ,...,i N .

Algorithm 9 Parallel CP-ALS with nonzero factors.

Input: X : An N -mode tensor, X 2 R I1,...,I N R: The rank of CP decomposition U (1) ,...,U (N ) : Initial factor matrices with nonzero entries Output: [[ λ; U (1) ,...,U (N ) ]] : The rank-R CP decomposition of X 1: Initialize Z and W (n) for all n 2 N N .

2: repeat 3:

for n =1,...,N do for i n 2I

(n) p do 6:

M (n) (i n , :) 0 7:
for z i1,...,i N 2 rl (n) (i n ) do 8:

M (n) (i n , :) += z i1,...,i N /U (n) (i n , :) 9: z i1,...,i N = z i1,...,i N ↵ U (n) (i n , :)
10:

H (n) ⇤ i6 =n W (i) I Matrix Hadamard product 11: U (n) M (n) H (n) † I Row-parallel GEMM 12: λ NONZERO-COLUMN-NORMALIZE(U (n) )
13:

W (n) U (i) T U (i) I Row-parallel SYRK 14:
parallel for p =1...P do I Update Z 15:

for i n 2I (n) p do 16: for z i1,...,i N 2 rl (n) (i n ) do 17: z i1,...,i N = z i1,...,i N ⇤ U (n) (i n , :) 18:
until convergence or reaching maximum number of iterations 19: return [[ λ; U (1) ,...,U (N ) ]]

The parallel algorithm for computing the CP decomposition is shown in Algorithm 9. In the main subiteration loop, we first compute the MTTKRP result M (n) in parallel using P processes at Lines 4-9. This step assumes a partition I p of rows of the matrix M (n) independently thanks to the precomputed reduction lists. Immediately after the process uses the entry z i 1 ,...,i N in MTTKRP, it divides it by the old value of U (n) (i, :). Once M (n) is formed, at Line 10 we compute H (n) by performing the Hadamard product of R ⇥ R matrices W (k) for k 6 = n, whose cost is negligible as R is a small constant in practice. Next, at Line 11 we update the factor U (n) in a parallel dense matrix multiplication step, in which each process p performs the multiplication M (n) (I (n) p , :)H (n) † . Once U (n) is computed, we swap its small entries with ✏ or -✏ and normalize its columns in a combined step at Line 12 in which each process p works on the sub-matrix U (n) (I (n) p , :). Using the updated U (n) , we first compute the new W (n) in another parallel dense matrix multiplication step at line Line 13, where the process p similarly performs U (n) (I

(n) p , :) T U (n) (I (n)
p , :), then multiply the entries of Z with the corresponding rows of U (n) using the same parallelization scheme as in Line 4. The initialization of matrices Z as well as W (n) at Line 1 are done in parallel similar to the manner of updating these matrices in the iteration loop. At the end of each iteration, one has to check the convergence as well. This computation takes insignificant amount of time [START_REF]Parallel CP decomposition of sparse tensors using dimension trees[END_REF], hence we skip the details.

Reduction lists for a dimension n can be computed by making two passes over the tensor nonzero indices in nth dimension to form a very efficient compressed data structure consisting of rl (n) , which correspond to reduction list pointers, and rli (n) , which correspond to the elements in the reduction list, in O(nnz(X )) time. This yields O(N nnz(X )) cost for all dimensions, and we can process each dimension in parallel. In contrast, existing methods in the literature require sorting the tensor indices which takes O(N nnz(X ) log(nnz(X ))) time for SPLATT [START_REF]Tensor-matrix products with a compressed sparse tensor[END_REF], and O(N log N nnz(X ) log(nnz(X ))) for HYPERTENSOR [START_REF]Parallel CP decomposition of sparse tensors using dimension trees[END_REF]. We show this data structure for a small tensor X 2 R 5⇥5⇥5⇥5 in Figure 4.1, and skip the computational details.

Load balancing

For a P -way parallel execution of Algorithm 9, one needs to partition the row indices 1,...,I n into

P sets I (n) 1 ,...,I (n)
P for each dimension n. There are two types of computational costs imposed on each process p by a such partition. First, the process p performs O( P

i2I (n) p |rl (n) (i)|)
p | should be balanced. Though these rows can be partitioned arbitrarily, in practice we desire to assign a contiguous set of rows to each thread to preserve the data locality. This not only helps improve the memory footprint of the MTTKRP step, but also increases the efficiency of BLAS routines used for dense matrix computations.

We define partitioning problem in this case as follows. For each row i n , we have an associated pair (|rl (n) (i n , :)|, 1) of costs that corresponds to sparse tensor and dense matrix computations, respectively. We aim to partition this "chain" of rows into P contiguous parts so that both cost metrics are balanced across processes. The single-cost version of this problem corresponds to the chains-onchains partitioning (CCP) problem in the literature for which many fast optimal algorithms and effective heuristics exist [START_REF] Pinar | Fast optimal load balancing algorithms for 1D partitioning[END_REF]. We employ CCP algorithms by using only the first cost metric |rl (n) (i, :)| for partitioning, as we observe that in practice, balancing this metric also establishes good row-balance.

|

rl (1) rli (1) Figure 4.1: Performing MTTKRP for a 4-dimensional tensor X 2 R 5⇥5⇥5⇥5 in the first mode. Red (dark) and yellow (light) colors represent memory regions that are first touched by two different threads.

Optimizations for NUMA scalability

Performing MTTKRP for sparse tensors is an extremely memory bound operation as the tensor is very sparse in general, and the data accesses due to tensor nonzeros lack locality. Therefore, optimizing the memory footprint of the implementation plays a crucial role in obtaining high performance.

Particularly on a NUMA architecture, one has to carefully allocate memory pages in NUMA nodes to be able to utilize the available memory bandwidth at maximum, and distribute the memory pages equitably across NUMA nodes. In most systems, this can be ensured by properly using memory firsttouch policies after allocation, which in turn yields adequate memory page-to-socket bindings. In our implementation, after the allocation each thread performs a first-touch on the matrix rows as well as the rows of rl (1) and rli (1) that it owns. For the matrix Z, each thread initializes a block of nnz(X )/P vectors. This way, we not only maximize the NUMA bandwidth utilization, but also aim to reduce the inter-NUMA node memory requests as much as possible. This allocation scheme together with a balanced partitioning is shown in Figure 4.1 for two threads.

Experiments

We compared our algorithm with two state-of-the-art implementations, SPLATT and HYPERTENSOR, and ran them on the Grunch computing environment detailed in Section 2.7. We performed two types of experiments to measure the parallel scalability and the accuracy of our method. The first experiment compares the runtime of three methods using synthetically generated high dimensional sparse tensors. The second experiment uses a real-world tensor to compare the quality of approximation of the standard CP-ALS computation with our method. In all experiments, we use ✏ = 10 -6 as the threshold parameter.

Scalability

We compare the runtime of three methods using 4, 8, 16, and 32-dimensional randomly (uniform) generated tensors of size 10M at each dimension, and having 100M nonzero elements. We employ synthetic data instead of real-world tensors for two reasons. First, with random data we are able to control the dimensionality of the tensor while fixing other tensor parameters, e.g., dimension sizes, the number of nonzeros, and the distribution of nonzero indices; thereby, observe the performance of the algorithms with the increasing tensor dimensionality in a controlled manner. Second, there is a lack of available big high dimensional sparse tensors in the literature in parallel to the lack of efficient computational tools to handle such tensors. We run all three implementations using 1, 14 (single socket), and 28 cores/threads (two sockets) for 20 CP-ALS iterations using R = 16, and report the average time spent per iteration in Table 4.1 with labels splatt, hypertensor, and cp-eps corresponding to SPLATT,HYPERTENSOR, and our algorithm provided in Algorithm 9, respectively.

In Table 4.1, we observe that using ten-4D, hypertensor runs the fastest using 1 and 14 cores, yet cp-eps surpasses hypertensor using 28-cores owing to better NUMA optimizations described in Section 4.1.4. In all other instances, cp-eps stays the fastest among all three methods, and the performance gains increase steadily as the tensor dimensionality grows. Using 4-dimensional to 32dimensional tensors, we observe that the speedup of cp-eps over splatt consistently increases from 1.39x to 10.47x using 28 threads, and from 1.34x to 16.67x using a single thread, which conforms with the algorithm complexities provided in Section 4.1.1. Similarly, the speedup of cp-eps over hypertensor varies from 1.56x to 3.29x using 28 threads, and from 0.75x to 2.40x using a single thread. Figure 4.2 demonstrates the speedup results of cp-eps and hypertensor over splatt for ten-32D.

In Table 4.2, we compare the time spent on setting up data structures for MTTKRP using 28 threads on all datasets. We observe that cp-eps performs the preprocessing step up to 15x faster than splatt, and up to 24x faster than hypertensor. This significant improvement is possible owing to the smaller asymptotic complexity described in Section 4.1.2, as well as the parallelization using cp-eps.

Accuracy

In computing Algorithm 9, we impose the constraint on factor matrices that they do not contain very small elements, which perturbs the decomposition slightly and can potentially affect the quality of approximation. To assess this, we compare the accuracy of this method with that of the original CP-ALS algorithm. We employ the 3-dimensional NELL2 tensor described in Section 2.6. We run hypertensor and cp-eps 100 times (each with a random initialization of factor matrices) with the rank of approximation R 2{25, 50, 75, 100}, and compute the geometric mean of the approximation quality in each case. In Figure 4.3 we detail all these results. We observe that both methods produce equally good approximations to the original tensor up to a small margin of error for R = 75 mostly due to randomization in factor matrix initialization. This shows that the nonzero constraint on factor matrices indeed has a negligible effect on the accuracy.

Conclusion

In this work, we propose an efficient parallel algorithm for computing a CP decomposition in which the factor matrices are assumed to not contain any zero elements. This constraint enables a computational scheme that provides O(log k) and O(log N log k) faster preprocessing times for a tensor with k nonzero entries, and performs O(N ) and O(log N ) less MTTKRP work over two efficient state-ofthe-art implementations splatt and hypertensor with a negligible effect on the accuracy. We achieve up to 16.7x speedup in sequential and 10.5x speedup in parallel executions over these methods, with up to 24x less data preprocessing time, using up to O(log N ) less memory for storing intermediate computations. With these advancements, our approach renders the analysis of higher order big sparse Chapter 5

Parallel Sparse Tucker Decomposition

This chapter concerns an efficient computation as well as shared and distributed memory parallelizations of HOOI algorithm for computing the Tucker decomposition of sparse tensors. Recall that the two main operations in an iteration of HOOI are the TTMc step, which involves the multiplication of the tensor with N -1 matrices, and the TRSVD step, which computes the truncated SVD of the matricization of the tensor obtained from the TTMc step. For an efficient algorithm using sparse tensors, both these steps should be performed in a way that effectively exploits the sparsity of the input tensor.

Towards this end, we present the following contributions in this chapter. We design efficient shared and distributed memory parallel algorithms for the TTMc operation on sparse tensors. To this end, we first introduce a particular nonzero-based reformulation of TTMc. Using this formulation, we then introduce a preprocessing step called symbolic TTMc to identify data dependencies and perform all index computations before the HOOI iterations for efficiency. Then, we provide a sharedmemory parallel algorithm for the main iteration of HOOI which makes use of the symbolic TTMc step. Next, we introduce a coarse and a fine-grain task definition for TTMc and TRSVD steps within the HOOI algorithm, and propose a hybrid shared-distributed memory parallel algorithm based on the distribution of these tasks. We discuss the computational and the communication requirements of the algorithm for a given task distribution, and make use of the hypergraph models from our earlier work on CP-ALS [START_REF] Kaya A N D B | Scalable sparse tensor decompositions in distributed memory systems[END_REF] for reducing communication and achieving load balance during each HOOI iteration. Third, we stress how to efficiently perform the TRSVD step in a distributed memory setting, and make use of the PETSc [START_REF] Balay | PETSc web page[END_REF] and SLEPc [START_REF] Roman | SLEPc users manual[END_REF] libraries in this step. We carefully designed this step so that the communication requirements in parallel iterative algorithms used for computing the singular vectors are reduced, and the load balance is achieved by making use of the data decomposition of the TTMc step. Finally, we propose an efficient OpenMP-MPI hybrid parallel implementation of the HOOI algorithm in C++, and present scalability results on a high-end parallel system using up to 4096 cores on real world tensors. To the best of our knowledge, this is the first high performance parallel implementation of the HOOI algorithm for sparse tensors in shared/distributed memory environments using OpenMP/MPI. The organization of this chapter is as follows. We provide a reformulation of the TTMc operation for sparse tensors in the next section. Then, in Section 5.2, we propose shared and distributed memory parallel HOOI algorithms based on this formulation, and discuss in detail the TTMc and the TRSVD steps. Next, we give a brief summary of recent related work in Section 5.3. Finally, we provide experimental results in Section 5.4, and conclude the chapter in Section 5.5. For the sake of simplicity

PERFORMING MULTIPLE TTMS FOR SPARSE TENSORS

57

of the notation and the discussion, we occasionally discuss the case of a 3-dimensional tensor, even though our algorithms and implementations have no such restriction. We adequately generalize the discussion to N -dimensional tensors whenever necessary.

Performing multiple TTMs for sparse tensors

We recall that the TTM of a tensor

X 2 R I 1 ⇥•••⇥I N with a matrix U 2 R J⇥In is denoted by Y X ⇥ n U. and the result Y is a tensor of size I 1 ⇥•••⇥I n-1 ⇥ J ⇥ I n+1 ⇥•••⇥I N . A particular element of Y is given by y i 1 ,...,i n-1 ,j,i n+1 ,...,i N = In X in=1 x i 1 ,i 2 ,...,i N u jin .
(5.1)

A tensor can likewise be multiplied by a set of matrices along a given set S of modes. We use the notation TTMc(X ,S,{U (n) : n 2 S}) to refer to the tensor n-dimensional product of X with matrices U (n) T for n 2 S. We use TTMc(S) for clarity, as the tensor X and the matrices U (n) T should be clear from the context. The operation Y X ⇥ i6 =n U (i) T stands for the TTMc of a tensor in all modes except n, which we also denote as TTMc({1 ...N}\{n}).

The way to perform the TTMc operation at Line 4 of Algorithm 3 is especially important. These TTM's can be performed in any order [START_REF]Tensor decompositions and applications[END_REF] using various schemes [START_REF] Kolda | Scalable tensor decompositions for multi-aspect data mining[END_REF] that formulate TTMc in terms of multiple tensor-times-vector (TTV) operations. As we are concerned with sparse tensors, we formulate TTMc in a way that specifies the computational task to be carried out for each nonzero x i 1 ,...,i N 2 X , which in turn enables expressing parallelism with different task granularities. Let Z = X ⇥ 2 U T 2 for X 2 R I 1 ⇥I 2 ⇥I 3 and U (2) 2 R I 2 ⇥R 2 . By considering (5.1) for Z, and performing the summation over nonzeros we obtain:

z i,t,k = X x i,j,k 2X x i,j,k U (2) (j, t) . (5.2)
We can vectorize this as

Z(i, :,k)= X x i,j,k 2X
x i,j,k U (2) (j, :) .

Similarly, Y = Z ⇥ 3 U T 3 for U (3) 2 R I 3 ⇥R 3 can be written as:

y i,m,t = X z i,m,k 2Z z i,m,k U (3) (k, t) .
Rewriting z i,m,k as in (5.2) we get

y i,m,t = X x i,j,k 2Z x i,j,k U (2) (j, m)U (3) (k, t),
and finally by applying the vectorization in (5.3) twice we obtain Y(i, :, :) = X

x i,j,k 2X

x i,j,k U (2) (j, :) • U (3) (k, :) .

Here, for each nonzero x i,j,k 2 X , we perform the outer product U (2) (j, :) • U (3) (k, :), scale it with x i,j,k , then add the result to the R 2 ⇥ R 3 dense matrix Y(i, :, :). Using the matricization of Y in the first mode, this results in the following formula where a Kronecker product replaces the outer product:

Y (1) (i, :) = X x i,j,k 2X
x i,j,k U (2) (j, :) ⌦ U (3) (k, :) .

(5.4)

For N -dimensional case, the formulation (5.4) trivially generalizes to

Y (n) (i n , :) = X x i 1 ,...,in,...,i N 2X
x i 1 ,i 2 ,...,i N ⌦ t6 =n U (t) (i t , :).

(5.5)

This formulation specifies the operations performed for each nonzero element of a tensor in TTMc.

The resulting computation of TTMc is called nonzero-based, and is given in Algorithm 10.

Algorithm 10 TTMc operation Y = X ⇥ i6 =n U (i) T

Input: X : An N -dimensional tensor, X 2 R I1,...,I N U (1) ,...,U (N ) : Matrices for TTMc,

U (n) 2 R In⇥Rn n: The mode for performing TTMc Output: Y X ⇥ i6 =n U (i) T Y (n) 0 for all x i1,...,i N 2 X do Y (n) (i n , :) Y (n) (i n , :) + x i1,...,i N ⌦ t6 =n U (t) (i t , :)

Parallel HOOI for sparse tensors

An efficient shared memory parallelization of the TTMc operation given in Algorithm 10 should avoid expensive lock mechanisms to resolve data dependencies. We perform a preprocessing step to organize the computations in a way that the subsequent numeric computations can be performed in parallel without any write conflicts. Following the TTMc step, computing the TRSVD of the matricized tensor Y (n) requires special attention. Direct SVD methods that are employed to compute TRSVD in dense Tucker decomposition algorithms are not feasible for sparse Tucker decomposition due to computational and memory constraints. For this reason, we resort to iterative methods for TRSVD, which not only reduces the computational cost by exploiting the low rank of approximation, but also renders TRSVD computation feasible in terms of memory utilization.

For the distributed memory parallelism, we employ coarse-and fine-grain task definitions. A coarse-grain task corresponds to computing a particular row Y (n) (i, :) of the TTMc result, as well as the corresponding row U (n) (i, :) of the factor matrix using TRSVD. In this scenario, the owner of this task possesses all the tensor nonzeros x i 1 ,...,i N where i n = i. Also, each such nonzero implies data dependencies to the tasks corresponding to rows U (1) (i 1 , :),...,U (N ) (i N , :) to be able to perform the computation of Y (n) (i, :) using Algorithm 10. Fine-grain task definition relaxes this constraint by allowing nonzeros to be distributed freely. It associates each nonzero x i 1 ,...,i N with a task which is responsible to compute x i 1 ,...,i N ⌦ t6 =n U (t) (i t , :) and generate a partial result for Y (n) (i n , :) of size Q t6 =n R t . This size is exponential in the ranks of approximation; therefore, merging the partial results can get very expensive in terms of communication, hence should be avoided. For this reason, we propose a method to effectively handle this communication within the TRSVD step.

Shared memory parallelism Parallel TTMc

As shown in Algorithm 10, each nonzero x i 1 ,...,i N contributes an outer product to Y (n) (i n , :) while performing TTMc in the first mode. For shared memory parallelism, this poses a write conflict whenever two threads simultaneously process nonzeros whose first index are i. To resolve this, we make a pass over the data to compute an update list ul n (i n ) that holds the list of nonzeros x i 1 ,...,i N that contribute to Y (n) (i n , :). In the actual implementation, we only store the index t of the nonzero x(t)=x i 1 ,...,i N to avoid duplicating the nonzero within ul n (i n ). This way, we untangle the write conflicts for each row of Y (n) and avoid using lock mechanisms. We also store the set J n of all indices i 2 I n such that ul n (i) 6 = ;. We repeat this computation in all dimensions, and name this step as symbolic TTMc, as it resolves all the index computations and dependencies once and for all outside the main loop of HOOI (shown at Lines 1-2 of Algorithm 11). This symbolic data can be reused many times for faster numeric TTMcs within the main loop of HOOI. Finally, symbolic TTMc of each dimension can be performed independently; hence, we perform this computation in parallel for each dimension.

After the symbolic TTMc, each row i of Y (n) can be updated independently in parallel by using ul n (i), which composes the parallel numeric TTMc step at Lines 5-8 of Algorithm 11. In our implementation, we use OpenMP parallel loop with dynamic scheduling to distribute the tasks to threads.

Algorithm 11 Shared memory parallel HOOI

Input: X : An N -dimensional tensor U (1) ,...,U (N ) : Initial factor matrices R 1 ,...,R N : Ranks of approximation Output: [[ G; U (1) ,...,U (N ) ]] : Tucker approximation of X 1: parfor n =1to N do 2:

{ul n , J n } SymbolicTTMc(X , {1,...,N}\{n}) 3: repeat 4:

for n =1to N do

Y (n) (i, :) 0 7:
for all x i1,...,i N 2 ul n (i) do 8:

Y (n) (i, :) += x i1,...,i N [⌦ t6 =n U (t) (i t , :)] 9: U (n) TRSVD(Y (n) ,R n ) 10:
G Y ⇥ N U (N ) 11: until convergence or maximum number of iterations 12: return [[ G; U (1) ,...,U (N ) ]]

Dimension tree-based parallel TTMc

Note that each TTMc call involves the multiplication of X with the set {U (1) ,...,U (N ) }\U (n) of matrices to eventually update U (n) at the end of the subiteration for mode n. Similar to the techniques used in Chapter 3, this opens the possibility to store and reuse partial TTM results across TTMc calls in different dimensions as follows. For a 4-dimensional tensor X 2 R I 1 ⇥I 2 ⇥I 3 ⇥I 4 , in the first two ALS subiterations, one needs to compute X ⇥ 2 U (2) ⇥ 3 U (3) ⇥ 4 U (4) and X ⇥ 1 U (1) ⇥ 3 U (3) ⇥ 4 U (4) and subsequently update the corresponding matrices U (1) and U (2) , respectively. Note that in doing so, U (3) and U (4) remains unchanged, which brings up the possibility of computing the intermediate tensor Z = X ⇥ 3 U (3) ⇥ 4 U (4) , then reusing this partial result in the first and second subiterations as Z ⇥ 2 U (2) and Z ⇥ 1 U (1) , respectively. This computation stays valid as one can perform TTM in a set of distinct modes in any order. This way, however, the cost of TTMs reduces significantly.

As we did in Chapter 3, we use BBDTs to minimize the number of TTMs and the number of intermediate tensors held in the dimension tree. Note also that setting the ranks of approximation to 1 equates a TTM to a TTV, in which case the nonzero structure of tree tensors are obtained by the symbolic TTV step described in Section 3.1.2. Increasing the ranks of approximation only increments the size of the value array per tensor element proportional to these ranks according to (5.5), while the sparsity structure in the modes in which TTM is not performed stays still. Therefore, we use the same symbolic TTV step described in Section 3.1.2 to construct the nonzero set of the tree tensors. Once the tree is constructed, we employ Algorithm 12, which is analogous to Algorithm 6, to carry out TTM operations for a node of a dimension tree constructed using a sparse tensor. The only major difference here is the outer product replacing the Hadamard product at Line 9, which is due to (5.5). We can then invoke SMP-DTREE-TTM(L n ) to compute TTMc for a dimension n. V t (i 1 ,...,i N ) zeros(1,R) I Initialize with a zero vector of size 1 ⇥ R.

6:

for all (j 1 ,...,j N ) 2R(i 1 ,...,i N ) do I Perform updates using elements in R.

7:

r V P (t) (j 1 ,...,j N )

8:

for all d 2 δ(t) do I Multiply the vector r with corresponding matrix rows. V t (i 1 ,...,i N ) V t (i 1 ,...,i N )+r I Do the update due to element (j 1 ,...,j N ).

Parallel truncated SVD

Following the TTMc in a mode n, HOOI requires finding the leading R n singular vectors of the matricized tensor Y (n) to update the matrix U (n) . Here,

Y (n) is of size I n ⇥ Q t6 =n R t .
In Algorithm 11, we directly compute this matricized tensor, and thereby avoid the cost of matricization. In a recent work on parallel Tucker decomposition of dense tensors [5], leading singular vectors of Y (n) are extracted by computing the eigenvalues of the Gram matrix

Y (n) Y (n)
T of size R In⇥In , where I n is typically in the order of thousands. However, I n can easily be in the order of millions for sparse tensors, rendering this method impractical for our purpose. Also, direct methods for eigenvalue and singular value problems typically compute all eigen/singular values at once, whereas in HOOI we need only R n left leading singular value/vector pairs out of min(I n , Q i6 =n R i ). For these reasons, we resort to using matrix-free iterative methods to compute the TRSVD of Y (n) [START_REF] Roman | SLEPc users manual[END_REF]. This way, we avoid forming the Gram matrix, and compute only the required singular value/vector pairs. In a shared memory context, this TRSVD can be parallelized by using optimized BLAS2 gemv kernel for the matrix-vector (MxV) and matrix transpose-vector (MTxV) multiplications, which dominate the computational cost of the TRSVD due to the matrix being dense. As we demonstrate in the next section, this approach also enables us to reduce the communication requirements in the distributed memory setting. One can also employ randomized linear algebra [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF] and cross-approximation techniques [START_REF] Goreinov | A theory of pseudoskeleton approximations[END_REF] to further accelerate this step, yet the proposed parallelization and the communication techniques stay pertinent.

After all factor matrices are updated, the core tensor is formed at Line 10 to check the convergence. Here, since the preceding subiteration of HOOI computes Y (N ) , which already holds the tensor Y = X ⇥ i6 =N U (i) T in the matricized form, we can multiply Y with U (N ) in mode N to obtain G. Both Y and G are dense tensors while

G 2 R R 1 ⇥•••⇥R N being significantly smaller than Y 2 R R 1 ⇥•••⇥R N -1 ⇥I N .
The parallel computation of dense G can be performed efficiently using BLAS3, and in practice its cost should be negligible compared to the cost of sparse irregular operations and TRSVDs carried out in the preceding N subiterations of HOOI. We skip the details of the parallelization of Line 10, and refer the reader to a recent work by Li et al. [START_REF] Li | An input-adaptive and in-place approach to dense tensor-times-matrix multiply[END_REF].

Distributed memory parallelism

Coarse-grain parallel HOOI

Recall that there are two main operations for each mode n in an iteration of HOOI: a TTMc step to obtain a matricized tensor Y (n) , and a TRSVD step to obtain the matrix U (n) from Y (n) . In the coarse-grain task decomposition, we define computing each row U (n) (i, :) as an atomic task in the TTMc step, and hold the owner of this task responsible for computing Y (n) (i, :) as well in the TRSVD step. We denote this task by t (n) i , and make the owner of t (n) i own the corresponding data elements U (n) (i, :) and Y (n) (i, :). For each mode n, we partition these tasks. As a result, the process p k owns the index set I (n) k of tasks, i.e., t

(n) i is owned by p k for each i 2I (n) k .
For a 3-dimensional X , Y (1) (i, :) receives a contribution x i,j,k U (2) (j, :) ⌦ U (3) (k, :) for each nonzero x i,j,k 2 X in Algorithm 10. Therefore, t

i needs all nonzeros in the tensor slice X (i, :, :) as well as the corresponding rows U (2) (j, :) and U (3) (k, :) to perform the Kronecker product. To compute Y (1) (i, :), t

i involves |X (i, :, :)| Kronecker products due to TTMc. During this computation, for each nonzero x i,j,k , t

i needs the data owned by t for each mode n, as well as the set of nonzeros X k that are needed to perform the local computations associated with these tasks. At Lines 1-6, p k performs the symbolic TTMc with its local tensor X k . Next, at Lines 9-12 the local TTMc operation Y (n)

X k ⇥ i6 =n U (i) is performed. Note that at Line 4, we set K n I (n)
k to make sure that the coarse-grain algorithm only computes TTMc results for the owned set of rows I (n) k . Also, p k does not need to store the whole matrix U (n) ; instead, it stores the set of owned rows U (n) (I

k : Set of task indices owned by p k in each mode U (1) ,...,U (N ) : Initial factor matrices R 1 ,...,R N : Ranks of approximation Output: [[ G; U (1) ,...,U (N ) ]] : Tucker approximation of X 1: parfor n =1to N do 2:

{ul n , J n } SymbolicTTMc(X k , {1,...,N}\{n})

3:
if type = 'coarse-grain' then 4:

K n I (n) k 5:
else 6:

K n J n 7: repeat 8:
for n =1to N do

Y (n) (i, :) 0 11:
for all x i1,...,i N 2 ul n (i) do 12:

Y (n) (i, :) += x i1,...,i N [⌦ t6 =n U (t) (i t , :)] 13: U (n) TRSVD(Y (n) ,R n ) 14:
Send/receive the updated rows of U (n) 15:

G k Y ⇥ N U (N )
16:

G ALLREDUCE(G k ) 17: until convergence or maximum number of iterations 18: return [[ G; U (1) ,...,U (N ) ]] we partition the vector y according to I (n) k . This way, after gathering all entries of x at all processes, the MxV operation y Y (n) x can be performed locally without any communication on y entries. Also, after the TRSVD solver converges, the computed left singular vectors have the same partition as y. As a result, p k ends up having all rows U (n) (I (n) k , :) in place, avoiding any post-communication. For x T y T Y (n) , we compute the local result y T Y (n) , then perform an all-to-all reduction to sum up the final results in the owner processes.

In our recent work [START_REF] Kaya A N D B | Scalable sparse tensor decompositions in distributed memory systems[END_REF], we used a similar coarse-grain task definition in the context of the parallel CP-ALS algorithm, and proposed a hypergraph model for representing the computational and communication requirements of the parallel algorithm. Here, we adopt the same hypergraph model to reduce the total communication volume and to balance the computational load during the HOOI iterations. In this model, we represent tasks with vertices and their interdependence using hyperedges. For each task t k , respectively. This becomes an exact measure for the computational cost of TTMc, as for each nonzero we perform an outer product with N -1 vectors. One can potentially associate a secondary computational cost of 1 with each task to balance the number of rows owned per process for the MxV and MTxV multiplication steps of the TRSVD step; yet in practice we observed that not adding this constraint still provides balanced number of rows. For each nonzero x i,j,k that imposes an update on tasks t k , we connect the vertices of each such task with the corresponding 3 hyperedges to model their data dependency. The standard partitioning problem of this hypergraph corresponds to reducing the total communication volume while establishing the load balance in the TTMc step. This could also ensure load balance in the TRSVD step if the processes have almost equal number of tasks (we investigate this in our experiments). Finally, as mentioned above, only the MTxV operation in the TRSVD step requires communication, which is regular and has a cost independent from the task distribution.

Fine-grain parallel HOOI

The coarse-grain approach has two main limitations. First, the number of tasks for a dimension n is limited by I n . In case the tensor is very small in one of the dimensions, this poses a granularity problem by not having enough tasks for parallelism, which hinders scalability. Also, coarse-grain tasks tend to be heavily interdependent with their data. As a result, there is typically little room in finding a good partition for parallelization. To address both limitations, we propose a fine-grain variant of the parallel HOOI which enables the partial computation of the rows of Y n . The fine-grain approach in Algorithm 13 differs from the coarse-grain one at Line 6, which enables computing partial results for the rows that are not owned by p k .

Similar to the coarse-grain algorithm, we define the task t (n) i to denote the ownership of U (n) (i, :). In the fine-grain case, the owner of t (n) i does not necessarily perform all the computations associated with Y (n) (i, :) nor U (n) (i, :). For each nonzero x i 1 ,...,i N 2 X , we define an associated task z i 1 ,...,i N , and let the process p k having the set of nonzeros X k also own the corresponding z-type tasks. For each mode n, the owner of z i 1 ,...,i N is responsible for performing the operation x i 1 ,...,i N [⌦ t6 =n U (t) (i t , :)] and generating a partial result for Y (n) (i n , :). One may consider merging these partial results in a way that the owner of the task t (n) i gets the final result Y (n) (i n , :). This way, one can proceed with the TRSVD computation of the row-wise distributed matrix Y (n) just as in the coarse-grain case. However, the problem in this scenario is that each partial result Y (n) (i, :) to be communicated is of size Q i6 =n R i , which is exponential in the ranks of approximation, and can easily get very large. In contrast, each message for communicating the rows of U (n) at Line 14 of Algorithm 13 is of size R n .

We make use of the following observation to asymptotically reduce the communication cost due to partial results of TTMc. In a P -way parallel execution, performing the local TTMc with a fine-grain task distribution produces the matrix

Y (n) in the sum-distributed form Y (n) = Y 1 (n) + ••• + Y P (n)
where Y k (n) is the partial TTMc result generated by the process p k . One should avoid assembling Y (n) , otherwise a high communication overhead incurs. Fortunately, we only need to provide MxV and MTxV operations associated with the matrix Y (n) in the subsequent TRSVD step. We can perform these multiplications without assembling Y (n) as follows. For the MxV operation, we perform

y Y k (n)
x at each process p k and generate a partial result on y. Then, we execute a point to point communication on the entries of y so that the owner of the task t (n) i sums all the partial results to obtain the final value of y i . In this way, instead of communicating a partial result Y k (n) (i, :) of size Q i6 =n R i , we communicate a single vector entry y i for each i requiring communication in each iteration of the TRSVD solver. The number of TRSVD iterations is typically in the order of R n or less, so it makes this communication cost conformal with the cost at Line 14 of Algorithm 13. We can easily perform the MTxV operation by computing x T y T Y k (n) at each process p k , and then by performing an all-to-all communication on x T just as in the coarse-grain algorithm.

We model the tasks of the fine-grain algorithm and their dependencies using the same hypergraph from our previous study [START_REF] Kaya A N D B | Scalable sparse tensor decompositions in distributed memory systems[END_REF], where we model tasks with vertices and the dependencies among tasks with hyperedges. With this model, the communication cost at Lines 13-14 of Algorithm 13 is equal to the cutsize of a partition of the corresponding hypergraph and can be effectively reduced by existing hypergraph partitioning tools. We refer the reader to [START_REF] Kaya A N D B | Scalable sparse tensor decompositions in distributed memory systems[END_REF] for the detailed analysis. We note also that by not combining the partial results of Y (n) , we increase the total computational load of matrixvector multiplications in the TRSVD solver, as we end up having more rows than |I n | to multiply in total. Fortunately, this increase is also equal to the cutsize of the hypergraph, and is significantly reduced with a good partition. Therefore, minimizing the cutsize is beneficial for reducing both the communication cost of the parallel HOOI and the redundant computation in its TRSVD step. Finally, as in the coarse-grain algorithm, the load balance in the local MxV and MTxV operations can be achieved if the MPI-ranks have equal number of rows in Y (n) . This is one of the complex partitioning problems where the total computational load can only be determined after a partition [START_REF] Kaya | On partitioning problems with complex objectives[END_REF][START_REF] Pinar | Partitioning for complex objectives[END_REF]. We do not explicitly address this problem and hope that assigning equal amount of n-dimensional indices will lead to load balance.

Related work

We give a brief overview of the recent progress on efficient tensor decomposition (CP and Tucker) algorithms. These can be categorized into four classes: (i) toolboxes for Matlab and similar environments [3,6,[START_REF] Bader | Matlab tensor toolbox version 2[END_REF][START_REF] Kolda | Scalable tensor decompositions for multi-aspect data mining[END_REF][START_REF] Papalexakis | ParCube: Sparse parallelizable CANDECOMP-PARAFAC tensor decomposition[END_REF]; (ii) implementations for shared memory systems [START_REF] Baskaran | Efficient and scalable computations with sparse tensors[END_REF][START_REF] Baskaran | Low-overhead load-balanced scheduling for sparse tensor computations[END_REF][START_REF] Li | An input-adaptive and in-place approach to dense tensor-times-matrix multiply[END_REF][START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF]; (iii) implementations based on MapReduce paradigm [START_REF] Jeon | Haten2: Billion-scale tensor decompositions[END_REF][START_REF] Kang | GigaTensor: Scaling tensor analysis up by 100 times -Algorithms and discoveries[END_REF]; (iv) implementations for distributed memory systems [5,[START_REF] Choi | DFacTo: Distributed factorization of tensors[END_REF][START_REF] Kaya A N D B | Scalable sparse tensor decompositions in distributed memory systems[END_REF][START_REF] Smith | DMS: Distributed sparse tensor factorization with alternating least squares[END_REF]. The implementations in the first group are very useful tools that enable fast prototyping. Those in the second group and similar work are helpful when data fits into the memory of a single machine, which is nowadays large enough to accommodate tensors from many applications. Those in the third and fourth groups enable computations on tensors that do not fit into the memory of a single machine. The ones in the third group are not designed for high performance, as MapReduce paradigm is meant to perform multiple passes over out-of-core data and perform global communication shuffling the input data.

To the best of our knowledge, there is no high performance distributed memory implementation of algorithms for the sparse Tucker decomposition. Among the cited references above, HaTen2 [START_REF] Jeon | Haten2: Billion-scale tensor decompositions[END_REF] is a MapReduce based HOOI implementation. Li et al. [START_REF] Li | An input-adaptive and in-place approach to dense tensor-times-matrix multiply[END_REF] investigate efficient shared memory execution of tensor times matrix products and as a future work mention how this can be used to perform intra-node TTM computations in a distributed memory setting. This work does not discuss other components of an HOOI implementation. Austin et al. [5] propose a distributed memory parallel implementation of HOOI for dense tensors. The challenges that are faced are very different from those faced in the sparse case, essentially due to all communications involving all processes, and memory accesses being regular in the dense case.

Experiments

We conducted our experiments on two parallel computing platforms. The first platform is the Blue Gene/Q described in Section 2.7. On this cluster, we ran our experiments up to 256 nodes (4096 cores) where we achieved the maximum scalability. Each core of PowerPC A2 can handle one arithmetic and memory operation simultaneously; therefore we assigned 32 threads per node (2 threads per core) to benefit from this. All codes we used in our benchmarks were compiled using the Clang C++ compiler (version 3.6.0) with IBM MPI wrapper using -O3 option for compiler optimizations, and linked against IBM ESSL library for LAPACK and BLAS routines. Our code depends on PETSc and SLEPc (version 3.6.2) libraries for the distributed truncated SVD computations. The second one is the the Ada cluster of IDRIS, which consists of 304 nodes each having 128 GBs of memory and four 8core Intel Sandy Bridge E5-4650 processors running at 2.7 GHz. We ran our experiments using up to 2048 cores (64 nodes), which is the maximum allowed in the cluster. We tried different configurations of 32 cores within a node for the MPI/OpenMP rank/thread assignments, and found the assignment of 4 MPI ranks (and 8 OpenMP threads per MPI-rank) per node to give the best results. 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, -openmp flag to enable OpenMP, and -mkl option to use the Intel MKL library (version 11.1.1) for LAPACK, and BLAS routines. To obtain the sequential runtimes, we allocated 1 MPI rank with 1 threads, disabled the multithreading within MKL, and used a highmemory nodes in the cluster having 256 GB memory. Finally, we conducted our experiments using four tensors, namely Netflix, NELL, Flickr, and Delicious, described in Section 2.6.

Our parallel algorithms are independent from the partitioning method. Therefore, we use two partitioning methods to test their performance. The first partitioning method assigns the tasks uniformly at random to processes (for coarse-grain tasks we use a blocked variant), and the second one uses hypergraph partitioning tool PaToH [START_REF] Çatalyürek A N D C | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF]. The first method is fast, promises load balance, but it does not pay attention to the communication overhead. The second method achieves load balance, reduces communication, but is time consuming. Speedups using these two partitioning methods show the worst case behavior and the potential of the parallel algorithms, if one is willing to pay the preprocessing cost. In general, a good parallel algorithm should deliver good performance with the second partitioning method; but it should also enjoy acceptable speed up with the first partitioning method.

We used PaToH (version 3.2) with default options to partition the hypergraphs. We created all partitions offline, and ran our experiments on these partitioned tensors on the cluster. We do not report timings for partitioning hypergraphs with PaToH, which is costly. Yet in most applications, the tensors from the real-world data are built incrementally and analyzed repetitively. In this scenario, 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 [START_REF] Kiers | A fast method for choosing the numbers of components in Tucker3 analysis[END_REF]. In these cases, the time spent in partitioning can be amortized across multiple runs.

To the best of our knowledge, HaTen2 [START_REF] Jeon | Haten2: Billion-scale tensor decompositions[END_REF] is the only parallel implementation of HOOI (see Section 5.3). All reported parallel runtimes of HaTen2 [START_REF] Jeon | Haten2: Billion-scale tensor decompositions[END_REF] are larger than that of the sequential execution of MET algorithm [START_REF] Bader | Matlab tensor toolbox version 2[END_REF][START_REF] Kolda | Scalable tensor decompositions for multi-aspect data mining[END_REF], and the sequential runtime of our method is less than that of MET. For example, on a random tensor of size 10K ⇥ 10K ⇥ 10K with 1M nonzeros, Tucker decomposition with five HOOI iterations took 87.2 seconds in MET and 11.3 seconds in our method (on a single core), including all preprocessing. This difference is expected as neither HaTen2 (which uses MapReduce) nor MET (which is a Matlab tool) are made for high performance; thus, we do not report further comparisons.

We set the ranks of approximation R 1 = R 2 = R 3 = 10 for the 3-dimensional tensors, and R 1 = R 2 = R 3 = R 4 =5for the 4-dimensional tensors, which is a viable choice in data analysis applications [START_REF] Ymeonidis | Tag recommendations based on tensor dimensionality reduction[END_REF][START_REF] Zhang | Cubic analysis of social bookmarking for personalized recommendation[END_REF]. We then run the parallel HOOI for 5 iterations, and report the average time spent per HOOI iteration. The symbolic TTMc is expected take much less time than the HOOI iterations. For instance, in 256-way parallel execution of Algorithm 13 using the fine-grain hypergraph partitioning for 5 iterations, symbolic TTMc took 14%, 12%, 19%, and 5% of the total execution time for Delicious, Flickr, Netflix, and NELL tensors, respectively. In general, HOOI is expected to run for more iterations; so this cost is expected to become less important. In addition, finding a

Distributed memory results

In Table 5.1, we give the strong scalability results for our distributed memory parallel HOOI algorithm. We report the average time spent per iteration in Algorithm 13. For each tensor in our dataset, we report two results for fine-grain and coarse-grain parallel algorithms (with two different partitioning methods). Those with the suffix "-hp" uses PaToH's partitionings; "fine-rd" refers to a random partitioning, whereas and "coarse-bl" corresponds to a contiguous block partitioning of tasks. We evaluate the algorithms up to 256 compute nodes where each node executes Algorithm 13 using 32 threads. Since the amount of memory per node of Blue Gene/Q is only 16GBs, some runs were not feasible. These missing runs are shown with "-" in Table 5.1. In two tensors, using the fine-rd partition was feasible only after 16 nodes due to higher memory requirements for storing partial results (which is proportional to the communication cost).

We first observe in Table 5.1 that in all test cases our algorithm graciously scale up to 256 MPI ranks (or 4096 cores), except with the fine-hp partition on Netflix tensor, where we observe a slowdown at 256 nodes. On Delicious tensor, our parallel algorithm achieves 13.5x speedup using 256 nodes over the run on 8 nodes with the fine-hp partition. With the fine-rd partition, the algorithm also scales to 256 nodes, although it runs almost two times slower than fine-hp. This is expected, as fine-rd targets load balance but incurs more communication overhead. Similarly on Flickr, our parallel algorithm achieves 10.3x speedup using 256 nodes over its execution on 8 nodes with the fine-hp partition, and runs roughly twice as fast as the configuration with the fine-rd partition. Unfortunately, with Delicious, Flickr, and Netflix datasets, we were not able to get the sequential and shared memory parallel timings on a single node due to data not fitting into the memory; hence we cannot provide overall speedup results over the sequential execution. For those results, we refer the reader to a technical report [START_REF] Kaya A N D B | High-performance parallel algorithms for the Tucker decomposition of higher order sparse tensors[END_REF], where we present speedups up to 742x on a different architecture with more memory. On NELL tensor, we managed to get runs on a single node. Using the fine-rd partition using 256 nodes (4096 cores), we obtained 280x speedup over the sequential execution. This translates into 29x speedup over the execution on a single node. Yet, unlike other three tensors, on this data the fine-hp partition lead to slower execution than the fine-rd partition. We analyzed the underlying reason for this result on 256 MPI ranks and saw that the maximum communication volume per process for with the fine-hp partition in the dominant dimension was 543K in contrast to 366K in the fine-rd partition. Here, this entailed a large overhead that could not be compensated by the reduction in the total communication volume (20M vs 94M).

In Table 5.2, we give the computation and communication requirements of one HOOI iteration on the Flickr tensor with all partitionings for 256 MPI ranks (4096 cores). In this table, W TTMc and W TRSVD correspond to the computational load of the TTMc and TRSVD steps of Algorithm 13 and Comm. vol. corresponds to the volume of send/receive communication incurred by different partitioning methods. We give both the average and the maximum values across all processes for all three metrics. We observe in the W TTMc columns of Table 5.2 that TTMc work per MPI-rank is always well balanced with the fine-grain partitions. This is owing to the finer granularity of tasks which allows perfect balance. On the other hand, with the coarse-grain formulation, the TTMc tasks are not well balanced, as some tasks might be significantly more costly than others. Particularly in the TTMc computation of the 4th dimension, we observe some computational imbalance of 436% and 471% using the coarse-hp and coarse-block partitions.

We realize in the W TRSVD columns of Table 5.2 that the average TRSVD work W TRSVD given by the fine-hp partition is only slightly higher than that given by the coarse-grain partitions, which introduce no overhead to the TRSVD computation. Particularly, W TRSVD is dominant in the third dimension, and the fine-hp partition results in the same total/average work (110K) as in the coarsegrain partitions. Using the random partition (fine-rd), however, the average W TRSVD is drastically increased to 435K in the same dimension. Moreover, even though none of the partitioning methods explicitly try to establish load balance for W TRSVD , we observe that the obtained load balance is generally acceptable. In the computationally dominant third mode, the fine-hp partition leads to 100% load imbalance, whereas the fine-rd, coarse-hp, and coarse-block partitions lead to 25%, 79%, and 18%.

The last two columns of Table 5.2 show the maximum and the average communication volume per process. This involves the send and receive volumes at Lines 13 and 14 of Algorithm 13. Using the fine-hp and fine-rd partitions, the average communication volumes are 11K and 1735K, respectively. Recall that the average communication volume is proportional to the cutsize of the corresponding hypergraph partition, and the cutsize equals to the total redundancy in MxV and MTxV operations. That is why we observe higher W TRSVD value using the fine-rd partition. The maximum communication volume per process also decreases to 166K with fine-hp partition, in comparison to 1744K using fine-rd partition. Our final observation is that fine-hp partitions are more effective in reducing the communication volume than the coarse-hp partitions.

In Table 5.3, we provide the relative timings of TTMc, TRSVD, and the computation of the core tensor within an iteration of HOOI using the 256-way fine-hp partition. The TRSVD timings include the time spent in the communication of the vector entries in the MxV and MTxV operations, as the communication takes place within the PETSc calls. We see in Table 5.1 that on Netflix tensor with the fine-hp partition, we lose scalability at 256 nodes. Table 5.3 shows that for this instance TRSVD step begins to dominate the timings, and the communication cost starts to prevent further scalability. Note that in this case, the increase in the communication cost and imbalance also affects the computational cost of the TRSVD step. This is because of the fact that for each communicated vector entry there is an associated unit of redundant work in the MxV and MTxV computations. We also realize in these results that the cost of forming the core tensor G with a TTM followed by an ALLREDUCE communication is negligible, as we expected. Finally, we inform that in all instances, TRSVD (as provided by SLEPc) converged in less than 5 iterations.

Shared memory results

We evaluate the shared memory scalability of the distributed memory HOOI algorithm by varying the number of threads from 1 to 32, and using the minimum number of nodes possible. We needed 8, 8, 4, and 1 nodes to be able to execute the code on Delicious, Flickr, Netflix, and NELL, respectively, using the fine-hp partitions.

TTMc is a memory latency-bound operation; for each nonzero x i,j,k , the access to U (1) (i, :), U (2) (j, :), and U (3) (k, :) likely results in a cache miss, due to the irregular pattern of nonzeros. Multi-threading is an effective way of hiding this latency; therefore it offers a great opportunity of acceleration through parallelism. However, the MxV and MTxV operations in the TRSVD step are memory bandwidth-bound due to the matrices being dense. Once the memory bandwidth is saturated, one may not expect a notable speedup with multi-threading (except in a NUMA architecture where each socket has an independent memory bandwidth; yet in this case the parallelism within a socket has the same issue). The TTMc operation count is proportional to the number of nonzeros, whereas the TRSVD cost is proportional to the number of rows of the matrix (or equivalently, the size of a dimension of the tensor).

As seen in Table 5.4, we manage to improve the runtime using 32 threads for all tensors except Flickr. Using 32 threads, the speedups we obtain for Delicious and Flickr tensors are 7.2x and 5.1x, whereas on NELL and Netflix we get 9.8x and 20x, respectively. Delicious and Flickr tensors have very large third dimension of size 17M and 28M, whereas the largest dimensions of NELL and Netflix are of size 3.2M and 480K. Therefore, NELL and particularly Netflix have more latency-bound computations and provide more room for speedup, which explains the better speedup results in comparison to Delicious and Flickr tensors. Another interesting point is that on the Netflix tensor, using 16 threads we achieve 13.8x speedup over the single threaded execution. Increasing to 32 threads results in a superlinear speedup of 20x on 16 cores. We believe that this is mostly due to each core being able to execute two threads (one for memory and one for compute operations) simultaneously, which is particularly advantageous for latency-bound sparse irregular operations.

The effect of dimension trees

We now compare the sequential and parallel execution timings of HOOI with and without dimension trees to evaluate its effect. These timings are provided in Table 5.5 on Ada cluster for sequential as well as parallel executions using 256 ranks and 8 threads. Here, "flat" refers to the standard nonzerobased TTMc computation given in Algorithm 10, whereas "dtree" represents the dimensional treebased computational scheme.

First, we note that in all instances, using a dimension tree yields faster execution. In the sequential execution, dimension tree provides 35%, 36%, 19%, and 14% speedup over the traditional scheme on Netflix, NELL, Delicious, and Flickr tensors, respectively, wheres in parallel runs it results in 13%, 8%, 215%, and 16% faster executions. This clearly shows that storing and reusing partial TTMc results indeed reduces the computational cost notably both in sequential and parallel runs.

Conclusion

We discussed an efficient parallelization of the alternating least squares-based Tucker decomposition algorithm (HOOI, also called Tucker-ALS) for sparse tensors in shared and distributed memory systems. We introduced a nonzero-based TTMc formulation, and proposed a shared-memory parallel HOOI with a preceding symbolic computation step that uses this formulation. We proposed a coarse and a fine-grain parallel algorithm with their corresponding task definitions, and investigated the issues of load balance and communication cost reduction on different components of the parallel Chapter 6

Parallel Sparse NMF

Though the main focus of this thesis is mainly on tensor decompositions, this chapter investigates an efficient computation of a widely used tool in data analysis and machine learning applications, namely non-negative sparse matrix factorization, whose parallelization involves striking similarities with parallel sparse tensor decompositions presented in the preceding chapters. This opens up to possibility of properly adopting these techniques for obtaining a fast distributed sparse non-negative matrix factorization algorithm, which is the ultimate goal of the work presented in this chapter.

where kXk F =(

P ij x 2 ij ) 1/2
is the Frobenius norm. The functions φ(W) and (H) are called regularization functions that prevent the model from overfitting the data, and these functions are chosen depending on the characteristics of the data. In our case, we are considering φ(W)=↵kWk 2 F , called `2 regularizer, and (H)=β P n i=1 kh i k 2 1 , called `1 regularizer, for addressing the inherent sparsity in the input matrix. The lack of `1 regularization on the matrix H can result in numerical instability.

NMF is widely used in data mining and machine learning as a dimensionality reduction and factor analysis method. It is a natural fit for many real world problems as the non-negativity is inherent in many representations of real-world data, in which case the resulting low rank factors are expected to have a natural interpretation. The applications of NMF range from text mining [START_REF] Pauca | Text mining using nonnegative matrix factorizations[END_REF], computer vision [START_REF] Oyer | Non-negative matrix factorization with sparseness constraints[END_REF], and bioinformatics [START_REF] Kim | Sparse non-negative matrix factorizations via alternating nonnegativity-constrained least squares for microarray data analysis[END_REF] to blind source separation [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-Way Data Analysis and Blind Source Separation[END_REF], unsupervised clustering [START_REF] Kuang | Symmetric nonnegative matrix factorization for graph clustering[END_REF][START_REF] Kuang | SymNMF: Nonnegative low-rank approximation of a similarity matrix for graph clustering[END_REF], and many others. In most real-world applications, m and n can be on the order of millions or more while k being much smaller in the order of tens to thousands.

We would like to highlight that "non-negative" matrix factorization is not the matrix factorization in collaborative filtering for recommender systems, for which many implementations exist. The collaborative filtering problem is different than NMF -the focus of this chapter-because it interprets the "zero" entries of the input matrix as missing data, while the NMF problem is defined for a completely known input matrix and does not handle missing values. This leads to different optimization problems, algorithms, and computational steps. Specifically, computing NMF involves three main types of operations; multiplication of A with a factor matrix W or H, computing the Gram matrices W T W and HH T , and solving a non-linear least squares (NNLS) problem using these two resulting dense matrices to update factor matrices, which can be performed directly (all columns at once) or column-by-column within the ALS framework [2]. Also, while ALS framework is popularly used for NMF, direct techniques updating all factors at once are known to be more robust particularly in ill-conditioned cases albeit requiring higher per-iteration computational cost [START_REF] Debals | Nonnegative matrix factorization using nonnegative polynomial approximations[END_REF]. Nevertheless, the approaches share similarities from a computational perspective, hence the contributions in this chapter remain pertinent in both cases.

To the best of our knowledge, the only high performance software available for parallelizing this computation is MPIFAUN [START_REF] Kannan | MPI-FAUN: an MPI-based framework for alternating-updating nonnegative matrix factorization[END_REF] which operates on both sparse and dense input matrices. For parallelism, it employs a fixed 2D uniform partitioning on A, and partitions factor matrices W and H conformally with this 2D partition. Each process works on its local matrix block of A, and communicates the rows and columns of W and H corresponding to its process row and column in the 2D topology. This communication is optimal in the dense case in the sense that no process receives a data element not used in its local computations. In the sparse case, however, the actual communication requirements of processes depend on the sparsity of A and its partitioning, and is in general significantly less than the dense case. Therefore, employing the same collective communication scheme in the sparse case brings about a major redundancy in the communication cost. Also, MPIFAUN similarly uses only fixed 2D uniform partitioning on sparse matrices; this creates significant imbalance both in computation and communication as typically the data is not homogenously distributed in A. In overall, the lack of an efficient sparse communication scheme that allows flexibility in partitioning hinders the scalability of this approach.

We summarize the contributions in this chapter as follows:

• We propose an efficient parallel NMF algorithm called DIST-SPNMF for sparse matrices that employs a point-to-point communication scheme to leverage the sparsity of the input matrix, and eliminate any redundant communication existing in MPIFAUN. The algorithm is flexible to work with any partition of the input matrix A as well as factor matrices W and H.

• We employ a 2D cartesion process topology in a way that eliminates the communication cost of parallel NMF in one dimension while bounding the maximum number of messages sent per process in our MPI implementation.

• We introduce effective partitioning strategies to further reduce the communication cost and enhance the parallel scalability of our algorithm.

• We introduce regularization to NMF computations to prevent potential numerical instabilities using large matrices.

• We compare our implementation with MPIFAUN, and report scalability results up to 32678 processors on two different parallel computing platforms using two real-world datasets.

Table 6.1 summarizes the notation we use throughout this chapter. We use subscripts to refer to subblocks of matrices. For example, A ij refers to the sub-block (i, j) of A in a 2D partition. We use m and n to denote the numbers of rows and columns of A, respectively, and assume w.l.g that m > n throughout.

A

Sparse non-negative matrix factorization

The Alternating-Updating NMF (AU-NMF) algorithms are those that alternate between updating one of W and H using the given input matrix A and other 'fixed' factor -H for updating W, and W for updating H. This update is performed using the Gram matrix associated with the fixed factor matrix (HH T for W and W T W for H), and the product of the input matrix A with the fixed factor matrix (AH T or W T A). We show this framework in Algorithm 14.

After computing the Gram matrix and the multiplication of A with the fixed factor matrix, the specifics of the updates at Lines 3 and 4 depend on the NMF algorithm, and we refer to the computation associated with these lines as the Local Update Computations (LUC).

Algorithm 14 AU-NMF: Alternating-Updating NMF algorithm

Input: A: An m ⇥ n matrix k: The rank of approximation Output: [W, H]: Factor matrices 1: Initialize H with a non-negative matrix in R n⇥k + . 2: while stopping criteria not met do 3:

Update W using HH T and AH T

4:

Update H using W T W and W T A We note that AU-NMF is very similar to a two-block, block coordinate descent (BCD) framework as explained by Bertsekas [START_REF] Bertsekas | Nonlinear programming[END_REF]. The BCD framework expresses solving optimization variables in complex non-linear optimization problem as one block at a time, while keeping the others fixed. In NMF, the two blocks are the unknown factors W and H, and we solve the following subproblems,

W argmin W>0 A -WH F + φ( W)+ (H), H argmin H>0 A -W H F + φ(W)+ ( H). (6.2)
Since each subproblem involves non-negative least squares, this two-block BCD method is also called the Alternating Non-negative Least Squares (ANLS) method [START_REF] Kim | Algorithms for nonnegative matrix and tensor factorizations: A unified view based on block coordinate descent framework[END_REF]. Block Principal Pivoting (ABPP) is one algorithm that solves these NLS subproblems. In the context of the AU-NMF algorithm, an ANLS method maximally reduces the overall NMF objective function value by finding the optimal solution for given H and W in Lines 3 and 4, respectively.

From time to time, these updates do not necessarily solve each of the subproblems (6.2) to optimality but simply improve the overall objective function (6.3), such as in Multiplicative Update (MU) [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF] and Hierarchical Alternating Least Squares (HALS) [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-Way Data Analysis and Blind Source Separation[END_REF] methods.

The convergence properties of these different NMF algorithms are discussed in detail by Kim, He and Park [START_REF] Kim | Algorithms for nonnegative matrix and tensor factorizations: A unified view based on block coordinate descent framework[END_REF]. While we focus only on the MU algorithms in this paper, we highlight that our algorithm is not restricted to this, and is seamlessly extensible to other NMF algorithms as well, including HALS, ABPP, Alternating Direction Method of Multipliers (ADMM) [START_REF] Sun | Alternating direction method of multipliers for non-negative matrix factorization with the beta-divergence[END_REF], and Nesterovbased methods [START_REF] Guan | Nenmf: An optimal gradient method for nonnegative matrix factorization[END_REF].

Multiplicative update

As we are considering `2 regularization on the matrix W and `1 regularization on the matrix H to address the sparsity of the input matrix, the NMF problem becomes

min W>0,H>0 kA -WHk F + ↵(kWk) 2 F + β n X i=1 kh i k 2 1 . (6.3)
The values ↵ and β are fixed for the experiments. In the case of MU [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], individual entries of W and H are updated with all other entries fixed. In this case, the update rules are

w ij w ij (AH T ) ij (W(HH T +2β1 k )) ij , and
h ij h ij (W T A) ij ((W T W +2↵I k )H) ij . (6.4)
where 1 k is a matrix of k ⇥ k with all one's and I k is an identity matrix of size k ⇥ k.

After computing the Gram matrices HH T and W T W and adding the appropriate regularizers and the products AH T and W T A, the extra cost of computing W(HH T +2β1 k ) and

(W T W + 2↵I k ) is F (m, n, k) = 2(m + n)k 2 flops
to perform updates for all entries of W and H, as the other elementwise operations affect only lower-order terms. The details about using AU-NMF in Algorithm 14 for other algorithms HALS and ABPP are explained in [START_REF] Kannan | A high-performance parallel algorithm for nonnegative matrix factorization[END_REF][START_REF] Kannan | MPI-FAUN: an MPI-based framework for alternating-updating nonnegative matrix factorization[END_REF]. This update in (6.4) can be easily parallelized using dense matrix kernels.

It is important to observe that the update function of W is element-wise normalization of the sparse matrix-dense matrix multiplication AH T with the denominator. Given that all the processes own the k ⇥ k gram of the factor matrix H, computing the entire denominator (W(HH T +2β1 k )) does not require any communication, hence can be done locally. One can argue that the same holds for updating H as well. Therefore, for row and column index sets I p and J p , one can update these rows of the factor matrices as follows:

W(I p , :) W(I p , :) ~( W(I p , :)) ↵ (W(I p , :)(G H +2β1 k )), and
H(J p , :) H(J p , :) ~( H(J p , :) ↵ (G W +2↵I k )H(J p , :)). (6.5)
where ~and ↵ correspond to element-wise multiplication and division of matrices or vectors. This scheme provides row-wise parallelism in NNLS computation. HALS and ABPP can similarly be expressed in this row-parallel form.

Parallel sparse NMF

Here, we first introduce our parallel NMF algorithm that operates on a partition of the matrices A, W, and H. For a given partition, we describe how parallel computations and communications take place within the algorithm, and illustrate computational and communication costs associated with a partition. In doing so, we elaborate the advantages of our parallelization scheme compared to the state of the art. We then discuss efficient partitioning strategies to better establish computational load balance and reduce communication.

Distributed parallel sparse NMF algorithm

Parallelizing NMF involves partitioning data and computational tasks to processes, and communicating results wherever necessary, just as any parallel algorithm. Specifically for NMF, one needs to partition the sparse matrix A as well as the factor matrices W and H, where the former partitioning distributes the computational load of sparse matrix-dense matrix multiplications AH and W T A, whereas the latter divides the workload of NNLS computations to processes. We provide the execution of our parallel algorithm for computing a rank-k NMF of a sparse matrix A 2 R m⇥n in Algorithm 15, which is executed by each process p for 1  p  P . The algorithm starts with an arbitrary partition of the input matrix and the factor matrices; process p owns the submatrices W(I p , :) and H(:, J p ) as well as the nonzero elements of the sparse matrix A p where A = S P i=1 A i , i.e., A 1 ,...,A p partitions the nonzeros of A. The sets F p and G p denote the "footprints" of the process p on the rows and columns of matrices W and H, respectively, in the course of the algorithm. That is, we have i 2F p or j 2G p if only if i 2I p or j 2J p (row/column is owned), or there is a nonzero element a i,j 2 A p (row/column is used in local computations). At each iteration, the process p is responsible for gathering the new value of submatrices W(I p , :) and H(:, J p ), and sending these values to processes in need.

In an iteration of Algorithm 15, each process p possesses three types of computational tasks as well as associated pre-and post-communication steps. The first task involves performing sparse matrixdense matrix multiplications A p H(:, G p ) and W(F p , :) T A p , whose results are stored in distributed matrices W and H, which follow the same row-/column-wise data distribution as W and H. Note that carrying out these multiplications must be preceded by a communication step where each process COMM-EXPAND(H(:, G p ))

3:

W(F p , :) A p H(:, G p ) 4:
COMM-FOLD( W(F p , :))

G H ALL-REDUCE(H(:, J p )H(:, J p ) T ) 6:

W(I p , :) NNLS(G H , W(I p , :)) 7:
COMM-EXPAND(W(F p , :))

8:

H(:, G p ) W(F p , :) T A p 9:
COMM-FOLD( H(:, G p ))

10:

G W ALL-REDUCE(W(I p , :) T W(I p , :))

H(J p , :) NNLS(G W , H(J p , :)) 12: until convergence or maximum number of iterations p gets the submatrices H(:, G p \J p ) and W(F p \I p , :) that are accessed by entries of A p , and this step is performed at Lines 2 and 7. These multiplications performed by each process p generate partial results for the set F p and G p of the rows of W and the columns of H, respectively, which is highlighted at Lines 3 and 8. Indeed, partial results for the submatrices W(F p \I p , :) and H(:, G p \J p ) correspond to rows and columns owned by other processes; hence, they need to be communicated. The results for W(I p , :) and H(:, J p ), however, should be kept locally, and all partial results for these matrix rows and columns generated by other processes should similarly be received and accumulated in order to obtain the final value for these owned portions. The second task is to compute the Gram matrices G H = HH T and G W = W T W of size k ⇥ k, and to make these matrices available to all processes, which is performed at Lines 5 and 10. This is done in a row-parallel dense matrix multiplication step, in which the process p computes H(:, J p )H T (:, J p ) and W T (I p , :)W(I p , :), followed by an ALL-REDUCE communication of these partial multiplications. The third task pertains to updating the factor matrices W and H using matrices W and G H , or H and G W , which takes place at Lines 6 and 11. This corresponds to Lines 3 and 4 of Algorithm 14, and can be computed locally at each process p by executing NNLS algorithm on dense matrices W(I p , :) and G H , or H(:, J p ) and G W , to obtain new W(I p , :) or H(I p , :), respectively, as described in Section 6.1.1.

The first type of communication in Algorithm 15 pertains to an ALL-REDUCE of a dense matrix of fixed size k ⇥ k at Lines 5 and 10, and the cost of this step is typically negligible in compare to the rest. The other two communication types involve (i) transferring the partial row results of W and H to their owner processes at Lines 4 and 9 to accumulate at the owners, (ii) sending the updated rows of W and H to processes in need at Lines 2 and 7. We respectively call these steps fold and expand communications, following the convention used by the sparse matrix community. The way these two communications are carried out plays a vital role in obtaining parallel scalability as they dominate the communication cost of the algorithm. Here, we employ an efficient point-to-point communication framework to perform fold and expand communications. In this scheme, after determining a partition of the matrix, each process p determines the sets F p and G p using the row and column indices in its local matrix A p and its sets I p and J p of owned row indices, then identifies the owners of these rows in order to request these rows from the corresponding processes in the expand step. Next, these requests are exchanged in a communication step so that each process p determines the unique row indices that it needs to send to or receive from each other process in expand and fold steps. Once this communication structure is formed, it is reused throughout all iterations in Algorithm 15. This is possible thanks to the data partition staying fixed during the algorithm, hence data dependencies do not change.

The existing work in the literature for parallel sparse NMF [START_REF] Kannan | A high-performance parallel algorithm for nonnegative matrix factorization[END_REF][START_REF] Kannan | MPI-FAUN: an MPI-based framework for alternating-updating nonnegative matrix factorization[END_REF] amounts to employing a collective communication strategy for both expand and fold steps of Algorithm 15 for dense as well as sparse A. In this strategy, the rows and columns indices 1,...,mand 1,...,nare divided into P r and P c (P = P r P c ) sets I 1 ,...,I Pr and J 1 ,...,J Pc of equal size and having contiguous indices. This effectively defines a 2D uniform checkerboard partition of A with matrix blocks A (r,c) = A(I r ,J c ) for all 1  r  P r and 1  c  P c , as well as a row partition of W and H. Here, process p owns the matrix subblock A (r,c) , where r = bP/P c c +1and c =( P mod P c )+1, as well as m/P and n/P rows of W(I r , :) and H(J c , :). The communication is performed within each process row and column using ALL-GATHER and REDUCE-SCATTER routines, in which the process p receives all matrix rows W(I r , :) and H(I c , :) belonging to processes in the same row and column of the process grid. Despite being favorable due to small number of exchanged messages in collective routines, in this strategy processes receive many rows that they do not need in their local sparse matrix dense matrix multiplication, and this redundancy dramatically increases the communication volume, thus preventing scalability. To illustrate, in Figure 6.1 we provide a sparse matrix partitioned using a 5 ⇥ 5 checkerboard scheme where each process owns a subblock of the matrix, and show only the nonzero elements belonging to three columns with indices j 1 , j 2 , and j 3 . Here, in performing the expand communication of Line 2, the collective communication strategy sends these column vectors to all processes in the same column, whereas point-to-point communication identifies the exact set of processes in need beforehand, indicated with the same color, and communicates the data only with these processes. In this example, the collective communication scheme incurs 4, 3, and 2 units of redundant communication for columns j 3 , j 2 , and j 1 , respectively, and this redundancy only worsens with the increasing number of processes.

Partitioning

Algorithm 15 requires a partitioning of the nonzeros of A as well as the rows and columns of W and H, and these three partitions completely determine its computational and communication costs. Here, we compare different partitioning strategies, and argue how they relate to these two performance metrics.

Our algorithm is fine-grained; meaning that it can work any fine-grain, 1D, or 2D/checkerboard partition of the sparse matrix A. Here, however, we exclusively focus on P r ⇥ P c (where P = P r P c ) checkerboard partitionings, and always set P c to the number of cores available per compute node in the cluster. This helps us reduce the communication cost in two ways. The first one corresponds to that we do not employ shared-memory parallelism in our implementation, hence assign one MPI rank per core. In this scenario, by assigning each group of P c processes owning the same row block of A to the same node, and similarly assigning the corresponding row blocks of W and W to these processes, we effectively restrain the communication to stay within the same node in the row dimen- sion of the process grid (which pertains to communicating matrices W and W), whose cost become negligible (as it is performed by memory copies). This enables us to focus solely on reducing the communication due to columns of H and H. Second, this topology bounds the maximum number of messages sent per process by P r -1 in the column dimension, thereby significantly reduces the communication latency by a factor of P c in compare to a fine-grain or a 1D partition, which incur up to P -1 messages per process.

Partitioning sparse input matrix A

Hypergraph partitioning (P2PHP) We aim to partition the matrix A into P r row slices first, and P c column slices next, to obtain an P r ⇥ P c checkerboard partition. To this end, we adopt hypergraph partitioning techniques used in parallel sparse matrix literature [START_REF] Çatalyürek | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF], and modify them appropriately to match the requirements of sparse NMF computations.

The typical checkerboard hypergraph partitioning approach used for parallel sparse matrix-vector multiplication (SpMV) proceeds in two steps as in what follows. First, a column-net hypergraph model is formed; for each row i and column j of A, a vertex v i and a hyperedge h j are added to the hypergraph, and for each nonzero (i, j) 2 A v i is connected to h j . Note that in the sparse NMF case, this connection also models the dependency to H(:,j) in computing W(i, :) at Line 3. Next, the weight |A(i, :)| is assigned to each vertex v i , which represents the computational load for the SpMV. In partitioning this hypergraph into P r parts, minimizing the cutsize of the hypergraph corresponds to minimizing the total communication volume due to column vector entries, while establishing the balance for vertex weights corresponds to balancing the computational load of SpMV. After obtaining the row partition this way, a second multi-constraint hypergraph partitioning step is performed on the row-net hypergraph model; for each column j and row i of A, a vertex v j and a hyperedge h i are added to the hypergraph, and for each nonzero (i, j) 2 A, v j is connected to h i . Similarly, this connection models the dependency to W(i, ) in computing H(:,j) at Line 8. This time, each vertex v j is assigned P r weights, each corresponding to the number of nonzeros corresponding to different row partitions in the jth column of A. Partitioning this hypergraph into P c parts finally defines a P r ⇥ P c checkerboard topology. Minimizing the cutsize of this hypergraph corresponds to minimizing the total row communication volume, whereas balancing vertex weights corresponds to balancing the number of nonzeros of A assigned to each part (SpMV load).

In the context of parallel sparse NMF, we modify both row and column partitioning steps to better grasp its computational and communication requirements. Note that in the first step, the traditional scheme balances the number of nonzeros per row slice, which in turn balances the workload for multiplying A with dense matrices W and H. However, this can still yield a significant imbalance in the number of rows assigned to each row slice. In the SpMV case, this incurs an imbalance in the number of vector entries assigned to each row slice, which does not cause a significant problem as the cost of vector operations is typically negligible in an iterative solver. In the NMF case, however, this causes a computational imbalance in the NNLS step involving dense matrices W and W, which is a major computational step with a non-negligible cost, particularly when the rank of approximation gets higher. To achieve load balance in the number of rows per row slice, one can similarly employ multi-constraint partitioning in this step as well, by assigning a second fixed weight of 1 to each vertex, then balancing both costs in the partitioning. This formulation is expected to achieve balance in both SpMM and NNLS steps. In practice, however, multi-constraint hypergraph partitioning is known to not work too good in practice; therefore, we aim to achieve the same effect by assigning a single weight (nnz(A)/m + nnz(A(i, :)) to each vertex v i . This weighting is preferred to find a tradeoff between balancing the number of vertices per row slice, which corresponds to the cost of dense matrix operations, and balancing the number of nonzeros per part, which corresponds to the cost of SpMM. We establish these two load balance goals by assigning a fixed cost to each vertex in the first summand, and increasing the cost of the vertices corresponding to denser rows in the second summand. In Figure 6.1 we show this hypergraph construction for a 3 ⇥ 3 subset of the matrix A, in which the rows corresponding to each vertex are highlighted, and the vertices are connected to related hyperedges. The vertex weights are omitted in the figure for simplicity.

In the second partitioning step, the traditional scheme balances the number of nonzeros per part, which similarly balances the workload for multiplying A with dense matrices W and H, while minimizing the total communication volume due to vector rows. This approach, however, is not ideal for the sparse NMF for multiple reasons. First, while balancing the sparse matrix multiplication load, it similarly disregards the NNLS work proportional to the number of columns assigned per column slice. Second, it employs multi-constraint partitioning with P r constraints which can be quite high, in which case multi-constraint partitioners are known to yield poor load balance even for the SpMM step. Finally, it tries to minimize the total communication volume due to rows, whose cost is expected to be negligible as it is restrained to stay within a compute node. The important communication objective here is balancing the column communication volume per column slice, which is neglected in this approach.

In overall, in the second partitioning phase we can discuss 2P r +1weights corresponding to the computation and communication induced by each column: one fixed weight of value 1 to balance the number of columns per column slice for NNLS; P r weights for the number of nonzeros in the column correspond to each row slice; P r weights for the volume communication incurred to each row slice by the column. Our objective here is to partition these columns in a way that all 2P r +1weights are independently balanced. This is essentially a knapsack problem with multiple weights, which is NP-hard. One can employ a specific heuristic for this partitioning; however, the computatational and memory requirements might prevent this approach from being practical for a high n and P r value. We instead partition the columns randomly, which is expected to give acceptably good balance for all these weights in practice. Indeed, this approach disregards the row communication volume; yet this cost is expected to be unimportant.

Randomized checkerboard partitioning(P2PRP)

This scheme corresponds to partitioning both the rows and columns of A into R segments randomly. It is expected to provide good load balance both in sparse and dense matrix operations, but it overlooks the communication cost.

Uniform checkerboard partitioning(FAUN)

This partitioning variant forms an R ⇥ C partition of A by putting a contiguous set of m/R and n/C rows and columns in each slice. W and H are partitioned conformally with this topology; each process is assigned a contiguous set of m/P and n/P rows and columns, respectively, belonging to the corresponding row and column blocks of H. This is the only partitioning scheme employed by MPIFAUN [START_REF] Kannan | A high-performance parallel algorithm for nonnegative matrix factorization[END_REF][START_REF] Kannan | MPI-FAUN: an MPI-based framework for alternating-updating nonnegative matrix factorization[END_REF], hence we only use it for this implementation. We also use a randomized variant (FAUNRP) of this scheme in which the rows and columns of A are permuted randomly before executing MPIFAUN to better balance its nonzeros.

Partitioning factors W and H

Once A is partitioned with an R⇥C checkerboard topology, one has to partition the rows and columns of factor matrices conformally with this topology to form the sets I p and J p in Algorithm 15. In doing so, we are interested in assigning rows and columns to processes in the corresponding process row or column equitably, as this balances NNLS work as well as the memory cost for storing a part of W and H. For this purpose, we specify imbalance parameters ↵ r and ↵ c that correspond to the maximum imbalance we allow in this partitioning; i.e., |I p |↵ r m/P and |J p |↵ c n/P for each process p. We normally set these imbalance parameters to 1.1. However, when m>>n, we increase ↵ c proportionally in order to give more flexibility in column partitioning which results in smaller communication volume.

Next, for each row and column of W and H, we create a list of processes that have a dependency to that row or column, which corresponds to processes owning the matrix blocks of the same color in Figure 6.1. Then, we are to assign each row/column to, preferably, one of these dependent parts. Ideally, we would like to avoid a communication imbalance due to a significantly higher send or receive volume in either of fold and expand communications. We argue that randomly assigning a row or column to one of its dependent parts establishes this communication balance for the following reason. Consider any part dependent to column j, and let h j be the hyperedge connecting the parts for this column. In the random assignment, this part owns column j with a probability of 1/λ(h j ), which in turn incurs a receive/send volume of λ(h j ) -1 in the fold/expand communication. With (λ(h j )-1)/λ(h j ) probability, the part does not own the column j, which yields a send/receive volume of 1 in the fold/expand communication. In this case, both the expected send and receive volumes of fold and expand communications due to column j are (λ(h j )-1)/λ(h j ); therefore, this assignment is expected to balance the overall send and receive volumes of all parts in both communication steps. To satisfy the balance constraints, however, we do this random assignment only to one of the dependent parts satisfying the imbalance constraint. If all processes in the dependence list are overloaded, we randomly assign it to a process that is not overloaded in the same process row or column. Note that the latter assignment increases the communication volume due to that row or column by 1; hence, larger imbalance parameters yield smaller communication volume by reducing this type of assignment in general.

Related work

In the data mining and machine learning literature, there is an overlap between low rank approximations and matrix factorizations due to the nature of applications. Despite its name, non-negative matrix "factorization" is in fact a low rank approximation. Recently, there has been a growing interest in collaborative filtering based recommender systems. One of the popular techniques for collaborative filtering is matrix factorization, often with nonnegativity constraints, and its implementation is widely available in many off-the-shelf distributed machine learning libraries such as GraphLab [START_REF] Ow | Distributed GraphLab: A framework for machine learning and data mining in the cloud[END_REF], MLLib [START_REF] Meng | MLlib: Machine Learning in Apache Spark[END_REF], and many others [START_REF] Satish | Navigating the maze of graph analytics frameworks using massive graph datasets[END_REF][START_REF] Yun | Nomad: Nonlocking, stochastic multi-machine algorithm for asynchronous and decentralized matrix completion[END_REF]. However, we would like to emphasize that collaborative filtering using matrix factorization is a different problem than NMF: In the case of collaborative filtering, nonzeros in the matrix are considered to be observed ratings, and zeros are treated as missing entries, while in the case of NMF, there is no missing entries and all zeros are considered as observed entries.

There are several recent distributed NMF algorithms in the literature [START_REF] Faloutsos | Flexi-FaCT: Scalable flexible factorization of coupled tensors on Hadoop[END_REF][START_REF] Liao | CloudNMF: A MapReduce implementation of nonnegative matrix factorization for large-scale biological datasets[END_REF][START_REF] Liu | Distributed nonnegative matrix factorization for web-scale dyadic data analysis on MapReduce[END_REF][START_REF] Yin | Scalable nonnegative matrix factorization with block-wise updates[END_REF]. Liu et al. propose running Multiplicative Update (MU) for KL divergence, squared loss, and "exponential" loss functions [START_REF] Liu | Distributed nonnegative matrix factorization for web-scale dyadic data analysis on MapReduce[END_REF]. Matrix multiplication, element-wise multiplication, and element-wise division are the building blocks of the MU algorithm. The authors discuss performing these matrix operations efficiently on Hadoop for sparse matrices. Using similar approaches, Liao et al. implement an open source Hadoop-based MU algorithm and study its scalability on large-scale biological data sets [START_REF] Liao | CloudNMF: A MapReduce implementation of nonnegative matrix factorization for large-scale biological datasets[END_REF]. Also, Yin, Gao, and Zhang present a scalable NMF that can perform frequent updates, which aim to use the most recently updated data [START_REF] Yin | Scalable nonnegative matrix factorization with block-wise updates[END_REF]. Similarly, Faloutsos et al. propose a distributed, scalable method for decomposing matrices, tensors, and coupled data sets through stochastic gradient descent on a variety of objective functions [START_REF] Faloutsos | Flexi-FaCT: Scalable flexible factorization of coupled tensors on Hadoop[END_REF]. The authors also provide an implementation that can enforce non-negative constraints on the factor matrices. All of these works use Hadoop framework to implement their algorithms, hence are not very efficient.

Spark [START_REF] Zaharia | Spark: Cluster computing with working sets[END_REF] is a popular big-data processing infrastructure that is generally more efficient for iterative algorithms such as NMF than Hadoop, as it maintains data in memory and avoids file system I/O. Even with a Spark implementation of previously proposed Hadoop-based NMF algorithms, the performance still suffers from expensive communication of input matrix entries, and Spark does not have innate mechanisms to overcome this shortcoming. Spark has collaborative filtering libraries such as MLlib [START_REF] Meng | MLlib: Machine Learning in Apache Spark[END_REF] which use matrix factorization and can impose non-negativity constraints.

In parallel with the Hadoop and Spark implementations, there has been a growing interest in the HPC community towards efficiently computing these algorithms with tuned high performance implementations. Kannan, Ballard and Park [START_REF] Kannan | A high-performance parallel algorithm for nonnegative matrix factorization[END_REF][START_REF] Kannan | MPI-FAUN: an MPI-based framework for alternating-updating nonnegative matrix factorization[END_REF] proposed MPIFAUN framework to implement various NMF algorithms such as multiplicative update (MU), Hierarchical Alternating Least Squares (HALS), and Alternating Non-negative Least Squares using Block Principal Pivoting (ANLS-BPP). We choose this work as a baseline, as it is the only available high performance implementation of NMF, and it performs significantly faster than Hadoop and Spark-based approaches. To elaborate this, Gittens et.al. [START_REF] Gittens | Matrix factorization at scale: A comparison of scientific data analytics in Spark and C+MPI using three case studies[END_REF] recently benchmarked the implementations of different matrix factorization algorithms (such as NMF and Principal Component Analysis (PCA)) in Spark and in MPI/C. They claim that native MPI implementations on HPC platforms outperform Spark implementation by a speedup factor of 44x. Similar observations have been made by Sukumar, Kannan, Matheson and Lim [START_REF] Sukumar | Kernels for scalable data analysis in science: Towards an architecture-portable future[END_REF][START_REF] Sukumar | Mini-apps for high performance data analysis[END_REF] on supercomputers at Oak Ridge Leadership Computing Facility. Finally, there are implementations of the MU algorithm in a distributed memory setting using X10 [START_REF] Grove | Supporting array programming in X10[END_REF], and on a GPU [START_REF] Mejía-Roa | NMF-mGPU: Non-negative matrix factorization on multi-GPU systems[END_REF].

Experiments

In this section, we compare our algorithm DIST-SPNMF against MPIFAUN on two big sparse matrices formed from real world datasets. We analyze and compare the computation and the communication timings of these algorithms on a smaller cluster, then test the scalability limits of our method on a larger supercomputing environment.

Experimental setup

Datasets

We use two matrices formed from Flickr.com and Delicious.com that involve images tagged with different labels by users. These are essentially matrices obtained from Delicious and Flickr tensors described in Section 2.6 by discarding their first and second dimensions. Flickr and Delicious matrices are of size 28Mx1.6M and 17Mx2.5M, and have 112M and 72M nonzero elements, respectively. The current implementation of MPIFAUN can only operate when R and C can divide m and n, hence we trimmed the matrices slightly.

Parallel platform

We conducted our experiments on two different parallel computing platforms. The first platform is the "Rhea" cluster at the Oak Ridge Leadership Computing Facility (OLCF), which is a commodity-type Linux cluster with a total of 512 nodes and a 4X FDR Infiniband interconnect. Each node contains dual-socket 8-core Intel Sandy Bridge-EP processors operating at 2GHz clock frequency, and 128 GB of memory. Each socket has a shared 20MB L3 cache, and each core has a private 256K L2 cache. There, we ran our experiments up to 3072 cores, which is the maximum allowed in the cluster. The second platform is Blue Gene/Q described in Section 2.7. We ran both algorithms using 16 MPI ranks per node, and set P c = 16 in all partitionings.

Our code for local matrix operations is developed using the matrix library Armadillo [START_REF] Sanderson | Armadillo: An open source C++ linear algebra library for fast prototyping and computationally intensive experiments[END_REF]. We use BLAS and LAPACK for dense matrix operations by linking Armadillo with Intel MKL, Open-BLAS [109], or any other BLAS and LAPACK implementation. Both codes are compiled using the default GNU C++ Compiler (g++ (GCC) 5.3.0) and MPI library (Open MPI 1.8.4) on RHEA, and Clang compiler (3.5.0) with IBM MPI library on Blue Gene/Q.

Algorithms

In our experiments, we considered the following algorithms and partitionings:

• FAUN:M PIFAUN algorithm [START_REF] Kannan | A high-performance parallel algorithm for nonnegative matrix factorization[END_REF][START_REF] Kannan | MPI-FAUN: an MPI-based framework for alternating-updating nonnegative matrix factorization[END_REF] with uniform natural partitioning (FAUN) where each process holds an input matrix of size m/P r ⇥ n/P c .

• FAUNRP: The partitioning strategy in FAUN could result in a significant computational load imbalance with a skewed nonzero distribution of A. We alleviate this by randomly permuting the rows and columns of the matrix before executing MPIFAUN, and call this scheme FAUNRP.

• P2PHP:DIST-SPNMF (Algorithm 15) with hypergraph checkerboard partitioning explained in Section 6.2.2.

• P2PRP:DIST-SPNMF (Algorithm 15) with randomized checkerboard partioning explained in Section 6.2.2.

Strong scaling

In Figures 6.2a and 6.2b, we show the speedup results of all four instances on the Rhea cluster using up to 3072 MPI ranks/cores on Flickr data. The speedup values are with respect to the slowest running time among all four instances using 16 cores (single node). We observe in Figures 6.2a and 6.2b that all algorithms scale up to 1536 cores, yet MPIFAUN instances achieve this with significantly lower parallel efficiency. This mostly is due to higher communication costs involved in the collective communication scheme for both instances. We also realize that FAUNRP significantly improves the runtime with respect to FAUN, meaning that FAUN indeed causes load imbalance in partitioning the nonzeros of A. At 3072 processes, both FAUNRP and FAUN lose scalability and slow down, whereas P2PHP and P2PRP scale to 3072 processors. One point to note is that for k = 16, P2PHP is considerably slower than P2PRP, whereas for k = 48, P2PHP starts to get faster than P2PRP. The former result is due to that P2PHP typically causes more imbalance in factor matrix rows and the input matrix nonzeros, whereas in the latter result, the size of each communication unit is tripled; therefore, the reduction in the communication cost using P2PHP begins to compensate for the performance loss due to load imbalance.

In Figures 6.3a and 6.3b, we provide the same speedup results for Delicious matrix. We observe a similar trend in the comparisons of different methods, except that FAUNRP and FAUN scale even worse in this case. Our algorithm also loses scalability after 1536 processes, and similarly to the previous case P2PHP starts slower than P2PRP due to load imbalance, and catches up for k = 48. These two test cases clearly show that employing a point-to-point communication with effective nonzero and matrix row partitionings is essential for obtaining high performance in the NMF algorithm.

To better test the scalability limits of our algorithms, we ran them on an IBM Blue Gene/Q supercomputer up to 32768 processors using the same two matrices. The results of these two experiments are provided in Figures 6.2c and 6.3c. Our algorithm graciously scales up to 16384 cores in all four instances, and P2PHP manages to slightly improve the runtime using 32768 cores on Flickr, while other three instances witness a slowdown using 32768 cores. Again, P2PHP is slower than P2PRP using lower number of processors as the communication cost is negligible in these instances, and P2PHP introduces worse load balance than P2PRP. However, using 32768 processors P2PHP manages to outrun P2PRP by incurring less communication.

Runtime dissection per iteration

In this section, we provide the time spent on each individual operation type and communication within an NMF iteration. We report the averages over 30 iterations on Rhea, and 10 iterations on Blue Gene/Q for each run. As provided in Algorithm 15, there are three types computations and two types of communications within an NMF iteration, and we present timings for these steps with the following labels:

• Gram: Computing local contribution to Gram matrix, and performing an ALL-REDUCE to gather the final result.

• MM: Computing SpMM using A p and one of the factor matrices.

• LUC : Local NNLS computation to obtain the new factor matrix.

• Comm: Total expand and fold communication time in the case of P2PRP and P2PHP, and the total time spent on ALL-GATHER and REDUCE-SCATTER steps for FAUN and FAUNRP.

In our results, we do not distinguish the costs of these tasks for W and H separately; we instead report their sum.

We report the runtime dissection for Flickr and Delicious matrices in Figure 6.4, Figure 6.5 for Rhea and Figure 6.6 for Blue Gene/Q. For each cluster and data set, we show the timings for the smallest and the largest number of processors used. Our objective in this experiment is to better analyze the speedup results and by comparing the computational and communication costs of different communication schemes and partitionings.

Flickr on Rhea: We observe in Figure 6.4 that in one node configuration with p = 16, the FAUN and FAUNRP performs similar to P2PRP and P2PHP in terms of computation, and the communication time takes a small portion of the execution in all instances. As the number of processes increases to 3072, the communication time of P2PRP and P2PHP stays reasonably low, whereas in the case of FAUN and FAUNRP, we clearly observe that the communication cost dominates the execution time using both low rank values k = 16 and k = 48. Randomization offers load balance to FAUNRP which gives it a slight edge over FAUN, yet both instances suffer from the high communication cost associated with the collective communication strategy, which explains the drop in the scalability results.

Delicious on Rhea: In Figure 6.5, we observe that P2PRP and P2PHP perform better than FAUN even in single node configuration. Figure 6.5 shows that FAUN takes twice more time than P2PRP and P2PHP in the sparse matrix multiplication step. This highlights the skewed distribution of the matrix nonzeros, which is alleviated to a certain extend by randomly permuting the matrix. Similar to Flickr data, using 3072 processors, P2PRP and P2PHP perform significantly better than FAUN and FAUNRP, whose iteration times are dominated by communication cost.

Flickr and Delicious on Blue Gene/Q: In Figure 6.6, we give the timings for computation and communication steps using our methods with two different partitionings of matrices on Blue Gene/Q. We observe that using 512 processors, communication cost is negligible, and P2PRP beats P2PHP thanks to better load balance. Using 16384 processors, however, P2PHP gets faster than P2PRP on Flickr matrix due to significant reduction in the communication volume. On Delicious matrix, P2PHP similarly better reduces the communication, yet this is outweighed by the load imbalance in matrix multiplications.

Conclusion

In this chapter, we proposed a novel algorithm parallel non-negative matrix factorization of sparse matrices in distributed memory environments. To the best of our knowledge, this is the first high performance implementation that takes the sparsity of the input matrix into consideration in communication to reduce the communication cost, and employs effective partitionings to further enhance parallel scalability. We introduce effective partitioning strategies for both the sparse input matrix A as well as factor matrices W and H that helps us achieve good computational load balance while minimizing the communication costs. With all the algorithmic contributions and an efficient parallel implementation, our method outperforms the state of the art by a significant margin, and gracefully scales up to 32768 cores on an IBM Blue Gene/Q supercomputer. Our immediate next steps for extending our work involve adding shared memory parallelism to obtain further speedup.

Chapter 7 On Computing Dense Tensor Decompositions Using Dimension Trees

As already mentioned in the introduction, tensors have been increasingly used in many application domains in the recent past, and in particular, methods for dense tensors have proven to be among the most powerful tools in many signal processing applications [START_REF] Nion | Batch and adaptive PARAFAC-based blind separation of convolutive speech mixtures[END_REF][START_REF] Sidiropoulos | Tensor algebra and multidimensional harmonic retrieval in signal processing for mimo radar[END_REF][START_REF] Sidiropoulos | Parallel factor analysis in sensor array processing[END_REF]. Among these applications are analyzing speech signals for source separation [START_REF] Nion | Batch and adaptive PARAFAC-based blind separation of convolutive speech mixtures[END_REF], finding the localization of the signal source from radar signals [START_REF] Sidiropoulos | Tensor algebra and multidimensional harmonic retrieval in signal processing for mimo radar[END_REF], and other communication applications [START_REF] Sidiropoulos | Parallel factor analysis in sensor array processing[END_REF]. The computational core of CP and Tucker decomposition algorithms involve special operations called the matricized tensor-times Khatri-Rao product (MTTKRP), and tensor-times-matrix multiplication (TTM). Performing these operations on high dimensional big tensor datasets gets computationally expensive in these applications; hence, efficiently carrying out these steps is indispensable to be able to address the computational needs of such applications. Our aim in this chapter is to investigate a fast computation of MTTKRP and TTM steps in CP-ALS and HOOI algorithms for dense tensors using the dimension tree-based computational scheme we introduced in Chapter 3 to achieve this performance goal.

The contributions in this chapter are as follows. We investigate an efficient computation of MT-TKRP and TTM operations performed in a repetitive manner, as performed by CP-ALS and HOOI algorithms and others variants for computing CP and Tucker decompositions. To achieve this, we introduce a novel computational scheme based on the dimension tree data structure we introduced in Section 2.4, and significantly reduce the computational costs of these operations by reusing common partial results with the help of the tree structure. We discuss optimal dimension tree structures minimizing the computational cost in both these algorithms, and show that finding optimal dimension trees is NP-hard in both cases. Finally, we introduce an optimal greedy algorithm that finds the cheapest ordering for performing a sequence of TTM operations, which is found at the core of HOOI algorithm.

The organization of the chapter is as follows. In the next section, we show how dimension trees can be used for computing TTM and MTTKRP steps in HOOI and CP-ALS algorithms. In Section 7.1, we provide the dimension tree-based HOOI algorithm followed by all our theoretical findings regarding an optimal tree structure. Similarly, in Section 7.2, we give the CP-ALS algorithm using dimension trees followed by related theorems for optimality.

CHAPTER 7. ON COMPUTING DENSE TENSOR DECOMPOSITIONS USING DIMENSION TREES

Computing dense Tucker decomposition using dimension trees

In computing Tucker decomposition using Algorithm 3, performing successive TTMs to compute the resulting tensor Y constitutes the most expensive step in each ALS subiteration. This step involves the multiplication of X with the set {U (1) ,...,U (N ) }\U (n) of matrices to eventually update U (n) at the end of the subiteration for mode n. Here, we propose an efficient algorithm to perform this step faster with the following observation, which is similar to that discussed in Chapter 3. For a 4-dimensional tensor X 2 R I 1 ⇥I 2 ⇥I 3 ⇥I 4 , in the first two ALS subiterations, one needs to compute 4) and subsequently update the corresponding matrices U (1) and U (2) , respectively. Note that in doing so, U (3) and U (4) remains unchanged, which brings about the possibility of computing the intermediate tensor Z 4) , then reusing this partial result in the first and second subiterations as Z ⇥ 2 U (2) and Z ⇥ 1 U (1) , respectively. This computation stays valid as one can perform TTM in a set of distinct modes in any order. This way, however, the cost of TTMs reduces significantly.

X ⇥ 2 U (2) ⇥ 3 U (3) ⇥ 4 U (4) and X ⇥ 1 U (1) ⇥ 3 U (3) ⇥ 4 U (
= X ⇥ 3 U (3) ⇥ 4 U (
For an N -dimensional tensor, we are interested in identifying and reusing such common partial results as much as possible in a systematical way. This can be achieved by using a dimension tree T as follows. With each node t 2 T , we associate a tensor T (t) corresponding to the multiplication

X ⇥ m 1 U (m 1 ) ⇥ m 2 ••• ⇥ m |µ 0 t| U (m |µ 0 t| )
where µ 0 (t)={m 1 ,...,|µ 0 t|}. Note that this implies T (ROOT(T )) = X , and T (L n )=X ⇥ -n U (n) for all n 2 N N . With this tree structure, the tensor of any non-root tree node t can be computed from that of its parent as T (t)=T (P (t))

⇥ d 1 U (d 1 ) ⇥ d 2 •••⇥ d k U (d k )
where µ(P (t)) \ µ(t)={d 1 ,...,d k }. Indeed, if T (P (t)) does not exist, it needs to be similarly computed from its parent's tensor. We provide the algorithm for computing the tensor of any tree node t in Algorithm 16.

T (t) T (t) ⇥ d U (d) 6: return T (t)
One indeterminacy we have in this algorithm is the order of TTMs performed at Line 5. The order of these multiplications can have a significant impact on the overall computational cost. For minimizing this cost, we use an optimal greedy ordering method that minimizes the cost of a series of TTMs in distinct modes, which is demonstrated in the following theorem.

Theorem 4 (Optimal TTM order). Let X 2 R I 1 ⇥•••⇥I N be a tensor and D = {d 1 ,...,d |D| },D ✓ N N represent the set of dimensions in which X is to be multiplied with a set {U 1 ,...,U |D| } of matrices of corresponding sizes, i.e., U i 2 R

I d i ⇥R i . Let β i = R i /(1 -R i /I d i ), and suppose w.l.g that β 1  β 2 • • •β |D| . Then, the total cost of TTM is minimized with the multiplication order X ⇥ d 1 U 1 ⇥ d 2 U 2 ⇥ d 3 •••⇥ d |D| U |D| . Proof. Let ↵ n = R n /I n for n 2 N N and T = Q n2N N I n be the number of elements in X . Then, the cost of the multiplication X ⇥ 1 U (1) ⇥ 2 •••⇥ N U (N ) becomes T (R 1 + ↵ 1 R 2 + ↵ 1 ↵ 2 R 3 + •••+ TREES

Algorithm 17 DTREE-HOOI: Dimension tree-based HOOI algorithm

Input: X : An N -dimensional tensor R 1 ,...,R N : The rank of the decomposition for each mode T : A dimension tree for N dimensions Output: Tucker decomposition [[ G; U (1) ,...,U (N ) ]]

1: Initialize the matrix U (n) 2 R In⇥Rn for n =1,...,N 2: repeat 3:

for n =1,...,N do

4:

for all t 2 T do 5:

if n 2 µ 0 (t) then 6:

DESTROY(T (t))

I Destroy all tensors that are multiplied by U (n) .

7:

Y DTREE-TTM(L n ) I Perform TTM for the leaf node.

8:

U (n) R n leading left singular vectors of Y (n) 9: G Y ⇥ N U (N )
T I Form the core tensor. 10: until convergence or the maximum number of iterations 11: return [[ G; U (1) ,...,U (N ) ]] k children requires k -1 TTMs for each dimension in µ(t), which is minimized for k =2using a binary tree. Nevertheless, in some extreme cases, using a binary tree may not be optimal in terms of actual computational cost, despite avoiding extra TTMs, and we provide such a counter-example using a 6-dimensional tensor in Section 7.1.1. In most practical scenarios, however, binary dimension tree (BDT) is expected to provide optimal or close to optimal results; hence, we conduct the analysis using BDT in the rest of the discussion in this section.

We now investigate the computational complexity of finding an optimal BDT that minimizes the TTM cost in DTREE-HOOI. We will show that finding an optimal "balanced" binary dimension tree (BBDT) in which the dimensions µ(t) of each tree node t are equitably partitioned to its children t 1 and t 2 , i.e., |µ(t 1 )|, |µ(t 2 )|≥b | µ(t)|/2c, is NP-hard, and conjecture that the problem using arbitrary dimension trees still remains NP-hard. In doing so, we consider the PRODUCT-PARTITION problem, which is NP-hard [START_REF] Ng | Product partition and related problems of scheduling and systems reliability: Computational complexity and approximation[END_REF]: Given a set S = {s 1 ,...,s K } of positive integers, find a subset S 0 of S such that max(

Q i2S 0 i, Q i2S\S 0 i) is minimized.
We use a variant of this problem called BALANCED-PRODUCT-PARTITION where |S| =2 K is a multiset for K>0, and |S 0 | = K is another multiset containing exactly a half of the elements of S. This problem still remains NP-hard as one can trivially solve any instance of PRODUCT-PARTITION using an instance of BALANCED-PRODUCT-PARTITION with the following polynomial-time reduction. For the given set S of size K for PRODUCT-PARTITION, we create a multiset Q of size 2K having all elements of S, and K 1s. Let Q 0 be the optimal solution for BALANCED-PRODUCT-PARTITION of Q, and S 0 = Q 0 \ S be the set of elements of the optimal solution Q 0 belonging to S. Since the multiset Q 0 \ S can only contain 1s, we get

Q i2S 0 i = Q i2Q 0 i as well as Q i2S\S 0 i = Q i2Q\Q 0 . Therefore, max( Q i2Q 0 i, Q i2Q\Q 0 i) is minimized if and only if max( Q i2S 0 i, Q i2S\S 0 i
) is minimized, which makes S 0 an optimal solution for PRODUCT-PARTITION of S. Note that the objective of BALANCED-PRODUCT-PARTITION is also equivalent to minimizing the sum Q i2Q 0 i + Q i2Q\Q 0 i, and the same for PRODUCT-PARTITION.

Theorem 5 (Optimal BBDT for Tucker decomposition). Finding a BBDT that minimizes the cost of TTMs in computing Tucker decomposition is NP-hard.

Proof. We perform a reduction from BALANCED-PRODUCT-PARTITION using a set S = {s 1 ,...,s 2N }.

The intuition of the proof is as follows. We aim to find an optimal BBDT for computing the Tucker decomposition of a 2N +2 -dimensional tensor X 2 R I 1 ⇥•••⇥I 2N +2 using ranks of approximation R 1 ,...,R 2N +2 . The reduction associates the first 2N dimensions of the tensor with the corresponding elements of S, and uses two more auxiliary dimensions. Forming a BBDT for this tensor yields N +1 dimensions in both sets µ(t 1 ) and µ(t 2 ) of the children t 1 and t 2 of the root. The goal is to have exactly N dimensions corresponding to elements of S in both µ(t 1 ) and µ(t 2 ), and to show that these two subsets of S provide an optimal solution for BALANCED-PRODUCT-PARTITION for an optimal BBDT. We achieve this in two steps. First, we analyze the case where both µ(t 1 ) and µ(t 2 ) have N dimensions associated with S, and provide lower and upper bounds for the TTM cost of a such BBDT. We show in this case that the cost of the BBDT is minimized when the multiplications of the elements of two associated subsets of S (with µ(t 1 ) and µ(t 2 )) are balanced (i.e., the sum of two multiplications is minimized), which is effectively the objective of BALANCED-PRODUCT-PARTITION.

Next, we argue that if a child of the root has N -1 (or equivalently, N +1 ) dimensions associated with S, then the BBDT cannot be optimal, exceeding the provided upper bound in the previous case. As both children of the root must have N +1dimensions in total, there would be no other partitioning possibilities regarding these 2N dimensions associated with S, and the two subsets of S of size N associated with µ(t 1 ) and µ(t 2 ) in the former case would provide an optimal solution for BALANCED-PRODUCT-PARTITION of S.

We set R i = s i and

I i =3 Q s2S s for i 2 N 2N . We use R 2N +1 = R 2N +2 = K 1 (I 1 ) N = K 1 (3 Q s2S s i ) N and I 2N +1 = I 2N +2 = K 2 R 2N +1
where the coefficients K 1 ≥ 1 and K 2 ≥ 1 are left to be determined appropriately in the course of the proof. This yields

↵ i = s i /(3 Q s2S s)  1/3 for i 2 N 2N which also implies β i = R i /(1 -↵ i )  3s i /2. For dimensions 2N +1and 2N +2, we get ↵ 2N +1 = ↵ 2N +2 = K -1
2 , and

β 2N +1 = β 2N +2 = R 2N +1 /(1 -↵ 2N +1 ) >R 2N +1 = K 1 (3 Y s2S s) N > 3s max /2 ≥ β i ,i 2 N 2N ,
where s max = max s2S s. Hence, these two dimensions are to be multiplied after dimensions 1 ...2N in an optimal solution according to Theorem 4.

We start by investigating the cost of a BBDT where dimensions 2N +1and 2N +2reside in µ(t 1 ) and µ(t 2 ), respectively. Without loss of generality, suppose that µ(t 1 )={1,...,N,2N +1} and µ(t 2 )={N +1,...,2N, 2N +2}, and that

β 1 •••β N and β N +1 •••β 2N (
which can otherwise be obtained by a proper permutation of dimensions). Therefore, we perform TTMs with X in the sequence of dimensions 1,...,N,2N +1and N +1,...,2N, 2N +2for nodes t 2 and t 1 , respectively, as this corresponds to the optimal dimension ordering stated in Theorem 4. Let C represent the TTM cost of an optimal BBDT with the given µ(t 1 ) and µ(t 2 ), and Z = Q i2N 2N +2 I i be the number of elements in X . Then, we can express the TTM cost as

C = Z(s 1 + ↵ 1 s 2 + •••+ ↵ 1 ...↵ N -1 s N + ↵ 1 ...↵ N K 1 (3 Y s2S s) N + s N +1 + ↵ N +1 s N +2 + •••+ ↵ N +1 ...↵ 2N -1 s 2N + ↵ N +1 ...↵ 2N K 1 (3 Y s2S s) N )+ C(1,...,N,2N + 1) + C(N +1,...,2N, 2N + 2)), (7.1)
where the first two summands correspond to the TTM cost due to nodes t 2 and t 1 , whereas C(1,...,N,2N + 1) and C(N +1,...,2N, 2N + 2) denote the total TTM cost of all remaining nodes in the subtrees rooted at t 1 and t 2 , respectively, in an optimal BBDT. We rewrite (7.1) as follows

C = ZK 1 [K -1 1 (s 1 + ↵ 1 s 2 + •••+ ↵ 1 ...↵ N -1 s N + s N +1 + ↵ N +1 s N +2 + •••+ ↵ N +1 ...↵ 2N -1 s 2N )+ s 1 ...s N + s N +1 ...s 2N ]+ C(1,...,N,2N + 1) + C(N +1,...,2N, 2N + 2), (7.2)
and obtain a trivial lower bound

C > ZK 1 (s 1 ...s N + s N +1 ...s 2N ), (7.3)
as the rest of the summands are positive. Next, we aim find an upper bound for C by bounding the costs C(1,...,N,2N +1) and C(N +1,...,2N, 2N +2) as following. C(1,...,N,2N +1) corresponds to the TTM cost using a BBDT having N +1leaves (and at most 2N nodes in total excluding its root t 1 ). The number of elements in the tensor at the root of this sub-tree is Z 1 = Z↵ N +1 ...↵ 2N K -1 2 as the tensor is obtained from the multiplication of X in all modes in µ(t 2 ). Each node can require at most one TTM per each dimension in µ(t 1 ), and the cost of the TTM in a dimension d cannot exceed Z 1 R d since the tensor cannot grow larger after multiplications (R d  I d ). Therefore, we obtain the following upper bound on the total TTM cost C(1,...,N,2N + 1):

C(1,...,N,2N + 1)  2N Z 1 ( N X i=1 R i + R 2N +1 ) < Z↵ N +1 ...↵ 2N K -1 2 2N ( 2N +2 X i=1 R i ). Finally, setting K 2 =2N ( P 2N +2 i=1 R i ) yields C(1,...,N,2N + 1) < Z↵ N +1 ...↵ 2N .
We similarly obtain the upper bound Z↵ 1 ...↵ N for C(N +1,...,2N, 2N + 2). Using these upper bounds in (7.2) gives (7.4)

C < ZK 1 K -1 1 (s 1 + ↵ 1 s 2 + •••+ ↵ 1 ...↵ N -1 s N + s N +1 + ↵ N +1 s N +2 + •••+ ↵ N +1 ...↵ 2N -1 s 2N + ↵ N +1 ...↵ 2N + ↵ 1 ...↵ N )+ s 1 ...s N + s N +1 ...s 2N .
Next, we analyze the case where dimensions 2N +1 and 2N +2 reside in the same child of the root, namely t 1 . Suppose w.l.g that µ(t 1 )={1,...,N -1, 2N +1, 2N +2} and µ(t 2 )={N, . . . , 2N }. Then, the TTM cost for t 2 becomes

C t 2 = Z(s 1 + ↵ 1 s 2 + •••+ ↵ 1 ...↵ N -2 s N -1 + ↵ 1 ...↵ N -1 K 1 (3 Y s2S s) N + ↵ 1 ...↵ N -1 K -1 2 K 1 (3 Y s2S s) N ) >Z↵ 1 ...↵ N -1 K 1 (3 Y s2S s) N = ZK 1 (s 1 ...s N -1 3 Y s2S s) ≥ZK 1 (s 1 ...s N + s N +1 ...s 2N + 1),
which already exceeds the upper bound in (7.4). Therefore, we conclude that dimensions 2N +1 and 2N +2cannot reside in the same child of the root in an optimal solution. As a result, we obtain exactly N dimensions associated with S in µ(t 1 ) and µ(t 2 ) in an optimal solution.

Using this result and the bounds in (7.3) and (7.4) we can finally perform a reduction from the decision version of BALANCED-PRODUCT-PARTITION: Given a set S of 2N positive integers, is there a S 0 ⇢ S, |S 0 | = N such that Q s2S 0 s + Q s2S\S 0 s Cfor some C≥1? We claim that such S 0 exists if and only if there exists a BBDT constructed in the aforementioned manner whose cost is smaller than ZK 1 (C + 1). If a such S 0 exists, then (7.4) suggests that one can construct a corresponding BBDT whose cost is less than ZK 1 (C + 1). On the contrary, if there exists no such S 0 , then any S 0 yields Q s2S 0 s + Q s2S\S 0 s ≥C+1, in which case the cost of the associated BBDT exceeds ZK 1 (C + 1) due to (7.3), which concludes the proof.

Counter-example: Tensor having no optimal BDT

Here, we provide a counter-example using a 6-dimensional tensor X 2 R I 1 ⇥•••⇥I 6 to show that using a BDT is not necessarily optimal to compute Tucker decomposition. The first three and the last three dimensions of X have identical sizes and ranks of approximation. Specifically, we let I 1 = I 2 = I 3 = k and R 1 = R 2 = R 3 = R<kin the first three dimensions, and I 4 = I 5 = I 6 = R 4 = R 5 = R 6 = c in the last three dimensions. We call the first three and the last tree dimensions type-1 and type-2 dimensions, respectively. Note that ↵ = ↵ 1 = ↵ 2 = ↵ 3 <↵ 4 = ↵ 5 = ↵ 6 =1 , and similarly β 1 = β 2 = β 3 <β 4 = β 5 = β 6 = 1, therefore, type-1 dimensions are to be multiplied before type-2 dimensions according to Theorem 4. In Figure 7.1, we provide a ternary dimension tree with a total cost of (3 + 3↵ +3↵ 2 )k +9↵ 2 c. We can choose c arbitrarily large so that the term 9↵ 2 c dominates the cost, and ↵ small enough so that ↵ 0 c 9↵ 2 c and ↵ 1 c 9↵ 2 c. In this case, any BDT whose cost involves a term of order ↵ 0 c or ↵ 1 c cannot be optimal, having a greater cost than the provided ternary tree.

Figure 7.1: A ternary dimension tree for X . The µ set of the node is provided on the left, and the tensor size coefficient is provided on the right of each tree node. The TTM cost of each level is provided on the right.

{1,2,3,4,5,6} {1,4} {1} 1 αα α α α 2 α 2 α 2 α α 3(k+αk+2α 2 c) 3(α 2 k+α 2 c) Cost {2,5} {3,6} {4} {2} {5} {3} {6}
We now show that any BDT using X either has a cost with a term of order ↵ 0 c or ↵ 1 c, or with a term 9↵ 2 c + d↵ 3 c with d ≥ 1; hence, cannot be optimal with a sufficiently large c and a sufficiently small ↵. We do this by exhaustively considering all possible BDTs and analyzing their costs, while aggressively pruning tree configurations that cannot provide optimality. Luckily, most non-optimal configurations can easily be pruned due to symmetry (as we only have two types of dimensions), leaving us with only a handful of instances to consider. We begin by partitioning type-2 dimensions to the children t 1 and t 2 of the root. There are only two possibilities: 4 2 µ(t 1 ) and 5, 6 2 µ(t 2 ), or 4, 5, 6 2 µ(t 2 ). Since we have only three type-1 dimensions, in the former case either µ(t 1 ) or µ(t 2 ) will have one or zero type-1 dimension, while the other set having two or three of them. In this case, the TTM cost of the other vertex involves a term ↵ 0 c or ↵ 1 c, which already renders this configuration non-optimal. Therefore, we only consider the partition 4, 5, 6 2 µ(t 2 ), which is the only one that can possibly provide an optimal solution. In this case, there are three possible configurations after partitioning type-1 dimensions: µ(t 1 )={1} and µ(t 2 )={2, 3, 4, 5, 6}, µ(t 1 )={1, 2} and µ(t 2 )={3, 4, 5, 6}, or µ(t 1 )={1, 2, 3} and µ(t 2 )={4, 5, 6}. Note that in the second and third configurations, computing TTM for t 1 involves a term ↵ 1 c and ↵ 0 c, respectively, which prevents optimality. Hence, in the rest of the discussion we only consider the first configuration, in which t 1 incurs the TTM cost 3↵ 2 c. We focus on the cost of the sub-tree rooted at t 2 , and count only the cost due to terms with the coefficient c. Note that if the remaining type-1 dimensions, namely 2 and 3, reside in the same child of t 2 , the other child of t 2 incurs a cost of at least ↵ 1 c; therefore, 2 and 3 must reside in different children of t 2 . In Figure 7.2 we provide six such possibilities, all of which incur a cost of 9↵ 2 c + c↵ 3 c with c ≥ 1. Therefore, we conclude that a BDT cannot be optimal for the given X for sufficiently small ↵ and large c.

Computing dense CP decomposition using dimension trees

In computing Tucker decomposition, we have shown how dimension trees reduce the computational cost of the TTM step by storing and reusing common partial TTM results used across different subiterations of HOOI. In (2.4), we have a TTV formulation of the expensive MTTKRP step in CP-ALS, which similarly enables precomputing and reusing partial TTV results, e.g., for a 4-dimensional tensor X , one can similarly compute Z = X ⇥ 3 u

(3) r ⇥ 4 u (4)
r , then use this partial result to obtain both M (1)

(:,r)=Z ⇥ 2 u (2)
r and M (2) (:,r)=Z ⇥ 1 u

(1) r . One difference in this case is that a series of R Figure 7.2: All possible configurations and associated TTM costs at each level of BDTs rooted at t 2 . µ(t) is given on the left of each tree node t, On the right of a tree node, the size coefficient of its tensor is provided, i.e., t 2 has a tensor of size (I 1 I 2 I 3 I 4 I 5 I 6 )↵. TTVs is needed to compute all columns of M (1) and M (2) . For this reason, each non-root tree node t holds tensors T r (t) for each r 2 N R corresponding to the partial result for the rth TTV. The root node, however, possesses only the original tensor X .

Algorithm 18 DTREE-TTV: Performing TTV on a dimension tree for a dense tensor Input: t: A dimension tree node Output: T : (t) are computed. We provide the method for performing TTVs using a dimension tree in Algorithm 18 for a tree node t. The algorithm similarly returns the tensors T (t) if they are already computed. Otherwise, DTREE-TTV is called at Line 4 to obtain the parent P (t)'s tensors. If P (t) is the root, then we perform traditional MTTKRP by properly matricizing T (P (t)) (which is X ), and multiplying it with the Khatri-Rao product of all matrices corresponding to all dimensions except those in µ(t). This yields a matrix of size Q d2µ(t) I d ⇥ R, whose rth column corresponds to T r (t) in vectorized form. This step costs Q n2µ(P (t)) I n = Q n2N N I n operations for matricizing the tensor, R Q n2µ(t) I n operations for forming the Khatri-Rao product, and finally R Q n2µ(P (t)) I n = R Q n2N N I n operations for multiplying the matricized tensor with the Khatri-Rao product. If P (t) is not the root, then each T r (t) is computed using the corresponding parent tensor T r (P (t)). Here, we matricize the tensor, compute the Kronecker product of vectors to be multiplied, and finally execute the matrix-vector multiplication. For computing each tensor T r (t), this incurs Q n2µ(P (t)) I n operations for matricizing the tensor T r (P (t)), Q n2µ(t) I n operations for performing the Kronecker product of vectors, and finally Q n2µ(P (t)) I n operations for multiplying the matricized tensor with the Khatri-Rao product, resulting in the overall cost R(2 Q n2µ(P (t)) I n + Q n2µ(t) I n ) for computing all R tensors, which we denote as T : (t). This technique is more efficient than multiplying the tensor with vectors one by one, which requires a matricization and a multiplication of the tensor for each vector. This is because TTV is a memory-bound operation, and this way we access the tensor only twice (once for matricization, once for multiplication), and the cost of Kronecker product is negligible compared with the cost of multiplication.

The dimension tree-based Algorithm 5 provided in Chapter 3 equally applies to both sparse and dense tensors. In the dense case, however, Algorithm 18 is to be employed to carry out TTVs at Line 9.

Indeed, the structure of the dimension tree plays a crucial role in the computational cost of MT-TKRPs, and we are interested in finding the optimal tree structure minimizing this cost. In the following part, we provide two theorems pertaining to the computation of an optimal dimension tree for CP-ALS. Theorem 6. For any N -dimensional tensor X 2 R I 1 ⇥•••⇥I N , except X 2 R 2⇥2⇥2 , an optimal dimen-sion tree that minimizes the TTV cost in the execution of CP-ALS must be binary.

Proof. The proof is by construction of a BDT in the following manner. We take any dimension tree T which is not binary, focus on its any subtree rooted at a node t having K>2 children, and finally replace this subtree with another having K -1 children and a smaller cost. This ensures that a nonbinary dimension tree cannot be optimal, as any such tree can be transformed into a BDT having less cost using a sequence of such transformations.

Let t 1 ,...,t K be the children of the root t, and ⇡(m)= Q n2µ(m) I n denote the size of the tensor of a tree node m. Let also C(t) be the TTV cost of an optimal dimension tree rooted at t, excluding the cost of the root. We first consider the case where t is not the root of T . We can express the cost C k of the original subtree with K children as

C k = R 2K⇡(t)+ K X k=1 ⇡(t) ⇡(t k ) + K X k=1 C(t k ), (7.5)
where, for each dimension k 2 N K , R ⇡(t) ⇡(t k ) is the Kronecker product cost, and 2R⇡(t) corresponds to the total cost of matricization and the multiplication of all R tensors with this product.

We now consider a variant of this tree where the nodes t 1 and t 2 are joint using an auxiliary node t 12 , which is made a child of t with µ(t 12 )=µ(t 1 ) [ µ(t 2 ) (thus ⇡(t 12 )=⇡(t 1 )⇡(t 2 )). The rest of the children t 3 ,...,t K of t are kept the same to obtain a subtree whose root t has K -1 children and the following cost: ). Hence, the modified subtree whose root has K -1 children always provides a smaller cost in this case. TREES Next, we consider the case where t is the root of T . In this case, the matricization costs only ⇡(t) for the root's tensor, hence the cost of the original tree having K children becomes

C k-1 =R 2(K -1)⇡(t)+ K X i=3 ⇡(
C k = K⇡(t)+R K⇡(t)+ K X k=1 ⇡(t) ⇡(t k ) + K X k=1 C(t k ). (7.8)
This time, we assume w.l.g that t 1 and t 2 have the two minimum tensor sizes among all children, i.e., ⇡(t 1 ),⇡(t 2 )  ⇡(t k ) for all k =3 ,...,K, which also implies that ⇡(t 1 ),⇡(t 2 )  L. For the tree having K -1 children, we obtain the following cost: where the first summand is the cost of the matricization and the multiplication of X with the Khatri-Rao products of corresponding matrices for t 1 and t 2 , the second and the third summands correspond to the cost of forming the Khatri-Rao product for t 2 and t 1 , respectively. Here, C(1,...,K,N+1,N+ 2) and C(K +1,...,N,N+3,N+ 4) denote the total MTTKRP cost of subtrees rooted at t 1 and t 2 , respectively, excluding the cost of t 1 and t 2 . Both t 1 and t 2 have two children as |µ(t 1 )|, |µ(t 2 )|≥2.

C k-1 =(K -
Each child of t 1 incurs a cost 2Rk 2 Q K n=1 s n for matricizing the tensor of t 1 , then multiplying it with the Kronecker products. Similarly, each of two children of t 2 has a cost of at least 2Rk 2 Q N n=K+1 s n . With these costs at hand (which exclude the cost of forming the Kronecker products for the children of t 1 and t 2 , and further costs down the tree), we obtain the following lower bound from (7.11):

C > (2R + 2)k 4 Y n2N N s n +5Rk 2 ( K Y n=1 s n + N Y n=K+1 s n ). (7.12)
Next, we aim to find an upper bound for an optimal C(1,...,K,N +1,N + 2) and C(K + 1,...,N,N +3,N + 4), which in turn would yield an upper bound for C. Let t 11 and t 12 be the children of t 1 . We consider the cost for the case where N +1 2 µ(t 11 ) and N +2 2 µ(t 12 ). In this case, the size of the tensor of t 1 is Rk 2 Q K n=1 ; therefore, cost of both t 11 and t 12 are upper bounded by 2Rk 2 Q K n=1 +Rk Q K n=1 as the cost of forming the Kronecker products cannot exceed Rk Q K n=1 with the given partition of type-2 dimensions. In addition, note that t 1 is the root of a subtree having K +2leaf nodes, and having up to K +1non-leaf nodes each of which has two children; hence, it cannot have more than 2K nodes excluding t 1 , t 11 , and t 12 . These nodes are descendants of either t 11 or t 12 ; therefore, each of these nodes incurs a cost which cannot exceed 2Rk

Q K i=1 s i + Rk Q K i=1 s i =3 Rk Q K i=1 s i
, where the first and the second summands are upper bounds for the cost of the matricization of the tensor and its multiplication with the Kronecker product, and the cost of forming the Kronecker product, respectively. This results in the upper bound C(1,...,k,N +1,N+ 2) < 4Rk 2 Q K n=1 s n +(6K + 2)Rk Q K n=1 s n . We similarly obtain the upper bound 4Rk in both cases, which exceeds the upper bound provided in (7.13). Therefore, we conclude that these two partitionings cannot provide an optimal BDT.

We can now perform the reduction from the decision version of PRODUCT-PARTITION, knowing that an optimal BDT assigns exactly two dimensions of type-2 to each children t 1 and t 2 of the root t, and has the lower and the upper bounds provided in (7.12) and (7.13) with respect to the partitioning of type-1 dimensions to µ(t 1 ) and µ(t 2 ). We claim that a S 0 ✓ S with Q s2S 0 s + Q s2S\S 0  C exists for some C ≥ 1 if and only if there is a BDT for X constructed in the aforementioned manner whose MTTKRP cost is smaller than (2R + 2)k 4 Q n2N N s n +5Rk 2 (C + 1) for any positive integer R. If there exists a such S 0 , then (7.13) suggests that we can construct a BDT whose cost is inferior to (2R + 2)k 4 Q n2N N s n +5Rk 2 (C + 1). If Q s2S 0 s + Q s2S\S 0 ≥ C +1of all subsets S 0 ✓ S, then (7.12) implies that the cost of all BDTs exceed (2R + 2)k 4 Q n2N N s n +5Rk 2 (C + 1), which concludes the proof.

Conclusion

In this chapter, we introduced efficient algorithms to compute Tucker and CP decomposition of dense tensors. The algorithms utilize an novel computational scheme based on a dimension tree structure, and thereby enable re-using common partial TTM and MTTKRP results across different subiterations of the tensor decomposition algorithms. We also provided the first complexity results for this problem regarding an optimal dimension tree structure minimizing the associated computational costs. In particular, we proved that finding an optimal binary tree to minimize the cost of TTM and MTTKRP operations respectively in DTREE-HOOI and DTREE-CP-ALS algorithms is NP-hard. We further showed that the optimal tree must be binary for TTV, except in one degenerate instance. On the contrary, we provided a counterexample using an optimal ternary tree for TTM. All these results lay the foundation for a complete analysis of the use of dimension trees for tensor decomposition. Our further work will involve finding fast exact algorithms as well as effective heuristics for finding an optimal/good dimension tree, and a thorough experimentation with the implementations of these two algorithms using dimension trees.

Chapter 8

Sparse Tensor Software

An overview of recent sparse tensor software

At the beginning of this thesis, the only available parallel tensor factorization software was DFACTO; A C++ implementation for computing the CP decomposition of 3-dimensional tensors using MPI parallelism [START_REF] Choi | DFacTo: Distributed factorization of tensors[END_REF]. Even though this library significantly accelerated computing CP decomposition compared with the previous Hadoop implementations, it a major drawback preventing it from being scalable for analyzing massive tensor datasets using thousands of processors. The algorithm was not fully distributed in the sense that the factor matrices were entirely replicated in all MPI processes. This not only prevented scalability due to memory requirements, but also incurred a significant amount of redundant communication to transfer factor matrices to all processes after each update. In the meantime, SPLATT (C, OpenMP) [START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF] software came up with an efficient implementation and parallelization in shared memory, outperforming DFACTO by a large margin.

HYPERTENSOR [START_REF] Kaya A N D B | Scalable sparse tensor decompositions in distributed memory systems[END_REF] first came up in 2015 as a distributed-memory parallel implementation (MPI, C++) for computing the CP decomposition of sparse tensors. Its main advantage over DFACTO was to be able to distribute the tensor as well as factor matrices, and to employ a point-to-point messaging scheme which communicates only the required rows of factor matrices. This helped overcome memory limitations as well as reduce the communication cost dramatically. On top of this, employing effective hypergraph partitioning routines further enhanced scalability of HYPERTENSOR by incurring even less communication. A similar MPI parallelism approach using point-to-point messaging was also employed by SPLATT later on in its distributed memory parallelization [START_REF] Smith | A medium-grained algorithm for sparse tensor factorization[END_REF]. Next, HY-PERTENSOR added a parallel algorithm for computing the Tucker decomposition of sparse tensors with distributed as well as shared memory parallelism, and it is still the only high performance parallel software computing the Tucker decomposition of sparse tensors [START_REF]High performance parallel algorithms for the Tucker decomposition of sparse tensors[END_REF]. Finally, a dimension treebased computational scheme for computing CP decomposition and an associated sparse tensor storage scheme with an effective shared memory parallelization was introduced to HYPERTENSOR [START_REF]Parallel CP decomposition of sparse tensors using dimension trees[END_REF]. This rendered HYPERTENSOR the fastest tensor factorization software both in shared and distributed memory, even though a very recent extension, called ADATM [START_REF] Li | Model-driven sparse CP decomposition for higher-order tensors[END_REF], also introduced the idea of storing and reusing partial results to SPLATT to reduce the computational cost.

In the following sections, we provide an overview of HYPERTENSOR and PACOS, which drives the partitioning and communication of HYPERTENSOR, alongside with some code samples.

HYPERTENSOR: A HIGH PERFORMANCE PARALLEL SPARSE TENSOR FACTORIZATION LIBRARY 109

8.2 HYPERTENSOR: A high performance parallel sparse tensor factorization library HYPERTENSOR is the parallel sparse tensor library we developed during this thesis as the software outcome of the research work. It is a C++ implementation using OpenMP and MPI for shared and distributed memory parallelism. The communication and partitioning routines driving the distributed memory parallelization of HYPERTENSOR are later on extracted in a generic software package called PACOS in order to simplify the code as well as to enable developing likewise parallel applications.

HYPERTENSOR is fully templated to support single and double floating point numbers as well as 32-and 64-bit integers for tensor indexes. It can be compiled and run with and without MPI and OpenMP support. It does, however, require a compatible C++11 compiler.

The two principal routines provided by HYPERTENSOR are parallel CP-ALS and HOOI algorithms for sparse tensors. HOOI involves a parallel truncated SVD step, which is realized with the help of PETSC and SLEPC software packages. HYPERTENSOR employs an object-driven paradigm where all tensors and matrices are represented as objects with associated functions. Following are the main data structures used by HYPERTENSOR:

• SpTensor: Sparse tensor data structure, in coordinate or dimension tree format. It has member functions for computing CP and Tucker decompositions.

• KTensor: Kruskal tensor, obtained from a CP decomposition algorithm. It involves N dense factor matrices as well as a vector representing the column norms.

• TTensor: Tucker tensor, obtained from a Tucker decomposition algorithm. It consists of an N -dimensional dense core tensor and N dense factor matrices.

• DMatrix: Dense matrix data structure providing a wrapper for many BLAS, LAPACK and other matrix operations.

Tensor data structures are also accompanied by some utilities such as norm calculation, print functions, I/O routines for loading from and saving to disk, and many others. DMatrix provides a lightweight wrapper on top of dense matrix routines, with the principle goal of providing templated matrix routines to be employed by other templated functions.

In Figure 8.1, we provide a sample code snippet computing the CP and Tucker decomposition of a sparse tensor read from a file. Each MPI rank reads a part of a partitioned tensor along with a vector partition table for all dimensions. Both the partitioned tensor and the partition table are obtained in a preceding step involving the partitioner of PACOS. Next, the sparse tensor is loaded from the file, and the parameters of the tensor decomposition algorithms are set. Finally, executing the member CP-ALS and HOOI algorithms output Kruskal and Tucker tensors representing the CP and Tucker decompositions of the tensor with the provided rank of approximation.

Chapter 9 Conclusion

In this thesis, we have investigated efficient algorithms to compute tensor decompositions for both sparse and dense tensors. We have achieved this efficiency through shared and distributed memory parallelism, effective partitioning and load balancing methods for a scalable parallel execution, memory efficient sparse tensor representations, and a new computational scheme asymptotically reducing the computational cost of tensor decomposition algorithms.

We employed most of these techniques for computing the CP decomposition of high dimensional sparse tensors in Chapter 3, which serves as the "porte-parole" for the subsequent chapters with the diversity of techniques involved in the chapter. First, in this chapter, we introduce a novel computational scheme using a dimension tree structure to reduce the number of TTV operations within a CP-ALS iteration by a factor of O(N/ log N ) for an N -dimension tensor. We show that this approach provides significant computational gains even for 3-and 4-dimensional tensors, and these gains only improve as the tensor dimensionality grows. To effectively carry out sparse tensor operations, we provided a sparse tensor data structure associated with a dimension tree, which similarly reduces the number of index arrays used to store the tensor to O(N log N ) instead of O(N 2 ), which is found in some state-of-the-art implementations. We give a fast shared memory parallel algorithm using this data structure and computational scheme. Finally, we perform an effective parallelization of this algorithm in distributed memory, which also benefits from novel hypergraph partitioning techniques for establishing load balance while minimizing communication requirements. The work in this chapter is pioneer in establishing the bridge between sparse linear and multilinear algebra by adopting and generalizing fundamental parallelization and partitioning techniques in the existing sparse matrix literature to parallel sparse tensor factorization, and thereby sets the bar for the state-of-the-art in the field.

In Chapter 4, we present a shared memory parallel algorithm for computing a nonzero CP decomposition. This introduces a constraint on factor matrices to not involve any numeric zeros, which in turn enables a very fast computation of MTTKRP (or equivalently, TTV) operations for high dimensional sparse tensors. Specifically, the new method amounts to performing O(N ) TTVs in a CP-ALS iteration, whereas the existing methods require O(N 2 ) TTVs, and even our dimension treebased approach takes O(N log N ) TTVs per iteration. Also, this method performs significantly less precomputation to set up tensor data structures, which provides, for an N -dimensional tensor having k nonzero entries, O(log k) and O(log N log k) faster preprocessing in compare to SPLATT and our dimension tree-based method. We also discuss load balancing strategies and memory initialization

Chapitre 4 :

 4 Décomposition CP Parallèle Nonzéro des Tenseurs Creux. Pour traiter efficacement les grands tenseurs à grande dimension, nous présentons dans ce chapitre un algorithme parallèle en mémoire partagée pour calculer une décomposition de CP dans laquelle les matrices de facteurs sont supposées ne contenir aucun élément zéro. Cette contrainte supplémentaire permet un calcul très efficace de l'étape MTTKRP dans CP-ALS, et les gains de performance augmentent avec la dimensionnalité du tenseur. Pour un tenseur N -dimensionnel ayant k entrées nonzéro, notre méthode fournit O(log k) et O(log N log k) des temps de prétraitement plus rapides et effectue O(N ) et O(log N ) moins de travail dans le calcul d'une décomposition CP sur deux méthodes courants, ce dernier étant l'algorithme introduit au Chapitre 3. Avec ces contributions algorithmiques et une implémentation parallèle efficace, nous atteignons jusqu'à 16,7x d'accélération séquentielle et jusqu'à 10,5x d'accelération dans des exécutions parallèles sur ces bibliothèques sur une machine ayant à 28 coeurs. Ce faisant, nous avons un temps de prétraitement de jusqu'à 24x moins, et utilisons jusqu'à O(log N ) moins de mémoire pour stocker des résultats intermédiaires. Nous montrons en utilisant un tenseur réel que la précision de notre méthode est comparable à la décomposition CP standard.

Figure 2 .

 2 2 illustrates the Tucker decomposition of a 3-dimensional tensor.

  r . Hence, we can compute Y r = X ⇥ 3 u(3) r ⇥ 4 u (4)r , then reuse it in the first and the second subiterations as Y r ⇥ 2 u(2)

Figure 3 . 1 :

 31 Figure 3.1: BBDT of a 4-dimensional sparse tensor X 2 R 4⇥4⇥4⇥4 having 7 nonzeros. Each closed box refers to a tree node. Within each node, the index array and the mode set corresponding to that node are given. The reduction sets of two nodes in the tree are indicated with the dashed lines.

p

  to minimize the cost of communication. In Algorithm 7, we describe the fine-grain parallel algorithm that operates in this manner at process p. The elements X p , I (n) p , and F (n) p are determined in the partitioning phase, and are provided as input to the algorithm. Each process starts with the subset F (n)

Algorithm 8

 8 COMM-FACTOR-MATRIX: Communication routine for factor matricesInput: M: Distributed factor matrix to be communicated comm : The type of communication. "fold" sums up the partial results in owner processes, whereas "expand" communicates the final results from owners to all others. F p : The rows used by process p I p : The rows owned by process p σ: The ownership of each row of M Output: Rows of M are properly communicated.1: if comm = "fold" then 2:

Figure 3 . 2 :

 32 Figure 3.2: Fine-grain hypergraph model for the 3 ⇥ 3 ⇥ 3 tensor X = {(1, 2, 3), (2, 3, 1), (3, 1, 2)} and a 3-way partition of the hypergraph. The objective is to minimize the cutsize of the partition while maintaining a balance on the total part weights corresponding to each vertex type (shown with different colors).

( 3 )

 3 k to model this dependency. In Figure3.2, we demonstrate this fine-grain hypergraph model on a sample tensor X = {(1, 2, 3), (2, 3, 1), (3, 1, 2)}, yet we exclude the vertex weights for simplicity. Each vertex type and hyperedge type is shown using a different color in the figure.

Figure 3 . 3 :

 33 Figure 3.3: Medium-grain hypergraph models for partitioning the 3 ⇥ 3 ⇥ 3 tensor X = {(1, 2, 3), (2, 3, 1), (3, 1, 2)} in all three dimensions. Vertex weights are not shown. In (a), we have a single constraint partitioning, whereas in (b) and (c), we employ multi-constraint partitioning with P 1 and P 1 P 2 constraints, respectively.

CHAPTER 3 .

 3 PARALLEL CP DECOMPOSITION OF SPARSE TENSORS USING DIMENSION TREES

CHAPTER 3 .

 3 PARALLEL CP DECOMPOSITION OF SPARSE TENSORS USING DIMENSION TREES

4. 1 .

 1 PARALLEL CP DECOMPOSITION USING NONZERO FACTORS 49 memory for intermediate results for an execution using P threads, yet it incurs the highest computational cost.

4 :

 4 parallel for p =1...P do I Compute M(n)

P

  of row indices 1,...,I n . Each process p computes the set I (n)

  vector Hadamard operations at Lines 8, 9 and 17. Second, it performs the multiplication of matrices of size |I (n) p |⇥R and R ⇥ R at Line 11, and of two matrices of size |I (n) p |⇥R at Line 13. To balance the first cost pertaining to sparse tensor computations, one has to make sure that the associated cost P i2I (n) p |rl (n) (i)| is partitioned equitably to processes. Regarding the second cost for dense matrix operations, each process should have equal number of rows, i.e., |I

Figure 4 . 2 :

 42 Figure 4.2: Speedup over splatt on ten32D.

Figure 4 . 3 :

 43 Figure 4.3: Accuracy comparison of hypertensor and cp-eps on NELL2-2. We show the geometric mean of approximation values of 100 CP-ALS executions with random initial factor matrices using both methods.

5 :

 5 parfor i 2J n do I TTMc for mode n 6:

  (d) (j d , :) 10:

( 3 )

 3 k to perform the Kronecker product. This specifies the data to be exchanged in the communication step.Algorithm 13 gives the distributed memory parallel HOOI executed by the process p k . Initially, we assume a partition of task indices I (n) k

k

  , :), and the rows i n that are accessed in the local TTMc computations due to x i 1 ,...,in,...,i N 2 X k .After the local TTMc step, we obtain the row-wise distributed matrix Y (n) where p k owns the setI (n)k of rows. For the subsequent TRSVD step, we need to perform the MxV and MTxV multiplications y Y (n) x and x T y T Y (n) . We use a block partitioning for the vector x. On the other hand, CHAPTER 5. PARALLEL SPARSE TUCKER DECOMPOSITION Algorithm 13 Distributed memory parallel HOOI executed at process p k Input: X k : Partition of X owned by the process p k type: 'coarse-grain' or 'fine-grain' I(1)k ,...,I

9 :

 9 parfor all i 2 K n do I TTMc for mode n 10:

i

  , we add a vertex and a hyperedge to the hypergraph. For N =3 , we set the computational weights |X (i, :, :)|, |X (:,j,:)|, and |X (:, :,k)| to the tasks t

  Non-negative matrix factorization (NMF) is the problem of finding two low rank factors W 2 R m⇥k + and H 2 R k⇥n + for a given input matrix A 2 R m⇥n + such that A ⇡ WH. Here, R m⇥n + denotes the set of m ⇥ n matrices with non-negative real values. Formally, the NMF problem [90] can be defined as min W>0,H>0 kA -WHk F + φ(W)+ (H),

6. 1 .

 1 SPARSE NON-NEGATIVE MATRIX FACTORIZATION 75 which have unique solutions for a full rank H and W:

Algorithm 15

 15 DIST-SPNMF: Distributed memory Sparse NMF algorithm Input: A p : An m ⇥ n sparse matrix I p ,J p : Set of rows/columns of W and H owned by process p F p ,G p : Footprints of process p on W and H W(I p , :),H(:, J p ): Owned rows and columns of initial W and H k: The NMF rank Output: Process p gets the final value of W(I p , :) and H(:, J p ) 1: repeat 2:

Figure 6 . 1 :

 61 Figure 6.1: A 5x5 checkerboard partition of a sparse matrix.

  Flickr matrix with k = 16 on Rhea. Flickr matrix with k = 48 on Blue Gene/Q.

Figure 6 . 2 :

 62 Figure 6.2: Strong scaling on Flickr matrix.

Figure 6 . 3 :

 63 Figure 6.3: Strong scaling on Delicious matrix.

  Rhea p = 3072, k = 48.

Figure 6 . 4 :

 64 Figure 6.4: Flickr matrix runtime dissection for k = 16 and k = 48 on Rhea.

  Rhea p = 16, k = 16. Rhea p = 3072, k = 16. Rhea p = 3072, k = 48.

Figure 6 . 5 :

 65 Figure 6.5: Delicious matrix runtime dissection for k = 16 and k = 48 on Rhea.

Figure 6 . 6 :

 66 Figure 6.6: Runtime dissection of Flickr and Delicious for k = 48 on Blue Gene/Q. The left two bars are for the Flickr matrix, and the right two bars are for the Delicious matrix.

7. 1 .

 1 COMPUTING DENSE TUCKER DECOMPOSITION USING DIMENSION TREES 97 Since ↵ i  1/3 for i 2 N 2N , setting K 1 =(4s max + 2) yields the final upper bound C < ZK 1 (s 1 ...s N + s N +1 ...s 2N + 1).

100 CHAPTER 7 .

 1007 ON COMPUTING DENSE TENSOR DECOMPOSITIONS USING DIMENSION TREES

CHAPTER 7 .+

 7 ON COMPUTING DENSE TENSOR DECOMPOSITIONS USING DIMENSION TREES C(1,...,K,N +1,N + 2) + C(K +1,...,N,N +3,N + 4),(7.11)

  2 Q N n=K+1 s n +(6N -6K + 2)Rk Q N n=K+1s n for C(k +1,...,N,N +3,N + 4) for the case where N +32 µ(t 21 ) and N +42 µ(t 22 ). Finally, setting k = l(6N + 4) Q N i=1 s i for any l ≥ 1, and replacing these two upper bounds in (7.11) yieldsC <(2R + 2)k 4Y analyze the MTTKRP cost using two other partitionings of type-2 dimensions. Without loss of generality, we only consider the cases where µ 1 = {1,...,K} and µ 2 = {K +1,...,N + 1,N+2,N+3,N+4}, and µ 1 = {1,...,K,N+1} and µ 2 = {K +1,...,N,N+2,N+3,N+4}. Considering only the cost of t 1 and t 2 , the former case yieldsC >(2R + 2)k 4 Y

#Figure 8 . 2 :

 82 Figure 8.2: A PACOS code snippet for performing TTMVs in parallel in an iterative manner.

Table 2 .

 2 1: Real-world tensors used in the experiments.

	Tensor	I 1	I 2	I 3	I 4 #nonzeros
	Amazon	6.6M 2.4M	23K	-	1.3B
	Delicious 1.4K 532K	17M 2.4M	140M
	Flickr	731 319K	28M 1.6M	112M
	Netflix	480K	17K	2K	-	100M
	NELL	3.2M	301 638K	-	78M
	NELL2	1.6M	297 338K	-	2M

  Algorithm 4 DTREE-TTV: Dimension tree-based TTV with R vectors in each mode

	Input: t: A node of the dimension tree
	Output: Tensors T : (t) of t are computed. 1: if EXISTS(T : (t)) then 2: return 3: DTREE-TTV(P (t)) 4: for r =1,...,Rdo	I Tensors T : (t) are already computed. I Compute the parent's tensors T : (P (t)) first. I Now update all tensors of t using parent's.

Table 3 .

 3 1: Sequential CP-ALS run time per iteration. Timings are in seconds for splatt, whereas we report the relative speedup with respect to splatt for other methods.

		splatt ht-tree2 ht-tree3 ht-tree
	Delicious	66.6	1.11	*	1.98
	Flickr	43.6	1.23	*	1.98
	Netflix	8.2	0.60	*	0.98
	NELL	8.3	0.87	*	1.28
	Amazon	214.6	0.99	*	1.02
	Random4D	224.7	1.43	*	2.05
	Random8D	1527.1	2.08	2.70	3.97
	Random16D	4401.6	1.31	2.02	3.94
	Random32D 19919.9	-	2.38	5.96

Table 3 .

 3

		splatt ht-tree2 ht-tree3 ht-tree
	Delicious	8.3	0.93	*	2.00
	Flickr	5.8	1.09	*	1.81
	Netflix	0.7	0.55	*	0.87
	NELL	1.3	1.11	*	1.46
	Amazon	24.4	0.86	*	0.95
	Random4D	20.1	0.91	*	1.47
	Random8D	86.9	0.81	1.69	2.18
	Random16D	349.2	0.82	1.73	3.47
	Random32D 1601.8	-	2.18	5.65

2: Shared memory parallel CP-ALS run time per iteration (in seconds). Timings are in seconds for splatt, whereas we report the relative speedup with respect to splatt for other methods.

Table 3 .

 3 4: Memory usage of different methods (in GBs).

		factors	splatt ht-tree2 index index index buffer index buffer ht-tree3 ht-tree
	Delicious	3	7	17.1	*	*	8.5	5.7
	Flickr	4.5	4.8	11	*	*	6.5	4.3
	Netflix	0.1	3.4	6.3	*	*	4.4	1.4
	NELL	0.6	2.5	4.7	*	*	3.2	0.5
	Amazon	1.4	49.6	91.2	*	*	65.4	65.4
	Random4D	6	9.1	20.1	*	*	10.2	10.2
	Random8D	11.9	21	89.1	22.7	15.3	22.7	30.5
	Random16D	23.8	44.9	373.5	66.1	15.3	55.4	45.7
	Random32D	47.7	94.3	OOM	226	15.3 123.8	61

Table 3 .

 3 5: Time spent per iteration (in seconds) for our distributed memory parallel CP-ALS using two threads per core with different partitions. R and C correspond to number of nodes (MPI ranks) and cores used in each instance, respectively.

	R⇥C	Delicious med-gd fine-rd fine-lb fine-hp med-gd fine-rd fine-lb fine-hp Flickr
	8 ⇥ 1	45.078 69.350 46.741 45.256	29.870	-25.674 26.021
	8 ⇥ 16	3.192	9.717	3.603	2.487	2.394	-	2.011	1.512
	16 ⇥ 16	2.589	7.156	2.297	1.454	1.603	9.476	1.091	0.862
	32 ⇥ 16	1.612	4.995	1.458	0.873	1.120	6.365	0.694	0.479
	64 ⇥ 16	1.315	2.907	0.899	0.539	1.013	3.627	0.353	0.286
	128 ⇥ 16	0.822	1.746	0.533	0.351	0.652	1.911	0.222	0.179
	256 ⇥ 16	0.603	1.115	0.318	0.246	0.554	1.093	0.173	0.144
	R⇥C	Netflix med-gd fine-rd fine-lb fine-hp med-gd fine-rd fine-lb fine-hp NELL
	4 ⇥ 1	27.656 32.233 28.118 29.336		19.540 22.528 20.215 19.967
	4 ⇥ 16	1.095	1.685	1.191	1.185		1.191	1.927	1.205	1.109
	8 ⇥ 16	0.617	1.247	0.697	0.688		0.749	1.455	0.700	0.681
	16 ⇥ 16	0.360	0.988	0.421	0.420		0.448	0.998	0.433	0.444
	32 ⇥ 16	0.222	0.804	0.267	0.261		0.287	0.733	0.282	0.318
	64 ⇥ 16	0.138	0.660	0.178	0.173		0.179	0.537	0.201	0.239
	128 ⇥ 16	0.097	0.525	0.117	0.102		0.127	0.398	0.153	0.155
	256 ⇥ 16	0.086	0.410	0.117	0.106		0.099	0.316	0.124	0.119
			R⇥C	Amazon med-gd fine-rd fine-lb
			64 ⇥ 1		-65.545 36.303
			64 ⇥ 16		-	8.803	1.874
			128 ⇥ 16	1.141	6.192	1.007
			256 ⇥ 16	1.056	4.229	0.570

Table 3 .

 3 6: Load balance and communication statistics for 256-way partitioning.

	Partitioning	Nnz Max.	Avg.	Comp. Load Max. Avg.	Comm. Vol. Max. Avg. Max. Avg. Num. Msg
				Delicious
	fine-hp	547K	547K	1947K	1807K	199K	137K 2039
	fine-lb	564K	547K	1849K	1737K	309K	265K 2040
	fine-rd	550K	547K	2747K	2737K 1083K 1080K 2040
	medium-gd	598K	547K	2353K	2214K	624K	571K 1096
				Flickr
	fine-hp	441K	441K	1334K	1221K	54K	38K 1827
	fine-lb	454K	441K	1335K	1257K	107K	87K 2040
	fine-rd	443K	441K	2318K	2308K 1042K 1038K 2040
	medium-gd	443K	441K	1826K	1806K	576K	558K 1152
				Netflix
	fine-hp	392K	392K	1439K	1211K	49K	19K 1380
	fine-lb	404K	392K	1252K	1208K	48K	39K 1530
	fine-rd	394K	392K	1681K	1674K	412K	411K 1530
	medium-gd	394K	393K	1184K	1177K	26K	24K	642
				NELL
	fine-hp	307K	307K	1219K	787K	46K	23K 1513
	fine-lb	316K	307K	920K	888K	89K	84K 1522
	fine-rd	309K	307K	1211K	1205K	271K	269K 1530
	medium-gd	310K	307K	871K	855K	58K	51K	583
				Amazon
	fine-lb	5104K 4955K 16871K 16241K	211K	203K 1503
	fine-rd	4962K 4955K 20964K 20940K 4656K 4651K 1530
	medium-gd 19984K 4955K 50023K 16255K	230K	170K	550

Table 4 .

 4 1: Per-iteration CP-ALS runtime results (in seconds) for sequential, single-socket (14 threads) and dual-socket parallel executions of all methods.

	Method	ten-4D ten-8D ten-16D ten-32D
			P =1
	splatt	175.1	791.2	3766.9 19994.8
	hypertensor	98.3	306.3	929.3	2873.8
	cp-eps	130.3	284.9	586.9	1199.4
			P = 14
	splatt	15.1	68.4	280.3	1292.2
	hypertensor	11.8	33.3	88.5	354.5
	cp-eps	13.1	27.5	56.12	111.31
			P = 28
	splatt	11.5	47.2	190.7	683.8
	hypertensor	13.0	36.6	69.1	215.0
	cp-eps	8.3	17.5	36.2	65.3

Table 4 .

 4 2: Initial setup time (in seconds) for parallel CP-ALS.

	Method | Data ten-4D ten-8D ten-16D ten-32D
	splatt	60	84	232	617
	hypertensor	70	167	418	983
	cp-eps	16	18	23	41

Table 5 .

 5 2: Statistics for the computation and communication requirements with different partitionings of Flickr in one HOOI iteration for 4096-way parallel run with 256 MPI ranks.

	Mode	W TTMc Max Avg	W TRSVD Max Avg	Comm. vol. Max Avg
			fine-hp
		441K 441K	590	507	2218	2029
		441K 441K 8778 5656	24K	17K
		441K 441K 221K 110K	166K	11K
		441K 441K	32K	19K	77K	53K
			fine-rd
		443K 441K	668	648	5884	2597
		443K 441K	98K	96K	409K	385K
		443K 441K 545K 435K 1744K 1735K
		443K 441K 110K 100K	432K	413K
			coarse-hp
		718K 441K	22	3	4910	1213
		810K 441K 2700	248	118K	66K
		798K 441K 197K 110K 3187K	810K
		2368K 441K	13K 6250	170K	102K
		coarse-block
		958K 441K	252	3	18K	907
		756K 441K 5401	248	126K	80K
		441K 441K 130K 110K 3324K 1250K
		2518K 441K	60K 6250	296K	138K

Table 5 .

 5 3: Relative timings of TTMc, TRSVD, and core tensor computation steps within HOOI using 256-way fine-hp partition (in percentage).

	Step	Delicious Flickr NELL Netflix
	TTMc	75.6	64.6	71.2	27.7
	TRSVD+comm	19.2	32.6	24.8	71.6
	core+comm	5.2	2.8	4.0	0.7

Table 5 .

 5 4: Time spent per iteration (in seconds) for shared memory parallel HOOI on node(s) with a 16-core processor. The number of MPI ranks used for each tensor is shown in the parentheses.

	#threads Delicious (8) Flickr (8) NELL (1) Netflix (4)
	1	1182.7	1055.8	2173.6	660.1
	2	634.5	583.2	1146.3	330.8
	4	361.1	354.5	616.4	167.6
	8	227.5	241.6	354.1	87.3
	16	173.2	201.0	252.7	48.7
	32	164.9	206.2	222.7	33.7

Table 5 .

 5 5: Comparison of the sequential and 2048-way parallel run times (in seconds) using fine-hp partitioning, with (dtree) and without (flat) dimension trees

	TTMc scheme/#procs Netflix NELL Delicious Flickr
	flat/1	166	176	669	772
	dtree/1	123	129	560	678
	flat/2048	0.44	1.3	2.09	1.04
	dtree/2048	0.39	1.2	0.97	0.90

Table 6 .

 6 Submatrix of A owned by process p I p ,J p Set of rows/columns of W/H owned by process p F p ,G p Set of unique row and column indices of A p W(I p , :) Owned rows of W by process p H(:, J p ) Owned columns of H by process p

		Input matrix
	W	Left low rank factor
	H	Right low rank factor
	m	Number of rows of A
	n	Number of columns of A
	k	Low rank parameter
	P	Number of parallel processes
	P r P c	Number of rows in the processor grid Number of columns in the processor grid
	A p

1: Notation.

  Algorithm 16 DTREE-TTM: Performing TTM on a dimension tree

	Input: t: A dimension tree node
	Output: T (t): The tensor of t
	1: if EXISTS(T (t)) then
	2:	return T (t)
	3: T (t)	DTREE-TTM(P (t))
	4: for d 2 µ(t) \ µ(P (t)) do
	5:

  1: if EXISTS(T : (t)) then

	2:	return
	3: δ(t)	µ(P (t)) \ µ(t)
	4: DTREE-TTV(P (t))
	5: if IS-ROOT(P (t)) then
	6: 7:	M for r=1...Rdo T (P (t)) (µ(t)) ( (n2δ(t)) U (n) )
	8:		vec(T r (t))	M(:,r)
	9: else
	10: 11:	for r=1...Rdo vec(T r (t))	T r (P (t)) (µ(t)) (⌦ (n2δ(t)) u	(n)

r )

  )⇡(t 2 ) . By taking the difference of (7.5) and (7.6) we finally obtain)⇡(t 2 ) -4⇡(t 1 )⇡(t 2 ) -⇡(t 2 ) -⇡(t 1 )

				t) ⇡(t i )	+	⇡(t) ⇡(t 12 )
	+4⇡(t 12 )+	⇡(t 12 ) ⇡(t 1 )	+	⇡(t 12 ) ⇡(t 2 )	+	K i=1 X	C(t i ).	(7.6)
	Let L = ⇡(t) ⇡(t 1 C k -C k-1 =R 2⇡(t)+ ⇡(t 1 ≥R ⇡(t) ⇡(t 1 ) + ⇡(t) ⇡(t 2 ) ⇡(t) -⇡(t) ⇡(t 1 ) + ⇡(t) ⇡(t 2 ) -⇡(t) ⇡(t 1 )⇡(t 2 ) -⇡(t 2 ) -⇡(t 1 )
	=R ⇡(t

2 )L + ⇡(t 1 )L -L -⇡(t 2 ) -⇡(t 1 ) =R (L -1) ⇡(t 2 )+⇡(t 1 ) -L ≥R 4(L -1) -L ≥ 2R,

(7.7)

as Q (t k ) ≥ 2 for all k 2 N K , and L ≥ 2 (since I n ≥ 2 for all n 2 N N

  )⇡(t 2 ) -4⇡(t 1 )⇡(t 2 ) -⇡(t 2 ) -⇡(t 1 ) =⇡(t)+R L⇡(t 1 )⇡(t 2 )+L⇡(t 2 )+L⇡(t 1 ) -L -4⇡(t 1 )⇡(t 2 ) -⇡(t 2 ) -⇡(t 1 ) =⇡(t)+R (L -4)⇡(t 1 )⇡(t 2 )+(L -1)(⇡(t 1 )+⇡(t 2 )) -L .For L =3, we obtainC k -C k-1 =⇡(t)+R -⇡(t 1 )⇡(t 2 ) + 2(⇡(t 1 )+⇡(t 2 )) -3 .As ⇡(t 1 ),⇡(t 1 )  L =3, we only have four possibilities for ⇡(t 1 ) and ⇡(t 2 ). For ⇡(t 1 )=⇡(t 2 )=3, we get C

									K X i=3	⇡(t) ⇡(t i )	+	⇡(t) ⇡(t 12 )
	+4⇡(t 12 )+	⇡(t 12 ) ⇡(t 1 )	+	⇡(t 12 ) ⇡(t 2 )	+	K i=1 X	C(t i ).	(7.9)
	By taking the difference of (7.8) and (7.9) we obtain
	C k -C k-1 =⇡(t)+R ⇡(t)+	⇡(t) ⇡(t 1 )	+	⇡(t) ⇡(t 2 )	-	⇡(t) ⇡(t 1

1)⇡(t)+R (K -1)⇡(t)+ (7.10) For L ≥ 4,(7.10) yields

C k -C k-1 >R (L -1)(⇡(t 1 )+⇡(t 2 )) -L ≥R 4(L -1) -L ≥ 8R. k -C k-1 =⇡(t)+R -9 + 2(6) -3 = ⇡(t) > 0. ⇡(t 1 )=2,⇡(t 2 )=3and ⇡(t 1 )=3,⇡(t 2 )=2yields

datasets amenable both in terms of computational and memory requirements for real world applications.

Authors are listed alphabetically, except for 'in submission' 1.

Acknowledgements

The work in this thesis is performed using the compute resources generously provided by ENS de Lyon (workstations and PSMN cluster), IDRIS Grants (i2016067501) for Ada and Turing supercomputers, and Rhea cluster at the Oak Ridge Leadership Computing Facility (OLCF).

Part II

Dense Tensor Decompositions

Consider the permutation σ of N N that minimizes the cost of the multiplication X ⇥ σ(1) U (σ(1)) ⇥ σ(2) (N )) . This cost is C σ = T (R σ(1) +↵ σ(1) R σ(2) +↵ σ(1) ↵ σ(2 R σ(3) +•••+↵ σ(1) ...↵ σ(N -1) R σ (N ) ). Assume by contradiction that σ is not the identity permutation, and let i be the first index such that β σ(i) >β σ(i+1) . We rewrite the cost C σ as C σ = C 1 + C 2 + C 3 where

) corresponds to terms i and i +1;

)R σ(j) corresponds to the last ni -1 terms. Now let us permute dimensions σ(i) and σ(i + 1). Formally, this amounts to replacing σ by σ ⇤ = τ i,i+1 • σ, where τ i,i+1 is the transposition that exchanges elements in position i and i +1. The costs C 1 and C 3 are not modified, so that

Hence, permuting dimensions σ(i) and σ(i + 1) does decrease the optimal cost, the desired contradiction. This concludes the proof.

In Algorithm 17, we provide the HOOI algorithm that performs TTMs using a dimension tree T and the TTM routine in Algorithm 16. In the ALS subiteration for mode n, it starts at Line 6 with destroying all tensors in the tree that involves a multiplication with U (n) , as U (n) will subsequently be updated in that subiteration, invalidating these tensors. Note that every tree node t that does not lie in the path from L n to ROOT(T ) has n 2 µ(t) 0 , hence their tensors are destroyed. Next, DTREE-TTM is performed for the leaf node L n at Line 7 to compute Y = X ⇥ -n U (n) . Algorithm 16 only computes tensors of nodes in the path from L n to ROOT(T ), and all other tensors not lying on this path are already destroyed at Line 6; therefore, at any instant of Algorithm 17, only the tensors of nodes lying on a path from a leaf to the root are stored. Following the computation of TTM for L n , U (n) is similarly updated by performing a truncated SVD on Y (n) at Line 8. After all N ALS subiterations, the core tensor G is formed using U (n) and the last Y corresponding to X ⇥ -n U (n) at Line 9.

Since L 1 ,...,L N corresponds to a post-order traversal of the tree, the tensor of each tree node is computed when DTREE-TTM is called for its first leaf descendant node, stays valid and reused throughout the subiterations of all its leaf descendants, and is destroyed in the subiteration following its last leaf descendant. This implies the tensor of every tree node is computed and destroyed exactly once per HOOI iteration.

An important point in Algorithm 17 is that T can be of any shape (balanced, unbalanced, kary), and can employ arbitrary partitions of dimensions at each level. Indeed, each tree configuration yields a different computational cost for TTMs; hence, we are interested in finding a tree topology that minimizes this cost. Intuitively, we argue that using a binary tree is a good choice because of the following observation. For a non-leaf node t 2 T with children t 1 ,...,t k (k ≥ 2), µ(t) is partitioned into k disjoint sets, namely µ(t 1 ),...,µ(t k ). Performing DTREE-TTM on each child t l for l 2 N k requires TTM in dimensions µ 0 (t l )=µ(t) \ µ(t l ). Overall, calling DTREE-TTM on all

The last possibility ⇡(t 1 )=2and ⇡(t 2 )=2gives

Finally, for L =2,(7.10) becomes

Note that in this case, the only option is to have ⇡(t 1 )=⇡(t 2 )=2 , which also implies that X 2 R 2⇥2⇥2 , hence ⇡(t)=8. As a result, we obtain

for all R ≥ 1. Therefore, we conclude that for all tensors except those in R 2⇥2⇥2 , the optimal dimension tree is binary.

In the next theorem, we show that finding an optimal dimension tree is NP-hard.

Theorem 7 (Optimal dimension tree for CP decomposition). Finding an optimal dimension tree that minimizes the cost of TTVs in computing CP decomposition is NP-hard.

Proof. We perform a reduction from PRODUCT-PARTITION using a set S = {s 1 ,...,s N } of positive integers. We construct an N +4-dimensional tensor X 2 R I 1 ⇥•••⇥I N +4 whose first N dimensions correspond to the elements of S, i.e., I n = s n for n 2 N N . We let

where k is to be determined appropriately in the course of the proof. We similarly name the first N and the last 4 dimensions as type-1 and type-2, respectively.

We now analyze the cost of an optimal dimension tree that minimizes the MTTKRP cost in executing CP-ALS for X . This tree must be a BDT as suggested by Theorem 6. Let t be the root this with children t 1 and t 2 . In an optimal BDT, we expect µ(t 1 ) to consist of dimensions corresponding to elements in an optimal solution S 0 ⇢ S of PRODUCT-PARTITION. All dimensions corresponding to S \ S 0 are similarly expected to belong to µ(t 2 ). In this scenario, we would only have three configurations for the partitioning of type-2 dimensions to µ(t 1 ) and µ(t 2 ) due to symmetry, namely N +1,N +2,N +3,N +4 2 µ(t 1 ), N +1,N +2,N +3 2 µ(t 1 ) and N +4 2 µ(t 2 ), and N +1,N +22 µ(t 1 ) and N +3,N +42 µ(t 4 ).

We first consider the case where N +1,N+2 2 µ(t 1 ) and N +3,N+4 2 µ(t 4 ), and analyze the cost of a BDT for a given partition of type-1 dimensions to µ(t 1 ) and µ(t 2 ). Without loss of generality, let µ(t 1 )={1,...,K,N +1,N +2} and µ(t 2 )={K +1,...,N,N +3,N +4} represent a such partition for some K, 0  K  N (all possible partitions can be obtained by a proper K and a permutation of the elements of S). We can express the cost of this BDT as

Part III

Sparse Tensor Software

PACOS: A partitioning and communication framework for sparse irregular applications

PACOS is a library providing a unified partitioning and communication framework for developing parallel sparse (multi)linear algebra routines, including but not limited to sparse tensor factorization (CP, Tucker, tensor completion, etc.), sparse matrix factorization variants (PCA, non-negative), and algorithms involving sparse matrix-vector multiplication (iterative solvers, page rank, etc.).

The execution of a parallel algorithm using PACOS involves two steps. The first one is the partitioning step, in which nonzero elements of a sparse struct (tensor or matrix) as well as vector entries (or matrix rows) in each dimension are partitioned. It supports many partitioning routines for nonzero elements such as randomized checkerboard partitioning, checkerboard hypergraph partitioning, finegrain hypergraph partitioning, and fine-grain random partitioning. It provides effective heuristics for partitioning vector entries in a way that balances the number of owned vector entries by processes as well as the amount of communication volume per process. Finally, it partitions the sparse struct accordingly, and creates per-process sparse struct files along with vector partition maps for all dimensions.

Once the sparse struct is partitioned, each process can read its data and vector mappings, and create a communicator for each dimension using PACOS, which is formed conformally with the partition created in the partitioning phase. Finally, each process can execute its computations on its local data, and call PACOS communication routines to carry out sends and receives at the synchronization points. PACOS provides fold and expand communication routines enabling a fine-grain algorithm, whereas coarse-grain variants can similarly be developed by skipping one of the two communication routines. PACOS employs non-blocking MPI routines underneath, giving the possibility of overlapping communication with local computations. It also uses a lightweight MPI wrapper that provides templating and 64-bit integer compatibility with MPI. Indeed, there are other software packages providing similar functionalities, ZOLTAN and PETSC to name a few, yet the main convenience of PACOS is providing a unified partitioning and communication framework geared towards matrix and tensor factorization variants, and giving a drop-in support for effective nonzero and vector partitioning routines through hypergraph partitioning. It is owing to this functionality that we were able to replace the entire set of partitioning and communication routines of NMFLIBRARY with PACOS using less than 100 lines of code modifications, and thereby achieve scalability up to 32K cores on an IBM BlueGene/Q supercomputer, which were simply impossible with its existing parallelization scheme.

To illustrate the usage, in Figure 8.2, we provide a PACOS code snippet performing iterative TTMVs in an alternating manner over the dimensions. Note that this concise code is the backbone of most iterative tensor factorization algorithms. PACOS effectively handles all the communications underneath, and also provides the potential to overlap computations with communications.

CHAPTER 9. CONCLUSION

techniques to obtain high performance on a NUMA parallel architecture. In overall, in comparison to the state of the art, this yields a parallel algorithm achieving up to 16.7x speedup in sequential and 10.5x speedup in parallel executions, requiring up to 24x less data preprocessing time, and using up to O(log N ) less memory for storing intermediate computations.

Chapter 5 adopts our parallelization techniques for computing CP decomposition to the computation of the Tucker decomposition for sparse tensors. In this chapter, we introduce a coarse-and a fine-grain algorithm for computing the Tucker of sparse tensors. The coarse-grain variant avoids an expensive fold communication step whereas limiting the possibility for a good partition, whereas the fine-grain variant provides the most flexibility in partitioning for establishing load balance and reducing communication costs, but incurs the fold communication step. To render the fine-grain algorithm amenable, we couple the fold communication with an iterative truncated SVD algorithm, which effectively reduces the unit communication cost of the fold communication from R N -1 to K where K is the number of iterations of the truncated SVD algorithm in computing in a rank-R Tucker decomposition. Finally, we similarly employ hypergraph partitioning to significantly reduce the number of unit communications required to execute the parallel algorithm. All these techniques yield scalability up to 4096 cores on an IBM BlueGene/Q supercomputer using real-world sparse tensors.

In Chapter 6, we propose a distributed memory parallelization of the non-negative matrix factorization problem, which shares a great amount of similarities with the CP decomposition in 2 dimensions. Our algorithm is the first to use a fine-grain parallelism and a point-to-point messaging scheme. We perform an MPI-only parallelization of the existing implementation in the NMFLIBRARY, and use 2D partitionings that eliminate the need for intra-node communication in one dimension. We similarly used appropriate hypergraph models to reduce communication with an effective partitioning. All these contributions yield a parallel algorithm achieving scalability up to 32768 cores on an IBM BlueGene/Q supercomputer.

We investigate the use of dimension trees in computing CP and Tucker decompositions of dense tensors Chapter 7. We provide the dimension tree-based CP-ALS and HOOI algorithms for dense tensors, and show many theoretical results related to NP-completeness of finding an optimal binary tree in both cases. We also introduce a greedy algorithm that finds the optimal ordering for a series of TTMs, which constitutes the most expensive step in HOOI algorithm.

In the course of this thesis, we developed two software packages, namely HYPERTENSOR and PACOS, that address the challenges in high performance parallel computations of sparse CP and Tucker decompositions, which we detailed in Chapter 8.H YPERTENSOR uses efficient tree-based sparse tensor data structures and computational schemes to execute CP-ALS and HOOI algorithms. It uses OpenMP for shared memory parallelization, and depends on PACOS for distributed memory parallelization and partitioning routines. PACOS provides a generic communication and partitioning framework for sparse (multi)linear algebra and machine learning applications, including but not limited to tensor factorization, matrix factorization, iterative solvers, and iterative graph algorithms. With these two software packages, direct outcomes of this thesis research, we provide the fastest available sparse tensor factorization software as of today.

Altogether, with all algorithmic contributions, parallelization techniques, and optimized software implementations, the work in this thesis enables performing tensor-and matrix-based analysis of big data problems having billions of entries in a matter of seconds on modern HPC platforms, and sets a new bar for the state of the art in the domain of high performance parallel tensor computations.

Perspectives

Even though we have not focused on it in this thesis, there has been recent work on the parallelization of tensor factorization on accelerators, including GPUs and Xeon Phis. These methods use the traditional representations of sparse tensors and employ tiling strategies to better utilize the available high-bandwidth memory of the accelerator. Our tree-based tensor storage and computational schemes provide significant computational advantages over the traditional techniques, and require non-trivial techniques to be adapted for the accelerator use. This extension would be a natural future direction for our work.

In the longer term, there are still computational challenges left to be addressed in tensor decompositions. HPC techniques for computing non-negative tensor factorization and coupled tensor-matrix and -tensor factorizations are still to be tackled, which might find significant value in modern big data applications. From an application perspective, the work in this thesis renders tensor decomposition methods amenable to use on datasets of massive scale, which was simply impractical before. As a result, we hope the outcome of this thesis to significantly expand the applications of tensor decompositions on big data analysis and machine learning problems including but not limited to graph analysis, web link analysis, healthcare data analytics, and recommender systems, in the years to come.