# A Sparse Tensor Generator with Efficient Feature Extraction

Tugba Torun<EMAIL_ADDRESS>0000-0001-6790-7094 ,  Eren Yenigul
<EMAIL_ADDRESS>0009-0000-8209-8037 ,  Ameer Taweel<EMAIL_ADDRESS>0009-0006-3444-1439 and  Didem Unat<EMAIL_ADDRESS>0000-0002-2351-0770 Koç
UniversityİstanbulTurkey

###### Abstract.

Sparse tensor operations are gaining attention in emerging applications such
as social networks, deep learning, diagnosis, crime, and review analysis.
However, a major obstacle for research in sparse tensor operations is the
deficiency of a broad-scale sparse tensor dataset. Another challenge in sparse
tensor operations is examining the sparse tensor features, which are not only
important for revealing its nonzero pattern but also have a significant impact
on determining the best-suited storage format, the decomposition algorithm,
and the reordering methods. However, due to the large sizes of real tensors,
even extracting these features becomes costly without caution. To address
these gaps in the literature, we have developed a smart sparse tensor
generator that mimics the substantial features of real sparse tensors.
Moreover, we propose various methods for efficiently extracting an extensive
set of features for sparse tensors. The effectiveness of our generator is
validated through the quality of features and the performance of decomposition
in the generated tensors. Both the sparse tensor feature extractor and the
tensor generator are open source with all the artifacts available at
https://github.com/sparcityeu/feaTen and https://github.com/sparcityeu/genTen,
respectively.

sparse tensors, feature extraction, synthetic data generation

##  1\. Introduction

Several applications such as social networks, deep learning, diagnosis, crime,
and review analysis require the data to be processed in the form of multi-
dimensional arrays, namely tensors (Kuang et al., 2014; Leclercq and Savonnet,
2018; Acar et al., 2007; Mu et al., 2011; Jukić et al., 2013; Yelundur et al.,
2019). Tensors are extensions of matrices that can have three or more
dimensions, or so-called modes. Tensor decomposition techniques such as CPD
and Tucker are widely used to analyze and reveal the latent features of such
real-world data in the form of sparse tensors (Kolda and Bader, 2009; Goulart
et al., 2015). For making the tensor decomposition methods faster or more
memory-efficient, numerous storage formats and partitioning or reordering
schemes are introduced in the literature (Li et al., 2019; Acer et al., 2018;
Kaya and Uçar, 2018; Li et al., 2018b; Shivakumar et al., 2021; Helal et al.,
2021). However, the performance of these schemes highly depends on the
sparsity pattern of the tensor.

Efficiently extracting the sparsity pattern of a given tensor is essential for
optimizing various aspects of tensor analysis and manipulation. The structural
features can inform format selection, aiding in determining the most suitable
storage format. It can influence the algorithm selection, with different
tensor operations exhibiting varying performance based on the sparsity
characteristics. Understanding the sparsity pattern also provides insights
into the performance of decomposition techniques. In a recent work, (Sun et
al., 2021), some tensor features are used to automatically predict the best
storage format for a sparse tensor via machine learning techniques. However,
the work considers the features for only a single mode, which may cause lack
of some critical intuition from other dimensions.

Sparse tensors in real-world scenarios often exhibit extreme sparsity, with
densities as low as $10^{-15}$. Unlike the sparse matrices often containing at
least one nonzero in their rows and columns, sparse tensors contain numerous
empty fibers and slices, which are one- and two-dimensional fragments of
tensors. A naive approach for extracting sparse tensor features involves
traversing the tensor nonzeros in coordinate (COO) format and updating the
nonzeros of respective slices and fibers. However, this approach becomes
impractical for large tensors with increasing dimensions. Another approach
(Sun et al., 2021) is to assume that the tensor is already sorted to extract
the features, yet it reveals the features from only a single mode order. Since
the sizes of a sparse tensor diverge a lot along different modes, focusing the
features solely on one perspective might lead to the loss of crucial
structural information inherent in the tensor.

To tackle these challenges, we develop a sparse tensor feature extraction
framework, FeaTen, which extracts a detailed and exclusive set of sparse
tensor features, encapsulating the features along all modes. It extends the
feature set by including important size-independent features such as
coefficient of variation and imbalance to gain a deeper insight into the
nonzero distribution. Additionally, FeaTen offers four alternative feature
extraction methods for efficiency, providing flexibility to select the most
suitable method based on machine and tensor characteristics.

By utilizing the generated features, machine learning tools can be used to
reveal the most suitable storage format, partitioning, or reordering method
that fits best with that tensor. Nevertheless, the primary challenge facing
this research stems from the necessity of having thousands of samples to train
machine learning models. Conversely, the majority of multi-dimensional real-
world data require manual cleaning to become readily usable in research.
Meanwhile, efforts to gather publicly available real-world data as sparse
tensor collections yield only a few instances (Jeon et al., 2015; Smith et
al., 2017). Moreover, tensor sizes can be large, making it inconvenient to
download, read, and use them in computation. A tensor generator for
performance analysis purposes can be handy, enabling the study of algorithms
by generating tensors on the fly and discarding them if necessary.

To fill these gaps in the literature, we propose a smart sparse tensor
generator, GenTen, that takes the significant features of tensors into
consideration. The proposed generator can be utilized to produce sparse
tensors with features resembling those of real tensors, thereby forming a
large suite of sparse tensor dataset. Moreover, an important benefit of this
generator is its use of size-independent features like coefficient of
variation, imbalance, and density, enabling easy generation of instances with
different sizes, thereby enhancing its versatility across various contexts.
The effectiveness of the generator is validated by comparing the features and
CPD performance of the generated tensors with the naive random and the
original real-world tensors.

The main contributions of this work are:

  * •

A sparse tensor feature extraction framework, FeaTen, is developed. It
includes four different feature extraction methods, which can be used
alternatively depending on the computation needs and the characteristics of
input tensors. All methods in FeaTen are parallelized using OpenMP.

  * •

We develop a smart sparse tensor generator, GenTen, which considers
significant tensor features. It can be used to generate an artificial tensor
dataset in which the properties and characteristics of the generated tensors
are similar to the real-word tensors. GenTen is parallelized with OpenMP for
faster generation.

  * •

Both tools are open source 111https://github.com/sparcityeu/feaTen
222https://github.com/sparcityeu/genTen, accompanied by comprehensive
documentation and illustrative examples tailored for the community’s usage.

This manuscript is organized as follows. Section 2 provides the background
information on sparse tensors and tensor decomposition. The proposed feature
extraction tool is introduced in Section 3. In Section 4, we present our
sparse tensor generator, GenTen. Experimental results are reported in Section
5. We discuss the related works in Section 6 and conclude in Section 7.

##  2\. Background

A tensor with $M$ dimensions is called an $M$-$mode$ or $M$th $order$ tensor,
where the mode count, $M$, is referred to as the $order$ of tensor. $Mode$ $m$
of a tensor refers to its $m$th dimension. Fibers are defined as the one-
dimensional sections of a tensor obtained by fixing all but one index. Slices
are two-dimensional sections of a tensor obtained by fixing every index but
two. Figure 1 depicts sample slice and fibers of a 3-mode tensor. The numbers
in the naming of slices and fibers derive from the mode indices that are not
fixed while forming them. For instance, $\mathcal{X}(i,:,:)$ is a mode-(2,3)
slice, and $\mathcal{X}(i,j,:)$ is a mode-3 fiber of a 3-mode tensor
$\mathcal{X}$, where a tensor element with indices $i,j,k$ is denoted by
$\mathcal{X}(i,j,k)$. A fiber (slice) is said to be a $nonzero$ $fiber$
($slice$) if it contains at least one nonzero element; and an $empty$ $fiber$
($slice$), otherwise.

Figure 1. Sample slice and fibers of a 3-mode tensor.

To reveal the relationship of data across different modes, tensor
decomposition techniques are widely used. Canonical Polyadic Decomposition
(CPD) and Tucker decomposition are the most two popular ones among them. In
Tucker decomposition, a tensor is decomposed into a much smaller core tensor
and a set of matrices; whereas in CPD, a tensor is factorized as a set of
rank-1 tensors, which can be considered as a generalization of matrix singular
value decomposition (SVD) method for tensors.

There are various implementations for CPD and each has different storage
schemes and reordering algorithms proposed specifically for that scheme.
SPLATT (Smith et al., 2015) and ParTI (Li et al., 2018a) are the two most
commonly-used libraries for CPD. ParTI is a parallel tensor infrastructure
that supports essential sparse tensor operations and tensor decompositions on
multicore CPU and GPU architectures. SPLATT provides a shared-memory
implementation for CPD while adopting a medium-grain partitioning for sparse
tensors for parallel execution of CPD.

##  3\. Feature Extraction

The features of a sparse tensor have the capacity to reveal the most suitable
storage format, partitioning or reordering method that fits best with that
tensor if well-examined. The main challenge in extracting the sparse tensor
features is to determine which and how many fibers and slices are nonzero.
This is because the real-world sparse tensors are highly sparse and contain
many empty fibers and slices. One naive approach for sparse tensor feature
extraction is to traverse the tensor nonzeros in coordinate (COO) format and
update the number of nonzeros of the respective slices and fibers. However,
this approach is not practical for large tensors since the real-world tensors
often have huge numbers of fibers and slices. For instance, some real tensors
have more than $10^{19}$ fibers, and even storing a boolean array of such a
large size requires a space of 10 million Terabytes, which is impractical in
modern machines (See the properties of some real tensors in Table 3). To
overcome these challenges, we propose and implement four alternative methods
for sparse tensor feature extraction.

