Jiajia Li
email: jiajia.li@pnnl.gov

Bora Uçar

Ümit V Çatalyürek

Jimeng Sun
email: jsun@cc.gatech.edu

Kevin Barker
email: kevin.barker@pnnl.gov

Richard Vuduc

Efficient and effective sparse tensor reordering

Keywords: Sparse tensor, tensor decomposition, reordering, hierarchical coordinate, HiCOO

published or not. The documents may come

INTRODUCTION

We consider the problem of how to improve memory reference locality for sparse tensor computations [START_REF] Bader | Efficient MATLAB computations with sparse and factored tensors[END_REF][START_REF] Kolda | Tensor Decompositions and Applications[END_REF]. A tensor is a multiway (N -way) array, with a vector and a matrix being the 1-way and 2-way (rows and columns) specializations thereof. It is sparse if most of its entries are zero and, therefore, need not be stored nor explicitly computed upon. A sparse tensor is often a natural way to represent a multifactor or multirelational dataset, and has found numerous applications in data analysis and mining [START_REF] Bader | Efficient MATLAB computations with sparse and factored tensors[END_REF][START_REF] Baskaran | Efficient and scalable computations with sparse tensors[END_REF][START_REF] Choi | Blocking Optimization Techniques for Sparse Tensor Computation[END_REF][START_REF] Cichocki | Era of Big Data Processing: A New Approach via Tensor Networks and Tensor Decompositions[END_REF][START_REF] Lathauwer | Decompositions of a Higher-Order Tensor in Block Terms-Part III: Alternating Least Squares Algorithms[END_REF][START_REF] Ho | Marble: High-throughput Phenotyping from Electronic Health Records via Sparse Nonnegative Tensor Factorization[END_REF][START_REF] Kaya | Scalable Sparse Tensor Decompositions in Distributed Memory Systems[END_REF][START_REF] Kolda | Tensor Decompositions and Applications[END_REF][START_REF] Li | An input-adaptive and in-place approach to dense tensor-times-matrix multiply[END_REF][START_REF] Liu | A Unified Optimization Approach for Sparse Tensor Operations on GPUs[END_REF][START_REF] Ma | Optimizing sparse tensor times matrix on GPUs[END_REF][START_REF] Perros | Sparse Hierarchical Tucker Factorization and its Application to Healthcare[END_REF][START_REF] Perros | SPARTan: Scalable PARAFAC2 for Large & Sparse Data[END_REF][START_REF] Ravindran | Memory-Efficient Parallel Computation of Tensor and Matrix Products for Big Tensor Decompositions[END_REF][START_REF] Smith | Streaming Tensor Factorization for Infinite Data Sources[END_REF][START_REF] Smith | An Exploration of Optimization Algorithms for High Performance Tensor Completion[END_REF][START_REF] Smith | SPLATT: Efficient and Parallel Sparse Tensor-Matrix Multiplication[END_REF] for health care [START_REF] Ho | Marble: High-throughput Phenotyping from Electronic Health Records via Sparse Nonnegative Tensor Factorization[END_REF][START_REF] Wang | Rubik: Knowledge Guided Tensor Factorization and Completion for Health Data Analytics[END_REF], natural language processing [START_REF] Kang | GigaTensor: Scaling Tensor Analysis Up by 100 Times -Algorithms and Discoveries[END_REF], machine learning [START_REF] Abadi | TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems[END_REF][START_REF] Anandkumar | Tensor Decompositions for Learning Latent Variable Models[END_REF][START_REF] Novikov | Tensorizing Neural Networks[END_REF], and social network analytics [START_REF] Papalexakis | ParCube: Sparse Parallelizable Tensor Decompositions[END_REF], among many others.

Like sparse matrices, there are a variety of data structures and techniques for storing and operating efficiently with sparse tensors [START_REF] Bader | Efficient MATLAB computations with sparse and factored tensors[END_REF][START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF][START_REF] Liu | A Unified Optimization Approach for Sparse Tensor Operations on GPUs[END_REF][START_REF] Smith | Tensor-Matrix Products with a Compressed Sparse Tensor[END_REF]. Moreover, one technique is to reorder the tensor, which means relabeling the indices-and, therefore, reorganizing the nonzero structure of the tensor-in the hopes of improving spatial or temporal locality when operating with it. One form of reordering explored in prior work is to "sort" the input tensor [START_REF] Kaya | Scalable Sparse Tensor Decompositions in Distributed Memory Systems[END_REF][START_REF] Liu | A Unified Optimization Approach for Sparse Tensor Operations on GPUs[END_REF][START_REF] Smith | SPLATT: Efficient and Parallel Sparse Tensor-Matrix Multiplication[END_REF]. In this paper, we formalize the reordering problem and propose and evaluate new heuristic algorithms for it.

While tensor reordering can improve locality, it can also aggravate load balance as measured by nonzeros assigned to each thread. We observe that reordering increases imbalance by as much as 6.7× in practice. In this paper, we propose improvements to the data structure and corresponding algorithms that address this problem.

Our evaluation is in the context of a widely used tensor computation known as the CANDECOMP/PARAFAC decomposition (CPD), a kind of generalization of the (truncated) singular value decomposition (SVD) from matrices to tensors [START_REF] Kolda | Tensor Decompositions and Applications[END_REF]. Within a CPD, the most costly computational kernel is the matricized tensor-times-Khatri-Rao product, or Mttkrp [START_REF] Kolda | Tensor Decompositions and Applications[END_REF][START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF][START_REF] Liu | A Unified Optimization Approach for Sparse Tensor Operations on GPUs[END_REF][START_REF] Smith | SPLATT: Efficient and Parallel Sparse Tensor-Matrix Multiplication[END_REF]. Thus, we examine the effects of reordering on the locality of sparse Mttkrp operations, which we then evaluate within a CPD on a modern multicore CPU platform. We also consider this evaluation when using three state-of-the-art tensor storage formats, known as coordinate (COO), compressed sparse fiber (CSF), and hierarchical COO (HiCOO). Thus, this evaluation gives practitioners a sense of when reordering pays off in commonly occurring scenarios.

The main claimed contributions of this work are as follows.

• We propose two heuristic tensor reordering schemes, BFS-MCS and Lexi-Order, to enhance data locality by relabeling mode indices. ( § 3) • We improve the state-of-the-art parallel strategy, superblock scheduling [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF], by suggesting a partition-based version thereof. ( § 4) • Lexi-Order achieves up to 4.14× speedup on sequential HiCOO-Mttkrp and 11.88× speedup on its parallel counterpart. BFS-MCS obtains up to 1.88× speedup on sequential HiCOO-Mttkrp and 1.94× speedup on its parallel case. ( § 5)

• For COO and CSF formats, Lexi-Order improves the performance of sequential COO-and CSF-Mttkrps by up to 4.29× and 2.33× and parallel COO-and CSF-Mttkrps by up to 1.48× and 1.86× respectively. BFS-MCS and Lexi-Order behave better than the state-of-the-art reordering based on graph and hypergraph partitionings [START_REF] Smith | SPLATT: Efficient and Parallel Sparse Tensor-Matrix Multiplication[END_REF] on CSF-Mttkrp. ( § 5) • We further parallelize the reordering and HiCOO conversion to reduce the pre-processing overhead. When CPD iterates for tens of times, which is common, then the reordering is worthwhile. ( § 5) Table 1: List of symbols and notation.

BACKGROUND

We summarize the related symbols and notation of tensors, CPD, and the three state-of-the-art tensor formats (COO, CSF, HiCOO) in Table 1.

i = 1,…,I j = 1, …, J k = 1 , … , K Figure 1: A third-order tensor X ∈ R I ×J ×K .

Tensors and CPD

The order of a tensor, N , is the number of its dimensions or modes.

We follow the notation in Kolda and Bader's survey [START_REF] Kolda | Tensor Decompositions and Applications[END_REF]. A firstorder tensor (N = 1) is a vector, denoted by a boldface lowercase letter, e.g., v; A second-order tensor (N = 2) is a matrix, denoted by a boldface capital letter, e.g., A. Higher-order tensors (N ≥ 3) are denoted by bold capital calligraphic letters, e.g., X. We show an example of a sparse third-order tensor, X ∈ R I ×J ×K , in Figure 1. In this example, a scalar entry of X at position (i, j, k) is x i jk . We assume an example sparse tensor X ∈ R I 1 ×I 2 ו••×I N with M nonzeros in the following context. CPD decomposes a tensor into a sum of component rank-one tensors [START_REF] Kolda | Tensor Decompositions and Applications[END_REF]. It approximates an N th-order tensor

X ∈ R I 1 ו••×I N as X ≈ R r =1 λ r a (1) r • • • • • a (N ) r ≡ λ; A (1) , . . . , A (N ) , ( 1
)
where R is the canonical rank of tensor X, taken as the number of component rank-one tensors [START_REF] Kolda | Tensor Decompositions and Applications[END_REF]. In a low-rank approximation, R is usually chosen to be a small number less than 100. The outer product of the vectors a produces R rank-one tensors, and A (n) ∈ R I n ×R , n = 1, . . . , N are the factor matrices, each one formed by taking the corresponding vectors as its columns, i.e.,

A (n) = [a (n) 1 a (n) 2 . . . a (n) R ].
We normalize these vectors to unit magnitude and store the factor weights in the vector λ = {λ 1 , . . . , λ r }.

The bottleneck of CPD is the matricized tensor-times-Khatri-Rao product (Mttkrp). Given an N th-order tensor X and matrices A (1) , . . . , A (N ) , the mode-n Mttkrp is

Ã(n) ← X (n) A (N ) ⊙ • • • ⊙ A (n+1) ⊙ A (n-1) ⊙ • • • ⊙ A (1) , (2)
where X (n) is the mode-n matricization (or unfolding) of tensor X, ⊙ is the Khatri-Rao product. Mode-n matricization reshapes a tensor into an equivalent matrix [START_REF] Kolda | Tensor Decompositions and Applications[END_REF]. The Khatri-Rao product is a "matching column-wise" Kronecker product between two matrices. Given matrices A ∈ R I ×R and B ∈ R J ×R , their Khatri-Rao product is denoted by

C = A ⊙ B where C ∈ R (I J )×R , C = A ⊙ B = [a 1 • b 1 , a 2 • b 2 , . . . , a R • b R ] , (3)
where a r and b r , r = 1, . . . , R are columns of A and B, • is the outer product of vectors, a special case of Kronecker product [START_REF] Kolda | Tensor Decompositions and Applications[END_REF]. In particular for a third-order tensor, Mttkrp multiples each nonzero entry x i, j,k with the R-vector formed by the entry-wise product of the jth row of A (2) and kth row of A (3) when computing Ã(1) . Detailed description can be found in Kolda and Bader's survey [START_REF] Kolda | Tensor Decompositions and Applications[END_REF]. In this paper, sparse Mttkrp will mean an Mttkrp between a sparse tensor and dense matrices.

Sparse Tensor Formats