###  3.1. Feature Set

We consider three kinds of statistics for sparse tensors: (i) nonzeros per
nonzero slices, (ii) nonzeros per nonzero fibers, (iii) nonzero fibers per
nonzero slices. Table 2 depicts the global features of the tensor that are
independent of these kinds. Table 2 shows the mode- and kind-dependent
features that we have extracted for all of these three kinds and for all
modes. Here by referring to all modes, we mean all possible angles that a
slice or fiber can have through fixing different modes. For example,
mode-(1,2), mode-(2,3), and mode-(3,1) slices represent all modes for slices
of a 3-mode tensor. For ease of expression, we refer to the set of slices
(fibers) along different modes as $slice$-$modes$ ($fiber$-$modes$).

A tensor $\mathcal{X}$ of size $I_{1}\times I_{2}\times\dots\times I_{M}$ is
assumed to be given as an input in Coordinate format, i.e. the extension of
matrix-market format for tensors. In the formulas, nz($\cdot$) is used as a
function returning the number of nonzeros. $\mathcal{F}_{m}$ refers to the set
of fibers in mode-$m$, whereas $\mathcal{S}_{k,\ell}$ refers to the set of
slices obtained by fixing indices in modes $k$ and $\ell$. In the feature
names, $\bullet_{nz}$ refers to considering only nonzero kind-entries in the
computation, i.e. nonzero fibers or slices, by ignoring the empty fibers or
slices. $\bullet_{all}$ refers to considering all kind-entries in the
computation, including the empty ones. For instance, for the nonzeros-per-
fiber kind in mode 1, $n_{all}$ gives the number of all mode-1 fibers
including empty, whereas $n_{nz}$ gives the number of only nonzero fibers
along mode-1.

Table 1. Extracted global features of tensors.

Feature | Description | Formula
---|---|---
size_m | Size of tensor in mode m | $I_{m}$
nnz | Number of nonzeros | $\texttt{nz}(\mathcal{X})$
$\mathrm{d}_{\mathrm{nz}}$ | Density of nonzeros | $\texttt{nz}(\mathcal{X})/\prod I_{m}$
$\mathrm{nfib}_{\mathrm{all}}$ | Number of all fibers | $\sum|\mathcal{F}_{m}|$
$\mathrm{nslc}_{\mathrm{all}}$ | Number of all slices | $\sum|\mathcal{S}_{k,\ell}|$
$\mathrm{nfib}_{\mathrm{nz}}$ | Number of nz fibers | $\sum|\texttt{nz}(\mathcal{F}_{m})|$
$\mathrm{nslc}_{\mathrm{nz}}$ | Number of nz slices | $\sum|\texttt{nz}(\mathcal{S}_{k,\ell})|$
$\mathrm{d}_{\mathrm{fib}}$ | Density of nz fibers | $\mathrm{nfib}_{\mathrm{nz}}/\mathrm{nfib}_{\mathrm{all}}$
$\mathrm{d}_{\mathrm{slc}}$ | Density of nz slices | $\mathrm{nslc}_{\mathrm{nz}}/\mathrm{nslc}_{\mathrm{all}}$

Table 2. Mode- and kind- dependent features that are extracted for each mode
and for each kind: (i) nonzeros per slice, (ii) nonzeros per fiber, (iii)
fibers per slice.

Feature | Description | Formula
---|---|---
$\mathrm{n}_{\mathrm{all}}$ | All count including empty | *
$\mathrm{n}_{\mathrm{nz}}$ | Nonzero count | *
nz_density | Nonzero sparsity | $\mathrm{n}_{\mathrm{nz}}/\mathrm{n}_{\mathrm{all}}$
max | Maximum value | *
min | Minimum value | *
dev | Deviation | max - min
sum | Summation of values | *
$\mathrm{avg}_{\mathrm{all}}$ | Average value | sum /$\mathrm{n}_{\mathrm{all}}$
$\mathrm{imbal}_{\mathrm{all}}$ | Imbalance | (max-$\mathrm{avg}_{\mathrm{all}}$)/max
$\mathrm{stDev}_{\mathrm{all}}$ | Standard deviation | *
$\mathrm{cv}_{\mathrm{all}}$ | Coefficient of variation | $\mathrm{stDev}_{\mathrm{all}}/\mathrm{avg}_{\mathrm{all}}$
$\mathrm{avg}_{\mathrm{nz}}$ | Average by excluding empty | sum /$\mathrm{n}_{\mathrm{nz}}$
$\mathrm{imbal}_{\mathrm{nz}}$ | Imbalance by excluding empty | (max-$\mathrm{avg}_{\mathrm{nz}}$)/max
$\mathrm{stDev}_{\mathrm{nz}}$ | stDev by excluding empty | *
$\mathrm{cv}_{\mathrm{nz}}$ | cv by excluding empty | $\mathrm{stDev}_{\mathrm{nz}}/\mathrm{avg}_{\mathrm{nz}}$

*Omitted due to the complex or trivial formula.

In addition to the tensor features utilized in the related work (Sun et al.,
2021), we also include the features of load imbalance, standard deviation, and
coefficient of variation. In (Nisa et al., 2019), it is shown that a high
standard deviation of fiber length causes inter-warp load imbalance and low
occupancy; whereas high standard deviation of the slice volume is related to
significant inter-thread-block load imbalance. The coefficient of variation is
another commonly-used metric for analysis which allows to compare between data
sets with widely different means. For example, it is used to evaluate the
dispersion of the number of nonzero elements per row for sparse matrix
computations (Abu-Sufah and Abdel Karim, 2013).

###  3.2. Extraction Workflow

The feature extraction process consists of two main stages, namely array
construction and the final reduction phases. In the array construction phase,
we construct three arrays for all modes: (i) $n^{nz}_{slc}$ (number of
nonzeros per slice), (ii) $n^{nz}_{fib}$ (number of nonzeros per fiber), (iii)
$n^{fib}_{slc}$ (number of fibers per slice). Here and hereafter, only the
nonzero slices and nonzero fibers are considered when referring to slices and
fibers, if not stated otherwise. Then in the final reduction phase, we extract
all the sparse tensor features by traversing or applying a reduction on these
smaller arrays. Both the array construction and the final reduction phases of
FeaTen are parallelized using OpenMP.

The arrays of $n^{nz}_{slc}$ and $n^{fib}_{slc}$ are constructed for each
slice-mode, so that there will be ${M\choose 2}$ of them. The array
$n^{nz}_{fib}$ is constructed for each fiber-mode, hence yielding $M$
different $n^{nz}_{fib}$ arrays. For instance, for a 3-mode tensor, there are
3 different types of statistics for $n^{nz}_{slc}$ and $n^{fib}_{slc}$ which
correspond to mode-(1,2), (2,3), and (3,1) slices; whereas for $n^{nz}_{fib}$,
there are 3 different types of statistics corresponding to mode-1, 2, and 3
fibers. Since there are 11 global features and 15 mode- and kind-dependent
features, the total number of features extracted for a 3-mode tensor is
$15\times 3\times 3+11=146$ in total. For a 4-mode tensor, slices are along
${4\choose 2}=6$ different mode pairs, and fibers are along 4 different modes,
so in total $15\times 6\times 2+15\times 4+11=251$ features are considered.
For 5-mode tensors, slices are along ${5\choose 2}=10$ different mode pairs,
and fibers are along 5 different modes, yielding a total of $15\times 10\times
2+15\times 5+11=386$ features.

One approach in FeaTen is to consider the slice and fiber modes in relation to
mode execution orders. By saying $mode$-$order$ $\langle i,j,k\rangle$, we can
think that the COO-based index arrays of tensor $\mathcal{X}$ are virtually
rearranged to $\langle\mathcal{X}_{i},\mathcal{X}_{j},\mathcal{X}_{k}\rangle$
instead of the original order
$\langle\mathcal{X}_{1},\mathcal{X}_{2},\mathcal{X}_{3}\rangle$. In this
context, we consider the slices are obtained by fixing all indices except the
last two index, and the fibers are obtained by fixing all indices except the
last one. For instance in mode order $\langle 1,2,3\rangle$, we consider
mode-(2,3) slices and mode-3 fibers. The advantage of this approach is to
combine slices and fibers so that their corresponding computations may
overlap.

Note that especially when computing $n^{fib}_{slc}$ arrays, two different
fiber modes can be associated with each slice mode. To be more specific, to
extract the number of fibers per mode-$(i,j)$ slices, it is possible to
consider either mode-$i$ or mode-$j$ fibers. For this reason, to cover all
slice and fiber modes, one can use both the mode-order set ${\langle
1,2,3\rangle,\langle 2,3,1\rangle,\langle 3,1,2\rangle}$ or the mode-order set
${\langle 1,3,2\rangle,\langle 2,1,3\rangle,\langle 3,2,1\rangle}$. Fig. 2
depicts the workflow of the feature extraction process for a 3-mode sparse
tensor for the first case. In the final reduction phase, nine different sets
of mode- and kind-dependent features are extracted in parallel, corresponding
to three different versions of $n^{nz}_{slc}$, $n^{nz}_{fib}$, and
$n^{fib}_{slc}$ arrays from distinct mode-orders.

Figure 2. The workflow of feature extraction for a 3-mode tensor using the
mode-order approach.

The mode-order approach is especially practical for 3-mode tensors. It is
because the number of slice-modes and fiber-modes are both equal to three,
making a one-to-one mapping possible between slice- and fiber-modes, hence
execution on three distinct mode-orders is sufficient to cover all cases.
However for higher dimensions, it can cause some calculations for
$n^{nz}_{fib}$ to be repeated since the number of slice-modes becomes larger
than the number of fiber-modes. Therefore we obey different approaches in
different methods of FeaTen.

###  3.3. Extraction Methods

In order to overcome the memory and speed limitations in trivial feature
extraction approaches, we propose four different feature extraction methods
that can be used alternatively depending on the computation needs and the
characteristics of input tensors. These methods are heap-based, sorting-based,
grouping-based, which is first proposed in this work, and a hybrid method,
which is a combination of sorting- and grouping-based methods.

##### Hash-based method:

We implement a hash-table-based method to solve the memory issue of the COO-
based naive method for large sparse tensors. Instead of long arrays that
mostly contain zero values, we use a hash table that has keys as slice or
fiber indices and values as the number of nonzeros in the corresponding slice
or fiber. One main difference of the hash-based method from other methods in
FeaTen is that it does not follow the mode-order approach, rather it extracts
the slice and fiber features independently.

##### Sorting-based method:

Another approach for feature extraction is based on sorting the tensor. It is
the conventional approach for extracting the features or constructing
compressed storage formats for tensors in the literature (Sun et al., 2021;
Smith and Karypis, 2015; Li et al., 2018b). It naturally follows the mode-
order approach by sorting the tensor indices according to a given mode-order.
After sorting, the array construction phase becomes easier since the nonzeros
belonging to the same fiber or slice are positioned consecutively. The related
work (Sun et al., 2021) assumes the tensor to be given as sorted and computes
the features only for a single mode-order. However, this assumption is not
valid when it comes to executing in the upcoming mode-orders. In other words,
a tensor $\mathcal{X}$ sorted in $\langle 1,2,3\rangle$ order has to be re-
sorted for $\langle 2,3,1\rangle$ and $\langle 3,1,2\rangle$ orders to execute
the features for other modes. Instead, we find the tensor features along all
modes by sorting the related indices of tensor.

##### Grouping-based method:

Instead of sorting the tensor nonzeros fully, we group the slices and fibers
according to their indices and keep the last indices of the COO format in
their original order. We first find the number of nonzero fibers and slices by
traversing the tensor entries. We store the indices of these nonzero slices
and fibers in a compressed manner and track their nonzero counts in the
respective arrays. This algorithm takes inspiration from the construction
process of Compressed Sparse Row and Column (CSR and CSC) formats for sparse
matrices.

Algorithm 1 shows the pseudo-code of the grouping-based feature extraction
method for the mode-order $\langle 1,2,3\rangle$ of a 3-mode tensor. In this
case, the algorithm considers mode-(2,3) slices and mode-3 fibers. $cnt1$ is a
temporary array that keeps the nonzero count information for all slices (Lines
1-4). It is used to obtain the number of nonzero slices ($nslc$), the array of
nonzero counts for each nonzero slice ($n^{nz}_{slc}$), the array of indices
for nonzero slices ($ind_{slc}$), and the array of starting locations for
nonzero slices (($xn^{nz}_{slc}$) that is similar to the $row\\_ptr$ array in
CSR format for matrices (Lines 5-8). $loc$ and $order$ are temporary arrays
that are constructed (Lines 10-16) to obtain the $cnt2$ array, which keeps the
nonzero count information for all fibers in each nonzero slice (Lines 18-21).
Finally, $cnt2$ is used to obtain the number of nonzero fibers
($n^{fib}_{slc}$) and the array of nonzero counts per nonzero fiber
($n^{nz}_{fib}$) in each nonzero slice (Lines 22-23).

Algorithm 1 Grouping-Based Extraction Method in FeaTen

Input: 3-mode tensor with COO-arrays $\mathcal{X}_{1}$, $\mathcal{X}_{2}$,
$\mathcal{X}_{3}$

Output: $n^{nz}_{slc}$, $n^{nz}_{fib}$, $n^{fib}_{slc}$ arrays for mode-order
$\langle 1,2,3\rangle$

1:Create the array $cnt1$ of size $I_{1}$ and initialize to zero

2:for $i\leftarrow 0$ to NNZ do

3: $cnt1(\mathcal{X}_{1}(i)))+\\!+$ $\rhd$ Temporary array for nonzeros per
slices

4:end for

5:$nslc\leftarrow$ # nz values in $cnt1$ $\rhd$ Number of nonzero slices

6:$n^{nz}_{slc}\leftarrow$ nz values in $cnt1$ (size $I_{1}$) $\rhd$ Number of
nonzeros per nonzero slice

7:$ind_{slc}\leftarrow$ indices of nz values in $cnt1$ (size $I_{1}$)

8:$xn^{nz}_{slc}\leftarrow\texttt{PrefixSum}(n^{nz}_{slc})$

9:Create the array $loc$ of size $ind_{slc}(I_{1})\\!+\\!1$

10:for $i\leftarrow 0$ to $nslc$ do

11: $loc(ind^{nz}_{slc}(i))\leftarrow xn^{nz}_{slc}(i)$

12:end for

13:Create the array $order$ of size $NNZ$

14:for $i\leftarrow 0$ to NNZ do

15: $order(loc(\mathcal{X}_{1}(k))+\\!+)\leftarrow k$

16:end for

17:for $i\leftarrow 0$ to $nslc$ do

18: Create the array $cnt2$ of size $I_{2}$ and initialize to zero $\rhd$
Temporary array for nonzeros per fiber

19: for $j\leftarrow xn^{nz}_{slc}(i)$ to $xn^{nz}_{slc}(i+1)$ do

20: $cnt2(\mathcal{X}_{2}(order(j)))+\\!+$

21: end for

22: $n^{fib}_{slc}(i)\leftarrow$ # nz in $cnt2$ $\rhd$ Number of nonzero
fibers per nonzero slice

23: $n^{nz}_{fib}(i)\leftarrow$ nz values in $cnt2$ (size $n^{fib}_{slc}(i)$)
$\rhd$ Number of nonzeros per nonzero fiber

24:end for

The memory requirement of this algorithm is $I_{1}\\!+\\!I_{2}\\!+\\!NNZ$ for
executing the $(i_{1},i_{2},:)$ fibers of an
$I_{1}\\!\times\\!I_{2}\\!\times\\!I_{3}$ tensor with $NNZ$ nonzeros. The
worst case serial runtime is $\mathcal{O}(I_{1}+I_{2}\mathcal{S}_{1})$ where
$\mathcal{S}_{1}$ is the number of nonzero slices in mode 1. For the general
case, $size(m)\times size(m+1)$ is a loose upper bound for mode $m$.

##### Hybrid method:

We propose a hybrid method by combining the sorting- and grouping-based
methods. The idea is to utilize different methods for extracting features in
different modes within the same tensor. This is because the performance of
different methods varies depending on the size of the respective mode. Since
we observe that sorting- and grouping-based methods are the best-performing
ones at the mode level, we use them interchangeably according to the
respective mode size. For mode $m$, we apply the grouping-based method if
$size(m)\times size(m+1)<10^{9}$ holds, and we apply the sorting-based method
otherwise.

###  3.4. Extension to Higher Orders

Hash-based method is the only method in FeaTen that does not follow the mode-
order approach. Despite its slight disadvantage for the 3-mode case, this
difference makes it more applicable for higher-order tensors due to the
increased gap between the number of slice-modes and fiber-modes. Additionally,
the vast number of features required for extracting higher-order tensors
complicates the evaluation process. For these reasons, we only allow the hash-
based method to extract all features along all modes for tensors with an order
higher than 3. For other methods in FeaTen, we opt to extract the features
along the modes with the three largest sizes. For this, we simply sort the
tensor modes according to their sizes, select the three largest ones, and
extract the features corresponding solely to those modes. We refer to it as
the $only$-3-$mode$ option, which is also available for the hash-based method,
in addition to the choice to extract features along all $M$ modes.

##  4\. Sparse Tensor Generator

Despite the increased need for research in sparse tensors, publicly available
sparse tensor datasets remain limited, comprising only a few instances (Jeon
et al., 2015; Smith et al., 2017). To address this scarcity and facilitate the
study of machine learning models with a larger variety of sparse tensors, as
well as to expand the size of open datasets, we introduce a smart sparse
tensor generator, called GenTen. This generator also enables rapid evaluation
of proposed methods and algorithms without the necessity of storing the
tensor.

###  4.1. Overview