We consider the three state-of-the-art tensor formats (COO, CSF, HiCOO) which are all for general unstructured sparse tensors. Each level corresponds to a tensor mode, and each nonzero is a path from a root to a leaf. The indices of CSF are stored in cinds with the pointers stored in cptrs to indicate the locations of nonzeros at the next level. Since cptrs needs to show the range up to M, we use β int or β long accordingly for differently sized-tensors .

2.2.3

HiCOO. Hierarchical Coordinate (HiCOO) [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF] format derives from COO format, but improves upon it by compressing the indices in units of sparse tensor blocks. HiCOO stores a sparse tensor in a sparse-blocked pattern with a pre-specified block size B, meaning in B ו • •×B blocks. It represents every block by compactly storing its nonzero triples using fewer bits. Figure 2 (c) shows the example tensor given 2 × 2 × 2 blocks (B = 2). For a third-order tensor, bi, bj, bk are block indices in β int bits, indexing tensor blocks; ei, ej, ek are element indices in β byte bits, indexing nonzeros within a tensor block. A bptr array in β long bits stores the pointers of every block's beginning locations, and val saves all the nonzero values.

Compare the three tensor formats: HiCOO and COO treat every mode equally and do not assume any mode order, these preserve the mode-generic orientation [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF]. CSF has a strong mode-specificity, since the tree structure implies a mode ordering for enumeration of nonzeros. Besides, compared to COO format, HiCOO and CSF save storage space and memory footprints whereas achieve higher performance generally.

PROBLEM DEFINITION AND OUR SOLUTIONS

Our objective is to improve the performance of Mttkrp and therefore CPD algorithms based on the three tensor storage formats described above. We will achieve this objective by reordering (or relabeling) the indices in one or more modes of the input tensor so as to improve the data locality in the tensor and the factor matrices of an Mttkrp operation. Take for example two nonzeros (i 2 , j 2 , k 1 ) and (i 2 , j 2 , k 2 ) of a third-order tensor in Figure 3 (a). Relabeling k 1 and k 2 to other two indices close to each other will potentially improve cache hits for both tensor and their corresponding rows of factor matrices in the Mttkrp operation [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF][START_REF] Smith | SPLATT: Efficient and Parallel Sparse Tensor-Matrix Multiplication[END_REF]. Besides, it may also influence the tensor storage for some formats. (Analysis will be illustrated in § 3.3.)

We propose two heuristics for reordering indices. The aim of the heuristics is to arrange the nonzeros close to each other, in all modes. If we were to look at matrices, this would correspond to reordering the rows and columns so that all nonzeros are clustered around the diagonal. This way, nonzeros in a row or column would be close to each other, and any blocking (by imposing fixed sized blocks) would have nonzero blocks only around the diagonal. The proposed heuristics are based on these observations and try to obtain similar behavior for tensors. The output of a reordering algorithm is the permutations for all modes being used to relabel tensor indices of them.

BFS-MCS

BFS-MCS is a Breadth First Search (BFS)-like heuristic approach based on the maximum cardinality search family [START_REF] Tarjan | Simple Linear-Time Algorithms to Test Chordality of Graphs, Test Acyclicity of Hypergraphs, and Selectively Reduce Acyclic Hypergraphs[END_REF]. We first construct a hypergraph for a sparse tensor, where vertices are tensor indices in all modes and hyperedges represent its nonzero entries. For a third-order sparse tensor X ∈ R I 1 ×I 2 ×I 3 with M nonzeros, its hypergraph H = (V , E) consists of |V | = I 1 + I 2 + I 3 vertices and |E| = M hyperedges. A nonzero entry x i 1 i 2 i 3 connects the three vertices i 1 , i 2 , i 3 . Figure 3 (b) shows an example of the hypergraph for a sparse tensor. Vertices are blank circles, and hyperedges are represented by grouping vertices.

For an N th-order sparse tensor, we need to find the permutations for N modes. We determine a permutation for a given mode n (perm n ) of a tensor X ∈ R I 1 ו••×I N as follows. Suppose some of the indices of mode n are already ordered. Then, BFS-MCS picks the next to-be-ordered index as the one with the strongest connection to the currently ordered index set. In the case of ties, it selects the index with the smallest number of nonzeros in the corresponding sub-tensor. Intuitively, in the hypergraph of a sparse tensor, a stronger connection represents more common indices in the modes other than n among an unordered vertex and the already ordered ones. This means more data from factor matrices can be reused in Mttkrp, if found in cache.

This process is implemented by maintaining a max-heap (H n , shown in Line 3 in Algorithm 1) for all mode-n indices (e.g., i n ) with two keys. The primary key of an index is the number of connections to the currently ordered indices. The secondary key is the number of nonzeros in the corresponding (N -1)th-order sub-tensor, e.g., X(:, . . . , :, i n , :, . . . , :), where the smaller secondary key values signify higher priority. Note that the secondary keys are static. Initially, H n is constructed according to the secondary keys of mode-n indices. For each mode-n index (e.g., i n ) of this max-heap, BFS-MCS traverses all connected tensor indices in the modes except n, calculates the number of connections of mode-n indices, and then updates the heap (H n ) using the primary and secondary keys.

Algorithm 1 makes BFS-MCS more precise. Let m V denote a size-|V | array that tracks whether a vertex has been visited. They are initialized to zeros (unvisited) and changed to 1s once after being visited. We record a mode-n index v (0) n , obtained from the max-heap Algorithm 1 BFS-MCS ordering based on maximum cardinality search for a given mode.

Input: An N th-order sparse tensor X ∈ R I 1 ו••×I N , hypergraph G = (V , E)
with vertex weights w, mode n; Output: Permutation perm n ;

1: Initialize w p to zeros.

2: Initialize w s to the minus of number of nonzeros of each corresponding sub-tensor as a static secondary key. 3: Build max-heap H n for mode-n indices with w p and w s as the primary and secondary keys respectively. 4: Initialize m V to zeros.

5: for i n = 1, . . . , I n do 6: v (0) n = GetHeapMax(H n ); 7: perm n (v (0) n ) = i n ; 8: m V (v (0) n ) = 1; 9: for e n ∈ hyperedges of v (0) n do 10: for v ∈ e n , v is not in mode n and m V (v) == 0 do 11: m V (v) = 1; 12:
for e ∈ hyperedges of v and e e n do 13:

v n = vertex in mode n of e; 14: if inHeap(v n ) then 15: w p (v n ) + +; 16: heapUpdateKey(H n , w p (v n )); 17: return perm n ;
H n , to the permutation array perm n . The algorithm visits all its hyperedges (i.e., all nonzeros with i n , Line 9) and their connected and unvisited vertices from the other (N -1) modes except n (Line 10). For these vertices, it again visits their hyperedges (e in Line 12) and then checks if the connected vertices in mode n (v n ) are in the heap (Line 14). If so, the primary key (connectivity) of v n is increased by 1 and the max-heap (H n ) is updated. To summarize, this algorithm locates the sub-tensors of the neighbor vertices of v (0) n to increase the primary key values of the occurred mode-n indices in H n .

The BFS-MCS heuristic has low time complexity. We explain it using tensor terms. The innermost heap update operation (Line 16) costs O(log(I n )). Each sub-tensor is only visited once with the help of the marker array m V . For all the sub-tensors in one tensor mode, the heap update is only performed for M nonzeros (hyperedges). Overall, the time complexity of BFS-MCS is O(N M log(I n )) for an N th-order tensor with M nonzeros, when computing perm n for mode n.

BFS-MCS does not exactly catch the memory access pattern of an Mttkrp. It treats the contribution from the indices of all modes except n equally to the connectivity, thus the connectivity might not match the actual data reuse. Also, it uses a greedy strategy to determine the next-level vertices, which could miss the optimal global orderings, as is common to greedy heuristics for hard problems.

(1, 1, 0, 0)

(1, 0, 1, 0) 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 1 1 1 0 0 1 0 0 0 0 1 0 0 0 0 1 0 original ordered (a) Vector comparison (b) Matrix > Figure 4
: A doubly lexical ordering of a zero-one matrix and its vector comparison operation.

Lexi-Order

Lexi-Order is an extension of doubly lexical ordering of matrices [START_REF] Lubiw | Doubly Lexical Orderings of Matrices[END_REF][START_REF] Paige | Three Partition Refinement Algorithms[END_REF] to tensors. A lexicographic ordering of an integer vector is the standard dictionary ordering of its elements, defined as follows.

Given two equal-length vectors, x and y, we say x ≤ y iff either (i) all elements are the same, i.e., x = y; or (ii) there exists an index j such that x(j) < y(j) and x(i) = y(i) for all 0 ≤ i < j. For example, Figure 4 (a) shows the vector comparison of two zero-one vectors, x ≡ (1, 0, 1, 0) and y ≡ (1, 1, 0, 0). x ≤ y because x(1) < y(1) and

x(i) = y(i) for all 0 ≤ i < 1.
A doubly lexical ordering of a matrix [START_REF] Lubiw | Doubly Lexical Orderings of Matrices[END_REF] is an ordering of the rows and columns of the matrix so that both the row and column vectors are in non-increasing lexicographic order. A row vector is read from left to right, and a column vector is read from top to bottom. Every real-valued matrix has a doubly lexicographic ordering, and such an ordering is not unique (the proof can be found in [START_REF] Lubiw | Doubly Lexical Orderings of Matrices[END_REF]). Figure 4 (b) shows a doubly lexical ordering of a zero-one example matrix. The row vectors are in non-increasing lexicographic order from top to bottom, and the column vectors are in non-increasing lexicographic order from left to right. Assume the ordered matrix in Figure 4 (b) is A, then rows a(1, :) ≥ a(2, : ) ≥ a(3, :) ≥ a(4, :) and columns a(:, 1) ≥ a(:, 2) ≥ a(:, 3) ≥ a(:, 4). Even from this simple example, the ordered matrix shows better data locality than the original one. Doubly lexical ordering has a number of applications, including the efficient recognition of totally balanced, subtree, and plaid matrices and (strongly) chordal graphs. To the best of our knowledge, the merits of this ordering for high performance have not been investigated for sparse matrices.

We propose an iterative algorithm called matLexiOrder for doubly lexical ordering of matrices, where in an iteration either rows or columns are sorted. The known doubly lexical ordering algorithms for matrices, which were first explored by Lubiw [START_REF] Lubiw | Doubly Lexical Orderings of Matrices[END_REF] and then improved by Paige and Tarjan [START_REF] Paige | Three Partition Refinement Algorithms[END_REF], are "direct" ordering methods with a non-linear runtime of O(M log(I + J ) + J ) and O(M + I + J ) space, for an I × J matrix with M nonzeros. We find the time complexity of these algorithms to be too high for our purpose. Furthermore, the data structures are too complex to allow an efficient generalized implementation for tensors. By appealing to a result of Lubiw [START_REF] Lubiw | Doubly Lexical Orderings of Matrices[END_REF]Claim 2.2], one can show that our matLexiOrder algorithm finds a doubly lexical ordering in a finite number of iterations. We do not show the proof, since we do not aim to obtain an exact doubly lexical ordering; a close-by ordering will likely suffice to improve the Mttkrp performance.