Our generator considers the significant features of tensors. Consequently, the
artificial tensors produced by GenTen closely emulate real tensors with their
respective features and characteristics. A notable advantage of this developed
generator is its utilization of size-independent features, such as coefficient
of variation, imbalance, and density. This allows for easy generation of
instances with varying sizes, enhancing its versatility and applicability in
diverse contexts.

In GenTen, we provide the option to get the seed for pseudo-randomness from
the user. Through this, one can obtain the exact same tensor when providing
the same seed in different tests. Moreover, the user can create tensors with
almost the same properties by simply changing the seed. The effect of seed
selection will be discussed in Section 5.3.3. All levels of GenTen are
parallelized using OpenMP for faster execution.

The main idea of the proposed tensor generator is to first determine fiber
counts per slice, and the number of nonzeros per fiber according to the given
coefficient of variation values. The generator has the flexibility to employ
any distribution to determine these counts, yet we are currently utilizing
normal or log-normal distributions to determine the nonzero layout. The
positions of the nonzeros and nonzero slice or fiber indices are selected
uniformly. The nonzero values are selected uniformly from the interval (0,1).

The generator is general and can generate any $M$-mode tensor. For ease of
expression, we describe the algorithm for a 3-mode tensor. Since generating
random sparse tensors while simultaneously adhering to features along all
modes is nearly impossible, we opt to consider the features of mode-(2, 3)
slices and mode-2 fibers as inputs.

###  4.2. The Proposed Algorithm

The pseudocode of the proposed sparse tensor generator (GenTen) is shown in
Algorithm 2. The generator utilizes a set of given metrics to create a tensor
having these features: (i) sizes of the tensor $I_{1}$, $I_{2}$, $I_{3}$, (ii)
densities of slices, fibers, and nonzeros $d^{slc}$, $d^{fib}$, $d^{nz}$,
(iii) coefficient of variations for fibers per slice and nonzeros per fiber
$cv_{slc}^{fib}$, $cv_{fib}^{nz}$, and (iv) imbalance for fibers per slice and
nonzeros per fiber $imbal_{slc}^{fib}$, $imbal_{fib}^{nz}$.

Algorithm 2 GenTen: Sparse Tensor Generator

Input: $I_{1}$, $I_{2}$, $I_{3}$, $d^{slc}$, $d^{fib}$, $d^{nz}$,
$cv_{slc}^{fib}$, $cv_{fib}^{nz}$, $imbal_{slc}^{fib}$, $imbal_{fib}^{nz}$

Output: Tensor $\mathcal{X}$ stored in COO format

1:Calculate $nslc$, $avg_{slc}^{fib}$, $stDev_{slc}^{fib}$, $max_{slc}^{fib}$

2:$ind_{slc}\leftarrow\texttt{RandInds}(nslc,I_{1})$ $\rhd$ Array of nonzero
slice indices

3:$(n^{fib}_{slc},ind_{fib})\leftarrow\texttt{Distribute}(nslc,avg_{slc}^{fib},stDev_{slc}^{fib},max_{slc}^{fib},I_{2})$
$\rhd$ Number and indices of fibers per slice

4:$xn^{fib}_{slc}\leftarrow\texttt{PrefixSum}(n^{fib}_{slc})$

5:$nfib\leftarrow xn^{fib}_{slc}(nslc)$ $\rhd$ Number of nonzero fibers

6:Calculate $avg_{fib}^{nz}$, $stDev_{fib}^{nz}$, $max_{fib}^{nz}$

7:$(n^{nz}_{fib},ind_{nz})\leftarrow\texttt{Distribute}(nfib,avg_{fib}^{nz},stDev_{fib}^{nz},max_{fib}^{nz},I_{3})$
$\rhd$ Number and indices of nonzeros per fiber

8:$xn^{nz}_{fib}\leftarrow\texttt{PrefixSum}(n^{nz}_{fib})$

9:for $i\leftarrow 0$ to $nslc$ do

10: for $j\leftarrow 0$ to $n^{fib}_{slc}(i)$ do

11: $j^{\prime}\leftarrow xn^{fib}_{slc}(i)+j$ $\rhd$ Global nonzero fiber
index

12: for $k\leftarrow 0$ to $n^{nz}_{fib}(j^{\prime})$ do

13: $k^{\prime}\leftarrow xn^{nz}_{fib}(j^{\prime})+k$ $\rhd$ Global nonzero
index

14: $\mathcal{X}(k^{\prime})\leftarrow\langle
ind_{slc}(i),ind_{fib}(i)(j),ind_{nz}(j^{\prime})(k)\rangle$ $\rhd$ Fill
nonzeros of tensor $\mathcal{X}$

15: end for

16: end for

17:end for

The algorithm first calculates the requested nonzero slice count ($nslc$), and
the average, standard deviation, and maximum values for fibers per slice
($avg_{slc}^{fib}$, $stDev_{slc}^{fib}$, $max_{slc}^{fib}$) from the given
inputs. In line 2, the nonzero slice indices are determined by the RandInds
function, which returns $nslc$ many different indices uniformly distributed in
the interval $[1,I_{1}]$. In line 3, the number of fibers per slice
($n_{slc}^{fib}$) array is constructed and the indices of these nonzero fibers
($ind_{fib}$) are determined respecting the $avg_{slc}^{fib}$,
$stDev_{slc}^{fib}$, and $max_{slc}^{fib}$ values, using the function
Distribute. The details of the Distribute method will be discussed later. A
simple prefix sum is applied on $n_{slc}^{fib}$ and we obtain
$xn_{slc}^{fib}$, to be mainly used in future calculations; yet the number of
nonzero fibers ($nfib$) is derived by using $xn_{slc}^{fib}$ (Lines 4-5).

We calculate the average, standard deviation, and maximum values for nonzeros
per fiber ($avg_{fib}^{nz}$, $stDev_{fib}^{nz}$, $max_{fib}^{nz}$) according
to the determined $nfib$ value. In line 7, the number of nonzeros per fiber
($n_{fib}^{nz}$) array is constructed and the indices of these nonzeros
($ind_{nz}$) are determined respecting the $avg_{fib}^{nz}$,
$stDev_{fib}^{nz}$, $max_{fib}^{nz}$ values, again using the Distribute
method. Similarly, the array $xn_{fib}^{nz}$ is obtained by applying prefix
sum on $n_{fib}^{nz}$ (Line 8). In the last stage of the algorithm (Lines
9-17), the indices of tensor $\mathcal{X}$ are filled using the arrays
$ind_{slc}$, $ind_{fib}$, and $ind_{nz}$.

The pseudocode of the Distribute method is given in Algorithm 3. This method
takes five values as inputs, then returns a count array $cnt$ and an index
array $inds$. We utilize the Box-Muller method (Golder and Settle, 1976) which
generates random numbers with normal distribution obeying a given standard
deviation and mean. Since this method obeys a continuous distribution, the
values are real and sometimes might be negative. Hovewer, our aim here is to
construct a count array, e.g. keeping track of the nonzero count in nonzero
fibers, so the target values should be positive integers. Thus we round the
generated values to the nearest positive integers. To avoid negative values,
we switch to the log-normal distribution when needed, which guarantees
positive real values.

Algorithm 3 Distribute

Input: Values $n$, $avg$, $std$, $max$, $limit$

Output: Arrays $cnt$ and $inds$

1:for $i\leftarrow 0$ to $n$ do

2: if $avg>3\times std$ then $\rhd$ Apply normal distribution

3: $cnt(i)\leftarrow\texttt{BoxMuller}(avg,std)$

4: else$\rhd$ Apply log-normal distribution

5: $avg_{logNorm}\leftarrow\texttt{log}(avg^{2}/(avg^{2}+std^{2})^{1/2})$

6: $std_{logNorm}\leftarrow(\texttt{log}(1+std^{2}/avg^{2}))^{1/2}$

7:
$cnt(i)\leftarrow\texttt{exp}(\texttt{BoxMuller}(avg_{logNorm},std_{logNorm}))$

8: end if

9:end for

10:$avg_{result}\leftarrow\texttt{sum}(cnt)/n$

11:$ratio=avg/avg_{result}$

12:if $ratio<0.95$ or $ratio>1.05$ then $\rhd$ Apply scaling if needed

13: for $i\leftarrow 0$ to $n$ do

14: $cnt(i)\leftarrow cnt(i)\times ratio$

15: end for

16:end if

17:for $i\leftarrow 0$ to $n$ do

18: $cnt(i)\leftarrow\texttt{min}(cnt(i),max)$

19: $cnt(i)\leftarrow\texttt{max}(cnt(i),1)$

20: $ind(i)\leftarrow\texttt{RandInds}(cnt(i),limit)$ $\rhd$ Array of size
$cnt(i)$, with elements in range $[1,limit]$

21:end for

The algorithm applies normal distribution if most values are expected to be
positive ($avg-3\times std>0$), or log-normal distribution, otherwise (Lines
1-8). This choice stems from the fact that approximately 99.7% of values
sampled from a normal distribution fall within three standard deviations of
the mean. Then the algorithm scales the values of $cnt$ if the ratio of the
resulting over the expected average is outside of a predetermined range, which
is determined as (0.95, 1.05) in our case (Lines 10-16). Finally in Lines
17-21, the values of $cnt$ are adjusted to obey the minimum and maximum
values, and the indices are selected uniformly using the RandInds method,
which returns $cnt(i)$ many indices in the interval $[1,limit]$.

###  4.3. Extension to Higher Orders

The algorithm of GenTen for higher orders is similar to the one for the 3rd-
order case with the main lines. For an $M$-mode tensor, we consider the
features of mode-($M\\!-\\!1,M$) slices and mode-($M\\!-\\!1$) fibers as
inputs. The primary additional challenge in a higher-order setting is
determining the indices of nonzero slices. It is because for an $M$-mode
tensor, the nonzero slice indices themselves form a $(M-2)$-mode smaller
boolean tensor, while they were simply scalars for a 3-mode tensor. Therefore
we propose a different scheme to determine the indices of nonzero slices
efficiently.

We consider four distinct cases regarding the given slice density, $d^{slc}$,
which is the ratio of nonzero slice count ($nslc_{nz}$) over the total slice
count ($nslc_{all}$). If $d^{slc}>0.97$, we round it to 1.0 and assign the
slice indices in sorted order. If $d^{slc}\in(0.5,0.97]$, we create an array
of size $nslc_{all}$ to keep track of empty slice indices with the assumption
that it can fit into memory, since $nslc_{all}$ can be at most 2 times larger
than $nslc_{nz}$. If $d^{slc}\in[0.1,0.5]$, instead of creating such an array,
we traverse all possible $nslc_{all}$ indices and select a fraction of
$d^{slc}$ of these indices uniformly. Finally if $d^{slc}<0.1$, we simply
generate $nslc_{nz}$ many random indices for each of the respective $(M-2)$
modes.

##  5\. Evaluation

###  5.1. Experimental setup

The experiments are conducted on an AMD EPYC 7352 CPU of 3200MHz with 512 GB
of memory. It has Zen 2 microarchitecture, which includes 2 sockets, and each
socket has 24 cores with 2-way simultaneous multi-threading. Both tools are
implemented in C/C++ utilizing OpenMP for shared memory parallelism and
compiled with GCC using optimization level O2.

The dataset is taken from two real-world sparse tensor collections, namely
FROSTT (Smith et al., 2017) and HaTen2 (Jeon et al., 2015). We have excluded
the tensors whose nonzero count is more than 1 billion. We have also excluded
the delicious-4d tensor since its fiber count exceeds the maximum value
($1.8\times 10^{19}$) for an unsigned long long int in the C language. As a
result, the dataset consists of 16 real-world sparse tensors whose properties
are given in Table 3.

Table 3. Properties of the real sparse tensors from the FROSTT and HaTeN2
collections.

|  |  |  |  |  |  | Slice Count | Fiber Count |
---|---|---|---|---|---|---|---|---|---
Name | Order | $I_{1}$ | $I_{2}$ | $I_{3}$ | $I_{4}$ | $I_{5}$ | All | Nonzero | All | Nonzero | NNZ
LBNL-network | 5 | 1,605 | 4,198 | 1,631 | 4,209 | 868,131 | 4.1E+13 | 1.0E+07 | 6.9E+16 | 6.8E+06 | 1,698,825
NIPS | 4 | 2,482 | 2,862 | 14,036 | 17 | - | 8.2E+07 | 3.8E+06 | 1.0E+11 | 7.3E+06 | 3,101,609
uber | 4 | 183 | 24 | 1,140 | 1,717 | - | 2.6E+06 | 2.6E+05 | 4.2E+08 | 2.8E+06 | 3,309,490
chicago-crime-comm | 4 | 6,186 | 24 | 77 | 32 | - | 8.3E+05 | 7.1E+05 | 3.1E+07 | 7.8E+06 | 5,330,673
chicago-crime-geo | 5 | 6,185 | 24 | 380 | 395 | 32 | 1.2E+09 | 3.0E+07 | 5.6E+10 | 2.8E+07 | 6,327,013
vast-2015-mc1-3d | 3 | 165,427 | 11,374 | 2 | - | - | 1.8E+05 | 1.8E+05 | 1.9E+09 | 2.6E+07 | 26,021,854
vast-2015-mc1-5d | 5 | 165,427 | 11,374 | 2 | 100 | 89 | 3.6E+11 | 1.2E+08 | 1.7E+13 | 1.0E+08 | 26,021,945
DARPA1998 | 3 | 22,476 | 22,476 | 23,776,223 | - | - | 2.4E+07 | 2.4E+07 | 1.1E+12 | 5.5E+07 | 28,436,033
enron | 4 | 6,066 | 5,699 | 244,268 | 1,176 | - | 3.2E+09 | 2.2E+07 | 1.2E+13 | 9.0E+07 | 54,202,099
NELL-2 | 3 | 12,092 | 9,184 | 28,818 | - | - | 5.0E+04 | 5.0E+04 | 7.2E+08 | 3.8E+07 | 76,879,419
freebase_music | 3 | 23,343,790 | 23,344,784 | 166 | - | - | 4.7E+07 | 4.6E+07 | 5.4E+14 | 2.2E+08 | 99,546,551
flickr-3d | 3 | 319,686 | 28,153,045 | 1,607,191 | - | - | 3.0E+07 | 3.0E+07 | 5.5E+13 | 1.5E+08 | 112,890,310
flickr-4d | 4 | 319,686 | 28,153,045 | 1,607,191 | 731 | - | 5.5E+13 | 1.9E+08 | 1.5E+19 | 2.8E+08 | 112,890,310
freebase_sampled | 3 | 38,954,435 | 38,955,429 | 532 | - | - | 7.8E+07 | 7.4E+07 | 1.5E+15 | 3.2E+08 | 139,920,771
delicious-3d | 3 | 532,924 | 17,262,471 | 2,480,308 | - | - | 2.0E+07 | 2.0E+07 | 5.3E+13 | 1.6E+08 | 140,126,181
delicious-4d | 4 | 532,924 | 17,262,471 | 2,480,308 | 1,443 | - | 5.3E+13 | 2.3E+08 | 4.4E+18 | 4.0E+08 | 140,126,220
NELL-1 | 3 | 2,902,330 | 2,143,368 | 25,495,389 | - | - | 3.1E+07 | 3.1E+07 | 1.3E+14 | 2.5E+08 | 143,599,552

###  5.2. Performance of Feature Extraction Methods

The performance of four different feature extraction methods in FeaTen are
compared. For this experiment, we use the maximum number of available hardware
threads in the machine, which is 96. Since all feature extraction methods give
the same and exact results as feature sets, we do not present an accuracy
comparison among them. Since our dataset includes both 3rd order and higher
order tensors, for a fair comparison, we run all four methods in FeaTen with
the only-3-mode option, i.e. we extract the features corresponding to modes
with the largest three sizes.

Figure 3 shows the total feature extraction time for each method and each
tensor. The tensors are shown in increasing order of nonzero counts. Excluding
the hybrid method, we observe that the grouping-based method is superior for
smaller tensors, whereas the sort-based method is better for larger tensors.
Overall, the hybrid method is the best-performing one for all cases except
some small tensors, for which the grouping-based method is the best.

Figure 3. Runtime comparison for different feature extraction methods.

A more detailed mode-level comparison of the sorting- and grouping-based
methods is given in Table 4. It is seen that the superiority of a method
varies between different modes within the execution of each tensor. Although
the grouping-based method tends to be costly in total time due to some modes
with larger sizes, it is still the winner for some modes of even large
tensors. This is the main reason why integrating the grouping-based method for
specific modes within the hybrid approach yields improved performance.

Table 4. Comparison of sort- and group-based methods in terms of mode-level
runtime in seconds. Best results for each mode are shown in bold.

| Sort-based | Group-based
---|---|---
Tensor | mode-1 | mode-2 | mode-3 | mode-1 | mode-2 | mode-3
lbnl-network | 0.6 | 0.4 | 0.3 | 0.3 | 0.6 | 0.6
nips | 1.0 | 0.7 | 0.7 | 0.1 | 0.1 | 0.2
uber | 1.1 | 0.7 | 0.7 | 0.1 | 0.1 | 0.1
chicago-crime-comm | 1.8 | 1.2 | 1.2 | 0.1 | 0.1 | 0.1
chicago-crime-geo | 2.1 | 1.3 | 1.3 | 0.2 | 0.2 | 0.2
vast-2015-mc1-3d | 8.3 | 5.7 | 5.2 | 0.6 | 0.8 | 0.9
vast-2015-mc1-5d | 8.0 | 5.4 | 5.4 | 0.8 | 0.8 | 1.0
DARPA1998 | 8.9 | 6.2 | 6.7 | 1.3 | 103.1 | 138.0
enron | 16.7 | 14.6 | 12.4 | 0.7 | 0.8 | 1.2
NELL-2 | 26.0 | 20.8 | 19.1 | 0.8 | 0.9 | 1.5
freebase_music | 34.4 | 28.6 | 22.0 | 3.4 | 16.7 | 36,000.0
flickr-3d | 38.2 | 29.8 | 24.6 | 318.7 | 4,846.7 | 13,158.4
flickr-4d | 38.9 | 39.1 | 24.4 | 326.7 | 4,798.0 | 13,109.2
freebase_sampled | 50.6 | 42.6 | 32.0 | 8.3 | 28.4 | 36,000.0
delicious-3d | 51.1 | 39.6 | 31.4 | 1,209.3 | 4,840.3 | 13,517.2
NELL-1 | 47.4 | 38.3 | 33.0 | 4,346.7 | 21,478.6 | 25,526.0