We first describe matLexiOrder algorithm for a sparse matrix. This aims to show the efficiency of our iterative approach compared to others. A partition refinement technique is used to order the rows and columns alternatively. Given an ordering of the rows, the columns can be sorted lexically. This is achieved by an order preserving variant of the partition refinement method [START_REF] Paige | Three Partition Refinement Algorithms[END_REF], which is called orderlyRefine. We briefly explain the partition refinement technique for ordering the columns of a matrix. Given an I × J matrix A, all columns are initially put into a single part. Then, A's nonzeros are visited row by row. At a row i, each column part C is split into two parts C 1 = C ∩ a(i, :) and C 2 = C \ a(i, :), and these two parts replace C in the order C 1 ≻ C 2 (empty sets are discarded). Note that this algorithm keeps the parts in a particular order which generates an ordering of the columns. orderlyRefine is used to refine all parts that have at least one nonzero in row i in O(|a(i, :)|) time, where |a(i, :)| is the number of nonzeros of row i. Overall, matLexiOrder costs a linear total time of O(M + I + J ) (for rows and columns ordering) per iteration and O(J ) space. We also observe that only a small number of iterations will be enough (will be shown in § 5.8 for tensors), yielding a more storage-efficient algorithm compared to the prior doubly lexical ordering methods [START_REF] Lubiw | Doubly Lexical Orderings of Matrices[END_REF][START_REF] Paige | Three Partition Refinement Algorithms[END_REF]. matLexiOrder, in particular the use of orderlyRefine routine a few times, is sufficient for our needs. This is so, as we do not need a full lexicographic order, thereby making our approach faster and simpler to implement. Algorithm 2 Lexi-Order for a given mode.

Input: An N th-order sparse tensor X ∈ R I 1 ו••×I N , mode n; Output: Permutation perm n ;
▷ Sort all nonzeros along with all but mode n.

1: quickSort(X, coordCmp); ▷ Matricize X to X (n) . 2: r = compose (inds ([-n]), 1); 3: for m = 1, . . . , M do 4: c = inds(n, m); ▷ Column index of X (n) 5:
if coordCmp(X, m, m -1) == 1 then

6: r = compose (inds ([-n]), m); ▷ Row index of X (n) 7: X (n) (r, c) = val(m);
▷ Use a variation of partition refinement in [START_REF] Paige | Three Partition Refinement Algorithms[END_REF] 8: perm n = orderlyRefine (X (n) ); 9: return perm n ;

▷ Comparison function for two indices of X 10: Function: coordCmp(X, m 1 , m 2 ) 11: for n ′ = 1, . . . , N do 12:

if n ′ ! = n then 13: if m 1 (n ′ ) < m 2 (n ′ ) then 14: return -1; ▷ Entry m 1 < entry m 2 15:
if m 1 (n ′ ) > m 2 (n ′ ) then (1, 2, 3, 0) To order tensors, we propose the Lexi-Order function as an extension of matLexiOrder. The basic idea of Lexi-Order is to determine the permutation of each tensor mode independently, while considering the order in other modes fixed. Lexi-Order sets the indices of the mode to be ordered as the columns of a matrix, the other indices as the rows and sorts the columns as described for matrices (with the order preserving partition refinement method). The precise algorithm appears in Algorithm 2, which we also illustrate in Figure 5 when applied to mode 1. Given a mode n, Lexi-Order first builds a matricized tensor in Compressed Sparse Row (CSR) sparse matrix format by a call to quickSort with the comparison function coordCmp and then by partitioning the nonzeros into the row segments (Lines 3-7). This comparison function coordCmp does a lexicographic comparison of all-but-mode-n indices, which enables efficient matricization. In other words, sorting of the tensor X is the same as building the matricized tensor X (n) by rows in the fixed lexical ordering, where mode n is the column dimension and the remaining modes constitute the rows. In Figure 5, the sorting step orders the COO entries by (j, k) tuples, which then serve as the row indices of the matricized CSR representation of X [START_REF] Abadi | TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems[END_REF] . Once the matrix is built, we construct zero-one row vectors in Like BFS-MCS approach, Lexi-Order also finds the permutations for N modes of an N th-order sparse tensor. Figure 6 (a) illustrates the effect of Lexi-Order on an example 4 × 4 × 3 sparse tensor. The original tensor is converted to a HiCOO representation with block size B = 2 consisting of 5 blocks, with maximum 2 nonzeros per block. After reordering with Lexi-Order, the new HiCOO has 3 nonzero blocks with up to 4 nonzeros per block. Thus, the blocks are denser, which should exhibit better locality behavior. However, this reordering scheme is heuristic. Figure 6

Original COO i j k val 0 0 0 1 0 1 1 2 0 1 2 3 1 2 2 4 1 3 2 5 2 3 2 6 3 0 0 3 1 0 7 8 1 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 0 1 0 0 0 1 1 0 i (j,k) (0,0) (1,0) (1,1) (1,2) (2,2) (3,2) Zero-one Representation Sorted COO 0 3 3 0 0 1 i (j,k) (0,0) (1,0) (1,1) (1,2) (2,2) (3,2) 1 2

JK*I CSR Matrix

Analysis

We take Mttkrp operation to analyze reordering behavior for COO, CSF, and HiCOO formats. Recall that for a third-order sparse tensor X, Mttkrp multiples each nonzero entry x i, j,k with the R-vector formed by the entry-wise product of the jth row of A (2) and kth row of A (3) when computing Ã(1) . The arithmetic intensity of Mttkrp algorithms on these formats is approximately 1/4 [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF]. Thus, Mttkrp can be considered memory-bound for most computer architectures, especially CPU platforms.

We analyze the memory access of HiCOO-Mttkrp algorithm in bytes (Equation ( 4), details in Equation ( 16) of [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF]) to roughly reflect its real performance. β ? is different bit-lengths for HiCOO (Table 1). The block ratio (α b ) and the average slice size per tensor block (c b ) are two critical parameters of HiCOO. Smaller α b and larger c b are favorable for good HiCOO-Mttkrp performance. (Readers could refer to [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF] for more explanation.) Our reordering algorithms tend to closely pack nonzeros together and get denser blocks (larger c b ), which potentially generates less tensor blocks (smaller α b ), thus HiCOO-Mttkrp has less memory access, more cache hits, and better performance further. Take for example two nonzeros (i 2 , j 2 , k 1 ) and (i 2 , j 2 , k 2 ) again. A HiCOO representation after reordering could combine them in a single block. (A detailed example is shown in Figure 6.) Moreover, from the HiCOO storage equation (Equation ( 5), Equation ( 13) from [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF]), smaller α b can also reduce tensor memory requirement as a side benefit. That is, for HiCOO format a reordering scheme should increase the block density, reduce the number of blocks, and increase cache localitythree related performance metrics. Reordering is more beneficial for HiCOO format than COO and CSF formats.

Bytes hicoo ≈ M 8 [2α b • β int + β byte + R min{ 1 c b , 1} • β float ]N (4)
S hicoo ≈ M [α b • β long + α b N • β int + N • β byte ] (5)
COO stores all nonzero indices, so relabeling does not make a difference in its storage. The same is also true for CSF; relabeling does not change the its tree structure, while the order of nodes may change. For an Mttkrp with tensors stored in these two formats, the performance gain from reordering is only the improved data locality and therefore cache behavior. Though it is hard to do theoretical analysis for them, the potential better data locality could also brings performance advantages, but might be less than HiCOO's.

Parallelize Preprocessing

For an efficient parallel HiCOO-Mttkrp algorithm, superblocks, an extra blocking level above blocks, are proposed to increase the workload granularity of scheduling [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF]. A superblock is essentially a "logical" subtensor that can potentially consist of many small blocks. During HiCOO format conversion in Figure 7, we first sort all nonzeros only by i indices in mode 1 ("rowblock sorting"), to ensure that a mode-1 slice is not split between superblocks. Then L × • • • × L nonzero superblocks are partitioned with an additional array lptr to store the beginning pointers of nonzero superblocks in size n l , the number of superblocks. We treat them as independent subtensors to convert each to the physical HiCOO format with B × • • • × B blocks (L ≥ B) through Z-Morton order sorting, partitioning, and compression steps. In the first rowblock sorting we do a quicksort for a single mode where the quicksort is parallelized. Afterwards, the operations on each superblock (Z-order sorting, partitioning, and compression) are naturally parallelized.

HYBRID PARALLEL HICOO-MTTKRP

An important side-effect of the proposed reordering schemes is that it can inadvertently create load imbalance, if one assumes the original "superblock scheduling" scheme proposed for HiCOO [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF]. This section summarizes the issue and proposes a modification to HiCOO's superblock scheduler to mitigate such imbalance. The baseline HiCOO superblock scheduler operates schematically as shown in Figure 8 (a). Each box is a superblock, where the number of nonzeros inside the block is shown in parentheses. Superblocks in the same row all write to the same output area and so are dependent (arrows), while columns indicate independent superblocks. Two parallel strategies are employed in [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF]: direct and privatized parallelization. Simply speaking, direct parallelization parallelizes rows while privatization strategy parallelizes iterations (or columns). Direct parallelization uses an owner-computes strategy to assign one (or more rows) to each available thread and then execute each column ("iteration") of independent superblocks in a bulk-synchronous fashion. Privatized parallelization uses threadlocal buffers to keep the updates of the superblocks from distinct iterations, and then accumulates these buffers together in parallel at the end. The superblock scheduler also varies the superblocksto-threads assignment from row-to-row or iteration-to-iteration if necessary for the two strategies respectively to achieve better load balance. These two parallel strategies are proposed to obviate the need for locks or atomic operations and to explore sufficient parallel degree for irregular sparse tensors.

However, neither of thee two strategies can obtain good load balance due to the uneven nonzero distribution of the superblocks in rows or columns, even under dynamic thread scheduling policy. This situation becomes worse after applying tensor reordering. We use the ratio of maximum number of nonzeros assigned to a thread to the average number of nonzeros per thread to represent the load balance. Under a random ordering, which might be expected to have good load balance, at the price of poor locality, most ratios are under 1.3. While after Lexi-Order reordering, using the baseline HiCOO superblock scheduler, the ratios become worse on most of tensors and raise to as high as 6.72, although it would be expected to have better locality. (Refer to details in Table 4.) Therefore, we are motivated to try to mitigate this imbalance.