###  5.3. Performance of Tensor Generator

We evaluate the effectiveness of the proposed tensor generator in terms of
feature quality, tensor decomposition performance, and sensitivity to seed
selection. The generated tensors are created by GenTen using the features of
real tensors, which are obtained via our feature extractor FeaTen. Except for
the sensitivity analysis part, we use the generated tensors in which the seed
for pseudo-randomness is set to 0.

####  5.3.1. Feature Quality

Table 5 shows the comparison of the generated tensors with their original
versions, i.e. real tensors, in terms of some important features. We present
the features of the original and the respective generated tensors, as well as
the ratio of the resulting value of the generated tensor over that of the
original tensor. The ratio values are colored green if the value is between
0.9 and 1.1; red if it is less than 0.5 or more than 2; and orange, otherwise.
For zero or too small (less than 0.1) coefficient of variation values, the
generator often yields values with zero coefficient of variation. For those
cases, the ratio values are omitted from the table since the ratio will appear
as either zero or undefined (0/0).

As can be seen in the table, the resulting densities, i.e. the nonzero count
of the generated tensors, are at least 0.96 times smaller or at most 1.05
times larger than the ones of the respective original tensors. The success of
the generator in terms of obeying the given density is seen in both levels of
nonzero slice, nonzero fiber, and nonzero density.

In methods that generate values to obey a given density and variation, there
is a trade-off between strictly obeying the density or the variation. In other
words, if one prioritizes achieving the exact coefficient of variation, the
density might get far from the desired value. However, since the number of
nonzeros is the most significant feature that a generator must obey for
performance concerns, we opted to prioritize adhering to density. For this, in
GenTen we apply some scalings during calculations to catch the given density.
It is the reason why the resulting ratios in coefficient of variation seem
relatively low compared to the ratios in densities.

We also present the tensor generation times in the last column of Table 5.
These are the runtime of GenTen in seconds when working with the maximum
available thread count, which is 96 in our case. We observe that GenTen takes
only a few seconds when generating medium-size tensors. To be precise, it
takes less than 5 seconds for 9 out of 16 tensors. Only for some large tensors
(4 cases in our dataset), the runtime of GenTen can go up to a few hours. Note
that the execution of GenTen depends on the requested nonzero count as well as
the slice and fiber counts; thus the sizes of the tensor also affect the
execution time.

Table 5.  Feature comparison of the generated tensors with their original
versions; and time to generate tensors.

| Coefficient of Variation | Density |
---|---|---|---
| Fiber per Slice | Nonzero per Fiber | Nonzero Slice | Nonzero Fiber | Nonzero |
Name | Org | Gen | Ratio | Org | Gen | Ratio | Org | Gen | Ratio | Org | Gen | Ratio | Org | Gen | Ratio | Time (s)
LBNL-network | 8.0 | 4.5 |  0.56 | 26.2 | 12.1 | 0.46 | 2.1E-06 | 2.1E-06 | 1.00 | 8.4E-10 | 9.5E-10 | 1.13 | 4.2E-14 | 4.4E-14 | 1.05 | 1.0
NIPS | 0.2 | 0.2 | 0.98 | 0.0 | 0.0 | - | 8.3E-04 | 8.3E-04 |  1.00 | 3.1E-05 | 3.0E-05 | 0.98 | 1.8E-06 | 1.8E-06 | 0.98 | 0.3
uber | 0.2 | 0.2 | 0.94 | 1.0 | 1.0 | 0.94 | 1.0E+00 | 1.0E+00 | 1.00 | 1.4E-01 | 1.3E-01 | 0.93 | 3.8E-04 | 3.8E-04 | 0.99 | 0.1
chicago-crime-comm | 0.4 | 0.3 | 0.82 | 0.5 | 0.4 | 0.82 | 9.5E-01 | 9.5E-01 | 1.00 | 3.2E-01 | 2.7E-01 | 0.83 | 1.5E-02 | 1.4E-02 | 0.96 | 0.3
chicago-crime-geo | 0.3 | 0.3 | 0.94 | 0.1 | 0.0 | - | 1.0E-01 | 1.0E-01 |  1.00 | 2.8E-04 | 2.9E-04 | 1.01 | 8.9E-06 | 8.9E-06 | 1.01 | 1.6
vast-2015-mc1-3d | 0.5 | 0.5 | 1.00 | 0.0 | 0.0 | - | 1.0E+00 | 1.0E+00 | 1.00 | 1.4E-02 | 1.4E-02 | 0.99 | 6.9E-03 | 6.9E-03 | 0.99 | 2.2
vast-2015-mc1-5d | 0.0 | 0.0 | - | 0.0 | 0.0 | - | 6.9E-03 | 6.9E-03 | 1.00 | 6.9E-05 | 6.9E-05 | 1.00 | 7.8E-07 | 7.8E-07 | 1.00 | 4.8
DARPA1998 | 13.1 | 8.2 | 0.63 | 23.1 | 14.0 | 0.61 | 8.0E-01 | 8.1E-01 |  1.02 | 1.5E-04 | 1.6E-04 | 1.03 | 2.4E-09 | 2.4E-09 | 1.00 | 147.8
enron | 4.1 | 3.6 | 0.87 | 1.8 | 1.4 | 0.76 | 4.4E-03 | 4.4E-03 | 1.00 | 3.7E-06 | 3.7E-06 | 0.99 | 5.5E-09 | 5.7E-09 |  1.05 | 3.1
NELL-2 | 3.3 | 3.1 | 0.94 | 0.9 | 1.1 | 1.25 | 1.0E+00 | 1.0E+00 | 1.00 | 3.0E-03 | 3.1E-03 | 1.01 | 2.4E-05 | 2.4E-05 | 0.99 | 0.6
freebase_music | 24.4 | 20.2 | 0.83 | 0.1 | 0.0 |  - | 9.7E-01 | 1.0E+00 | 1.03 | 1.8E-07 | 1.9E-07 | 1.04 | 1.1E-09 | 1.1E-09 | 1.03 | 4,370.6
flickr-3d | 3.3 | 3.2 | 0.97 | 1.0 | 1.0 | 0.99 | 1.0E+00 | 1.0E+00 | 1.00 | 3.1E-06 | 3.1E-06 | 1.00 | 7.8E-12 | 7.9E-12 |  1.01 | 2,080.4
flickr-4d | 1.0 | 1.0 | 0.99 | 0.0 | 0.0 | - | 3.1E-06 | 3.1E-06 | 1.00 | 7.8E-12 | 7.9E-12 | 1.01 | 1.1E-14 | 1.1E-14 | 1.01 | 1,180.8
freebase_sampled | 24.0 | 19.4 | 0.81 | 0.1 | 0.0 |  - | 9.1E-01 | 9.1E-01 | 1.00 | 9.2E-08 | 9.6E-08 | 1.05 | 1.7E-10 | 1.8E-10 | 1.04 | 10,608.6
delicious-3d | 2.8 | 2.7 | 0.99 | 1.4 | 1.0 | 0.71 | 1.0E+00 | 1.0E+00 | 1.00 | 5.1E-06 | 5.1E-06 | 1.00 | 6.1E-12 | 6.1E-12 |  1.00 | 6,332.0
NELL-1 | 13.6 | 10.8 | 0.80 | 7.5 | 4.5 | 0.60 | 1.0E+00 | 1.0E+00 | 1.00 | 2.8E-06 | 2.8E-06 | 1.01 | 9.1E-13 | 9.2E-13 | 1.01 | 23,384.5

####  5.3.2. CPD Performance

We evaluate the effectiveness of the generated tensors in mimicking the
behavior of real tensors, particularly in terms of tensor decomposition
performance. For this, the performance of the generated tensors is compared
with the performance of the naive random tensors, which have the same sizes
and nonzero counts as the original tensors but the nonzero locations are
uniformly random. The SPLATT (Smith et al., 2015) and ParTI (Li et al., 2018a)
tools are used for applying the CPD decomposition for 50 iterations, by
ignoring the first 5 iterations for a warm cache. For each case, we take 5
independent runs and choose the minimum runtime to represent the peak
performance of the system and to be less susceptible to noise than the
average.

Figure 4 illustrates the comparison between the generated and naive random
tensors in their ability to resemble real tensors regarding the CPD
performance. The runtime results for both the naive random and the generated
tensors are normalized with respect to the runtime obtained for the respective
original tensor. Therefore, the normalized values closer to 1.0 are
interpreted as more successful in terms of resembling the original tensor
performance.

(a) ParTI performance comparison for a single thread

(b) SPLATT performance comparison for a single thread

(c) SPLATT performance comparison for 48 threads

Figure 4.  CPD performance comparison of the naive random and the generated
tensors. The values are normalized with respect to the runtime obtained for
the respective original tensor.

Figure 4(a) depicts the CPD performance comparison by using the ParTI tool
with a single thread. It is evident that the generated tensors emerge as the
clear winner in most cases, and in the remaining scenarios, they are
comparable to the naive random tensors but never inferior. The superiority of
the generated tensors over the naive random ones is especially higher for
larger tensors. The ParTI tool gave an error for the chicago-crime-geo tensor,
so this tensor is not presented in Figure 4(a) and the rest of the experiments
are conducted by using only the SPLATT tool.