Here is our scheme that tries heuristically to improve the balance ratio. The main idea behind our modified superblock scheduler is illustrated in Figure 8 (b). We simply aggregate superblocks in the same row of the figure, which effectively creates "jagged" partition boundaries (partitions 1 and 2 in Figure 8 (b)), with any remainder superblocks in their own final "tail" partition (the last partition). To decide how many superblocks within each row may be aggregated at a time, there is an additional tuning parameter which is a target maximum number of nonzeros to aggregate at a time. In this example, that target is set to 8 nonzeros, superblocks 0 and 4 are aggregated, as are the pairs [START_REF] Akbudak | Exploiting Locality in Sparse Matrix-Matrix Multiplication on Many-Core Architectures[END_REF][START_REF] Baskaran | Memory-efficient parallel tensor decompositions[END_REF] and [START_REF] Anandkumar | Tensor Decompositions for Learning Latent Variable Models[END_REF][START_REF] Baskaran | Efficient and scalable computations with sparse tensors[END_REF]. Correspondingly, our scheme again uses two parallel strategies: direct and privatized parallelization. However, direct parallelization parallelizes rows in units of partitions, and privatization strategy parallelizes partitions by saving thread-local buffers. Under this aggregation scheme, the two partitions of this example combined represent a more balanced unit to schedule than the original superblock scheduler's units of rows or columns (iterations). The number of nonzeros of our new scheduling unit is around 8 for both row and partition scheduling whereas that of the baseline row and column scheduling units varies in the range of 1 -8. Note that our partitioning scheme does not split any superblock, this may introduce some imbalance. Imbalance also occurs in the last partition for the remainders. We also observe that the number of partitions in Figure 8 (b) could be easily smaller than the number of iterations in Figure 8 (a). This could lead to insufficient parallel degree for privatization strategy.

We record these partition-based partitions as blschr, companied with plschr to record the partitioning positions. In practice, because the superblock sizes are relatively large, there are not many superblocks to schedule on our evaluation tensors. Therefore, it is not expected that these auxiliary metadata will increase the size of the overall HiCOO data structure.

EXPERIMENTS 5.1 Experimental Setup

Platform. We perform experiments on a Linux-based Intel Xeon E5-2698 v3 multicore server platform with 32 physical tributed on two sockets, each with 2.3 GHz frequency. The processor microarchitecture is Haswell, having 32 KiB L1 data cache and 128 GiB memory. The code artifact was written in the C language using OpenMP parallelization, and was compiled using icc 18.0.1.

Dataset. We use the sparse tensors, derived from real-world applications, that appear in Table 2, ordered by decreasing nonzero density separately for third-and fourth-order tensors. Most of these tensors are included in The Formidable Repository of Open Sparse Tensors and Tools (FROSTT) dataset [START_REF] Smith | FROSTT: The Formidable Repository of Open Sparse Tensors and Tools[END_REF]. The darpa (source IP-destination IP-time triples), fb-m, and fb-s (short for "freebasemusic" and "freebase-sampled", entity-entity-relation triples) are from the dataset of HaTen2 [START_REF] Jeon | HaTen2: Billion-scale Tensor Decompositions[END_REF], and choa is built from the electronic health records (EHRs) of pediatric patients at Children's Healthcare of Atlanta (CHOA) [START_REF] Perros | SPARTan: Scalable PARAFAC2 for Large & Sparse Data[END_REF].

Configurations. We report the results under the best configurations for the following parameters for the highest Mttkrp performance with the three formats: (i) the number of reordering iterations, the superblock size L, and the block size B for HiCOO format; (ii) privatization or not for COO; and (iii) tiling or not for CSF. For HiCOO, B = 128 achieves the best results in most cases, and we use five reordering iterations which will be analyzed in § 5.8. All experiments use approximate rank of R = 16. The parallel experiments are run with 32 threads under dynamic scheduling strategy in units of one superblock/partition. We use the total execution time of Mttkrps in all modes for every tensor to calculate the speedup which is the ratio of the total Mttkrp time on a randomly reordered tensor over that using a specific reordering scheme. All the execution time is averaged over five runs.

HiCOO-Mttkrp with Reordering

Figure 9 (a) shows the speedup of the proposed reordering methods on sequential and parallel HiCOO-Mttkrp respectively. Lexi-Order reordering obtains 0.99-4.14× speedup (2.12× on average); while BFS-MCS reordering gets 0.99-1.88× speedup (1.34× on average). Lexi-Order and BFS-MCS do not behave as well on fourthorder tensors as on third-order tensors. Tensor flickr4d is constructed from the same data with flickr, with an extra short mode (refer to Table 2). Lexi-Order obtains 4.14× speedup on flickr while only 3.02× speedup on flickr4d. The same phenomenon is also observed on tensors deli and deli4d. This phenomenon indicates that it is harder to get good data locality on higher-order tensors, which will be justified in Table 3. shows the speedup of the proposed reordering methods on multicore parallel HiCOO-Mttkrp. We set the same superblock size L before and after reordering for a fair comparison. Overall, Lexi-Order results in 0.70-11.81× speedup (2.12× on average) for parallel HiCOO-Mttkrp; while BFS-MCS reordering gets 0.25-1.94× speedup (0.98× on average). Thus, the benefit from reordering using either Lexi-Order or BFS-MCS is generally less than that in the sequential case. This observation may indicate that it is harder to pursue a good load balance after reordering. The 11.81× speedup achieved on tensor flick4d is because Lexi-Order changes the nonzero distribution and hence the optimal superblock size L. Thus, under the same L, reordering gives a big performance improvement. It indicates that automatically tuning the parameters of HiCOO will be very helpful.

HiCOO Parameters

We investigate two critical parameters of HiCOO [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF]: the block ratio (α b ) and the average slice size per tensor block (c b ). Smaller α b and larger c b are favorable for good HiCOO-Mttkrp performance. Table 3 lists the parameter values for all tensors before and after Lexi-Order, the HiCOO-Mttkrp speedup (as shown in Figure 9), and the storage ratio of HiCOO over random ordering. Generally, when α b is reduced and c b is increased using Lexi-Order for a tensor, we see a good performance speedup and storage ratio. (Sequential speedup reflects better by involving less factors than parallel cases.) For the same data in different orders, e.g., flickr4d and flickr, the α b and c b values are the same for random reordering, after Lexi-Order, these values of flickr are better than those of flickr4d. This fact justifies the phenomenon that getting good data locality is harder for higher-order tensors.

Partition-Based Scheduler

Table 4 shows the load balance under three scenarios and the corresponding performance speedup on parallel HiCOO-Mttkrp over random ordering. Load balance is represented by the ratio of maximum number of nonzeros assigned to a thread to the average number of nonzeros per thread 1 , where lower numbers are favored.

The target maximum number of nonzeros of a partition is set to five times of the maximum number of nonzeros per superblock.

The second column shows the balance ratio under a random ordering, which might be expected to have good load balance, at the price of poor locality. Most ratios in this column are under 1.3. The third column shows the balance ratios under Lexi-Order reordering scheme, which increase for most of tensors, some even reach 6.72. Despite these balance ratios, HiCOO-Mttkrp on reordered tensors obtains speedup over random ordering on most of tensors, because of its better data locality. The resulting balance ratios under the partition-based scheduling appear in the fourth column. Only tensors vast, fb-m, fb-s, and nips get better balance ratios compared 1 The same metric is used in [START_REF] Smith | A Medium-Grained Algorithm for Distributed Sparse Tensor Factorization[END_REF].