In Figure 4(b), we present the CPD performance comparison by using the SPLATT
tool with a single thread. As seen in the figure, the generated tensors show
significantly better performance than the naive random ones for 11 out of 16
cases; and yield similar performance with naive random tensors for 3 tensors
in the serial setting. For instance, for the NIPS tensor, while the naive
random tensor is 3.3 times slower than the original tensor, the tensor
generated with GenTen yields a CPD time of only 1.1 times more than the
runtime of the original real tensor.

Figure 4(c) shows the CPD performance comparison by using SPLATT with 48
threads. We observe that the generated tensors are superior to the naive
random ones for 8 cases; and inferior for only 1 tensor (chicago-crime-geo).
The generated tensors show almost the same performance as the naive random
ones for 3 cases but for those, their performance is either already the same
as the respective real tensor (chicago-crime-comm), or only 1.3
(vast-2015-mc1-5d) and 1.5 (vast-2015-mc1-3d) times far from the performance
of the original tensors.

####  5.3.3. Sensitivity Analysis

To evaluate the sensitivity of our generator, we produce 3 different versions
for each generated tensor by taking the seed for pseudo-randomness as 0, 1,
and 2. Table 6 shows the nonzero counts and the CPD (SPLATT) runtime for 1 and
48 threads (in seconds). We present the coefficient of variation (CV) for each
case. As can be seen in the table, the CV for nonzero count is at most
$5\times 10^{-3}$, whereas the CV for CPD performance is at most 0.1. That is,
the generator is robust in terms of maintaining a consistent sparsity pattern,
and the performance of the generated tensor is not highly sensitive to changes
in the seed value. It implies that the proposed generator produces reliable
and predictable results, providing a dependable tool for generating sparse
tensors for various applications.

Table 6.  Comparison of the generated tensors with different seeds in terms of
nonzero count and CPD time in seconds.

| Nonzero Count | CPD Time (1 Thread) | CPD Time (48 Threads)
---|---|---|---
Name | Seed=0 | Seed=1 | Seed=2 | CV | Seed=0 | Seed=1 | Seed=2 | CV | Seed=0 | Seed=1 | Seed=2 | CV
LBNL-network | 1,699,088 | 1,699,162 | 1,699,444 | 9.0E-05 | 21.3 | 21.3 | 22.1 | 1.8E-02 | 3.9 | 5.0 | 4.9 | 1.1E-01
NIPS | 3,041,941 | 3,042,054 | 3,041,816 | 3.2E-05 | 6.1 | 6.1 | 6.7 | 3.9E-02 | 0.5 | 0.5 | 0.6 | 4.2E-02
uber | 3,283,106 | 3,282,919 | 3,283,336 | 5.2E-05 | 6.4 | 6.3 | 6.3 | 6.1E-03 | 0.3 | 0.3 | 0.3 | 0.0E+00
chicago-crime-comm | 5,131,130 | 5,131,402 | 5,131,234 | 2.2E-05 | 8.0 | 7.7 | 7.9 | 1.5E-02 | 0.4 | 0.4 | 0.4 | 7.0E-03
chicago-crime-geo | 6,364,282 | 6,364,306 | 6,365,068 | 5.7E-05 | 75.3 | 100.8 | 90.9 | 1.2E-01 | 7.9 | 9.8 | 8.8 | 8.9E-02
vast-2015-mc1-3d | 25,800,284 | 25,802,581 | 25,801,506 | 3.6E-05 | 64.1 | 67.1 | 70.4 | 3.8E-02 | 29.1 | 31.3 | 29.7 | 3.1E-02
vast-2015-mc1-5d | 26,021,841 | 26,021,841 | 26,021,841 | 0.0E+00 | 268.9 | 269.6 | 268.9 | 1.2E-03 | 218.4 | 214.2 | 212.2 | 1.2E-02
DARPA1998 | 28,732,452 | 28,406,813 | 28,405,526 | 5.4E-03 | 250.7 | 251.5 | 269.8 | 3.4E-02 | 84.6 | 94.6 | 93.0 | 4.8E-02
enron | 56,967,347 | 57,055,466 | 57,032,145 | 6.5E-04 | 588.9 | 670.4 | 651.7 | 5.5E-02 | 29.8 | 32.7 | 32.7 | 4.3E-02
NELL-2 | 76,203,551 | 76,202,209 | 76,203,884 | 9.5E-06 | 89.8 | 89.5 | 90.1 | 2.6E-03 | 7.7 | 7.9 | 7.8 | 1.4E-02
freebase_music | 102,877,680 | 102,916,007 | 102,874,661 | 1.8E-04 | 1,459.2 | 1,495.5 | 1,480.7 | 1.0E-02 | 243.3 | 316.4 | 312.1 | 1.2E-01
flickr-3d | 113,727,282 | 113,736,842 | 113,674,701 | 2.4E-04 | 1,471.7 | 1,458.7 | 1,544.7 | 2.5E-02 | 261.1 | 292.3 | 308.5 | 6.8E-02
flickr-4d | 113,641,027 | 113,641,445 | 113,640,707 | 2.7E-06 | 4,714.5 | 4,786.8 | 4,700.2 | 8.0E-03 | 596.7 | 579.3 | 627.7 | 3.3E-02
freebase_sampled | 145,748,228 | 145,801,470 | 145,741,097 | 1.8E-04 | 2,409.9 | 2,424.5 | 2,493.8 | 1.5E-02 | 380.6 | 465.8 | 482.4 | 1.0E-01
delicious-3d | 139,996,876 | 139,993,400 | 139,999,596 | 1.8E-05 | 1,845.9 | 1,835.4 | 1,833.6 | 3.0E-03 | 290.3 | 350.6 | 328.5 | 7.7E-02
NELL-1 | 145,175,629 | 145,155,883 | 145,155,081 | 6.5E-05 | 1,513.9 | 1,599.3 | 1,609.8 | 2.7E-02 | 256.0 | 300.6 | 292.8 | 6.9E-02

##  6\. Related Work

To the best of our knowledge, there is only one study (and its extension) that
extracts tensor features for optimizing sparse tensor computations in the
literature: Sun et al.(Sun et al., 2020) proposed a framework, namely SpTFS,
that automatically predicts the optimal storage format for CPD. SpTFS lowers
the sparse tensor to fixed-sized matrices and gives them to convolutional
neural networks (CNN) as inputs along with tensor features. The authors
improve SpTFS by adopting both supervised and unsupervised learning-based
machine learning models in a recent work (Sun et al., 2021).

The previous works (Sun et al., 2020, 2021) have considered the features for
only one mode, assuming that the tensor is already sorted along that mode.
However, real-world tensor sizes diverge significantly so that sizes in some
modes reach millions while some are only orders of ten or even less.
Therefore, considering the global values might result in losing some important
information about the structure of the tensor.

Although several sparse matrix and graph generators are proposed in the
literature (Luszczek et al., 2020; Mpakos et al., 2023; Kolda et al., 2014;
Armstrong and Rendell, 2008; Bonifati et al., 2020), the studies on sparse
tensor generators are very limited. Smith and Karypis (Smith and Karypis,
2017) used the method of Chi and Kolda (Chi and Kolda, 2012) to generate
synthetic tensors whose values follow a Poisson distribution. Due to the
deficiency of publicly available sparse tensors, Sun et al. (Sun et al., 2020)
produced synthetic tensors by combining sparse matrices in Suite Sparse (Davis
and Hu, 2011) collection such that the elements of matrices form the modes of
tensors. Baskaran et al. (Baskaran et al., 2012) used MATLAB Tensor Toolkit
(Kolda and Bader, 2006) to generate synthetic sparse tensors but these are
rather small tensors with less than one million nonzeros.

##  7\. Conclusion

In this study, we introduce two tools, namely FeaTen and GenTen, that we
proposed and developed to facilitate the research in sparse tensor operations.
FeaTen is a feature extraction framework for sparse tensors that provides four
different methods, with a focus on extracting the tensor features in the most
efficient way. This contribution is valuable because even the feature
extraction is challenging and costly due to huge fiber and slice counts in
real sparse tensors. We provide an evaluation of the performance of different
feature extraction methods and observe that the methods whose algorithms are
proposed first in this work are preferable to the conventional methods.

GenTen is a smart sparse tensor generator that obeys a significant and large
set of tensor features. The experimental results validate the success of the
generator in mimicking the real tensors regarding both sparsity pattern and
tensor decomposition performance. GenTen holds the benefit of utilizing size-
independent features so that it can be used to generate tensors with smaller
or larger sizes, yet still have similar properties with real tensors. It paves
the way for obtaining a large dataset of synthetic sparse tensors that have
similar characteristics and behavior with real-world data.

###### Acknowledgements.

This project has received funding from the European High-Performance Computing
Joint Undertaking under grant agreement No 956213. Koç University is supported
by the Turkish Science and Technology Research Centre Grant No 120N003.

## References

  * (1)
  * Abu-Sufah and Abdel Karim (2013) Walid Abu-Sufah and Asma Abdel Karim. 2013.  Auto-tuning of sparse matrix-vector multiplication on graphics processors. In _International Supercomputing Conference_. Springer, 151–164.
  * Acar et al. (2007) Evrim Acar, Canan Aykut-Bingol, Haluk Bingol, Rasmus Bro, and Bülent Yener. 2007.  Multiway analysis of epilepsy tensors.  _Bioinformatics_ 23, 13 (2007), i10–i18.
  * Acer et al. (2018) Seher Acer, Tugba Torun, and Cevdet Aykanat. 2018.  Improving medium-grain partitioning for scalable sparse tensor decomposition.  _IEEE Transactions on Parallel and Distributed Systems_ 29, 12 (2018), 2814–2825.
  * Armstrong and Rendell (2008) Warren Armstrong and Alistair P Rendell. 2008.  Reinforcement learning for automated performance tuning: Initial evaluation for sparse matrix format selection. In _2008 IEEE International Conference on Cluster Computing_. IEEE, 411–420.
  * Baskaran et al. (2012) Muthu Baskaran, Benoît Meister, Nicolas Vasilache, and Richard Lethin. 2012.  Efficient and scalable computations with sparse tensors. In _2012 IEEE Conference on High Performance Extreme Computing_. IEEE, 1–6.
  * Bonifati et al. (2020) Angela Bonifati, Irena Holubová, Arnau Prat-Pérez, and Sherif Sakr. 2020.  Graph generators: State of the art and open challenges.  _ACM computing surveys (CSUR)_ 53, 2 (2020), 1–30.
  * Chi and Kolda (2012) Eric C Chi and Tamara G Kolda. 2012.  On tensors, sparsity, and nonnegative factorizations.  _SIAM J. Matrix Anal. Appl._ 33, 4 (2012), 1272–1299.
  * Davis and Hu (2011) Timothy Davis and Yifan Hu. 2011.  The University of Florida Sparse Matrix Collection.  _ACM Trans. Math. Software_ 38, 1 (2011), 1–25.  https://doi.org/10.1145/2049662.2049663
  * Golder and Settle (1976) ER Golder and JG Settle. 1976.  The Box-Müller Method for Generating Pseudo-Random Normal Deviates.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ 25, 1 (1976), 12–20.
  * Goulart et al. (2015) José Henrique de M Goulart, Maxime Boizard, Rémy Boyer, Gérard Favier, and Pierre Comon. 2015.  Tensor CP decomposition with structured factor matrices: Algorithms and performance.  _IEEE Journal of Selected Topics in Signal Processing_ 10, 4 (2015), 757–769.
  * Helal et al. (2021) Ahmed E Helal, Jan Laukemann, Fabio Checconi, Jesmin Jahan Tithi, Teresa Ranadive, Fabrizio Petrini, and Jeewhan Choi. 2021.  Alto: Adaptive linearized storage of sparse tensors. In _Proceedings of the ACM International Conference on Supercomputing_. 404–416.
  * Jeon et al. (2015) Inah Jeon, Evangelos E Papalexakis, Uksong Kang, and Christos Faloutsos. 2015.  Haten2: Billion-scale tensor decompositions. In _2015 IEEE 31st international conference on data engineering_. IEEE, 1047–1058.
  * Jukić et al. (2013) Ante Jukić, Ivica Kopriva, and Andrzej Cichocki. 2013.  Noninvasive diagnosis of melanoma with tensor decomposition-based feature extraction from clinical color image.  _Biomedical Signal Processing and Control_ 8, 6 (2013), 755–763.
  * Kaya and Uçar (2018) Oguz Kaya and Bora Uçar. 2018.  Parallel candecomp/parafac decomposition of sparse tensors using dimension trees.  _SIAM Journal on Scientific Computing_ 40, 1 (2018), C99–C130.
  * Kolda and Bader (2006) Tamara G Kolda and Brett W Bader. 2006.  _MATLAB tensor toolbox_.  Technical Report. Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA (United States).
  * Kolda and Bader (2009) Tamara G Kolda and Brett W Bader. 2009.  Tensor decompositions and applications.  _SIAM review_ 51, 3 (2009), 455–500.
  * Kolda et al. (2014) Tamara G Kolda, Ali Pinar, Todd Plantenga, and Comandur Seshadhri. 2014.  A scalable generative graph model with community structure.  _SIAM Journal on Scientific Computing_ 36, 5 (2014), C424–C452.
  * Kuang et al. (2014) Liwei Kuang, Fei Hao, Laurence T Yang, Man Lin, Changqing Luo, and Geyong Min. 2014.  A tensor-based approach for big data representation and dimensionality reduction.  _IEEE transactions on emerging topics in computing_ 2, 3 (2014), 280–291.
  * Leclercq and Savonnet (2018) Éric Leclercq and Marinette Savonnet. 2018.  A tensor based data model for polystore: an application to social networks data. In _Proceedings of the 22nd International Database Engineering & Applications Symposium_. 110–118.
  * Li et al. (2018a) Jiajia Li, Yuchen Ma, and Richard Vuduc. 2018a.  ParTI!: A parallel tensor infrastructure for multicore CPUs and GPUs.
  * Li et al. (2018b) Jiajia Li, Jimeng Sun, and Richard Vuduc. 2018b.  HiCOO: Hierarchical storage of sparse tensors. In _SC18: International Conference for High Performance Computing, Networking, Storage and Analysis_. IEEE, 238–252.
  * Li et al. (2019) Jiajia Li, Bora Uçar, Ümit V Çatalyürek, Jimeng Sun, Kevin Barker, and Richard Vuduc. 2019.  Efficient and effective sparse tensor reordering. In _Proceedings of the ACM International Conference on Supercomputing_. 227–237.
  * Luszczek et al. (2020) Piotr Luszczek, Yaohung Tsai, Neil Lindquist, Hartwig Anzt, and Jack Dongarra. 2020.  Scalable Data Generation for Evaluating Mixed-Precision Solvers. In _2020 IEEE High Performance Extreme Computing Conference (HPEC)_. IEEE, 1–6.
  * Mpakos et al. (2023) Panagiotis Mpakos, Dimitrios Galanopoulos, Petros Anastasiadis, Nikela Papadopoulou, Nectarios Koziris, and Georgios Goumas. 2023.  Feature-based spmv performance analysis on contemporary devices. In _2023 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 668–679.
  * Mu et al. (2011) Yang Mu, Wei Ding, Melissa Morabito, and Dacheng Tao. 2011.  Empirical discriminative tensor analysis for crime forecasting. In _Knowledge Science, Engineering and Management: 5th International Conference, KSEM 2011, Irvine, CA, USA, December 12-14, 2011. Proceedings 5_. Springer, 293–304.
  * Nisa et al. (2019) Israt Nisa, Jiajia Li, Aravind Sukumaran-Rajam, Richard Vuduc, and P Sadayappan. 2019.  Load-balanced sparse MTTKRP on GPUs. In _2019 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 123–133.
  * Shivakumar et al. (2021) Shruti Shivakumar, Jiajia Li, Ramakrishnan Kannan, and Srinivas Aluru. 2021.  Efficient parallel sparse symmetric tucker decomposition for high-order tensors. In _SIAM Conference on Applied and Computational Discrete Algorithms (ACDA21)_. SIAM, 193–204.
  * Smith et al. (2017) Shaden Smith, Jee W Choi, Jiajia Li, Richard Vuduc, Jongsoo Park, Xing Liu, and George Karypis. 2017.  FROSTT: The formidable repository of open sparse tensors and tools.
  * Smith and Karypis (2015) Shaden Smith and George Karypis. 2015.  Tensor-matrix products with a compressed sparse tensor. In _Proceedings of the 5th Workshop on Irregular Applications: Architectures and Algorithms_. 1–7.
  * Smith and Karypis (2017) Shaden Smith and George Karypis. 2017.  Accelerating the tucker decomposition with compressed sparse tensors. In _European Conference on Parallel Processing_. Springer, 653–668.
  * Smith et al. (2015) Shaden Smith, Niranjay Ravindran, Nicholas D Sidiropoulos, and George Karypis. 2015.  SPLATT: Efficient and parallel sparse tensor-matrix multiplication. In _2015 IEEE International Parallel and Distributed Processing Symposium_. IEEE, 61–70.
  * Sun et al. (2020) Qingxiao Sun, Yi Liu, Ming Dun, Hailong Yang, Zhongzhi Luan, Lin Gan, Guangwen Yang, and Depei Qian. 2020.  Sptfs: Sparse tensor format selection for mttkrp via deep learning. In _SC20: International Conference for High Performance Computing, Networking, Storage and Analysis_. IEEE, 1–14.
  * Sun et al. (2021) Qingxiao Sun, Yi Liu, Hailong Yang, Ming Dun, Zhongzhi Luan, Lin Gan, Guangwen Yang, and Depei Qian. 2021.  Input-aware Sparse Tensor Storage Format Selection for Optimizing MTTKRP.  _IEEE Trans. Comput._ (2021).
  * Yelundur et al. (2019) Anil R Yelundur, Vineet Chaoji, and Bamdev Mishra. 2019.  Detection of review abuse via semi-supervised binary multi-target tensor decomposition. In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_. 2134–2144.