to reordering-only, but better performance is achieved on about half of tensors including these three. The improvement of the Re-ordered+Balanced scheme over Reordered-only is up to 42%. Two factors may harm the load balance in our new scheduling. One is whether we can get sufficient amount of partitions for parallelization; the other is that the cubical superblocks limits the parallelism for highly irregular tensors. Since the balanced superblock scheme does not always improve performance, our final scheduler selects between these two options. Table 4: Load balance and its effect before and after Lexi-Order tensor reordering and partition-based scheduling. 5.5 Reordering Effect on Other Formats 5.5.1 COO-Mttkrp with Reordering. We show the effect of the two reordering approaches on sequential and parallel COO-Mttkrp from ParTI! [START_REF] Li | ParTI!: A Parallel Tensor Infrastructure for Multicore CPU and GPUs[END_REF] in Figure 10. This COO-Mttkrp is implemented in C and OpenMP parallelized with or without privatization determined for different modes, using the same algorithm with Tensor Toolbox [START_REF] Bader | MATLAB Tensor Toolbox (Version 3.0[END_REF]. For any reordering approach, after doing a BFS-MCS, Lexi-Order, or random reordering on the input tensor, we still sort the tensor in the mode order of 1 The performance is shown by the speedup of each implementation to reordered COO-Mttkrp for every tensor. HiCOO achieves the best performance for most third-order tensors and is better than COO and CSF with tiling for most cases. In these experiments, we found that Splatt without tiling option is even faster than with tiling especially on fourth-order tensors. After reordering, HiCOO is 0.74-7.21× (3.08× on average) faster than COO and 0.25-28.84× (5.29× on average) faster than CSF-tiling. Compared to CSF with no tiling, HiCOO obtains 0.19 -5.30× (1.40× on average) speedup. On flickr, deli, nell1, and enron data, HiCOO behaves worse than both the CSF settings because of the tensors' hypersparsity property. Although the reordering methods have improved the nonzero locality on them, due to their severe load imbalance in Table 4, their performance is still not comparable to CSF format. Overall, as also stated in [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF], HiCOO works the best on short tensor modes and with reasonable sparsity of tensor blocks; while CSF is the complimentary for hypersparse tensors. 3 2 Our tests using one CSF representation (ONEMODE) show quite similar results of reordering effect to the ALLMODE setting. Same for using tiling or no tiling options. 3 HiCOO shows more performance advantage with more CPU cores as illustrated in [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF].

Reordering Methods Comparison

As seen above, Lexi-Order improves performance more than BFS-MCS in most cases. Compared to the reordering method used in Splatt [START_REF] Smith | SPLATT: Efficient and Parallel Sparse Tensor-Matrix Multiplication[END_REF], by setting ALLMODE (identical to the work [START_REF] Smith | SPLATT: Efficient and Parallel Sparse Tensor-Matrix Multiplication[END_REF]) to CSF-Mttkrp, BFS-MCS gets 1.04, 1.64, and 1.61× speedups on tensors nell2, nell1, and deli respectively, and Lexi-Order obtains 1.04, 1.70, and 2.24× speedups. By contrast, the speedups using graph partitioning [START_REF] Smith | SPLATT: Efficient and Parallel Sparse Tensor-Matrix Multiplication[END_REF] on these three tensors are 1.06, 1.11, and 1.19× and 1.06, 1.12, and 1.24× by using hypergraph partitioning [START_REF] Smith | SPLATT: Efficient and Parallel Sparse Tensor-Matrix Multiplication[END_REF] respectively. Our BFS-MCS and Lexi-Order schemes both outperform graph and hypergraph partitionings [START_REF] Smith | SPLATT: Efficient and Parallel Sparse Tensor-Matrix Multiplication[END_REF].

The available methods in the state-of-the-art are based on graph and hypergraph partitioning. Partitioning is a successful approach when the number of partitions is known, while the number of blocks in HiCOO is unknown ahead of time. In our case, the partitions should also be ordered for better cache reuse. Additionally, partitioners are less effective for tensors than their usual applications for matrices, as some dimensions could be very short (creating very high degree vertices). Thus, our proposed ordering based methods deliver better results.

Reordering Overhead

Figure 13 shows the parallel speedup of Lexi-Order reordering and the HiCOO conversion steps respectively. Parallel Lexi-Order gains 3.61× speedup on average while parallel HiCOO conversion gets 5.84× average speedup on 32 physical cores. There is still room to further accelerate these parallel algorithms, which will be one of our future work. Figure 14 shows the overhead of Lexi-Order with 5 iterations. We specifically report the ratio of parallel Lexi-Order reordering time to parallel HiCOO construction time. (That is, we show how much more expensive it might be to reorder than to construct the HiCOO representation in the first place.) The results lie in the range of 1.97 to 12.91×. However, we could use less iterations for Lexi-Order to compromise some Mttkrp performance, which will be shown in Figure 15 (a). Thus, the reordering overhead can be further reduced.

Effect of the Number of Iterations in Lexi-Order

Since Lexi-Order improves the ordering iteratively, we evaluate the effect of the number of iterations on HiCOO-Mttkrp performance.

The results appear in Figure 15 (a), which is normalized to the runtime of 10 iterations. Mttkrp on most tensors does not vary a lot by setting different number of iterations, except vast, nell2, uber, and nips. We use 5 iterations to get good Mttkrp performance similar to that of 10 iterations, with about half of the overhead (shown in Figure 15 (b)). But 3 or fewer iterations will get an acceptable performance when users care much about the pre-processing time.

CPD Application

Figure 16 shows the Lexi-Order effect on the parallel HiCOO-CPD algorithm using alternating least squares [START_REF] Kolda | Tensor Decompositions and Applications[END_REF]. Similar to the numbers in Figure 9 (b), CPD after reordering achieves 0.65-10.44× speedup (1.81× on average). Reordering helps us to enhance the performance of a whole tensor application.

RELATED WORK

Plenty recent research studied the optimization of tensor algorithms [6, 8-10, 12, 13, 16, 19, 20, 23, 25, 27, 29, 38, 40]. Our work emphasizes on reordering schemes to get better nonzero structure. Various reordering methods have been considered for matrix algebra [START_REF] Akbudak | Exploiting Locality in Sparse Matrix-Matrix Multiplication on Many-Core Architectures[END_REF][START_REF] Karypis | A Parallel Algorithm for Multilevel Graph Partitioning and Sparse Matrix Ordering[END_REF][START_REF] Mellor-Crummey | Improving Memory Hierarchy Performance for Irregular Applications Using Data and Computation Reorderings[END_REF][START_REF] Pichel | Optimization of sparse matrix-vector multiplication using reordering techniques on GPUs[END_REF][START_REF] Pothen | Computing the Block Triangular Form of a Sparse Matrix[END_REF][START_REF] Yzelman | High-Level Strategies for Parallel Shared-Memory Sparse Matrix-Vector Multiplication[END_REF]. For tensor operations, Smith et al. [START_REF] Smith | SPLATT: Efficient and Parallel Sparse Tensor-Matrix Multiplication[END_REF] proposed two reordering methods based on partitioning: the first is the partitioning of a graph that models the interactions between slices; the second is the partitioning of a hypergraph that models the memory access of the Mttkrp operation. A tensor reordering is induced such that the vertices in the same partition are set to consecutive labels. Hypergraphs have also been used to efficiently parallelize Mttkrp in distributed memory systems [START_REF] Kaya | Scalable Sparse Tensor Decompositions in Distributed Memory Systems[END_REF]. These two hypergraphs are constructed differently from our BFS-MCS, they use fibers and tasks as vertices respectively, while our hypergraph takes tensor indices as vertices. Besides, we propose matLexiOrder, an alternative doubly lexical ordering method for sparse matrices, and design Lexi-Order for sparse tensors based on it.

CONCLUSION

Motivated by an interest in further improving the HiCOO implementation, we investigated the problem of reordering a tensor to improve block density for tensor computations. Inspired by algorithms for sparse matrices, we proposed two heuristics, BFS-MCS and Lexi-Order. BFS-MCS is based on the maximum cardinality search principle; Lexi-Order is based on matLexiOrder, which we propose here as an alternative to existing doubly lexical ordering method for sparse matrices. matLexiOrder has near linear runtime per tensor dimension, rendering the overall tensor reordering method Lexi-Order practical. Lastly, we improve the superblock scheduling strategy of HiCOO [START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF] with a partition-based scheme, which eases the side effect of increased load imbalance that can occur after reordering, thereby improving Mttkrp performance for several sample inputs. Overall, BFS-MCS is asymptotically faster while Lexi-Order is more effective. Lexi-Order obtains performance speedup on HiCOO-Mttkrp of 2.12× on average for both sequential and multicore parallel implementations. For COO and CSF-Mttkrps, the average sequential and parallel speedups are 1.79×, 1.21× and 1.50×, 1.27×, respectively. Users can control a parameter of Lexi-Order to reduce the reordering overhead without much performance drop and also tune the HiCOO parameters for the reordered tensors to pursue the highest Mttkrp performance. As future work, we plan to pursue the automatic performance tuning for different sparse tensor formats, HiCOO parameters, the tradeoff of reordering methods and parallel strategy parameters, and improve the runtime of Lexi-Order. Our proposed heuristics set a basis for improving Mttkrp performance through tensor reordering. Furthermore, we anticipate the development of additional heuristics for the general tensor reordering problem. Additionally, we expect reordered tensors will prove useful for any tensor operation in which the memory access behavior depends on the locality of tensor indices, as in Mttkrp. Examples include Tensor-times-matrix and tensortimes-vector operations in Tucker decompositions and higher-order power methods, respectively. Our approaches will be especially useful as local improvements for intra-node data locality in large HPC environment. contract FA8750-18-2-0108, under the DARPA MTO Software Defined Hardware program, and the Laboratory Directed Research and Development program at Sandia National Laboratories under contract DE-NA-0003525. Disclaimer: The views, opinions, and/or findings contained in this document are those solely of the author(s) and should not be interpreted as representing the official views or policies of any of its funding sources.

  Khatri-Rao product between two matrices a • b Outer product between two vectors Input parameters N Tensor order I, J, K, In Tensor mode sizes M #Nonzeros of the input tensor X R Approximate tensor rank (usually a small value) Word sizes βint Bit-length of an integer β long Bit-length of a long integer β byte Bit-length of a byte or character β float Bit-length of a single-precision floating point value Data structures inds Indices of COO, βint bits cinds Indices of CSF, βint bits cptrs Pointers of CSF, βint or β long bits einds Element indices of HiCOO, β byte bits binds Block indices of HiCOO, βint bits bptr Block pointers of HiCOO, β long bits lptr Superblock pointers of HiCOO, β long bits lschr Superblock scheduler of HiCOO, βint bits blschr Partition-based superblock scheduler of HiCOO, βint bits plschr Partition pointers of blschr, βint bits val Nonzero value array of COO, CSF, HiCOO, β float bits Performance parameters perm n Permutation for a given mode n L Tensor superblock size B Tensor block size, B ≪ L α b Block ratio, α b = n b M , n b is #Nonzero tensor blocks M b Geometric mean of #Nonzeros per tensor block c b Average slice size per tensor block, c b =

2. 2

 2 .1 COO. The coordinate (COO) format is the simplest yet arguably most popular format by far. It stores each nonzero value along with all of its position indices, shown in Figure 2 (a). i, j, k are indices (inds) of the nonzeros stored in the val array. Generally, 32-bit (β int ) integers are large enough to represent tensor indices.

Figure 2 :

 2 Figure2: COO, CSF, and HiCOO formats for an example third-order tensor (This example is originally from[START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF]). 2.2.2 CSF. Compressed Sparse Fiber (CSF) is a hierarchical, fibercentric format that effectively generalizes the CSR matrix format to tensors. An example of its representation appears in Figure2 (b). Conceptually, CSF organizes nonzero indices into trees (a forest). Each level corresponds to a tensor mode, and each nonzero is a path from a root to a leaf. The indices of CSF are stored in cinds with the pointers stored in cptrs to indicate the locations of nonzeros at the next level. Since cptrs needs to show the range up to M, we use β int or β long accordingly for differently sized-tensors .

Figure 3 :

 3 Figure 3: A hypergraph example of a sparse tensor.

▷ 2 i

 2 Entry m 1 > entry m 2 17: return 0; ▷ Entry m 1 = entry m

Figure 5 :

 5 Figure 5: The steps of Lexi-Order illustrated for mode 1.To order tensors, we propose the Lexi-Order function as an extension of matLexiOrder. The basic idea of Lexi-Order is to determine the permutation of each tensor mode independently, while considering the order in other modes fixed. Lexi-Order sets the indices of the mode to be ordered as the columns of a matrix, the other indices as the rows and sorts the columns as described for matrices (with the order preserving partition refinement method). The precise algorithm appears in Algorithm 2, which we also illustrate in Figure5when applied to mode 1. Given a mode n, Lexi-Order first builds a matricized tensor in Compressed Sparse Row (CSR) sparse matrix format by a call to quickSort with the comparison function coordCmp and then by partitioning the nonzeros into the row segments (Lines 3-7). This comparison function coordCmp

   to illustrate its nonzero distribution, which could seamlessly call orderlyRefine function. Apart from the quickSort, the other parts of Lexi-Order are of linear time complexity (linear in terms of tensor storage). We use OpenMP Tasks to parallelize quickSort to accelerate Lexi-Order.

Figure 6 :

 6 Figure 6: Comparison of HiCOO representations before and after Lexi-Order.Like BFS-MCS approach, Lexi-Order also finds the permutations for N modes of an N th-order sparse tensor. Figure6(a) illustrates the effect of Lexi-Order on an example 4 × 4 × 3 sparse tensor. The original tensor is converted to a HiCOO representation with block size B = 2 consisting of 5 blocks, with maximum 2 nonzeros per block. After reordering with Lexi-Order, the new HiCOO has 3 nonzero blocks with up to 4 nonzeros per block. Thus, the blocks are denser, which should exhibit better locality behavior. However, this reordering scheme is heuristic. Figure6 (b)shows HiCOO representation after reordering of the tensor shown in Figure2. The number of blocks is unchanged (4), although the maximum number of nonzeros per block increases to 4. For this tensor, Lexi-Order may not show a big advantage.

  Figure 6: Comparison of HiCOO representations before and after Lexi-Order.Like BFS-MCS approach, Lexi-Order also finds the permutations for N modes of an N th-order sparse tensor. Figure6(a) illustrates the effect of Lexi-Order on an example 4 × 4 × 3 sparse tensor. The original tensor is converted to a HiCOO representation with block size B = 2 consisting of 5 blocks, with maximum 2 nonzeros per block. After reordering with Lexi-Order, the new HiCOO has 3 nonzero blocks with up to 4 nonzeros per block. Thus, the blocks are denser, which should exhibit better locality behavior. However, this reordering scheme is heuristic. Figure6 (b)shows HiCOO representation after reordering of the tensor shown in Figure2. The number of blocks is unchanged (4), although the maximum number of nonzeros per block increases to 4. For this tensor, Lexi-Order may not show a big advantage.

Figure 7 :

 7 Figure 7: The conversion between COO and HiCOO formats for an example third-order tensor. HiCOO uses 2×2×2 blocks (B = 2) with word sizes marked above.

Figure 8 :

 8 Figure8: Baseline[START_REF] Li | HiCOO: Hierarchical Storage of Sparse Tensors[END_REF] and partition-based superblock scheduling tables for mode-1 Mttkrp. Each box represents a superblock identified by its number, along with its number of nonzeros in parentheses.

vast 3 165K × 11K × 2 26M 6 . 9 × 10 - 3 nell2 3 12K 4 6K × 24 × 77 × 32 5M 1 . 5 × 10 - 2 uber 4 183 × 24 × 1140 × 1717 3M 3 . 9 × 10 - 4 nips 4 2K × 3K × 14K × 17 3M 1 . 8 × 10 - 6 enron 4 6K× 4 320K× 4 533K×

 32691033424321510242417173910441718106444 6K × 244K × 1K 54M 5.5 × 10 -9 flickr4d 28M × 2M × 731 113M 1.1 × 10 -14 deli4d 17M × 3M × 1K 140M 4.3 × 10 -15

Figure 9 :

 9 Figure 9: Reordered HiCOO-Mttkrp speedup over a random ordering implementation.Figure9(b) shows the speedup of the proposed reordering methods on multicore parallel HiCOO-Mttkrp. We set the same superblock size L before and after reordering for a fair comparison. Overall, Lexi-Order results in 0.70-11.81× speedup (2.12× on average) for parallel HiCOO-Mttkrp; while BFS-MCS reordering gets 0.25-1.94× speedup (0.98× on average). Thus, the benefit from reordering using either Lexi-Order or BFS-MCS is generally less

Figure 9 (

 9 Figure 9: Reordered HiCOO-Mttkrp speedup over a random ordering implementation.Figure9(b) shows the speedup of the proposed reordering methods on multicore parallel HiCOO-Mttkrp. We set the same superblock size L before and after reordering for a fair comparison. Overall, Lexi-Order results in 0.70-11.81× speedup (2.12× on average) for parallel HiCOO-Mttkrp; while BFS-MCS reordering gets 0.25-1.94× speedup (0.98× on average). Thus, the benefit from reordering using either Lexi-Order or BFS-MCS is generally less

Figure 10 :

 10 Figure 10: Reordered COO-Mttkrp speedup over a random reordering implementation.

Figure 11 :

 11 Figure 11: Reordered CSF-Mttkrp speedup over a random reordering implementation.

Figure 12 :

 12 Figure 12: Format comparison with Lexi-Order on parallel Mttkrps.

Figure 13 :

 13 Figure 13: Parallel conversion and Lexi-Order speedup. Times of HiCOO Constrution 0 3 6 9 12 15

Figure 14 :

 14 Figure 14: The reordering overhead over HiCOO construction.

Figure 15 :

 15 Figure 15: The performance and overhead of different numbers of iterations.

Figure 16 :

 16 Figure 16: CPD performance on reordered tensors.

Table 2 :

 2 Description of sparse tensors.

	Tensors	Order	Dimensions	#Nnzs	Density

Table 3 :

 3 HiCOO parameters change before and after Lexi-Order reordering.

	Tensors	Random reordering α b c b	Lexi-Order α b c b	Speedup seq omp	Storage ratio
	vast	0.004	1.758	0.004	1.562	1.01	1.03	0.999
	nell2	0.020	0.314	0.008	0.074	1.04	1.04	0.966
	choa	0.089	0.057	0.016	0.056	1.16	1.07	0.833
	darpa	0.796	0.009	0.018	0.113	3.74	1.50	0.322
	fb-m	0.985	0.008	0.086	0.021	3.84	1.21	0.335
	fb-s	0.982	0.008	0.099	0.020	3.90	1.23	0.336
	flickr	0.999	0.008	0.097	0.025	4.14	3.66	0.277
	deli	0.988	0.008	0.501	0.010	2.24	0.83	0.634
	nell1	0.998	0.008	0.744	0.009	1.70	0.70	0.812
	crime	0.001	37.702	0.001	8.978	0.99	1.42	1.000
	uber	0.041	0.469	0.011	0.270	1.00	0.78	0.838
	nips	0.016	0.434	0.004	0.435	1.03	1.34	0.921
	enron	0.290	0.017	0.045	0.030	1.25	1.36	0.573
	flickr4d	0.999	0.008	0.148	0.020	3.02	11.81	0.214
	deli4d	0.998	0.008	0.596	0.010	1.76	1.26	0.697

  Mttkrp with Reordering. We show the effect of the two reordering approaches on sequential and parallel CSF-Mttkrp from Splatt v1.1.1[START_REF] Smith | SPLATT: The Surprisingly ParalleL spArse Tensor Toolkit[END_REF] in Figure11. CSF-Mttkrp is set to use all CSF representations (ALLMODE) for Mttkrps in all modes and with tiling option on.2 Lexi-Order improves sequential CSF-Mttkrp perfor-Mttkrp on average. While BFS-MCS is always helpful in the sequential case, Lexi-Order is not helpful on only one tensor crime. The improvements achieved by the two reordering approaches for CSF are less than those for HiCOO and COO formats. We conclude that both reordering methods are helpful for CSF-Mttkrp, but in a lesser extend than for HiCOO and COO based Mttkrp.

	5.5.2 CSF-mance by 0.65-2.33× (1.50× on average) and accelerates parallel
	CSF-Mttkrp by 0.86-1.88× (1.27× on average). BFS-MCS improves
	sequential CSF-Mttkrp performance by 1.00-1.86× (1.22× on aver-
	age) and accelerates parallel CSF-Mttkrp by 0.59-1.36× (1.04× on
	average). Both ordering approaches improves the performance of
	CSF-Speedup	0.5 1.0 1.5 2.5 2.0														Lexi-Order BFS-MCS
		0.0	vast	nell2	choa	darpa	fb-m	fb-s	ickr	deli	nell1	crime	uber	nips	enron	deli4d ickr4d
					3-D Tensors						4-D Tensors
																≻ • • • ≻ N . Observe that Lexi-
																Order improves sequential COO-Mttkrp performance by 1.00-
																4.29× (1.79× on average), while BFS-MCS gets 0.95-1.27× (1.10×
																on average). Also, Lexi-Order improves parallel COO-Mttkrp
																performance by 1.01-1.48× (1.21× on average), while BFS-MCS
																improves by 0.70-1.68× (1.11× on average). Note that Lexi-Order
																improves the performance of COO-Mttkrp for all tensors, both in
																sequential (79%) and parallel (21%). We conclude that this ordering
																is always helpful for COO-Mttkrp, while the improvements being
																less than what we saw for HiCOO-Mttkrp.

ACKNOWLEDGMENTS

This research was partially funded by the US Department of Energy, Office for Advanced Scientific Computing (ASCR) under Award No. 66150: "CENATE: The Center for Advanced Technology Evaluation". Pacific Northwest National Laboratory (PNNL) is a multiprogram national laboratory operated for DOE by Battelle Memorial Institute under Contract DE-AC05-76RL01830. This research was also funded in part under Defense Advanced Research Projects Agency (DARPA)


P Albrodt

M Niemeyer

P Crump

J Hamperl

F Moron

P Georges

And G Lucas-Leclin
email: gaelle.lucas-leclin@institutoptique.fr

Coherent beam combining of high power quasi continuous wave tapered amplifiers

We demonstrate coherent beam combining of four high brightness tapered amplifiers in pulsed, quasi continuous wave (QCW) operation, seeded by a 976 nm laser diode. The maximum power of 22.7 W was achieved with > 64 % combining efficiency in a close to diffraction limited beam. We discuss turn-on dynamics of tapered amplifiers operated in pulsed mode in detail.

Introduction

The power and brightness of laser diodes are continuously increasing and as the cost per watt is decreasing at the same time, they become more attractive for industrial application requiring high brightness [START_REF]High Power Diode Lasers: Technology and Applications[END_REF]. It is however necessary to make use of beam combining technologies in order to meet the power requirements in this regard [START_REF] Fan | Laser beam combining for high-power, high-radiance sources[END_REF]. Especially spectral beam combining is a widely used approach for development of high brightness diode laser systems [START_REF] Strohmaier | High-power, high-brightness direct-diode lasers[END_REF][START_REF] Daneu | Spectral beam combining of a broad-stripe diode laser array in an external cavity[END_REF], the ultimate limitation is however the power available within a certain spectral linewidth. Coherent beam combining (CBC) relies on constructive interference of several separate beams and allows to scale the power available within a very narrow spectral bandwidth in a beam with nearly diffraction limited beam quality. CBC with diode lasers has been demonstrated in many different external cavity configurations [START_REF] Bartelt-Berger | Power-scalable system of phase-locked single-mode diode lasers[END_REF][START_REF] Liu | Coherent addition of high power laser diode array with a V-shape external Talbot cavity[END_REF][START_REF] Paboeuf | Narrowline coherently combined tapered laser diodes in a Talbot external cavity with a volume Bragg grating[END_REF][START_REF] Montoya | External cavity beam combining of 21 semiconductor lasers using SPGD[END_REF][START_REF] Schimmel | Rear-side resonator architecture for the passive coherent combining of high-brightness laser diodes[END_REF] and in master oscillator power amplifier (MOPA) configurations [START_REF] Redmond | Active coherent beam combining of diode lasers[END_REF][START_REF] Creedon | High efficiency coherent beam combining of semiconductor optical amplifiers[END_REF]. The highest combined power has been achieved in a MOPA setup using large arrays of single mode amplifiers [START_REF] Creedon | High efficiency coherent beam combining of semiconductor optical amplifiers[END_REF]. Using a few high brightness tapered amplifiers is a promising approach to reduce the complexity of CBC setups and yields in a higher power per element [START_REF] Schimmel | Coherent beam combining architectures for high power tapered laser arrays[END_REF][START_REF] Albrodt | Coherent combining of high brightness tapered amplifiers for efficient non-linear conversion[END_REF]. Operating those amplifiers in pulsed, quasi continuous wave (QCW) mode, allows even higher powers per amplifier to be reached as the devices are no longer limited by thermal rollover. However CBC of diode lasers has commonly been demonstrated in continuous wave (CW) operation and there are only few demonstrations of CBC of pulsed diode lasers [START_REF] Osinski | Phased array of high power, coherent, monolithic flared amplifier master oscillator power amplifiers[END_REF]. Pulsed (QCW) semiconductor amplifiers are not in steady state but exhibit dynamic temperature changes. This consequently leads to a dynamic change of the accumulated phase and of the beam propagation parameters, which makes CBC with pulsed amplifiers more challenging. Additionally, high brightness pulsed sources with high spectral purity are in demand as NIR pump sources of solid state lasers and for nonlinear frequency conversion to address medical applications of modulated high power visible laser sources [15][16][START_REF] Christensen | Deep modulation of second-harmonic light by wavelength detuning of a laser diode[END_REF]. Furthermore CBC of pulsed tapered amplifiers may improve the performance of diode-laser-based portable differential absorption lidar systems used for atmospheric measurements [START_REF] Nehrir | Eye-Safe Diode-Laser-Based Micropulse Differential Absorption Lidar (DIAL) for Water Vapor Profiling in the Lower Troposphere[END_REF][START_REF] Nehrir | Micropulse water vapor differential absorption lidar: transmitter design and performance[END_REF].

We use a master oscillator power amplifier (MOPA) configuration with four parallel tapered amplifiers operated in QCW mode. Their coherent superposition is achieved in a filled aperture approach using standard beamsplitters as the combiner elements. The experimental setup is described in section 2. The performance of the tapered amplifiers is characterized in section 3 especially regarding the dynamic change in phase as the devices are turned on. We observe and analyze effects associated with the operation of the pulsed amplifiers at high bias, which induce strong peak power fluctuations and deteriorate the seed line stability. Finally we demonstrate CBC of four high brightness tapered amplifiers in QCW mode and compare the results to operation in CW mode (see section 4).

Experimental setup

We used 5 mm long amplifiers consisting of a 1 mm long 4 µm wide ridge waveguide (RW) followed by a 4 mm long tapered section (TP, 6° taper angle). The facets of the devices were passivated and anti-reflection coated (R 1 = 0.01 % for RW facet and R 2 = 0.05 % for TP facet) [20] The epitaxial design of the devices was the so called ELOD2 design [START_REF] Crump | Development of highpower diode lasers with beam parameter product below 2 mm mrad within the BRIDLE project[END_REF], having an extremely low vertical divergence, that has proven to maintain good beam quality at high drive currents in tapered laser format [START_REF] Crump | Development of highpower diode lasers with beam parameter product below 2 mm mrad within the BRIDLE project[END_REF], with good results also in related designs in truncated tapered amplifier format [START_REF] Wang | 17-W Near-Diffraction-Limited 970-nm Output From a Tapered Semiconductor Optical Amplifier[END_REF]. We used individually mounted emitters instead of an array of amplifiers to achieve a simple modular setup and benefit from the efficient heat removal in each heatsink. The devices were mounted p-side up on CuW heat sinks onto Cmounts with separate electrical contacts for the currents I rw into the RW and I tp into the tapered section of each device respectively. An additional CuW heat spreader was also attached to the p-side of the tapered amplifier section for improved cooling. The C-mounts were fixed on individual temperature-stabilized mechanical mounts. The current into the RW was not pulsed as it has low impact on the total heat load of the device (< 1 W). We used small variations (+/-50 mA) of the continuous current into the ridge section of each amplifier for phase-control [START_REF] Albrodt | Coherent combining of high brightness tapered lasers in master oscillator power amplifier configuration[END_REF]. We used Arroyo 4320-QCW laser diode drivers for the current into the tapered section of each device. The pulse durations were in the order of a few milliseconds (typically Δt QCW = 2 ms) and a duty cycle of 10 %. The amplifiers were arranged in a multiarm-interferometer as shown in Fig. 1 similar to the experimental setup used in previous work [START_REF] Albrodt | Coherent combining of high brightness tapered amplifiers for efficient non-linear conversion[END_REF]. The seed source was a 100 mW CW narrow linewidth (< 20 MHz) DFB laser at 976 nm. After optical isolation the beam was split using non-polarizing 50:50 beamsplitters. The seed power incident for each amplifier was 20 mW. The amplified beams were collimated in slow axis (SA) and fast axis (FA) individually and the astigmatism was corrected depending on the current into the tapered section. We used also 50:50 plate beamsplitters for recombining the beams in a filled aperture appr measured wit phase correcti

Investig

Optical p

The QCW op typically belo thermal rollov is clearly vis thermal rollov the case for Q 0.7 W/A resu for moderate imaged beam into the tapere mode were m time, the devi tapered ampli energy and lo side lobes gra content in the differences in currents but amplifiers are at 13 A. The the central lob quality also d atism that nge in the lementary ng that the beam is focusing rapidly at the beginning of the pulse and then converges to a stable shape after about 500 µs. The measured 2 nd moment beam diameter decreases from 1100 µm to 875 µm, which corresponds to a displacement of the waist in SA by about 80 % of the corresponding Rayleigh length at a current of 13 A into the tapered section. The spatial dynamics were similar for all four amplifiers in our setup. In order to ensure a good spatial overlap of the beams during the full pulse duration, it is therefore necessary to operate the different amplifiers at similar bias conditions in order to optimize the coherent beam combining efficiency.

Temperature -driven piston phase drift of pulsed amplifiers

The rise in temperature as the device is turned on does also increases the effective optical path length in the amplifier and leads to a drift of the accumulated phase in each amplifier. The piston phase drift is approximately proportional to the change in temperature. We measured the phase drift during the QCW pulses by overlapping the output beam from an amplifier with a reference beam as described in detail in [START_REF] Albrodt | Recent progress in brightness scaling by coherent beam combining of tapered amplifiers for efficient high power frequency doubling[END_REF]. Figure 3(a) shows the measurements for one amplifier at different currents into the tapered section. The amplitude of the piston phase drift increases with the current, and the experimentally measured phase drift was reproduced accurately at all bias points using the following fit function (Eq. 1), as seen in Fig. 3(a):

(

The time constants , τ 1 2 of the two exponential functions were found by fitting to be 0.1 ms and 1 ms respectively and these were independent of the current into the tapered section I tp .

The amplitudes of the exponential contributions , are proportional to I tp and roughly equal for the two exponential contributions. The linear term C t × describes a linear increase in the average temperature of the heatsink during the pulse and is also proportional to I tp . The amplitude of the overall change in piston is large, reaching 30 π at 13 A for 5 ms long pulses. The dynamic variation in phase is however highly reproducible, consistent between samples and we propose it to be defined by the geometry and efficiency of the amplifiers, the heatsink and the rise-time of the current driver. We thus achieved a very good overlap of the piston phase drifts in all four amplifiers of the interferometer as shown in Fig. 3(b). We used small adjustments (+/-0.5 A) of the pulsed current into the taper sections in order to achieve a piston phase drift of the same magnitude for all four amplifiers. The inset in Fig. 3(b) shows a zoom into the same data, where one can see that the relative differences in phase between the four amplifiers are in the range of π / 10 only.

Power st

We measured oscillations in with time an differences be significant bu the oscillation Devices with oscillations a currents (e.g. e narrow oth device w-finesse scanning Fabry-Perot cavity whose resonance peaks shift during the current pulse following the dynamic phase drift discussed earlier. Consequently, the period of the oscillations increases during the pulse (cf. Fig. 3 and Eq. 1). Furthermore the cavity finesse is enhanced with the amplifier gain, leading to an increase of the modulation depth. The slight disparities in the behavior of the different devices are therefore attributed to small differences in the coating reflectivity, coupling efficiency and gain in the device.

It would be desirable for further studies to use amplifiers with the best possible ARcoating on both facets in order to push the onset of these oscillations to even higher powers, which was not possible in the framework of this work. Besides the degradation of the beam quality at high bias, the poor pulse stability is another limitation for power scaling by QCW operation in QCW mode. Especially regarding coherent superposition, it is important to have a reasonable stability of the power within the pulse. We therefore limited the current in the tapered section to below 14.5 A for the subsequent CBC experiments, as the pulse stability was insufficient at higher currents.

Optical feedback from overdriven tapered amplifiers

It is well known that optical feedback from a distant reflector may perturb the lasing properties of a DFB laser, such as frequency, spectral linewidth and power [START_REF] Tkach | Regimes of feedback effects in 1.5-µm distributed feedback lasers[END_REF]. In particular, weak time-dependent feedback may lead to a rich variety of effects including frequency hopping. We started our experiment with one double-stage Faraday isolator (-55 dB optical isolation) between the seed DFB laser diode and the CBC interferometer to limit reflection effects. A small part of the amplified beam gets reflected at the output facet (R 2 = 0.5 %) of the amplifier and propagates backwards through the setup towards the seed laser. This optical feedback is spectrally identical to the seed laser with only minor ASE contribution. As several parallel overdriven tapered amplifiers lead to high levels of absolute optical feedback, it appears that this isolation was not sufficient to prevent perturbations of the seed laser. Indeed we measured the seed laser frequency using a Fabry-Perot interferometer (Thorlabs SA210-8B, FSR = 10 GHz, finesse > 150). Figure 5(a) shows the measured seed laser frequency perturbations (blue line) induced by optical feedback from one tapered amplifier during a 5 ms-long QCW pulse (black line). The feedback leads to a frequency shift and a frequency hopping up to 80 MHz. The decreasing period of the observed frequency hopping appears to be linked to the phase dynamics discussed in section 3.2.

In the presence of an optical feedback, the frequency shift of a single-mode laser is described by [START_REF] Tkach | Regimes of feedback effects in 1.5-µm distributed feedback lasers[END_REF] (2) with the coupling coefficient ( ) ,

s e s s R R R - κ = τ 1 (3)
where s R is the reflectivity of the laser output facet, e R is the effective reflectivity of the external optical setup, s τ is the roundtrip time in the laser, ( ) e t τ is the time dependent round trip time of the feedback and α is the linewidth enhancement factor. Frequency hopping occurs for higher levels of the product e κτ when there are multiple solutions for Eq. 2. We calculated the frequency change using the measurement of the phase drift at I tp = 13 A given in Fig. 3(a) for defining ( ) e t τ with 0 = 10 ns related to the distance from the seed laser to the amplifier (about 1.5 m in this laboratory scale experiment). We compared in Figs. 5(a) and 5(b) the experimental evolution of the seed frequency during the pulse with a numerical resolution of Eq. 2. Each black dot in Fig. 5(b) is one calculated solution to this equation. We can easily identify regimes with multiple solutions where frequency hopping is likely to occur. The red line represents the expected evolution of the seed laser frequency, assuming that the frequency hopping occurs when the change in the laser frequency has reached the local minimum. Fig. 5(a) shows that the measurement (blue) and simulation (red) are in good agreement during the whole pulse. The value for the external reflection of the extended cavity = 9 × 10 used in these simulations may appear low but corresponds to the small part of the optical feedback from the amplifier that couples back into the seed laser. Indeed feedback from a distant reflector is critical even at very low level (see [START_REF] Tkach | Regimes of feedback effects in 1.5-µm distributed feedback lasers[END_REF]). In comparison feedback from near passive optical elements (such as L 1 in Fig 1 .) may have a similar level, but would be less critical as it would be induced at short distance and would result in a stable piston phase. Equation ( 2) has multiple solutions, resulting in frequency hopping, if

e + α κτ > 2 1
1.

With our simulation parameters, the product e + α κτ we estimate that the optical isolation in our experimental conditions should be better than -62.7 dB to avoid those frequency hops. Furthermore the maximum change in frequency would be below 10 MHz for an optical isolation better than -73 dB considering that ∝ . We repeated our experiment after increasing the optical isolation from -55 dB to -75 dB by using a third stage of optical isolation and did no longer measure any effect on the laser frequency during the QCW pulses. We conclude that the requirements regarding the optical isolation of a narrow linewidth seed laser can be extremely high in CBC architectures, especially when the amplifiers are overdriven and generate therefore a high level of feedback. The sensitivity to optical feedback could be reduced in an integrated optical setup (to minimize e τ ), by increasing the reflectivity s R of the seed laser and by further improvement of the AR coatings of the tapered amplifiers (to reduce e R ). Small perturbations of the seed laser frequency in the range of 80 MHz (see Fig. 5(a)) are not necessarily a problem for successful CBC, as this corresponds still to a coherence length in the range of 1 m, but is critical for applications requiring a stable frequency.

Coherent beam combining in QCW operation

The phase matching between the four amplified beams in the coherent beam combining experiment was achieved by a sequential hill-climbing algorithm with stepwise adjustment of the current into the RW from pulse to pulse. The error signal was measured with a fast photodiode (see Fig. 1) measuring the power at the CBC output. The currents into the amplifiers were optimized separately: the current into the RW of A2 was optimized first from pulse to pulse, then the current into A3 similarly and finally the current into A4. A variation of ± 35 mA corresponded to a phase shift of ± π and we used steps of 0.75 mA for the hill- climbing algorithm. Phase control with this simple sequential pulse-to-pulse algorithm was sufficient as the phase noise in the interferometer was dominated by very low frequencies (<1 Hz). The feedback loop was triggered to maximize the power of the combined beam after 50 % of the pulse duration as the pulse stability was better in the second half of the pulse (cf. Fig. 4). The pulse duration was adjustable from 1 to 10 ms and the duty cycle from 0.01 to 10 % (we typically used 2 ms pulses and a duty cycle of 10 %). The achieved coherently combined power P CBC and the corresponding combining efficiency, defined as the ratio of the combined power over the sum of the individual powers ( CBC

CBC i P P η =  4 1
), are summarized in Table 1. Power and combining efficiency for operation in continuous wave mode are also given for comparison. The current values I tp into the tapered section correspond to the average of the four currents, which were slightly adjusted in order to optimize the phase matching during the pulse (see Fig. 3 (b)). For each measurement, the SA cylindrical lenses (L 4 see Fig. 1) were positioned to compensate for the amplifier astigmatism and optimize the overlap of the beams during the current pulse despite the wandering SA waist position mentioned in section 3.1.

At the maximum operating current I tp = 14.5 A per amplifier the QCW output power of the coherently combined output beam reached 22.7 W. This is the highest power reported for CBC of tapered lasers or amplifiers. In contrast we measured 16.9 W as the highest power in CW, which is about 25 % lower than the maximum QCW power. The combining efficiency at the highest bias was > 64% limited by the degradation of the beam quality at high currents.

The combining efficiency at moderate currents was measured to be > 70 % and almost identical for QCW and CW operation. The combining efficiency of the coherent superposition of only two amplifiers was slightly higher, in the range of 79-84 %, as each combining step adds combining losses. A large part of the combining losses were attributed to a spatial clean-up of the beam. The beam from each amplifier has a central lobe power content in the order of 75 % (cf. Fig. 2), which is significantly lower than the central lobe power content of the coherently combined beam of 92 % as shown in Fig. 6. The M 2 beam propagation factor of the combined beam was below 1.3 in both directions considering the 1/e 2 beam diameters. The more accurate and standard definition of the M 2 factor, using the 2 nd moment beam diameters, resulted in 1.3 in FA and 5 in SA : this is below the beam propagation factor of the individual amplifiers M 2 > higher central of tapered las phase of the h those individu efficiency con mismatches in were identifie shows the mea of the pulses o he peak to vall 14.5 A, mainl with (feedback ical isolation ( fects perturb th ned beam was individual amp lifiers [START_REF] Schimmel | Rear-side resonator architecture for the passive coherent combining of high-brightness laser diodes[END_REF][START_REF] Albrodt | Coherent combining of high brightness tapered amplifiers for efficient non-linear conversion[END_REF] be achieved even at high currents. Power instabilities in the QCW pulses occur at high bias as resonances in the device interfere with the seed laser beam. Further increase of the current into the taper was limited by the pulse stability and the decreasing beam quality. CBC of four tapered amplifiers in QCW lead to a total power of 22.7 W with 64 % combining efficiency. The advantage in using tapered amplifiers for CBC is the extremely high power per device (here 5.67 W per amplifier) that exceeds the power achievable with single-mode amplifiers by more than a factor of 5. Higher total powers have so far only been achieved by coherent beam combining of large arrays of single-mode amplifiers at the cost of an increased complexity (CBC of 47 single-mode slab-coupled optical waveguide amplifiers yielded 40 W [START_REF] Creedon | High efficiency coherent beam combining of semiconductor optical amplifiers[END_REF]). The achieved power in the central lobe of the combined beam is higher than 20 W (92 % of total power). This is the highest power in a diffraction-limited beam ever reported from a single-frequency system based on tapered lasers or amplifiers. The operation in QCW mode allowed to overcome thermal limitations and to reach about 25 % higher power in comparison to CW operation.

The power could be further increased by using more efficient semiconductor structures, either directly as here or fabricated as single-contact truncated tapered amplifiers, which provide the highest reported diffraction limited power (QCW power in the central lobe up to 17 W per device [START_REF] Wang | 17-W Near-Diffraction-Limited 970-nm Output From a Tapered Semiconductor Optical Amplifier[END_REF]). Higher powers are also likely to be possible by using lower reflectivity coatings and tilted facets the tapered amplifiers to suppress internal oscillations. Combining these approaches offers a potential path for the development of compact CBC modules exceeding 50 W output power with close to diffraction limited beam quality. Efforts to reduce the pulse durations to the microsecond and nanosecond regime are also promising for potential applications in biomedicine [15,16] and remote sensing [START_REF] Nehrir | Eye-Safe Diode-Laser-Based Micropulse Differential Absorption Lidar (DIAL) for Water Vapor Profiling in the Lower Troposphere[END_REF][START_REF] Nehrir | Micropulse water vapor differential absorption lidar: transmitter design and performance[END_REF].

Fig. 1 .

 1 Fig. 1. Experimental setup of the CBC-interferometer with four high brightness tapered amplifiers (A 1-4 ). The currents into the DFB seed laser and the RWs of the amplifiers were continuous and only the currents into the taper of each amplifier were pulsed. (L 1 and L 2 aspheric lens f = 8 mm NA = 0.5, L 3 aspheric lens f =2.35 mm NA=0.55, L 4 cylindrical lens f = 19 mm).

  Fig. 2 ampli differe Gauss condit duty c shows end of Tapered depends on th beam shape in Visualization

  in wn for one typi CW and QCW m ms after the sta e and beam sha close to diffra higher angle si ses at higher c , which is give hapes between erve a signifi n QCW operati of the four amp tent at equal bi nment toleranc tput power versus QCW operation. ns. The power co ntral lobe. The int urements: P in = 20 lization 1 for the e beam shape meas mplifiers have a the tapered sec art of each QC d 20 images of mall part of the he signal for t int . (b) Measuremen ontent in the centr tegration time of 0 mW, I rw = 350 m spatial turn-on dy sured with 50 µs i a temperature ction [24,25]. W pulse. This f the beam shap e output beam the microcontr ers gths in the msd the limitation d in QCW on th 2(a). While th n CW mode (d efficiency rem amplifiers at 1 CW operation. T (A3) in Fig. 2( am shapes show e with 50 µs in d. The observe central lobe c e slow axis (SA the taper leadi for each bias CW operation tion of the be a central lobe quite similar wi (+/-5 %), con pling into the n to the tapered sec nts of the beam ral lobe was calc the CCD camera mA, T = 25°C, Δ ynamics in QCW integration time fr induced wand This leads to a effect is show ape into a short (less than 0.5 roller based clo -range and du ns normally se he optical outp here is an ons dashed lines), t mains constant 16 A. The outp The beam shap (b) for differen wn in Fig. 2(b) ntegration time ed shapes are ty containing mo A). The intens ing to decrease point in Fig 2 are minor at eam quality w power content ith minor diffe nsidering that narrow RW. ction for the four shape at waist at culated with a 2D was 50 µs. (Bias Δt qcw = 1 ms, 10 % W mode. The video rom the start to the dering astigma a dynamic chan wn in the suppl t video showin

  Fig. 3 curren points operat in ord device lines. ms, 10

  Fig. 4 condit duty c These qua spectral line s facets are AR

Fig. 5 .

 5 Fig. 5. (a) Measurement of the seed laser frequency perturbations -(blue), and simulated laser frequency perturbations (red) induced by optical feedback from a single tapered amplifier (A1) operated at I tp = 13 A (black). (b) Zoom into the measured frequency change with the calculated solutions (black dots) for the seed laser frequency considering the phase change in the amplifier and the expected frequency hopping in the multi-stable regime (red solid line). The simulation parameters used were: = 5, = 0.01, = 9 × 10 , = 50 , = 10 + , where was taken from Fig 3(a)). The optical isolation was -55 dB.

  Fig. 6 A into (see t combi (Bias ms, 10

  Figure 6 the envelope o thin lines). Th to 12.5 % at section 3.3, w levels of opti resonance eff of the combin quality of an lasers or ampl %, with only W at 14.5 W beam from ta semiconducto5. ConclusIn this paper continuous op the beam shap the rising tem tapered sectio

  ut of the CBC inter s we measured the shape (thin line). nt and the measure P in = 20 mW, I rw = tability at thre e (see thick lin of the pulse en the resonance ability in the s C at even high ol and the oscil iffraction limit eam clean-up is ntent in the cen s in SA. The ev action

Table 1 .

 1 Measured combined power P CBC and combining efficiency for coherent superposition of two (A1+A2) and four (A1+A2+A3+A4) amplifiers at different currents into the tapered sections. Operating conditions for each individual amplifier: P in = 20 mW, I rw = 300-400 mA, T = 25°C, Δt qcw = 2 ms, 10 % duty cycle.

			QCW			CW
		2 amplifiers	4 amplifiers	4 amplifiers
	I tp [A]	P CBC [W]	η CBC [%]	P CBC [W]	η CBC [%]	P CBC [W]	η CBC [%]
	8	-		-		10.8	72
	10	9.77	84	16.8	71	14.4	71
	12.3	12.4	85	20.4	68	16.9	71
	13	12.7	81	21.7	69	-	-
	14.5	13.7	79	22.7	64	-	-

Acknowledgments

The authors are grateful to all colleagues involved into device fabrication and mounting at the Ferdinand-Braun-Institut in Berlin.

Parts of this work have been previously been presented at the CLEO®/Europe-EQEC 2019 conference in Munich, Germany (Paper CB-4.6).