Why PMI?: The intuition behind computing $\mathbf{S}$ using PMI comes from the observation that sometimes, in the co-clustering of user-product data, a user cluster is associated with another product cluster that does not exist. This leads to the idea of an update rule based on co-occurrence between users and products for all cluster pairs in the given data. PMI is an information-theoretic approach that measures how often two clusters ($\mathbf{A}_m,\mathbf{B}_n$) occur as compared with what we expect if they were independent. The numerator of PMI informs us how often we observed the two clusters together in user-product context consumption. The denominator informs us how often we would anticipate both to co-cluster, assuming they are independent clusters. Thus, the ratio provides us an estimation of how much more the two clusters co-cluster than we anticipate by chance. Our choice for PMI is encouraged by <cit.>. Using PMI, we encourage users of similar interactions with products to have closer representation in latent space. For example, when representing co-clusters, we can easily think of positive relationship (e.g. “Female” and “Nurse”) but find it much harder to relate negative ones (“Female” and “Carpenter”). Therefore, we focus on positive point-wise mutual information. Next, Equ. (<ref>) and (<ref>) can be re-written as:

\begin{equation}
\label{eq:lossfunction_A1}
\mathbf{A} = \mathbf{X}\mathbf{B} \log\Big( \frac{\mathbf{A}^T\mathbf{X}\mathbf{B}}{\sum_{i=1}^I \mathbf{A}^T\mathbf{X} \sum_{i=1}^I \mathbf{X}\mathbf{B}}\Big)^T
\end{equation}

\begin{equation}
\label{eq:lossfunction_B1}
\mathbf{B} = \mathbf{X}^T\mathbf{A} \log\Big( \frac{\mathbf{A}^T\mathbf{X}\mathbf{B}}{\sum_{j=1}^J \mathbf{A}^T\mathbf{X} \sum_{i=1}^I \mathbf{X}\mathbf{B}}\Big)
\end{equation}

Now, we assign each user/product to a single cluster by finding the cluster with maximum membership. This translates to finding the maximum column index for each row.

§.§ Step 2: Co-Cluster Explanation

After discovering outstanding co-clusters, we next aim to identify the user and product feature subsets that best explain the co-clusters. The process involves two steps: auxiliary feature matrix creation to infer the user preferences over product features to get insights about implicit user similarities, and feature selection using point wise mutual information to compute the association between features and user/product cluster centroids. Through experiments we show that the auxiliary feature matrix helps to improve the explainability of our proposed method (see Section <ref>). The quest for similarities plays an important role in co-cluster analysis <cit.>. In quadripartite graphs (see Figure <ref>), one can derive several semantics on similarity by considering different paths in a graph. Upon deriving these, we have all information to learn important co-cluster features: the final suitability score for each feature is captured as a linear combination of both steps. We briefly explain each step as follow:

§.§.§ Auxiliary Feature Matrix

We have the user and product feature vectors ($\mathbf{U}_f \in \mathbb{R}^{I \times F_1}$ and $\mathbf{P}_f \in \mathbb{R}^{J \times F_2} $) which are independent and do not infer any user preferences for the features that appear in products and vice-versa. This could lead to undermining the user/product similarities for feature learning. To overcome this issue, we compute the user's proximity to a product feature through a meta-path as an indication of the user’s possible “preference” towards the product feature and vice versa. The meta-path <cit.> is a powerful mechanism for a user to select an appropriate similarity semantics to learn from a set of examples of similar objects. Formally, a meta-path can be defined as:

A meta-path is a sequence of relations $\mathcal{R}$ between object types $\mathcal{O}$, which defines a new composite relation between its starting type and ending type. It is denoted in the form of $\mathcal{O}_1 \xrightarrow[]{\mathcal{R}_1} \mathcal{O}_2 \xrightarrow[]{\mathcal{R}_2}  \dots \xrightarrow[]{\mathcal{R}_{m-1}} \mathcal{O}_m$, which defines a composite relation  between $\mathcal{R}_1 \circ \mathcal{R}_2 \circ \dots  \circ \mathcal{R}_m $ between types $\mathcal{O}_1$ and $\mathcal{O}_m$, where $\circ$ denotes the composition operator on relation. For example, a meta-path (Figure <ref>, purple path) user (male) $\rightarrow$ product (movie) $\rightarrow$ feature (action) $\rightarrow$ product (movie) $\rightarrow$ feature (comedy) (denoted as UPFPF) indicates user preferences based on content similarities.

The process of auxiliary feature matrix creation is adapted from <cit.>. In <cit.>, movie preferences are learned by leveraging user similarities defined through different types of meta paths or relations to successfully design a new movie. Here, the auxiliary feature matrix helps to leverage both explicit features and user/product similarities via a graph-theoretic approach. For user auxiliary feature matrix as shown in Figure <ref>, the red path $\mathbf{U}_{p_f}^{UPF}$ (i.e., starting from a user and ending on a product feature via a product) finds the preferences for the product features for each user based on its interactions. The blue path $\mathbf{U}_{p_f}^{UPUPF}$ finds user preferences for the product features based on user-user similarity, and the purple path $\mathbf{U}_{p_f}^{UPFPF}$ finds user preferences based on product-product (content) similarity.

Quadripartite graph of the users $\mathbf{U} \in \mathbb{R}^{I}$, Product $\mathbf{P} \in \mathbb{R}^{J}$ and the user features $\mathbf{U}_f \in \mathbb{R}^{I \times F_1}$ and product features $\mathbf{P}_f \in \mathbb{R}^{J \times F_2}$. The resultant user  auxiliary feature matrices is $\mathbf{U}_{p_f} \in \mathbb{R}^{I \times F_2}$. V: can we make this figure larger? There is a lot of whitespace that probably can be covered without the figure taking up any extra space (but also keeping the aspect ratio)

For the final weighted  matrix $\mathbf{U}_{p_f} \in \mathbb{R}^{I \times F_2}$, which represents user $\mathbf{U}$ preference over product features $\mathbf{P}_f$ is a linear combination of the weighted $\mathbf{U}_{p_f}$ over the three predefined meta-path types:

\begin{equation}
\mathbf{U}_{p_f} = \alpha \mathbf{U}_{p_f}^{UPF} + \beta \mathbf{U}_{p_f}^{UPUPF} + \gamma \mathbf{U}_{p_f}^{UPFPF}
\label{eqn:metapath_upf}
\end{equation}

This linear combination helps to smooth the information in case of sparse direct user-product preferences. Similarly, we compute weighted auxiliary matrix for product features as:
to obtain $\mathbf{P}_{u_f} \in \mathbb{R}^{J \times F_1}$ (2nd figure), red path  $\mathbf{P}_{u_f}^{1}$ finds preferences for the user features for each product based its interaction, blue path $\mathbf{P}_{u_f}^{2}$ finds product preferences for the user features based product-product similarity and purple path $\mathbf{P}_{u_f}^{3}$ finds product preferences based on  user-user similarity. should explain these choices for parameters.

\begin{equation}
\mathbf{P}_{u_f} = \alpha \mathbf{P}_{u_f}^{PUF} + \beta \mathbf{P}_{u_f}^{PUPUF} + \gamma \mathbf{P}_{u_f}^{PUFUF}
\label{eqn:metapath_puf}
\end{equation}

where $\alpha, \beta, \gamma \in \mathbb{R}^{+}$ are combination parameters satisfying criterion $\alpha + \beta + \gamma = 1$. We set $\alpha = 0.5$, $\beta = 0.25$ and $\gamma = 0.25$ in both cases. We give higher importance to direct user/product preferences ($\alpha$) and lower importance to 4-step path ($\beta, \gamma$) because it could 'obscure' the individual preferences by depending on 'similar' users/products. Section <ref> shows sensitivity analysis for these parameters.
Should we add (maybe in the beginning of this subsection or even here, what is the disadvantage of using $U_f$ and $P_f$ directly instead of those new matrices. I think this currently doesn't come across as strongly as we would want.

§.§.§ Feature Selection

For feature selection, we propose a PMI-based approach that leverages the association between user and product cluster centroids. We compute information about how often we observed the certain features for each user cluster using $\mathbf{A}$ and $\mathbf{U}_f$ as $\frac{(\mathbf{A}^T\mathbf{U}_f)_{ij}}{\sum_{i=1}^I\sum_{j=1}^{J}(\mathbf{A}^T\mathbf{U}_f)_{ij} } \in \mathbb{R}^{M \times F_1} $. Next, we compute information about their association while independently drawn as $p(\mathbf{U}_f) = \frac{\sum_{j=1}^J U_f}{\sum_{i=1}^I\sum_{j=1}^{J}(\mathbf{A}^T\mathbf{U}_f)_{ij} }$ and $p(\mathbf{A}) = \frac{\sum_{i=1}^I A}{\sum_{i=1}^I\sum_{j=1}^{J}(\mathbf{A}^T\mathbf{U}_f)_{ij} }$.

We specifically for user feature learning, obtain matrix $\mathbf{T}$ by multiplying user clustering matrix $\mathbf{A}$ and user feature matrix $\mathbf{U}_f$ as $\mathbf{T} = \mathbf{A}^T\mathbf{U}_f \in \mathbb{R}^{M \times F_1}$. Now, we compute information about how often we observed the certain features for each user cluster by using $p(\mathbf{T}) = \frac{T_{ij}}{\sum_{i=1}^I\sum_{j=1}^{J}T_{ij}} $. Next, we compute information about their association while independently drawn as $p(\mathbf{U}_f) = \frac{\sum_{j=1}^J U_f}{\sum_{i=1}^I\sum_{j=1}^{J}T_{ij}} $ and $p(\mathbf{A}) = \frac{\sum_{i=1}^I A}{\sum_{i=1}^I\sum_{j=1}^{J}T_{ij}}$.

Finally, we compute PMI relation as:

\begin{equation}
\label{sec:pmi_feature}
\mathbf{U}_{PMI}^1 =\log\frac{p( \mathbf{A}^T\mathbf{U}_f)}{p(\mathbf{U}_f)p(\mathbf{A})} = \frac{ \mathbf{A}^T\mathbf{U}_f * \sum_{i=1}^I\sum_{j=1}^{J}(\mathbf{A}^T\mathbf{U}_f)_{ij}}{\sum_{j=1}^J \mathbf{U}_f^\mathbf{T} * \sum_{i=1}^I \mathbf{A}}
\end{equation}

where $\mathbf{U}_{PMI}^1 \in \mathbb{R}^{M \times F_1}$. Now, we compute PMI relation for user auxiliary features matrix $\mathbf{U}_{p_f}$ as:

\begin{equation}
\label{sec:pmi_aux}
\mathbf{U}_{PMI}^2 = \log\frac{p( \mathbf{A}^T\mathbf{U}_{p_f})}{p(\mathbf{U}_{p_f})p(\mathbf{A})} \in \mathbb{R}^{M \times F_2}
\end{equation}

Similarly, we compute both PMIs ($\mathbf{P}_{PMI}^1 \in \mathbb{R}^{N \times F_2}$ and $\mathbf{P}_{PMI}^2 \in \mathbb{R}^{N \times F_1}$) for product features also. Due to space limitations and symmetry, we do not include derivations for them.

To select the most relevant user and product features for the niche, we linearly combine the PMIs associated with co-clustering (i.e. via summary matrix $\mathbf{S}$), attribute matrices (via. Equ. <ref> ), and auxiliary matrices (via. Equ. <ref>) as following:

\begin{equation}
\label{user}
\textbf{For users:} \quad   \mathbf{e}^{u}  = \mathbf{U}_{PMI}^1  + \mathbf{S} * \mathbf{P}_{PMI}^2 \in  \mathbb{R}^{M \times F_1}
\end{equation}

\begin{equation}
\label{prod}
\textbf{For products:} \quad  \mathbf{e}^{p} = \mathbf{P}_{PMI}^1  + \mathbf{S}^T * \mathbf{U}_{PMI}^2 \in \mathbb{R}^{N \times F_2}
\end{equation}

Finally, we choose the top $N$ highest values for each cluster (or each row) from $\mathbf{e}^{u}$ (Equ. <ref>) and $\mathbf{e}^{p}$ (Equ. <ref>) to explain the niche.

§ EXPERIMENTS

In this section, we aim to answer the following research questions:

* RQ1 Accuracy: Can outperform state-of-the-art alternatives at effectively capturing co-clusters?

* RQ2 Explainability: Can help learn meaningful explanations of co-clusters, and thus better niches?

* RQ3 Scalability How efficient and scalable is with respect to the size of the input graphs?

We discuss these after detailing our experimental setup.

§.§ Datasets and Experiment Setup

The details of the synthetic and the real data used for experiments are given in Table <ref>.

§.§.§ Synthetic Data

In order to fully control and evaluate the niches in our experiments, we generate synthetic data i.e. user-product engagement graph $\mathbf{X} \in \mathbb{R}^{I \times J}$, user attributes $\mathbf{U}_f \in \mathbb{R}^{I \times F_1}$ (e.g. age, gender, country, app engagement etc.) and product attributes $\mathbf{P}_f \in \mathbb{R}^{J \times F_2}$   (e.g. name, publisher name, country, category,subcategory etc.) as discussed in Supplementary Material.

2*Dataset   2*#users    2*#products    #features       #clusters

(users, products)      (users, products)

Syn-I   $10K$   $1K$   $(22,43)$   ($14,20$)

Syn-II   $50K$   $5K$   $(22,55)$   ($70,35$)

Syn-III   $500K$   $10K$   $(22,63)$   ($140,50$)

Syn-IV    $1M$   $50K$   $(22,86)$   ($280,70$)

Cora     $3K$   $1.5K$   $-$   ($7,-$)

WebKB4      $4K$   $1K$   $-$   ($4,-$)

MovieLens     $6K$   $4K$   $(25,23)$   ($-,20$)

News20      $19K$   $61K$   $-$   ($20,-$)

Caltech     $2K$   $300$   $-$   ($3,-$)

Snapchat     $500K$   $2.5K$   $(22,238)$   ($-,-$)

Details for the datasets. “-” indicates unknown/unavailable public information.

§.§.§ Real Data

We employ several real-world public datasets from different domains:  Cora is publications dataset,  WebKB4 [<http://membres-lig.imag.fr/grimal/data.html>] consists of classified web page information, MovieLens [<https://grouplens.org/datasets/movielens/>]has data has 6000 users and 4000 movie rating information, News20 [<http://qwone.com/ jason/20Newsgroups/>] is a collection of approximately 20,000 newsgroup documents, Caltech <cit.> is an image dataset. We also use private Snapchat dataset containing full-duration views between users and public feed contents from the Snapchat platform; we assume full views to indicate persisted interest in the consumed content, as in <cit.>.

§.§ Baseline Methods

The two major components of are in (a) the discovering of coherent co-clusters, and (b) their concise explanation via nodal attributes. Thus, we conduct experiments with two categories of baseline to evaluate each contribution separately., demonstrating 's outperformance in both categories.

§.§.§ Co-clustering
We compare with recent, state-of-the-art co-clustering approaches:

* NEO-CC <cit.>: A non-exhaustive overlapping co-clustering method based on the minimum sum-squared residue objective.

* DeepCC <cit.>: A deep autoencoder based co-clustering method which employs a variant of Gaussian Mixture Model (GMM) to infer cluster assignments.

* CoClusInfo <cit.>: An information-theoretic approach which uses mutual information to define its objective function.

* FNMTF <cit.>: A fast, NMTF method based on projected gradients, coordinate descent, and alternating least squares optimization.

* WC-NMTF <cit.>: A word co-occurrence NMTF method that leverages mutual information for co-clustering the word and documents.

* SNCC <cit.>: a sparse neighbor constrained co-clustering with dual regularizers for learning category consistency.

Note that in <cit.>, WC-NMTF method performed better than original NMF <cit.>, orthogonal NMF  (ONMF) <cit.>,  projective NMF (PNMF) <cit.>,  graph regularized NMF (GNMF) <cit.>, NMTF <cit.>, orthogonal NMTF (ONMTF) <cit.> and graph regularized NMTF (GNMTF) <cit.>. Similarly, in <cit.>, paper evaluated performance of SNCC against K-means <cit.>, NMF <cit.>, SNMF <cit.>, graph regularized NMF (GNMF) <cit.>, dual regularization NMTF (DNMTF) <cit.>, dual local learning co-clustering (DLLC) <cit.> and structured optimal bipartite graph (SOBG). Therefore, to avoid the repetitive comparison, we chose to compare our proposed method's performance with SNCC and WC-NMTF.

§.§.§ Explainability
Although we could not find any explicit co-clustering explainability baselines, we adapted the recently proposed LightGBM <cit.> for our explainability baseline. It is a boosting decision tree-based method that employs feature bundling to deal with a large number of features. To use it as our baseline, we fed co-clustering outcomes from the above-discussed method as labels for the user and product clusters along with user/product features data matrices to discern feature importance per co-cluster. We also compared our method with recently proposed method BMGUFS  <cit.> which is rigorous approximation algorithms for block model guided unsupervised feature selection and helps in finding high-quality features for cluster explanation.

Dataset    	Cluster     	Metric    NEO-CC   DeepCC   CoClusInfo   FNMTF   WC-NMTF   SNCC

4*I    2*User    NMI   $ 0.672 \pm 0.052$   $0.771 \pm 0.031$   $0.879 \pm 0.021$   $0.601 \pm 0.023$   $0.824 \pm 0.022$   $0.831 \pm 0.016$   0.948 $\pm$ 0.017

Accuracy     $0.532\pm 0.029$     $0.597\pm 0.034$     $0.781\pm 0.049$     $0.703\pm 0.021$     $0.695 \pm 0.019$    $0.780 \pm 0.021$    0.865 $\pm$ 0.045

2*Product    NMI   $0.772 \pm 0.075$   $ 0.736 \pm 0.002$   $0.886 \pm 0.001$   $0.851 \pm 0.043$   $0.801 \pm 0.032$   $0.854 \pm 0.161$   0.949 $\pm$ 0.096

Accuracy     $0.776\pm 0.105$     $0.651\pm 0.005$     $0.789\pm 0.057$     $0.719\pm 0.001$     $0.683 \pm 0.045$    $0.741 \pm 0.122$    0.865 $\pm$ 0.032

4*II    2*User    NMI   $0.793 \pm 0.062$   $ 0.769 \pm 0.047 $   $0.901  \pm 0.008$    $0.883 \pm 0.041$   $0.746 \pm 0.052$   $0.902 \pm 0.021$   0.960 $\pm$ 0.011

Accuracy     $0.682\pm 0.062$     $0.568\pm 0.075$     $0.794\pm 0.052$     $0.703\pm 0.028$     $0.683\pm 0.119$    $0.761 \pm 0.061$    0.867 $\pm$ 0.038

2*Product    NMI   $0.673 \pm 0.024$   $ 0.743 \pm 0.064$   $0.924 \pm 0.064$   $0.904 \pm 0.058$   $0.839 \pm 0.072$   $0.911 \pm 0.141$   0.947 $\pm$ 0.014

Accuracy     $0.540\pm 0.086$     $0.694\pm 0.109$     $0.799\pm 0.034$     $0.834\pm 0.058$     $0.788\pm 0.001$    $0.801 \pm 0.014$    0.853 $\pm$ 0.042

4*III    2*User    NMI   4*-   4*-   $0.720 \pm 0.011$   $0.913 \pm 0.046$   $0.763 \pm 0.023$   $0.921 \pm 0.046$   0.973 $\pm$ 0.006

Accuracy           $0.694 \pm 0.017$     $0.523\pm 0.034$     $0.492\pm 0.002$    $0.632 \pm 0.021$    0.882 $\pm$ 0.023

2*Product    NMI         $0.840 \pm 0.001$   $0.842 \pm 0.074$   $0.763 \pm 0.028$   $0.856 \pm 0.121$   0.933 $\pm$ 0.019

Accuracy          $0.640 \pm 0.068$     $0.653\pm 0.104$     $0.535\pm 0.112 $   $0.716 \pm 0.214$    0.799 $\pm$ 0.023

4*IV    2*User    NMI   4*-   4*-   $0.879 \pm 0.021$   $0.879 \pm 0.021$   4*-   $0.881 \pm 0.012$   0.959 $\pm$ 0.017

Accuracy             $0.534\pm 0.011$    $0.495\pm 0.133$          $0.563 \pm 0.112$    0.587 $\pm$ 0.072

2*Product    NMI           $0.947 \pm 0.001$    $0.893 \pm 0.129$        $0.891 \pm 0.012$   0.971 $\pm$ 0.007

Accuracy               $0.832 \pm 0.064$     $0.793 \pm 0.101$          $0.801 \pm 0.102$    0.902 $\pm$ 0.032

Experimental results for NMI and Accuracy for synthetic data. Boldface indicates the best results.

The co-clustering result of on synthetic data. Left to right: (a) original synthetic data with $7$ users and $7$ product clusters, (b) shuffled synthetic data, (c) the second-best performing baseline (CoClusInfo) result (91% accurate), and (d)  's result (100% accurate). is successfully able to detect large as well as small sized co-clusters more accurately.

§.§ Evaluation Measures

§.§.§ Co-clustering
We evaluate and the baselines for co-clustering using three criteria: Normalized Mutual Information (NMI), Accuracy and CPU time (sec).

§.§.§ Explainability
We evaluate and the baselines for explainability using three criteria namely Average Precision, Stability score <cit.> and our proposed Compression Score. The MDL based compression score is given by:

\begin{equation}
\begin{aligned}
Score &= \log^{*}r_m+\log^{*}c_n+ \log^{*}uf_m + \log^{*}pf_n + E(B_{mn})  \\
& - r_m \log_2\frac{r_m}{I} -  c_n \log_2\frac{c_n}{J} +  \log_2(r_mc_n +1)  \\
& - uf_m \log_2\frac{uf_m}{F_1} + \log_2(r_muf_m +1) + E(B_{mm}) \\
& - pf_n \log_2\frac{pf_n}{F_2} + \log_2(c_npf_n +1) + E(B_{nn})
\end{aligned}
\end{equation}

where $\log^{*}$ is the universal code length for integers <cit.>, $r_m$ is $m^{th}$ the row cluster encoding bits, $c_n$ is $n^{th}$ the column cluster encoding bits, $uf_m$ is the $m^{th}$ feature cluster encoding bits corresponding to the $m^{th}$ row cluster, $pf_n$ is $n^{th}$ feature cluster encoding bits corresponding to the $n^{th}$ column cluster. $ E(B_{mn})$, $ E(B_{mm})$ and $ E(B_{nn})$ are the number of bits required to encode block of user-product, user-feature and product-feature, respectively. The lower the value, the better is the compression. For all metrics other than CPU time and compression score, higher evaluation measures are advantageous.

§.§ Quantitative Analysis

In this section, we provide quantitative evaluation and analysis of our method for co-clustering and explainability on synthetic and real-world datasets. Table <ref> and <ref> report results on co-clustering performance on synthetic and real data, respectively. Table <ref> reports results on explainability performance.
Recommendation is undoubtedly one key task for businesses. However, explainable and thematic engagement is useful for both individual content creators (influencers, etc.) and companies (Snapchat, Netflix, etc.) who aim to attract certain audiences with original content. The most concrete and evaluate-able components are discovering and describing the co-clusters and explanations found, upon which these creators can leverage findings. Thus, our quantitative evaluation is focused around these points, and we also show associated qualitative end-to-end findings on real data.

§.§.§ Co-clustering Performance

To answer the first experimental question, we report co-clustering  performance of all methods in Table <ref> and <ref>.  We cap the allowed run-time of all methods to 24 hours, indicating unfinished/non-converged results as “-”.

Synthetic Data:  As we can see in Table <ref>, achieves the best performance for all synthetic datasets. As expected, significantly out-performs two related NMTF-based methods, SNCC, FNMTF and WC-NMTF, due to its use of mutual information via the summary matrix as described in Section <ref>.  Moreover, outperforms CoClusInfo, NEO-CC and the neural method DeepCC, in both accuracy and NMI, surpassing them with substantial accuracy improvements ( average $\approx 14\%$). We visualize the co-clustering ability of and second highest performing baseline (CoClusInfo) on an additional small synthetic data (SmallSYN). The SmallSYN data has 500 users, 500 features and $7$ co-clusters, and is shown in Figure <ref>(a). The shuffled data fed into CoClusInfo and is shown in Figure <ref>(b).  Subsequently, the users and products are rearranged to show discovered co-clusters as shown in Figure <ref>(d), according to which detects the co-clusters precisely (100% accuracy).  CoClusInfo's results (achieving 91% accuracy) are shown in Figure <ref>(c). Note that the inferred co-cluster sequences in Figure <ref>(c-d) are not  the same as in Figure <ref>(a) – this is because the cluster assignment for the same user/product clusters may be arbitrarily permuted during inference.

Dataset    	Cluster     	Metric    NEO-CC   DeepCC   CoClusInfo   FNMTF   WC-NMTF   SNCC

3*Cora    3*User    NMI   $ 0.034 \pm 0.004$   $0.003 \pm 0.001$   $0.152 \pm 0.006$   $0.152 \pm 0.011$   $0.173 \pm 0.002$   $0.112 \pm 0.001$   0.181 $\pm$ 0.023

Accuracy     $0.206\pm 0.034$     $0.2943\pm 0.045$     $0.375\pm 0.032$     $0.373\pm 0.101$     $0.314\pm 0.133$    $0.213 \pm 0.022$    0.399 $\pm$ 0.034

Time (sec)     $1488.8\pm 126.6$     $1481.4\pm 256.3$    $6.6\pm 1.23$    $4.5\pm 1.2$     $121.3\pm 12.5$    4.1 $\pm$ 0.34    $7.9 \pm 0.7$

3*WebKB4    3*User    NMI   $0.136 \pm 0.001$   $ 0.379 \pm 0.084$   $0.383 \pm 0.095$   $0.255 \pm 0.054$   $0.241 \pm 0.058$   $0.149 \pm 0.003$   0.489 $\pm$ 0.043

Accuracy     $0.374\pm 0.075$     $0.568\pm 0.026$     $0.653\pm 0.082$     $0.549\pm 0.102$     $0.503\pm 0.121$    $0.461 \pm 0.101$    0.752 $\pm$ 0.029

Time (sec)     $4089.7\pm 118.8$     $1719.9\pm 493.5$     $4.7 \pm 1.6$     $ 3.7 \pm 0.63$     $286.3\pm 2.2$    3.2 $\pm$ 0.21    $3.8\pm 0.6$

3*MovieLens    3*Product    NMI   $0.356 \pm 0.020$   $ 0.636 \pm 0.102$   $0.742 \pm 0.112$   $0.394\pm 0.021$   $0.689 \pm 0.039$   $0.472 \pm 0.102$   0.783 $\pm$ 0.124

Accuracy     $0.413 \pm 0.020$      $0.550\pm 0.192$     $0.649\pm 0.102$     $0.382\pm 0.048$     $0.592\pm 0.124$    $0.564 \pm 0.116$    0.683 $\pm$ 0.017

Time (sec)     $8402.3 \pm 363.5$       $2067.78\pm 34.5$     $103.98 \pm 3.6$    $34.78 \pm 6.9$    $673.67\pm 5.7$    29.43 $\pm$ 6.4    $ 41.55\pm 3.6$

3*News20    3*User    NMI   3*-   $ 0.448 \pm 0.018$   $0.499 \pm 0.091$   $0.149 \pm 0.031$   $0.392 \pm 0.019$   $0.334 \pm 0.016$   0.541 $\pm$ 0.019

Accuracy           $0.452\pm 0.075$     $0.433\pm 0.012$     $0.108\pm 0.093$     $0.329\pm 0.121$    $0.316 \pm 0.021$    0.477 $\pm$ 0.029

Time (sec)          $4028.3\pm 42.6$     $114.4\pm 4.1$      $97.4 \pm 12.3$     $1143.1\pm 41.6$     92.6 $\pm$ 10.4   $111.1\pm 16.3$

3*Caltech    3*User    NMI   $0.324 \pm 0.086$   $ 0.636 \pm 0.087$   $0.753 \pm 0.059$   $0.211 \pm 0.101$   $0.593 \pm 0.109$   $0.462 \pm 0.011$   0.756 $\pm$ 0.019

Accuracy     $0.543\pm 0.001$     $0.778\pm 0.022$     $0.899\pm 0.012$     $0.531\pm 0.111$     $0.835\pm 0.291$    $0.791 \pm 0.091$    0.911 $\pm$ 0.102

Time (sec)     $2649.6\pm 132.4$     $281.4\pm 13.3$     $2.4 \pm 0.72$     $3.5\pm 0.41$     $103.5\pm 17.4$    2.9 $\pm$ 0.21   $3.9 \pm 0.24$

Experimental results for NMI, Accuracy and CPU Time in seconds  for real data. The boldface means the best results.

Real Data: Table <ref> presents ’s performance compared to state-of-the-art techniques on real datasets, using the NMI, accuracy and CPU Time (in seconds). For each dataset, we list the scores that correspond to a specific cluster type. For all real datasets except Movielens data, only user labels are available. For Movielens data, product labels are inferred from movie genre types. is the only method that shows consistently high performance on all different kinds of datasets, indicating its flexibility. Due to space limitations, here, we provide analysis of Movielens data below and analysis of other real datasets are provided in Supplementary Material.

Movielens Data: We observe that is able to detect $15$ clusters for movies. This is reasonable outcome for this data, as there are overlapping movie categories; for example, “Toy Story (1995)” can be categorized in Animation as well as Comedy. FNMTF detected only $5-6$ clusters as shown in Figure <ref>.  CoClusInfo, WC-NMTF, SNCC and DeepCC achieved better performance, but still lower than .

Visualized co-clustering result of on the Movielens data. From left to right: (a) Original data, (b) co-cluster detected by , (c) FNMTF, (d) WC-NMTF, and (e) DeepCC. Each method is run with $M = 50, N = 20$ (specifying a maximum of 50 user clusters and 20 product clusters). evidently produces the most coherent co-clustering structure.

We also visualize the co-clustering results on the Movielens dataset, to show a visual representation of improved co-clustering performance. We rearrange the original data matrix (Figure <ref>(a)) according to the user and product cluster assignments to show the co-clustering result.  We observe that the co-clusters for (Figure <ref>(b)) are more salient compared to those for baselines FNMTF, WC-NMTF and DeepCC (Figures <ref>(c-e)). Snapchat dataset is very interesting and detailed analysis is provided in Section <ref>.

§.§.§ Explainability Evaluation

To answer the second experimental question, we first analyze the explanations derived from , and compare the results with explanations from baseline methods. Table <ref> presents 's performance over LightGBM explanations on 2 synthetic and 2 real datasets, using the stability score, average precision, CPU time in sec. and our proposed compression score.

Setup: First, we obtain co-clusters using and other very closely related two matrix tri-factorization based methods, i.e FNMTF and WC-NMTF. Note that baselines  CoclusInfo, DeepCC and NEO-CC provide only user and product cluster factor matrices but do not provide any summary matrix $\mathbf{S}$, which is required for our computations (See Equ. <ref> and <ref>).  Next, co-clustering results are fed into explainability components for both LightGBM and .
For , first, we create auxiliary matrices for both users and products (see Section <ref>), and then compute PMI using co-clustering outcomes (see Section <ref>).  We then select the top-5 user and product features per co-cluster to explain it.

Results: We observe that drastically outperforms the baselines for all the datasets. For synthetic data SYN-I, each user cluster is created with combination of age and gender. Similarly, product clusters are focused on gender dominated viewership. For example, in our simulated data,  women associate highly with “makeup & cosmetic, weight loss, and pop music” and men associate highly “card games, driving and racing games, and body building.”  Our results suggest that LightGBM is not able to explain co-clusters well with the correct attributes. This is likely, partially because it trains a classifier using only co-cluster outcomes and does not leverage implicit similarities between users and products. Conversely, successfully captures these similarities in auxiliary matrices and is thus able to explain these relationships for each co-cluster. Similarly, is able to provide more complex explanations for SYN-II in which product clusters do not have any specific gender based dominance.

For Movielens, discovered an interesting co-cluster in which salesmen and programmers strongly associated with adventure movies, and in another co-cluster lawyers strongly associated with drama and fantasy movies.  Interestingly, all co-clustering baselines with combination of LightGBM for explanation underperformed in our experiments, compared to the solution adopted by .  consistently compressed the discovered co-clusters quite well, compared to LightGBM, suggesting that we can describe the co-clusters more concisely with better-quality attribute/feature explanations. Overall, outperforms state-of-the-art approaches by at least $\approx15\% $ stability and  $\approx20\%$ average precision improvement. Also, achieves at least $5\%$ compression bits and $20\%$ runtime reduction (see Table <ref>), In the table <ref>) boldface means the best results. The 'a' represents Stability, 'b' represents Avg. Precision, 'c' represents CPU Time (sec) and 'd' represents  Compression (Kb). Here 'U' = User, 'P' = Product and 'B' = Both..

Experimental results for explainability evaluation.

2*Data      2*Score     3cFNMTF     3cWC-NMTF     3c

LightGBM   BMGUFS      LightGBM   BMGUFS      LightGBM   BMGUFS

6*S-I    U-a    $0.306 \pm 0.04$   $0.294 \pm 0.01$   $0.445 \pm 0.06$   $0.381 \pm 0.06$   $0.264 \pm 0.02$   $0.447 \pm 0.07$   $0.501 \pm 0.06$   $0.416 \pm 0.08$   0.596 $\pm$ 0.06

U-b    $0.119 \pm 0.08$   $0.115 \pm 0.06$   $0.134 \pm 0.01$   $0.112 \pm 0.03$   $0.201 \pm 0.09$   $0.289 \pm 0.02$   $0.148 \pm 0.01$   $0.296 \pm 0.16$    0.496$\pm$0.07

P-a    $0.288 \pm 0.03$   $0.276 \pm 0.04$   $0.301 \pm 0.02$   $0.129 \pm 0.04$   $0.131 \pm 0.02$   $0.132 \pm 0.01$   $0.116 \pm 0.09$   $0.264 \pm 0.12$    0.306 $\pm$ 0.04

P-b    $0.218 \pm 0.01$   $0.201 \pm 0.06$   $0.311 \pm 0.02$   $0.109 \pm 0.06$   $0.146 \pm 0.14$1   $0.231 \pm 0.03$   $0.274\pm 0.01$   $0.336 \pm 0.17$    0.412 $\pm$ 0.02

B-c    $43.34 \pm 5.2$   $264.43 \pm 21.64$   $8.45 \pm 4.2$   $58.01 \pm 4.3$   $124.6 \pm 6.4$   $10.23 \pm 3.1$   $28.26 \pm 2.6$   $119.6 \pm 10.4$   5.86 $\pm$ 1.9

B-d    $518.55 \pm 0.55$   $598.65 \pm 0.84$   $499.69 \pm 0.4$   $659.9 \pm 0.16$   $652.6 \pm 0.26$   $636.4 \pm 0.07$   $340.26 \pm 0.06$   $387.4 \pm 0.04$   307.5$\pm$ 0.05

6*S-II    U-a    $0.241 \pm 0.05$   $0.224 \pm 0.07$   $0.356 \pm 0.1$   $0.202\pm 0.02$   $0.200 \pm 0.06$   $0.331 \pm 0.02$   $0.485 \pm 0.07$   $0.483 \pm 0.02$   0.8481$\pm$ 0.02

U-b    $0.105 \pm 0.01$   $0.106 \pm 0.03$   $0.321 \pm 0.09$   $0.09 \pm 0.01$   $0.121 \pm 0.04$   $0.298 \pm 0.04$   $0.185 \pm 0.08$   $0.321 \pm 0.02$   0.403 $\pm$ 0.02

P-a    $0.637 \pm 0.11$   $0.592 \pm 0.09$   $0.812 \pm 0.13$   $0.071 \pm 0.01$   $0.062 \pm 0.01$   $0.035 \pm 0.01$   $0.585 \pm 0.15$   $0.576 \pm 0.14$   0.836 $\pm$ 0.06

P-b    $0.121 \pm 0.06$   $0.122 \pm 0.08$   $0.224 \pm 0.01$   $0.101 \pm 0.03$   $0.102 \pm 0.04$   $0.101 \pm 0.02$   $0.100 \pm 0.04$   $0.101 \pm 0.03$   0.390$\pm$ 0.05
B-c    $1123.69 \pm 12.7$   $1202.1.7 \pm 0.04$   $102.454 \pm 11.9$   $1229.82 \pm 32.8$   $1364.7 \pm 10.4$   $116.42 \pm 18.4$   $1178.24 \pm 8.9$   $962.4 \pm 10.4$   64.1 $\pm$ 4.8

B-d    $6294.58 \pm 0.08$     $6746.24 \pm 0.02$   $5928.72 \pm 0.01$   $9761.52 \pm 0.1$   $9842 \pm 0.04$   $9285.6 \pm 0.02$   $4590 \pm 0.09$   $4656.8 \pm 0.02$   3920.18 $\pm$ 0.07

6*ML     U-a    $0.197 \pm 0.01$    $0.210 \pm 0.01$   $0.126 \pm 0.06$   $0.136 \pm 0.04$   $0.141 \pm 0.02$   $0.1679 \pm 0.05$   $0.129\pm 0.08$   $0.196 \pm 0.02$   0.215 $\pm$ 0.04

U-b    $0.124 \pm 0.02$   $0.132 \pm 0.11$   $0.246 \pm 0.05$   $0.136 \pm 0.01$   $0.142 \pm 0.06$   $0.389 \pm 0.08$   $0.197 \pm 0.01$   $0.271 \pm 0.06$   0.466 $\pm$ 0.11

P-a    $0.009 \pm 0.07$   $0.101 \pm 0.01$   $0.101 \pm 0.02$   $0.129 \pm 0.09$   $0.162 \pm 0.03$   $0.196 \pm 0.06$   $0.212 \pm 0.02$   $0.108 \pm 0.01$   0.228 v 0.06

P-b    $0.246 \pm 0.12$   $0.294 \pm 0.02$   $0.301 \pm 0.02$   $0.312 \pm 0.01$   $0.264 \pm 0.01$   $0.358 \pm 0.03$   $0.326 \pm 0.09$   $0.113 \pm 0.05$   0.424 $\pm$ 0.03

B-c    $38.6 \pm 7.4$   $124.6 \pm 2.5$   $30.5 \pm 1.4$   $35.1 \pm 2.9$   $112.2 \pm 3.6$   $28.4 \pm 6.4$   $21.24 \pm 1.4$   $98.7 \pm 3.3$   14.30 $\pm$ 2.9

B-d   $749.29 \pm 0.42$   $760.4 \pm 0.01$   $730.95 \pm 0.19$   $788.08 \pm 0.12$   $796.6 \pm 0.04$   $760.02 \pm 0.18$   $696.17 \pm 0.49$   $683.6 \pm 0.01$   671.7 $\pm$ 0.32

6*SC    U-a    $0.289 \pm 0.04$   $0.122 \pm 0.04$   $0.305 \pm 0.02$   $0.292 \pm 0.03$   $0.161 \pm 0.01$   $0.321 \pm 0.04$   $0.329 \pm 0.04$   $0.231 \pm 0.02$   0.351$\pm$ 0.08

U-b    $0.374 \pm 0.07$   $0.101 \pm 0.04$   $0.518 \pm 0.11$   $0.281 \pm 0.03$   $0.009 \pm 0.01$   $0.521 \pm 0.04$   $0.277 \pm 0.04$   $0.201 \pm 0.01$   0.527 $\pm$ 0.05

P-a    $0.298 \pm 0.01$   $0.106 \pm 0.07$   $0.344 \pm 0.09$   $0.312 \pm 0.01$   $0.112 \pm 0.01$   $0.241 \pm 0.09$   $0.311 \pm 0.09$   $0.261 \pm 0.09$   0.351 $\pm$ 0.01

P-b    $0.178 \pm 0.03$   $0.113 \pm 0.02$   $0.236 \pm 0.01$   $0.201 \pm 0.02$   $0.121 \pm 0.04$   $0.241 \pm 0.10$   $0.200 \pm 0.12$   $0.102 \pm 0.01$   0.257 $\pm$ 0.03

B-c    $1498.4 \pm 12.5$   $2642.6 \pm 23.7$   $149.5 \pm 38.8$   $1298.5 \pm 34.9$   $4364.36 \pm 12.4$   $246 \pm 13.9$   $1256.37 \pm 56.2$   $3641.78 \pm 34.7$   107.15 $\pm$11.3

B-d    $60457.47 \pm1.21$    $61362.2 \pm 0.04$   $54538.27 \pm1.14$   $58694.03 \pm3.21$   $58961.7 \pm 0.01$   $56022.87 \pm3.42$   $52889 \pm0.45$   $53316.8 \pm 0.14$   47765.7 $\pm$ 0.34

(a) Total running time  (averaged over 10 runs) of versus the total number of non-zeros for Snapchat data (b) Co-clustering running time for synthetic data.

§.§ Scalability

Finally, to answer the third question, we experimentally study the runtime of with respect to the input size on real graph. For runtime measurements, we use a large private viewer-publisher Snapchat data with $5$ million viewer and $7500$ publisher. To generate real graphs of growing size, we increasingly sample the Snapchat user-product user engagement matrix rows and report runtime averaged over $10$ runs in Figure <ref>.  We can observe that the run-time of scales approximately linearly; notably, can handle interaction matrices with many millions of interactions in mere minutes. We show only linear scaling for Snapchat dataset because (a) it is the largest dataset, and (b) the scaling trends are consistent with other datasets.

§.§ Effectiveness of Auxiliary Feature Matrix

To show the effectiveness of using meta path based auxiliary feature matrices (in short AFM), we computed average precision scores with and without them in Equ. <ref> and <ref> for Movielens data. We fed with varied number of user cluster $M$ (as this is unknown) but we fixed the number of movie clusters $N$ to 20. Table <ref> shows that performance improved by $\approx 7.5\%$ when using auxiliary feature matrices.

Given $M$      Predicted $M$    Without AFM       With AFM

$15$    $10$    $0.105\pm 0.11$   $0.113\pm 0.14$

$30$    $24$   $0.198\pm 0.02$     $0.215\pm 0.04$

$50$    $27$   $0.201\pm 0.01$    $0.214\pm 0.06$

Experiment evaluation of Auxiliary Feature Matrix on Movielens dataset.

§.§ Parameter Sensitivity Analysis

We evaluate the sensitivity of the interpolation parameters $\alpha,\beta$ and $\gamma$  in Equ. <ref>-<ref> which describe involvement of each meta-path matrix. We learn stability and average precision score for different combinations of $\{\alpha, \beta,\gamma\}$ on synthetic data SYN-I. Here, we kept $\beta = \gamma$ to give equal importance to 4-step meta path. Figure <ref>(left) shows that performs better when higher importance is given to direct paths i.e $\alpha \geq 0.5$. Beyond $\alpha = 0.5$ , there is no significant change in performance.  Moreover, Figure <ref> (right) shows that average precision achieves optimal values around $\alpha=0.5$, suggesting that incorporating indirect interactions via meta-paths beyond just direct paths does contribute to improved performance.  Hence, we choose $\alpha = 0.5$,$\beta = \gamma = 0.25$ in our experiments.

Sensitivity of interpolation parameters for co-cluster explanation; values of $\alpha = 0.5$ (interpolating direct path with indirect metapath matrices) produce best explanation results.  maybe 1 sentence take-home msg here?

§.§ at Work

We use on a private viewer-publisher interaction dataset from Snapchat, consisting of $500K$ viewers, $2500$ publishers, and $5$ million viewer-publisher interactions (here, viewers correspond to users, and publishers to products, for consistency with our prior discussion).  The dataset is naturally unlabeled, making NMI and accuracy analysis infeasible; therefore, we resort to qualitative discussion.  In this dataset, each viewer is described by $22$ associated features (age, gender, country, etc.) and each publisher has $238$ features (e.g. publisher demographics, publishing category etc.). We ran with the maximum number of user clusters $M = 100$ and product clusters $N = 25$. For explainability evaluation, we use top $N=5$ highest feature values computed with Equ. <ref> and <ref>. finds viewer-publisher co-clusters of various sizes in Snapchat data, as shown in Figure <ref>. Table <ref> provide three major niches based on high PMI values from the summary matrix $\mathbf{S}$.

Beauty & Lifestyle     Family & Hobbies       Music

Cosmetics   Parenting & child care   Urban hip-hop

Women Fashion   Home Improvement   Electronic Dance Music

Women Lifestyle   Shopping   Pop Music

Spa beauty     Gardening   Concerts

Nail Art    Decor Design   Classical Opera

Illustration of niches discovered by .

on Snapchat data finds clusters of viewers/publishers with similar attributes coherence. The interaction matrix is carefully arranged by ,  revealing patterns: e.g., the young women heavily view content regarding `Beauty' category published by women creators, and they focus on the same group of content.

The viewers clusters in the green square labeled as `1' mostly belong to groups of young women (age between $13-20$) associated with publisher who publish content regarding ‘Beauty & Lifestyles’ category. We observe that viewers in this dense group have similar content consumption. The viewers cluster `2' contains viewers with connections across many publishers clusters like 'Home & Family' and 'Hobbies & Interests'. All of the viewers are between two age groups $25-34$ and $>35$. This niche does not have any gender dominance in viewer cluster but we observe that most of the content creators or publishers are men. Niche represented by `3' is the most prominent: this is a dense group of male viewers, $>70\%$ of them are from North America and all of them are between age group $21-24$, highly associated with content related to 'World Music', published by group of male publishers between age group $25-34$. Most of the publishers are from Europe. Overall as shown in Table <ref>, outperforms state-of-the-art approaches by $\approx 6\% $ stability and $\approx 4\%$ average precision improvement. Also, achieves at least $10\%$ compression bits and $25\%$ run time reduction (See Table <ref>) for Snapchat dataset.  All in all, by using , we are able to understand and explain the user content consumption graph in a completely unsupervised fashion.

§ CONCLUSIONS

In this work, we tackle the problem of discovering market niches for strategic content creation to satisfy diverse audience groups.  We pose the niche detection problem as one which involves discovering coherent co-clusters in user-product (content) interaction graph data, as well as explaining the co-clusters using nodal attributes on user and product nodes.  To our knowledge, ours is the first work which tackles such an explainable co-clustering problem.  We proposed , the first niche detection framework for finding and explaining co-clusters in attributed interaction graphs.  utilizes principles from mutual information to (a) propose and solve a non-negative matrix tri-factorization oriented objective to discover cohesive co-clusters, and (b) select important user and product features associated with these co-clusters using meta-path driven feature selection. In doing so, we find that successfully discovers niches in user-content consumption data, and demonstrates improvements over both state-of-the-art co-clustering approaches ($\approx14\%$ accuracy) as well as candidate explanation approaches ($\approx20\%$ average precision) on multiple simulated and real-world datasets.  Finally, we show that provides an advantageous speed-quality tradeoff over alternative approaches, and scales effortlessly and discovers interesting insights on large-scale interaction data with over 60M interactions, via a private dataset from Snapchat.

The content of this chapter was under blind peer review at the time of thesis submission.

Mining Streaming Tensors

CHAPTER: SAMPLING-BASED BATCH INCREMENTAL TENSOR DECOMPOSITION

[backgroundcolor=Orange!20,linewidth=1pt,  topline=true,  rightline=true, leftline=true]
"How to summarize high-order data tensor? How to incrementally update those patterns over time?”

Tensor decompositions are invaluable tools in analyzing multimodal datasets. In many real-world scenarios, such datasets are far from being static, to the contrary they tend to grow over time. For instance, in an online social network setting, as we observe new interactions over time, our dataset gets updated in its "time" mode. How can we maintain a valid and accurate tensor decomposition of such a dynamically evolving multimodal dataset, without having to re-compute the entire decomposition after every single update? In this chapter we introduce , a Sampling-based Batch Incremental Tensor Decomposition algorithm, which incrementally maintains the decomposition given new updates to the tensor dataset. is able to scale to datasets that the state-of-the-art in incremental tensor decomposition is unable to operate on, due to its ability to effectively summarize the existing tensor and the incoming updates, and perform all computations in the reduced summary space. We extensively evaluate using synthetic and real datasets. Indicatively, achieves comparable accuracy to state-of-the-art incremental and non-incremental techniques, while being up to 25-30 times faster. Furthermore, scales to very large sparse and dense dynamically evolving tensors of dimensions up to $100K \times 100K \times 100K$ where state-of-the-art incremental approaches were not able to operate. The content of this chapter is adapted from the following published paper:

Gujral, Ekta, Ravdeep Pasricha, and Evangelos E. Papalexakis. "Sambaten: Sampling-based batch incremental tensor decomposition." In Proceedings of the 2018 SIAM International Conference on Data Mining, pp. 387-395. Society for Industrial and Applied Mathematics, 2018.

§ INTRODUCTION

Tensor decomposition is a very powerful tool for many problems in data mining <cit.>. , machine learning <cit.>, chemometrics <cit.>, signal processing <cit.> to name a few areas. The success of tensor decomposition lies in its  capability of finding complex patterns in multi-way settings, by leveraging higher-order structure and correlations within the data.  The dominant tensor decompositions are CP/PARAFAC (henceforth referred to as CP), which extracts interpretable latent factors from the data, and Tucker, which estimates the joint subspaces of the tensor. In this work we focus on the CP decomposition, which has been shown to be extremely effective in exploratory data mining time and time again <cit.>.

outperforms state-of-the-art baselines while maintaining competitive accuracy.

In a wide array of modern real-world applications, especially in the age of Big Data, data are far from being static. To the contrary, data get updated dynamically. For instance, in an online social network, new interactions occur every second and new friendships are formed at a similar pace. In the tensor realm, we may view a large proportion of these dynamic updates as an introduction of new “slices” in the tensor: in the social network example, new interactions that happen as time evolves imply the introduction of new snapshots of the network, which grow the tensor in the “time” mode. A tensor decomposition in that tensor can discover communities and their evolution over time. How can we handle such updates in the data without having to re-compute the decomposition whenever an update arrives, but incrementally update the existing results given the new data? In the community detection example, how can we track the evolution of the existing communities, and discover new ones, for the new updates that continuously arrive?

Computing the decomposition for a dynamically updated tensor is challenging, with the challenges lying, primarily, on two of the three V's in the traditional definition of Big Data: Volume and Velocity. As a tensor dataset is updated dynamically, its volume increases to the point that techniques which are not equipped to handle those updates incrementally, inevitably fail to execute due to the sheer size of the data. Furthermore, even though the applications that tensors have been successful so far  do not require real-time execution per se, the decomposition algorithm must, nevertheless, be able to ingest the updates to the data at a rate that will not result in the computation being “drowned” by the incoming updates.

The majority of prior work has focused on the Tucker Decomposition of incrementing tensors <cit.>, however very limited amount of work has been done on the CP. Nion and Sidiropoulos <cit.> proposed  two methods namely Simultaneous Diagonalization Tracking (SDT) and Recursive Least Squares Tracking (RLST) and most recently, <cit.> introduced the OnlineCP decomposition for higher order online tensors. Even though prior work in incremental CP decomposition, by virtue of allowing for incremental updates to the already computed model, is able to deal with Velocity, when compared to the naive approach of re-executing the entire decomposition on the updated data, every time a new update arrives, it falls short when the Volume of the data grows.

We show a snapshot of our results in Figure <ref>: is faster than all state-of-the-art methods on data that the baselines were able to operate on. Furthermore, was able to scale to, both dense and sparse, dynamically updated tensors, where none of the baselines was able to run. Finally, achieves comparable accuracy to  existing incremental and non-incremental methods.
Our contributions are summarized as follows:

* Novel scalable algorithm: We introduce  , a novel scalable  and  parallel incremental tensor decomposition algorithm for efficiently computing the CP decomposition of incremental tensors. The advantage of stems from the fact that it only operates on small summaries of the data at all times, thereby being able to maintain its efficiency regardless of the size of the full data.

To the best of our knowledge, this is the first incremental tensor decomposition which effectively leverages sparsity in the data.

* Quality control: As a tensor is dynamically updated, some of the incoming updates may contain rank-deficient structure, which, if not handled appropriately, can pollute the results. We equip with a quality control option, which effectively determines whether an update is rank-deficient, and subsequently handles the update in a way that it does not affect latent factors that are not present in that update.

* Extensive experimental evaluation: Through experimental evaluation on six real-world datasets with sizes that range up to $70$GB, and synthetic tensors that range up to $100K\times 100K \times 100K$, we show thatcan incrementally maintain very accurate decompositions, faster and in a more efficient and scalable manner than state-of-the-art methods.

* Stability and Scalability: Based on extensive experimental analysis on real and synthetic datasets, our algorithm produces more stable decompositions than existing online approaches, as well as better scalability.

Reproducibility: We make our Matlab implementation publicly available at link []. Furthermore, all the datasets we use for evaluation are publicly available.

§ RELATED WORK

Incremental tensor methods in the literature can be categorized into three main categories:  1) Tucker
decomposition, 2) CP decomposition, 3) Tensor completion

Tucker Decomposition:
Online tensor decomposition was first proposed by Sun el at.<cit.> as ITA (Incremental Tensor Analysis), describing the three variants of Incremental Tensor Analysis. First,  DTA i.e. Dynamic tensor analysis which is based on calculation of co-variance of matrices in traditional higher-order singular value decomposition in an incremental fashion. Second, with help of SPIRIT algorithm, they found approximation of DTA named as Stream Tensor Analysis (STA). Third, they proposed window-based tensor analysis (WTA). To improve the efficiency of DTA, it uses a sliding window strategy. Liu el at.<cit.> proposed an efficient method to diagonalize the core tensor to overcome this problem. Other approaches replace SVD with incremental SVD to improve the efficiency. Hadi el at. <cit.> proposed the multi-aspect-streaming tensor analysis (MASTA) method that relaxes  constraint and allows the tensor to simultaneously grow in all modes.

CP Decomposition:
There is very limited study on online CP decomposition methods. Phan el at. <cit.> had developed a theoretic approach GridTF to large-scale tensors processing based on an extension to CP's fundamental mathematics theory. They used divide and concur technique to get sub-tensors and fuses the output of all factorization to achieve final factor matrices which is proved to be same as decomposing the whole tensor using CP decomposition. Its potential of concurrent computing methods to adapt the engineering applications remains unclear. Nion el at. <cit.>, proposed two algorithms that focus on CP decomposition namely SDT (Simultaneous Diagonalization Tracking)  that incrementally perform the SVD of the unfolded tensor; and RLST (Recursive Least Squares Tracking) , which recursively updates the decomposition factors by minimizing the mean squared error. The latest related work is OnlineCP, proposed by Zhou, el at. <cit.>, is an online CP decomposition method, where the the latent factors are updated when there are new data.

Tensor Completion:
The main difference between completion and decomposition techniques is that in completion “zero” values are considered “missing” and are not part of the model, and the goal is to impute those missing values accurately, rather than extracting latent factors from the observed data.

The earliest work on incremental tensor completion traces back to  <cit.>, and recently,  Qingquan el at.<cit.>,  proposed streaming tensor completion based on block partitioning.

§ PROBLEM FORMULATION

In many real-world applications, data grow dynamically. In a time-evolving social network, we observe user interactions every few seconds, which can be translated to new tensor slices, after fixing the temporal granularity (a problem which, on its own merit, is very hard to solve optimally, and we do not address in this chapter). This incremental property of data  gives rise to the need for an on-the-fly update of the existing decomposition, which we name incremental tensor decomposition. Notice that the literature (and thereby this chapter) uses the terms “incremental”, “dynamic”, and “online” interchangeably.  In such scenarios, data updates happen very fast which make traditional (non-incremental) methods to collapse because they need to recompute the decomposition for the entire dataset.

Real life dynamic data examples (a) collecting new location update from satellite to GPS recommendations systems (b) growing Facebook interaction between people over the time.

We focus on a 3-mode tensor one of whose dimensions are growing with time. However, the problem definition (and our proposed method) extends to any number of modes. Let us consider $\tensor{X}(t)$=$\mathbb{R}^{I\times J \times K_1(t)}$ at time $t$. The CP decomposition of $\tensor{X}(t)$ is given as :

\[
\begin{array}{c}
\tensor{X}^{(1)}(t)\approx (\mathbf{A}(t) \odot \mathbf{B}(t))\mathbf{C}^{T}(t) \approx \mathbf{L}(t)\mathbf{C}^{T}(t)
\end{array}
\]

where $\mathbf{L}(t)=(\mathbf{A}(t) \odot \mathbf{B}(t))$ of dimension $IJ\times R$ and $\mathbf{C}^{T}(t)$ is of dimension $K_{1}\times R$. When new incoming slice $\tensor{X}(t^{'})$=$\mathbb{R}^{I\times J \times K_2(t^{'})}$ is added in mode 3, required decomposition at time $t^{'}$ is  :

\[
\begin{array}{c}
\tensor{X}^{(1)}(t+t^{'}) \approx \mathbf{L}(t+t^{'})\mathbf{C}^{T}(t+t^{'})
\end{array}
\]

where $\mathbf{L}(t+t^{'})=(\mathbf{A}(t+t^{'}) \odot \mathbf{B}(t+t^{'}))$ of dimension $IJ\times R$ and $\mathbf{C}^{T}(t+t^{'})$ is of dimension $(K_{1}+K_{2})\times R$.

The problem that we solve is the following:

[linecolor=red!60!black,backgroundcolor=gray!20,linewidth=1pt,    topline=true,rightline=true, leftline=true]
Given (a) an existing set of decomposition results $\mathbf{A}(t),\mathbf{B}(t)$ and $\mathbf{C}(t)$ of $R$ components, which approximate tensor $\tensor{X}_{old}$ of size $ I \times J \times K_1$ at time t , (b) new incoming batch of slices in form of tensor $\tensor{X}_{new}$ of size $ I \times J \times K_2$ at any time $t^{'}$, find updates of $\mathbf{A}(t^{'}),\mathbf{B}(t^{'})$ and $\mathbf{C}(t^{'})$ incrementally to approximate tensor $\tensor{X}$ of dimension $ I \times J \times K$, where K =$ K_1 + K_2$ after appending new slice or tensor to $3^{rd}$ mode while maintaining a comparable accuracy with running the full CP decompositon on the entire updated tensor $\tensor{X}$.

To simplify notation, we will interchangeably refer to $\mathbf{A}(t)$ as $\mathbf{A}_{old}$ (when we need to refer to specific indices of that matrix), and similarly for $\mathbf{A}(t')$ we shall refer to it as $\mathbf{A}^{'}$.

§ PROPOSED METHOD: SAMBATEN

As we mention in the introduction, there exists a body of work in the literature that is able to efficiently and incrementally update the CP decomposition in the presence of incoming tensor slices <cit.>. However, those methods fall short when the size of of the dynamically growing tensor increases, and eventually are not able to scale to very large dynamic tensors. The reason why this happens is because these methods operate on the full data, and thus, even though they incrementally update the decomposition (avoiding to re-compute it from scratch), inevitably, as the size of the full data grows, it takes a toll on the run-time and scalability.

In this chapter we propose , which takes a different view of the solution, where instead of operating on the full data, it operates on a summary of the data. Suppose that the “complete” tensor (i.e., the one that we will eventually get when we finish receiving updates) is denoted by $\tensor{X}$. Any given incoming slice (or even a batch of slice updates) can be, thus, seen as a sample of that tensor, $\tensor{X}$ where the sampled indices in the third mode (which we assume is the one receiving the updates) are the indices of the incoming slice(s).  Suppose, further, that given a set of sample tensors (which are drawn by randomly selecting indices from all the modes of the tensor) we can approximate the original tensor with high-accuracy (which, in fact, the literature has shown that it is possible <cit.>). Therefore, when we receive a new batch of slices as an update, if we update those samples with the new indices, then we should be able to compute a decomposition very efficiently which incorporates the slice updates, and approximates the updated tensor well. A visual summary of is shown in Figure <ref>.

: Sampling-based Batch Incremental Tensor Decomposition: 1) Sample incoming tensor into sub-tensors,  2) run parallel decompositions on the samples, 3) project back the results into the original space, and, finally, 4) update the incrementally growing factor matrix $\mathbf{C}$.

§.§ The heart of

The algorithmic framework we propose is shown in Figure <ref> and is described below:
We assume that we have an existing set of decomposition results, as well as a set of summaries of the tensor, before the update. Summaries are in the form of sampled sub-tensors, as described in the text below. For simplicity of description, we assume that we are receiving updated slices on the third mode, which in turn have to add new rows to the $\mathbf{C}$ matrix (that corresponds to the third mode). We, further, assume that the updates come in batches of new slices, which, in turn, ensures that we see a mature-enough update to the tensor, which contains useful structure. Trivially, however, can operate on singleton batches.
In the following lines, $\tensor{X}$ is the tensor prior to the update and $\tensor{X}_{new}$ is the batch of incoming slices. Given an update, performs the following steps:

Sample: The rationale behind is that each batch $\tensor{X}_{new}$ can be seen as a sample of third-mode indices of (what is going to be) the full tensor. In this step, we are going to merge these incoming indices with an already existing set of sampled tensors. In order to obtain those pre-existing samples, we follow a similar approach to <cit.>. Namely, we sample indices from the tensor $\tensor{X}$ based on a measure of importance.
To determine the importance for each mode $m$ and then sample the indices using this measure as a sampling weight divided by its probability. An appropriate measure of importance (MoI) is the sum-of-squares of the tensor for each mode.  For the first mode , MoI is defined as:
$x_a(i)=\sum_{j=1}^J \sum_{k=1}^K \tensor{X}(i,j,k)^2$
for  $i \in$ (1,I). Similarly, we can define the MoI for modes 2 and 3.

\begin{equation}\label{eq:2}
x_b(j)=\sum_{i=1}^I \sum_{k=1}^K \tensor{X}(i,j,k)^2, x_c(k)=\sum_{i=1}^I \sum_{j=1}^J \tensor{X}(i,j,k)^2
\end{equation}

for  $j \in$ (1, J) and $k \in$  (1,K).

We sample each mode of $\tensor{X}$ without replacement, using the above MoI to bias the sampling probabilities. With $s$ as sampling factor, i.e. if $\tensor{X}$ has size $I \times J \times K$, then $\tensor{X}_s$ will be of size $\frac{I}{s},\frac{J}{s},\frac{K}{s}$. Sampling rate for each mode is independent from each other, and in fact, different rates can be used for imbalanced modes. In the case of sparse tensors, the sample will focus on the dense regions of the tensor which contains most of the useful structure. In the case of dense tensors, the sample will give priority to the regions of the tensor with the highest variation.

After forming the sample summary $\tensor{X}_s$ for $\tensor{X}$, we merge it with the samples obtained from the intersection of the third-mode indices of $\tensor{X}_{new}$ and the already sampled indices in the remaining modes, so that the final sample is equal to
$\tensor{X}_s$=$\tensor{X}(I_s,J_s,K_s \cup [K+1 \cdots K_{new}])$, where $K+1 \cdots K_{new}$ are the third-mode indices of $\tensor{X}_{new}$.

Due to the randomized nature of this summarization, we need to draw multiple samples, in order to obtain a reliable set of summaries. Each such independent sample is denoted as $\tensor{X}_s^{(r)}$. In the case of dense tensors, obtaining multiple, independent random samples helps summarize as much of the useful variation as possible. In fact, we will see in the experimental evaluation that increasing the number of samples, especially for dense tensors, improves accuracy.

In <cit.> the authors note that in order for their method to work, a set of anchor indices must be common between all samples, so that, later on, we can establish the correct correspondence of latent factors. However, in we do not have to actively fix a set of indices across sampling repetitions. When we sample $I_s,J_s,K_s$ each time, those indices correspond to a portion of the decomposition that is already computed. Therefore, the entire set of indices $I_s,J_s,K_s$ can serve as the set of anchors. This is a major advantage compared to <cit.>, since 1) does not need to commit to a set of fixed indices for all samples a-priori, which, due to randomization may happen to represent a badly structured portion of the tensor, 2) does not need to be restricted in a “small” set of fixed common indices (which is required in <cit.> in order to ensure that sufficiently enough new indices are sampled across repetitions), but to the contrary, is able to use a larger number of anchor indices to establish correspondence more reliably, and 3) does not require any synchronization between different sampling repetitions, which results in higher parallelism potential.

Decompose: Having obtained $\tensor{X}_s$, from the previous step, decomposes the summary using any state-of-the-art algorithm, obtaining factor matrices $[\mathbf{A}_i^{'}, \mathbf{B}_i^{'}, \mathbf{C}_i^{'}]$. For the purposes of this chapter, we use the Alternating Least Squares (ALS) algorithm for the CP decomposition, which is probably the most well studied algorithm for CP.

Tensor $\tensor{X}_{new}$ of size $I \times J \times K_{new}$, Factor matrices $\mathbf{A}_{old},\mathbf{B}_{old}, \mathbf{C}_{old}$, sampling factor $s$ and number of repetitions $r$.
Factor matrices $\mathbf{A}, \mathbf{B}, \mathbf{C}$, $\lambda$.
$i = 1$ to $r$
Compute $\mathbf{x}_a, \mathbf{x}_b$ and $\mathbf{x}_c$.

Sample a set of indices $I_s,J_s,K_s$ from $\tensor{X}$ without replacement using $\mathbf{x}_a(i)/\sum\limits_{i=1}^I x_a(i) $ as probability (accordingly for the rest).

$\tensor{X}_s$=$\tensor{X}(I_s,J_s,K_s \cup [K+1 \cdots K_{new}])$

$[\mathbf{A}_i^{'}, \mathbf{B}_i^{'}, \mathbf{C}_i^{'}] = $ CP$\left( \tensor{X}_s, R \right)$.

Normalize $\mathbf{A}_i^{'}, \mathbf{B}_i^{'}, \mathbf{C}_i^{'}$ (as in the text) and absorb scaling in $\boldsymbol{\lambda}$.

Compute optimal matching between the columns of $\mathbf{A}_{old}, \mathbf{B}_{old}, \mathbf{C}_{old}$ and $\mathbf{A}_i^{'}, \mathbf{B}_i^{'}, \mathbf{C}_i^{'}$ as in the text

Update only zero entries of $\mathbf{A}, \mathbf{B}, \mathbf{C}$ that correspond to sampled entries of $\mathbf{A}_i^{'}, \mathbf{B}_i^{'}, \mathbf{C}_i^{'}$

Obtain $\mathbf{C}_{new}$ of size $K_{new} \times R$ by taking last $K_{new}$ rows of $\mathbf{C}^{'}$.

Use a shared copy of $\mathbf{C}_{new}$ and average out its entries in a column-wise fashion across different sampling repetitions.

Update $\mathbf{C}$ of size $(K_{new}+K_{old}) \times R$ as :
$C =\begin{bmatrix}
\mathbf{C}_{old} \\
\mathbf{C}_{new}\\
\end{bmatrix}$

Update scaling $\boldsymbol{\lambda}$ as the average of the previous and new value.

Project back: The CP decomposition is unique (under mild conditions) up to permutation and scaling of the components <cit.>. This means that, even though the existing decomposition $[\mathbf{A}_{old},\mathbf{B}_{old}, \mathbf{C}_{old}]$ may have established an order of the components, the decomposition $[\mathbf{A}_i^{'}, \mathbf{B}_i^{'}, \mathbf{C}_i^{'}]$ we obtained in the previous step is very likely to introduce a different ordering and scaling, as a result of the aforementioned permutation and scaling ambiguity. Formally, the sampled portion of the existing factors and the currently computed factors are connected through the following relation Equ. <ref>:

\begin{equation}
\label{sambateneq:rel}
[\mathbf{A}_{old}(I_s,:),\mathbf{B}_{old}(J_s,:), \mathbf{C}_{old}(K_s,:)] =
[\mathbf{A}_i^{'}(I_s,:) \boldsymbol{\Lambda \Pi}, \mathbf{B}_i^{'}(J_s,:) \boldsymbol{\Pi}, \mathbf{C}(K_s,:)_i^{'} \boldsymbol{\Pi}]
\end{equation}

where $\boldsymbol{\Lambda}$ is a diagonal scaling matrix (which without loss of generality we absorb on the first factor), and $\boldsymbol{\Pi}$ is a permutation matrix that permutes the order of the components (columns of the factors).

In order to tackle the scaling ambiguity, we need to normalize the results in a consistent manner. In particular, we normalize such that each column of the newly computed factors which spans the indices that are shared with $[\mathbf{A}_{old},\mathbf{B}_{old}, \mathbf{C}_{old}]$ has unit norm:
and accordingly for the remaining factors. Note that for $\mathbf{A}_i^{'}$,  trivially holds that $\mathbf{A}_i^{'}(I_s,f) = \mathbf{A}_i^{'}(:,f)$ and similarly for $\mathbf{B}_i^{'}$. After normalizing, the relation between the existing factors and the currently computed is $\mathbf{A}_{old}(I_s,:) = \mathbf{A}_i^{'} \boldsymbol{\Pi}$ (and similarly for the rest). Each iteration retains a copy of $[\mathbf{A}_{old}(I_s,:),\mathbf{B}_{old}(J_s,:), \mathbf{C}_{old}(K_s,:)]$ which will serve as the anchor for disambiguating the permutation of components. We normalize $[\mathbf{A}_{old}(I_s,:),\mathbf{B}_{old}(J_s,:), \mathbf{C}_{old}(K_s,:)]$ to unit norm as well, and the reason behind that lies in Lemma <ref>:

Consider $\mathbf{a} = \mathbf{A}_i^{'}(:,f_1)$ and $\mathbf{b} = \mathbf{A}_{old}(:,f_2)$. If $f_1$ and $f_2$ correspond to the same latent CP factor, in the noiseless case, then $\mathbf{a}^T \mathbf{b} = 1$ otherwise  $\mathbf{a}^T \mathbf{b} < 1$.

From Cauchy-Schwartz inequality  $\mathbf{a}^T \mathbf{b} \leq \|\mathbf{a}\|_2 \|\mathbf{b}\|_2$. The above inequality is maximized when $\mathbf{a} = \mathbf{b}$ and for unit norm $\mathbf{a}, \mathbf{b}$, $\mathbf{a}^T \mathbf{b} \leq 1$. Therefore, if $\mathbf{a} = \mathbf{b}$, which happens when $f_1$ and $f_2$ correspond to the same latent CP factor, $\mathbf{a}^T \mathbf{b} =  1$.

Lemma <ref> is a guide for identifying the permutation matrix $\boldsymbol{\Pi}$: For every column of $\mathbf{A}_i^{'}$ we compute the inner product with every column of $\mathbf{A}_{old}(I_s,:)$ and compute a matching when the inner product is equal (or close) to 1. Given a large-enough number of rows for  $\mathbf{A}_i^{'}$ (which is usually the case, since we require a large-enough sample of the tensor in order to augment it with the update and compute the factor updates accurately), this matching can be computed  reliably even in noisy real-world data, as we show in the experimental evaluation.

Update results:
After appropriately permuting the columns of $\mathbf{A}_i^{'}, \mathbf{B}_i^{'}, \mathbf{C}_i^{'}$, we have all the information needed to update our model. Returning to the problem definition of Section <ref>, $\mathbf{A}_i^{'}$ contains the updates to the rows within $I_s$ for $\mathbf{A}(t)$ (and similarly for $\mathbf{B}$ and $\mathbf{C}$). Even though $\mathbf{A, B}$ do not increase their number of rows over time, the incoming slices may contribute valuable new estimates to the already estimated factors. Thus, for the already existing portions of $\mathbf{A, B, C}$ we only update the zero entries that fall within the range of $I_s, J_s$, and $K_s$ respectively. Finally, $ \mathbf{C}_i^{'}([K+1 \cdots K_{new}],:)$ contains the factors for the newly arrived slices, which need to be merged to the already existing columns of $\mathbf{C}$. After properly permuting the columns of $\mathbf{C}_i^{'}$, we accumulate the lower portion of the  $\mathbf{C}_i^{'}$ (corresponding to the new slices) into $\mathbf{C}_{new}$ and we take the column-wise average of the rows to-be-appended to $\mathbf{C}$, across repetitions. Finally, we update
\mathbf{C}(t') =\begin{bmatrix}
\mathbf{C}_{old} \\
\mathbf{C}_{new}\\
\end{bmatrix}.

Guarantees for Correctness  Lemma <ref> is essentially providing a guarantee for the correctness of . In particular, Lemma <ref> ensures that, under mild conditions that make the CP decomposition unique and identifiable, will discover the correct latent factors, disambiguate the permutation and scaling ambiguity, and update the correct columns of the factor matrix that is to be augmented. As we will also empirically demonstrate in the experimental evaluation, indeed, is able to produce correct results that are on par with state-of-the-art methods.

§.§ Dealing with rank deficient updates

So far, the above algorithm description is based on the assumption that each one of the sampled tensors $\tensor{X}_s$ we obtain are full-rank, and the CP decomposition of that tensor is identifiable (assumption that is also central to previous works <cit.>). However, this assumption glosses over the fact that in reality, updates to our tensor may be rank-deficient. In other words, even though $\mathbf{A}(t),\mathbf{B}(t),\mathbf{C}(t)$ have $R$ components, the update may contain $R_{new}$ components, where $R_{new} < R$. If that happens, then the matching as described above is going to fail, and this inevitably leads to very low quality results. Here, we tackle this issue by adding an extra layer include graphics of quality control when we compute the CP decomposition, right before line 5 of Algorithm <ref>: we estimate the number of components $R_{new}$ in $\tensor{X}_s$ and instead of the “universal” rank $R$, which may result in low quality factors, we use $R_{new}$ and we accordingly match only those $R_{new}$ to their most likely matches within the existing $R$ components.  Estimating the number of components is a very hard problem, however, there exist efficient heuristics in the literature, such as the Core Consistency Diagnostic () <cit.> which gives a quality rating for a computed CP decomposition. By successively trying different candidate ranks for $\tensor{X}_s$, we estimate its actual rank, as shown in Algorithm  <ref> (), and use that instead. For efficiency we use a recent implementation of that is especially tailored for exploiting sparsity <cit.>. In the experimental evaluation we demonstrate that using indeed results in higher-quality latent factors

Tensor $\tensor{X}$, maximum rank $R$ , maximum no. of iterations $it$.
1 to $R$
$j=1$ to $it$
Run CP Decomposition on $\tensor{X}$ with rank $i$ and obtain $f_i$

Run ($\tensor{X}_s$,$f_i$) and obtain $p(i,j)$

Sort $p$ and get top 1 index $idx_1$.

$R_{new}$ =$idx_1$

§ EXPERIMENTS

In this section we extensively evaluate the performance of on multiple synthetic and real datasets, and compare its performance with state-of-the-art approaches. We experiment on the different parameters of and the baselines, and how that affects performance. We implemented in Matlab using the functionality of the Tensor Toolbox for Matlab <cit.>  which supports very efficient computations for sparse tensors. Our implementation is available at link []. We used Intel(R) Xeon(R), CPU E5-2680 v3 @ 2.50GHz machine with 48 CPU cores and 378GB RAM.

§.§ Data-set description

Below, we describe the process for generating synthetic data and the real datasets we used.

§.§.§ Synthetic data generation

We generate random tensors of dimension $I=J=K$ with increasing $I$. where I $\in$ [100, 500, 1000, 3000, 5000, 10000, 50000, 100000]. Those tensors are created from a known set of randomly generated factors, so that we have full control over the ground truth of the full decomposition. We dynamically calculate the size of batch or slice for our all experiments to fit the data into memory. The machine used in this experiments has $\approx$380GB of memory. For example, in case of $I=J=K=50000$, batch size is up to 5, the dense incoming slice requires memory space equivalent to $\approx$80 GB  and the rest of memory space is used for the computations.

The specifications of each synthetic dataset are given in Table <ref>.

Dimension      Density-     Density-    Batch     Sampling

($I=J=K$)    dense     sparse    size    factor

100     100%   65%   50    2

500     100%   65%   150   2

1000     100%   55%   150   2

3000     100%   55%   100   5

5000     100%   55%   100   5

10000    100%   55%   10   2

50000     100%   35%   5   2

100000     100%   35%   5   2

Table of Datasets analyzed

§.§.§ Real Data Description

In order to truly evaluate the effectiveness of , we test its performance against six real datasets that have been used in the literature. Those datasets  are summarized in Table <ref> and are publicly available at <cit.> <http://frostt.io/tensors>.

NIPS, NELL, Facebook-Wall, Facebook-links, Patents, and Amazon dataset. These datasets are publicly available at <cit.>. NIPS publication dataset is collected by Globerson et al.<cit.> and consists of papers published from 1987 to 2003 in NIPS. The modes are “paper”, “author” and “word relation”.
NELL <cit.> is collected as part of the Read the Web project at Carnegie Mellon University. It is an entity-relation-entity tuple snapshot of the Never Ending Language Learner knowledge base.
Facebook Wall posts dataset was first used in <cit.> and has modes “Wall owner”, “Poster”, and
“day”, where the Poster created a post on the Wall owner’s Wall on the specified timestamp.
Facebook Links posts <cit.> contains a list of all of the user-to-user links from the Facebook New Orleans sub-network.
Patents dataset is pairwise co-occurence of terms within window of 7 words in the US utility patents on a year basis. The modes of the tensor represents year-term-term, and the values are $\log(1+f_{i,j})$ , where $f_{i,j}$ is the frequency in which the words $i$ and $j$ appeared in the window.  Each slice of the tensor is symmetric.
Finally, Amazon dataset is collected by SNAP <cit.> and consists of product review from Amazon where mode are in form of user-product-word.

Name   	Description    	Dimensions   	NNZ      Batch     Sampling     Dataset

size     factor     File Size

NIPS <cit.>   (Paper,Author,Word)	   2.4K x 2.8K x 14Km	   3,101,609   500   10   57MB

NELL	 <cit.>   (Entity,Relation,Entity)	    12K x 9K x 28K   76,879,419   500   10   1.4GB

Facebook-wall <cit.>   (Wall owner, Poster, day)	   62K x 62K x 1,070   78,067,090   100   5   2.1GB

Facebook-links <cit.>    (User, Links, Day)	   62K	x  62K x	650	   263,544,295   50   2   3.8GB

Patents<cit.>	    (Term ,Term, Year)     239K x 239K x 46    3,596,640,708   10   2   73GB

Amazon<cit.>	    (User, Product ,Word)    4.8M x	1.7M x 1.8M   	1,741,809,018   50000   20   43GB

Real datasets analyzed

§.§ Evaluation Measures

We evaluate and the baselines using three criteria: Relative Error, Wall-Clock time and Fitness. These measures provide a quantitative way to compare the performance of our method.
More Specifically, Relative Error is effectiveness measurement and defined as :

\[
Relative Error=\frac{||\tensor{X}_{original}-\tensor{X}_{predicted}||}{||\tensor{X}_{original}||}
\]

where, the lower the value, the better.

CPU time (sec): indicates how much faster does the decomposition runs as compared to re-running the entire decomposition algorithm whenever we receive a new update on the existing tensor. The average running time denoted by $T_{tot}$ for processing all slices for given tensor, measured in seconds, and is used to validate the time efficiency of an algorithm.

Relative Fitness:
Relative Fitness is defined as:

\[
Relative Fitness=\frac{||\tensor{X}_{original}-\tensor{X}_{\sambaten}||}{||\tensor{X}_{original}-\tensor{X}_{BaseLine}||}
\]

where, again, lower is better.

§.§ Baselines for Comparison

Here we briefly present the state-of-the-art baselines we used for comparison. Note that for each baseline we use the reported parameters that yielded the best performance in the respective publications. For fairness, we compare against the parameter configuration for that yielded the best performance in terms of low wall-clock timing, low relative error and fitness.

Note that all comparisons were carried out over 10 iterations each, and each number reported is an average with a standard deviation attached to it.

CP_ALS <cit.>: is considered the most standard and well optimized algorithm for CP. We use the implementation of the Tensor Toolbox for Matlab <cit.>. Here, we simply re-compute CP using CP_ALS after every update.

SDT <cit.>: Simultaneous Diagonalization Tracking (SDT) is based on incrementally tracking the Singular Value Decomposition (SVD) of the unfolded tensor $\tensor{X}_{(3)}=U\sum V^T$.

RLST <cit.>: Recursive Least Squares Tracking (RLST) is another online approach in which recursive updates are computed to minimize the Mean Squared Error (MSE) on incoming  slice.

OnilneCP <cit.>: This is the most recent and state-of-the-art method in online computation of CP.

OnilneCP fixes $\mathbf{A}$ and $\mathbf{B}$ to solve for $\mathbf{C}$ and minimizes the cost as:

\[\mathbf{C}\leftarrow \argmin_c \frac{1}{2} ||\left(\begin{smallmatrix} \tensor{X}_{old(3)}\\  \tensor{X}_{new(3)}\end{smallmatrix}\right)-\left(\begin{smallmatrix} \mathbf{C}_{old}\\  \mathbf{C}_{new}\end{smallmatrix}\right)(\mathbf{B} \odot \mathbf{A})^T||^2
\]

Then $\mathbf{C}$ is updated as :

\[
\mathbf{C}=\left(\begin{smallmatrix} \mathbf{C}_{old}\\  \mathbf{C}_{new}\end{smallmatrix}\right)=\left[\left(\begin{smallmatrix} \mathbf{C}_{old}\\  \tensor{X}_{new(3)}\end{smallmatrix}\right)-((\mathbf{B} \odot \mathbf{A})^T)^\ddagger\right]
\]

Matrix $\mathbf{A}$ is updated as $\mathbf{A}=\mathbf{P}\mathbf{Q}^{-1}$ where $\mathbf{P}=\mathbf{P}_{old}+\tensor{X}_{new(1)}(\mathbf{C}_{new} \odot \mathbf{B})$ and Q is $(\mathbf{C}^T\mathbf{C} \otimes \mathbf{B}^T\mathbf{B})$. Matrix $\mathbf{B}$ is estimated as $\mathbf{B}=\mathbf{U}\mathbf{V}^{-1}$ where $\mathbf{U}=\tensor{X}_{(2)}(\mathbf{C}\odot \mathbf{A})$ and $\mathbf{V}=(\mathbf{C} \odot \mathbf{A})^T (\mathbf{C} \odot \mathbf{A})$.

We conduct our experiments on multiple synthetic datasets and six real-world tensors datasets. We set the tolerance rate for convergence between consecutive iterations to $10^{-5}$ and the maximum number of iteration to 1000 for all the algorithms. The batch size and sampling factor is selected based on dimensions of first mode i.e. $I$, provided in Table <ref> and <ref> for synthetic and real dataset respectively.
We use the publicly available implementations for the baselines, as provided by the authors. We only modified the interface of the baselines, so that it is consistent across all methods with respect to the way that they receive the incoming slices. No other functionality has been changed.

§.§ Experimental Results

The following major three aspects are analyzed.

Q1. Effectiveness and Accuracy: How effective is as compared to the baselines on different synthetic and real world datasets?

Q2. Speed & Scalability: How fast is when compared to the state-of-the-art methods on very large sized datasets?

Q3: Effectiveness and Accuracy What is the cost-benefit trade-off of computing the actual rank of the incoming batch?

Q3. Parameter Sensitivity: What is the influence of sampling factor $s$ and sampling repetitions $r$?

§.§.§ Baselines for Comparison

For all datasets we compute Relative Error,CPU time (sec) and Fitness. For , CP_ALS , OnlineCP, RSLT and SDT we use 10% of the data in each dataset as existing dataset. We experimented for both dense as well as sparse tensor to check the performance of our method. The results for the dense and sparse synthetic data are shown in Table <ref> - <ref>. For each of datasets , the best result is shown in bold. OnlineCP, SDT and RLST address the issue very well. Compared
with CP_ALS, SDT and RLST reduce the mean running time by up to 2x times and OnlineCP reduce mean time by up to $3\times$ for small dataset (I up to 3000). Performance of RLST was better than SDT algorithm on 8 out of 8 third-order synthetic tensor datasets. In fact, the efficiency (in terms of CPU time (sec)) of SDT is quite close to RLST. However, the main issue of SDT and RLST is their estimation of relative error and fitness.  For some datasets, such as $I=100$ and $I=3000$, they perform well, while for some others, they exhibit  poor fitness and relative error, achieving only nearly half of the fitness of other methods. For small size datasets, OnlineCP's efficiency and accuracy is better than all methods. As the dimension grows, however, the performance of OnlineCP method reduces,and particularly for datasets of dimension larger than  $5000 \times 5000 \times 5000$. Same behavior is observed for sparse tensors. is comparable to baselines for small dataset and outperformed the  baselines for large dataset. CP_ALS is the only baseline able to run on datasets up to size $3000 \times 3000 \times 3000$. These results answer Q1 as the have comparable accuracy to other baseline methods.

I=J=K   $CP_{ALS}$    OnlineCP   SDT   RLST

100   0.109 $\pm$ 0.01   0.107$\pm$ 0.02   	0.173$\pm$ 0.02   	0.151$\pm$ 0.02   0.115$\pm$ 0.02

500   	0.102 $\pm$ 0.09   	0.102$\pm$ 0.09   	0.217$\pm$ 0.06	   0.217$\pm$ 0.06   0.102$\pm$ 0.09

1000   0.103$\pm$ 0.01   	0.103$\pm$ 0.01   	0.287$\pm$ 0.01	   0.296$\pm$ 0.01   0.102$\pm$ 0.01

3000   0.119$\pm$ 0.01   	0.108$\pm$ 0.01   0.189$\pm$ 0.01	   0.206$\pm$ 0.01   0.109$\pm$ 0.01

5000	   N/A   	0.122$\pm$ 0.002   0.201$\pm$ 0.002	   0.196$\pm$ 0.04   0.115$\pm$ 0.009

10000   N/A   	0.173$\pm$ 0.04   0.225$\pm$ 0.04	   0.252$\pm$ 0.06   0.162$\pm$ 0.01

50000	   N/A   0.215$\pm$ 0.03   0.229$\pm$ 0.03   0.26$\pm$ 0.01   0.169$\pm$ 0.01

100000   N/A   N/A   	N/A   	N/A   0.275 $\pm$ 0.00

Experimental results for relative error for synthetic dense tensor. We see that gives comparable accuracy to baseline.

Experimental results for CPU time (sec) for (a) dense tensor (b) sparse tensor

Experimental results for relative fitness improvement for (a) dense tensor (b) sparse tensor

I=J=K   $CP_{ALS}$    OnlineCP   SDT   RSLT

100   0.169$\pm$ 0.01   0.154$\pm$ 0.02   0.306$\pm$ 0.01   0.313$\pm$ 0.01   0.178$\pm$ 0.01

500   0.175$\pm$ 0.01   0.188$\pm$ 0.01   0.43$\pm$ 0.01   0.421$\pm$ 0.01   0.184$\pm$ 0.01

1000   0.179$\pm$ 0.01   0.185$\pm$ 0.01   0.613$\pm$ 0.01   0.813$\pm$ 0.01   0.178$\pm$ 0.01

3000   0.177$\pm$ 0.02   0.171$\pm$ 0.04   0.446$\pm$ 0.03   0.513$\pm$ 0.02   0.176$\pm$ 0.03

5000   N/A   0.192$\pm$ 0.02   0.494$\pm$ 0.09   0.535$\pm$ 0.13   0.187$\pm$ 0.04

10000   N/A   0.173$\pm$ 0.01   0.212$\pm$ 0.01   0.224$\pm$ 0.11   0.198$\pm$ 0.12

50000   N/A   0.222$\pm$ 0.00   0.262$\pm$ 0.00   0.259$\pm$ 0.00   0.200$\pm$ 0.00

100000   N/A   N/A   	N/A   	N/A   0.283$\pm$ 0.00

Experimental results for relative error for synthetic sparse tensor. We see that works better in very large scale dataset such as 50000 $\times$ 50000 $\times$ 50000.

is efficiently able to compute $100K \times 100K \times 100K$ sized tensor with batch size of 5 and sampling factor 2. It took 58095.72s and 47232.2s to compute online decomposition for dense and sparse tensor, respectively. On other hand, state-of-art methods are unable to handle such large scaled incoming data.

Table <ref> shows the comparison between methods. outperforms other state-of-the-art approaches in most of  the real dataset. In the case of NIPS datset, gives better results compared to the baselines, specifically in terms of CPU Time (faster up to 20 times) and Fitness (better up to 15-20%). outperforms for NELL, Facebook-Wall and Facebook-Links dataset in terms of efficiency comparable to CP_ALS. For the NIPS dataset, is 25-30 times faster than OnlineCP method. Due to high dimensions of dataset, RSLT and SDT are unable to execute further. Note that all the real datasets we use are highly sparse, however, no baselines except CP_ALS actually take advantage of that sparsity, therefore, repeated CP_ALS tends to be faster because the baselines have to deal with dense computations which tend to be slower, when the data contain a lot of zeros.  Most importantly, however, performed very well on Amazon and Patent datasets, arguably the hardest of the six real datasets we examined and have been analyzed in the literature, where none of the baselines was able to run. These result answer Q1 and Q2 and show that is able to handle large dimensions and sparsity.

Dataset   5|c|CPU Time (sec)   4|c|Fitness w.r.t

$CP_{ALS}$    OnlineCP   SDT   RSLT       $CP_{ALS}$   OnlineCP   SDT   RSLT
NIPS   177.46$\pm$ 2.9   372.03$\pm$ 8.9   1608.23$\pm$ 37.5   1596.07$\pm$ 15.2   43.98$\pm$ 0.6   0.96$\pm$ 0.01   0.98$\pm$ 0.01   0.78$\pm$ 0.02   0.82$\pm$ 0.01

NELL   8783.27$\pm$ 11.2   42325.22$\pm$ 70.0   65325.22$\pm$ 25.2   63485.98$\pm$ 10.6   983.83$\pm$ 30.7   0.95$\pm$0.02   0.81$\pm$ 0.01   0.76$\pm$ 0.02   0.81$\pm$ 0.01

Facebook-wall   3041.98$\pm$3.8   N/A   N/A   N/A   736.07$\pm$4.1   0.97$\pm$0.01   N/A   N/A   N/A

Facebook-links   2689.69$\pm$7.9   N/A   N/A   N/A   343.32$\pm$6.3   0.96$\pm$0.06   N/A   N/A   N/A

Amazon   N/A   N/A   N/A   N/A   4892.07$\pm$61.8   N/A   N/A   N/A   N/A

Patent     N/A   N/A   N/A   N/A   8068.27$\pm$55.4    N/A   N/A   N/A   N/A

performance for real datasets. outperforms the baselines for all the large  tensors.

§.§.§ Sensitivity of Sampling Factor s

The sampling factor plays an important role in . We performed experiments to evaluate the impact of changing sampling factor on . For these experiments , we fixed batch size to 50 for all datasets. We see in figure <ref> that increasing sampling factor results in reduction of CPU time (as sparsity of sub sampled tensor increased) and it reduces the fitness of output up-to 2-3%. In sum, these observations demonstrate that: 1) a suitable sampling factor on sub-sampled tensor could improve the fitness and result in better tensor decomposition, and 2) the higher sampling factor is, the lower the CPU time. This result partially answers Q3.

CPU Time (sec)  and Relative Fitness vs. Sampling Factor $s$ on different datasets (lower is better).

§.§.§ Sensitivity of Repetition Factor r

We evaluate the performance for parameter setting $r$ i.e the number of parallel decompositions.  For these experiments, we choose batch size and sampling rate for synthetic $ 500 \times 500 \times 500$ dataset and  NIPS real world dataset as provided in table <ref> and <ref>, respectively. We can see that with higher values of the repetition factor $r$, FMS score  and  Relative Fitness (vs CP_ALS) is  improved as shown in Figure <ref>(a). We experiment on varying repetition factor r with Sampling factor s on NIPS real world dataset to check the performance of our method as shown in Figure <ref>(b). The lower the relative fitness score, the better the decomposition. This result completes the answer to Q3.

Relative Fitness vs. repetition factor $r$ on synthetic and NIPS datasets (lower is better).

§ CONCLUSIONS

We introduce , a novel sample-based incremental CP tensor decomposition. We show its effectiveness with respect to approximation quality, with its performance being on par with state-of-the-art incremental and non-incremental algorithms, and we demonstrate its efficiency and scalability by outperforming state-of-the-art approaches (25-30 times faster) and being able to run very large incremental tensors where none of the baselines was able to produce results.

Chapter based on material published in SDM 2018 <cit.>.

CHAPTER: ONLINE COMPRESSION-BASED TENSOR DECOMPOSITION

[backgroundcolor=Orange!20,linewidth=1pt,  topline=true,  rightline=true, leftline=true]
"How to leverage random compression for streaming high-order tensor decomposition? Can we guarantee the identifiability and speed-up with parallelism?”

Tensor decompositions are powerful tools for large data analytics as they jointly model multiple aspects of data into one framework and enable the discovery of the latent structure and higher-order correlations within the data. One of the most widely studied and used decompositions, especially in data mining and machine learning, is the Canonical Polyadic or CP decomposition. However, today's datasets are not static and these datasets often dynamically growing and changing with time. To operate on such large data, we present the first ever compression-based online parallel implementation for the CP decomposition.We conduct an extensive empirical analysis of the algorithms in terms of fitness, memory used and CPU time, and in order demonstrate the compression and scalability of the method, we apply to big tensor data. Indicatively, performs on-par or better than to state-of-the-art online and offline methods in terms of decomposition accuracy and efficiency, while saving up to 40-250 % memory space.	The content of this chapter is adapted from the following published paper:

Gujral, Ekta, Ravdeep Pasricha, Tianxiong Yang, and Evangelos E. Papalexakis. "Octen: Online compression-based tensor decomposition." In 2019 IEEE 8th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), pp. 455-459. IEEE, 2019.

§ INTRODUCTION

A Tensor is a multi-way array of elements that represents higher-order or multi-aspect data. In recent years, tensor decompositions have gained increasing popularity in big data analytics <cit.>. In higher-order structure, tensor decomposition are capable of finding complex patterns and  higher-order correlations within the data. Much like matrix factorization tools like the Singular Value Decomposition, there exist generalizations for the tensor domain, with the most widely used being CANDECOMP/PARAFAC or CP <cit.> which extracts interpretable factors form tensor data, and Tucker decomposition <cit.>, which is known for estimating the joint subspaces of tensor data. In this work we focus only on the CP decomposition, which is extremely effective in exploratory knowledge discovery on multi-aspect data.

In the era of information explosion, data are generated or modified at large volume. In such  environments, data may be increased or decreased in any of its dimensions with high velocity. When using tensors to represent this dynamically changing data, an instance of the problem is that of a “streaming”, “incremental”, or “online” tensors[Notice that the literature (and thereby this chapter) uses the above  terms as well as “dynamic” interchangeably.]. Considering an example of time evolving social network interactions, where a large number of users interact with each other every second (Facebook users update $\approx 684K$ piece of content and Twitter users send $\approx 100K$ tweets every single minute[http://mashable.com/2012/06/22/data-created-every-minute/]); every such snapshot of interactions is a new incoming slice(s) to the tensor on its “time” or mode, which is seen as a streaming update. Additionally, the tensor may be growing in all of its n-modes, especially in complex and evolving environments such as online social networks. In this chapter, our goal is, given an already computed CP decomposition, to track the CP decomposition of an online tensor, as it receives streaming updates, 1) efficiently,  being much faster than re-computing the entire decomposition from scratch every update, and utilizing small amount of memory, and 2) accurately, incurring an approximation error that is as close as possible to the decomposition of the full tensor. For exposition purposes, we focus on the  time-evolving/streaming scenario, where a three-mode tensor grows on the third (“time”) mode, however, our work extends to cases where more than one modes is online.

framework. Compressed tensor summaries $\tensor{Y_p}$ and $\tensor{Z_p}$  is obtained by applying compression matrices $(\mathbf{U}_p, \mathbf{V}_p,\mathbf{W}_p)$ and $(\mathbf{U}_p^{'}, \mathbf{V}_p^{'},\mathbf{W}_p^{'})$ to $\tensor{X_{old}}$ and $\tensor{X_{new}}$ or incoming slice(s) respectively. The updated summaries are computed by $\tensor{X_p}$ = $\tensor{Y_p}+\tensor{Z_p}$. Each $\tensor{X_p}$ is independently decomposed in parallel. The update step anchors all factor matrices to a single reference, and solves a linear equation for the overall A, B, and C.

As the volume and velocity of data grow, the need for time- and space-efficient online tensor decomposition is imperative. There already exists a modest amount of prior work in online tensor decomposition both for Tucker <cit.> and CP <cit.>. However, most of the existing online methods <cit.> We *may* wanna say here “with the exception of SamBaTen” but maybe not, because we may distract from the main point., model the data in the full space, which can become very memory taxing as the size of the data grows. There exist memory efficient tensor decompositions, indicatively MET for Tucker <cit.> and PARACOMP <cit.> for CP, neither of which are able to handle online tensors. In this chapter, we fill exactly that gap.

Online tensor decomposition is a challenging task due to the following reasons. First, maintaining high-accuracy (competitive to decomposing the full tensor) using significantly fewer computations and memory than the full decomposition calls for innovative and, ideally, sub-linear approaches. Second, operating on the full ambient space of the data, as the tensor is being updated online, leads to super-linear increase in time and space complexity, rendering such approaches hard to scale, and calling for efficient methods that work on spaces that are significantly smaller than the original ambient data dimensions. Third, in many real settings, more than one modes of the tensor may receive streaming updates at different points in time, and devising a flexible algorithm that can handle such updates seamlessly is a challenge.
To handle the above challenges, in this chapter, we propose to explore how to decompose online or incremental tensors based on CP decomposition. We specifically study: (1) How to  make parallel update method based on CP decomposition for online tensors? (2) How to identify latent component effectively and accurately after decomposition? Answering the above questions, we propose (Online Compression-based Tensor Decomposition) framework. Our contributions are summarized as follows:

* Novel Parallel Algorithm We introduce , a novel compression-based algorithm for online tensor decomposotion that admits an efficient parallel implementation. We do not limit to 3-mode tensors, our algorithm can easily handle higher-order tensor decompositions.

* Correctness guarantees By virtue of using random compression, can guarantee the identifiability of the underlying CP decomposition in the presence of streaming updates.

* Extensive Evaluation Through experimental evaluation on various datasets, we show that provides stable decompositions (with quality on par with state-of-the-art), while offering up to 40-250 % memory space savings.

Reproducibility: We make our Matlab implementation publicly available at link []. Furthermore, all the small size datasets we use for evaluation are publicly available on same link.

§ RELATED WORK

In this section, we provide review of the work related to our algorithm. At large,  incremental tensor methods in the literature can be categorized into two main categories as described below:

Tensor Decomposition: Phan el at. <cit.> had purposed a theoretic method namely GridTF to large-scale tensors decomposition based on CP's basic mathematical theory to get sub-tensors and join the output of all decompositions to achieve final factor matrices. Sidiropoulos el at.<cit.>, proposed algorithm that focus on CP decomposition namely RLST (Recursive Least Squares Tracking), which recursively updates the factors by minimizing the mean squared error. In 2014, Sidiropoulos el at. <cit.> , proposed a parallel algorithm for low-rank tensor decomposition that is suitable for large tensors. The Zhou, el at. <cit.> describes an online CP decomposition method, where the latent components are updated for incoming data. The most related work to ours was proposed by <cit.> which is sampling-based batch incremental tensor decomposition algorithm. These state of art techniques focus on only fast computation but not effective memory usage. Besides CP decomposition, tucker decomposition methods<cit.> were also introduced. Some of these methods were not only able to handle data increasing in one-mode, but also have solution for multiple-mode updates using methods such as incremental SVD <cit.>. Latest line of work is introduced in <cit.> i.e TuckerMPI to find inherent low-dimensional multi-linear structure, achieving high compression ratios. Tucker is mostly focused on recovering subspaces of the tensor, rather than latent factors, whereas our focus is on the CP decomposition which is more suitable for exploratory analysis.

Tensor Completion:
Another field of study is tensor completion, where real-world large-scale datasets are considered incomplete. In literature, wide range of methods have been proposed based on online tensor imputation<cit.> and  tensor completion with auxiliary information<cit.>. The most recent method in this line of work is by Qingquan el at.<cit.>,  who proposed a low-rank tensor completion with general multi-aspect streaming patterns, based on block partitioning of the tensor. However, these approaches can not directly applied when new batches of data arrived. This provides us a good starting point for further research.

§ PROBLEM FORMULATION

In many real-world applications, data grow dynamically and may do so in many modes. For example, given a dynamic
tensor in a movie-based recommendation system, organized as users $\times$ movie $\times$ rating $\times$ hours,the number of registered users, movies watched or rated, and hours may all increase over time. Another example is network monitoring sensor data where tons of information like source and target IP address, users ,ports etc., is collected every second. This nature of data gives rise to update existing decompositions on the fly or online and we call it incremental decomposition [We uses the terms “incremental”, “dynamic”, and “online” interchangeably in chapter].  In such conditions, the  update needs to process the new data very quickly, which makes non-incremental methods to fall short because they need to recompute the decomposition for the full dataset.
The problem that we solve is the following:

[linecolor=red!60!black,backgroundcolor=gray!20,linewidth=1pt,    topline=true,rightline=true, leftline=true]

Given (a) an existing set of summaries {$\tensor{Y}_{1},\tensor{Y}_{2} \dots \tensor{Y}_{p}$} of $R$ components, which approximate tensor $\tensor{X}_{old}$ of size { $ I^{(1)} \times I^{(2)} \times \dots I^{(N-1)} \times  t_{old}$} at time t , (b) new incoming batch of slice(s) in form of tensor $\tensor{X}_{new}$ of size {$ I^{(1)} \times I^{(2)} \times \dots I^{(N-1)} \times  t_{new}$}, find updates of ($\mathbf{A}^{(1)},\mathbf{A}^{(2)}$ , …, $\mathbf{A}^{(N-1)}$, $\mathbf{A}^{(N)}$) incrementally to approximate tensor $\tensor{X}$ of dimension {$ I^{(1)} \times I^{(2)} \times \dots I^{(N-1)} \times I^{(N)}$}, where $ I^{(N)} = (t_{old}+t_{new})= I_{1\dots n}^{(N)} + I_{(n+1)\dots m}^{(N)}$ after appending new slice or tensor to $N^{th}$ mode.

§ PROPOSED METHOD: OCTEN

As we mention in the previous section, to the best of our knowledge, there is no algorithm in the literature that is able to efficiently compress and incrementally update the CP decomposition in the presence of incoming tensor slices. However, there exists a method for static data <cit.>. Since this method considers the tensor in its entirety, it cannot handle streaming data and as the data size grows its efficiency may degrade, since it handles the full data in one shot. In this section, we introduce , a new method for parallel incremental decomposition designed with two main goals in mind: G1:  Compression, speed , simplicity, and parallelization; and G2: correctness in recovering compressed partial results for incoming data, under suitable conditions.

The algorithmic framework we propose is shown in Figure <ref> and is described below:

We assume that we have a pre-existing set of summaries of the $\tensor{X}$ before the update. Summaries are in the form of compressed tensors of dimension [$Q \times Q \times Q$]. For simplicity of description, we assume that we are receiving updated slices on the third mode. We, further, assume that the updates come in batches of new slices, which, in turn, ensures that we see a mature-enough update to the tensor, which contains useful structure. Trivially, however, can operate on singleton batches and for $>3$ modes also.

In the following, $\tensor{X}_{old}$ is the tensor prior to the update and $\tensor{X}_{new}$ is the batch of incoming slice(s). Considering $S=\prod_{i=1}^{[N-1]}I^{(i)}$ and $T=\sum_{i=1}^{[N-1]}I^{(i)}$ , we can write space and time complexity in terms of $S$ and $T$. Given an incoming batch, performs the following steps:

Parallel Compression and Decomposition:
When handling large tensors $\tensor{X}$ that are unable fit in main memory, we may compress the tensor $\tensor{X}$ to a smaller tensor that somehow apprehends most of the systematic variation in  $\tensor{X}$. Keeping this in mind, for incoming slice(s) $\tensor{X}_{new}$, during the parallel compression step, we first need to create 'p' parallel triplets of random compression matrices {$\mathbf{U_p}, \mathbf{V_p}, \mathbf{W_p}$} of $\tensor{X}$.  Thus, each worker (i.e. Matlab parpool) is responsible for creating and storing the three compression matrices of dimension $\mathbb{R}^{I \times Q},\ \ \mathbb{R}^{J \times Q}$ and $\mathbb{R}^{t_{new} \times Q}$. These matrices share at least 'shared' amount of columns among each other. Mathematically, we can describe it as follows:

\begin{equation}
\tensor{X}=  \begin{bmatrix}
\{\mathbf{U_1}, \mathbf{V_1},\mathbf{ W_1}\} \\
\{\mathbf{U_2}, \mathbf{V_2}, \mathbf{W_2}\}\\
\dots\\
\{\mathbf{U_p}, \mathbf{V_p}, \mathbf{W_p}\}
\end{bmatrix}
= \begin{bmatrix}
\{(\mathbf{u} \ \ \mathbf{U_{1'}}),(\mathbf{v} \ \ \mathbf{V_{1'}}),(\mathbf{w}  \ \ \mathbf{W_{1'}})\} \\
\{(\mathbf{u} \ \ \mathbf{U_{2'}}),(\mathbf{v} \ \ \mathbf{V_{2'}}),(\mathbf{w}  \ \ \mathbf{W_{2'}})\} \\
\dots\\
\{(\mathbf{u} \ \ \mathbf{U_{p'}}),(\mathbf{v} \ \ \mathbf{V_{p'}}),(\mathbf{w}  \ \ \mathbf{W_{p'}})\} \\
\end{bmatrix}
\end{equation}

where u, v and w are shared and have dimensions of $\mathbb{R}^{I \times Q_{shared}}, \mathbb{R}^{J \times Q_{shared}}$ and $\mathbb{R}^{t_{new} \times Q_{shared}}$.

For compression matrices, we choose to assign each worker to create a single row of each of the matrices to reduce the burden of creating an entire batch of {$\mathbf{U_p^{'}}, \mathbf{V_p^{'}}, \mathbf{W_p^{'}}$} of $\tensor{X}_{new}$. We see that for each worker is sufficient to hold {$\mathbf{U_p}, \mathbf{V_p}, \mathbf{W_p}$} in memory. Now, we created compressed tensors or replicas $\{\tensor{Z_1}, \tensor{Z_2} \dots \tensor{Z_p}\}$ from each triplets of 3-mode compression matrices generated from $\tensor{X}_{new}$;see Figure <ref>.  $\tensor{Z}_{p}$ is 3-mode tensor of size  $\mathbb{R}^{Q \times Q \times Q}$. Since $Q$ is considerably smaller than [I ,J,  K], we use $ O(Q^3)$ of memory on each worker.

For $\tensor{X}_{old}$, we already have replicas $\{\tensor{Y_1}, \tensor{Y_2} \dots \tensor{Y_p}\}$ from each triplets of 3-mode compression matrices {$\mathbf{U_p}, \mathbf{V_p}, \mathbf{W_p}$} generated from $\tensor{X}_{old}$;see Figure <ref>. In general, the compression comprises N-mode products, leads to complexity of $(Q_{(1)}St_{new} + Q_{(2)}St_{new} +Q_{(3)}St_{new} + \dots Q_{(N)}St_{new})$ overall for dense tensor $\tensor{X}$, if the first mode is compressed first, followed by the second, and then the third mode and so on. We choose to keep $Q_1, Q_2,Q_3 \dots Q_N$ of same order as well non-temporal dimensions are of same order in our algorithm, so time complexity of parallel compression step for N-mode data is $O(QSt_{new})$ for each worker. The summaries are always dense, because first mode product with tensor is dense, hence remaining mode products unable to exploit sparsity.

After appropriately computing summaries $\{\tensor{Z_1}, \tensor{Z_2} \dots \tensor{Z_p}\}$ for incoming slices, we need to update the old summaries $\{\tensor{Y_1}, \tensor{Y_2} \dots \tensor{Y_p}\}$. Each worker reads its segment and processing update parallel as given below.

\begin{equation}
\begin{bmatrix}
\tensor{X_1} \\
\tensor{X_2}\\
\vdots\\
\tensor{X_p}
\end{bmatrix}
\begin{bmatrix}
\tensor{Y_1} \\
\tensor{Y_2}\\
\vdots\\
\tensor{ Y_p}
\end{bmatrix}
\oplus
\begin{bmatrix}
\tensor{Z_1} \\
\tensor{Z_2}\\
\vdots\\
\tensor{Z_p}
\end{bmatrix} \odot  \begin{bmatrix}
\mathbf{W^{'}_{1}(k,q)} \\
\mathbf{W^{'}_{2}(k,q)}\\
\vdots\\
\mathbf{W^{'}_{p}(k,q)}
\end{bmatrix} \implies \begin{bmatrix}
(\mathbf{A_{s(1)}}, \mathbf{B_{s(1)}} , \mathbf{C_{s(1)}}) \\
(\mathbf{A_{s(2)}}, \mathbf{B_{s(2)}} , \mathbf{C_{s(2)}})\\
\vdots\\
(\mathbf{A_{s(p)}}, \mathbf{B_{s(p)}} , \mathbf{C_{s(p)}})
\end{bmatrix}
\end{equation}

Where $k$ is the number of slices of the incoming tensor and $q$ is the slice number for thecompressed tensor. Further, for the decomposition step, we processed 'p' summaries on different workers, each one fitting the decomposition to the respective compressed tensor $\{\tensor{X_1}, \tensor{X_2} \dots \tensor{X_p}\}$ created by the compression step.  We assume that the updated compressed tensor $\{\tensor{X_1}, \tensor{X_2} \dots \tensor{X_p}\}$ fits in the main memory, and performs in-memory computation.  We denote $p^{th}$ compressed tensor decompositions as $(\mathbf{A_{s(p)}}, \mathbf{B_{s(p)}} , \mathbf{C_{s(p)}})$ as discussed above. The data for each parallel worker $\tensor{X_p}$ can be uniquely decomposed, i.e. $(\mathbf{A_p,B_p,C_p})$ is unique up to scaling and column permutation.

Furthermore, parallel compression and decomposition is able to achieve Goal G1.

for incremental tensor decomposition

$\tensor{X}_{new} \in \mathbb{R}^{I \times J \times K_{(n+1)\dots m}}$,summary $ \tensor{Y}_i \in \mathbb{R}^{Q \times Q \times Q}$,  R,p,Q, shared S.
Factor matrices $\mathbf{A}, \mathbf{B}, \mathbf{C}$ of size $(I \times R)$, $(J \times R)$ and $(K_{1\dots n,(n+1) \dots m} \times R)$.

new slice(s) coming

$\tensor{Z}_i \leftarrow \{(\mathbf{U}_i^{'},\mathbf{V}_i^{'}, \mathbf{W}_i^{'})\}$ ,  $ \tensor{Z}_i \in \mathbb{R}^{Q \times Q \times Q}$,  $i \in (1,p)$

$\tensor{X}_i\leftarrow \tensor{Y}_i\oplus \tensor{Z}_i$ ,   $ \tensor{X}_i \in \mathbb{R}^{Q \times Q \times Q}$,  $i \in (1,p)$
$(\widetilde{\mathbf{A}}_{s(i)},\widetilde{\mathbf{B}}_{s(i)},\widetilde{\mathbf{C}}_{s(i)}) \leftarrow CP(\tensor{X}_i,R),\ \ i \in (1,p)$
$(\widetilde{\mathbf{P}}_{a(i)},\widetilde{\mathbf{P}}_{b(i)},\widetilde{\mathbf{P}}_{c(i)})\leftarrow \{(\mathbf{U}^{'}(i,[S,:],:)^T,\mathbf{V}^{'}(i,[S,:],:)^T,\mathbf{W}^{'}(i,[S,:],:)^T)\},\ \ i \in (1,p)$
$i \leftarrow 1$ to $p-1$
$(\widetilde{\mathbf{A}}_{s},\widetilde{\mathbf{B}}_{s},\widetilde{\mathbf{C}}_{s})\leftarrow \Pi\big[(\widetilde{\mathbf{A}}_{s(i)}, \widetilde{\mathbf{B}}_{s(i)} , \widetilde{\mathbf{C}}_{s(i)}) \ \  {;} \ \ (\mathbf{A}_{s(i+1)}, \mathbf{B}_{s(i+1)} , \mathbf{C}_{s(i+1)})\big]$

$(\widetilde{\mathbf{P}}_{a},\widetilde{\mathbf{P}}_{b},\widetilde{\mathbf{P}}_{c}) \leftarrow \big[(\widetilde{\mathbf{P}}_{a(i)},\widetilde{\mathbf{P}}_{b(i)},\widetilde{\mathbf{P}}_{c(i)}) \ \  {;} \ \ (\widetilde{\mathbf{P}}_{a(i+1)},\widetilde{\mathbf{P}}_{b(i+1)},\widetilde{\mathbf{P}}_{c(i+1)})\big]$

$\mathbf{A} \leftarrow \widetilde{\mathbf{P}}_a^{-1}*\widetilde{\mathbf{A}}_s$ ; $\mathbf{B} \leftarrow \widetilde{\mathbf{P}}_b^{-1}*\widetilde{\mathbf{B}}_s$ ;   $\mathbf{C} \leftarrow [\mathbf{C}_{old}; \widetilde{\mathbf{P}}_c^{-1}*\widetilde{\mathbf{C}}_s]$

Factor match for identifiability: According to Kruskal <cit.>, the CP decomposition is unique (under mild conditions) up to permutation and scaling of the components i.e. $\mathbf{A},\mathbf{B}$ and $\mathbf{C}$ factor matrices. Consider an 3-mode tensor $\tensor{X}$ of dimension $I$, $J$ and $K$ of rank $R$. If rank

\begin{equation}
\label{octen:uniq1}
r_c=F \implies K \geq R \  \&  \ I(I-1)(J-1)\geq 2R(R-1),
\end{equation}

then rank 1 factors of tensor $\tensor{X}$  can be uniquely computable<cit.>. With help of Kronecker product<cit.> property i.e. $(\mathbf{U}^T \otimes \mathbf{C}^T \otimes \mathbf{W}^T)(\mathbf{A} \odot\mathbf{B} \odot \mathbf{C}) = ((\mathbf{U}^T\mathbf{A})\odot(\mathbf{V}^T\mathbf{B})\odot(\mathbf{W}^T\mathbf{C})) \approx (\widetilde{\mathbf{A}} ,\widetilde{\mathbf{B}} ,\widetilde{\mathbf{C}})$. Now Kruskal's combining uniqueness and Kronecker product property, we can obtain identifiable factors from summaries if

\begin{equation}
\label{octen:uniq2}
\min(Q,r_A)+\min(Q,r_B)+\min(Q,r_C) \geq 2R+2
\end{equation}

where Kruskal-rank of A,  denoted $r_A$, is the maximum r such that any $r$ columns of $\mathbf{A}$  are linearly independent;see <cit.>. Hence, upon factorization of 3-mode $\tensor{X_p}$ into $R$ rank-1 components, we obtain $ \mathbf{A}= a^T_pA\Pi_{p}\lambda_p^{(1/N)}$ where $a$ is shared among summaries decompositions , $\Pi_p$ is a permutation matrix, and $\lambda_p$ is a diagonal scaling matrix obtained from CP decomposition. To match factor matrices after decomposition step, we first normalize the shared columns of factor matrices $(\mathbf{A_{s(i)}, B_{s(i)} , C_{s(i)}})$ and $(\mathbf{A_{s(i+1)}, B_{s(i+1)} , C_{s(i+1)}})$ to unit norm $||.||_1$ From V: I think you need $\| \|_2$ here, right?, From Ekta: here we need to normalize each column by dividing it with its max value. Basically converting values between 0-1.. Next, for each column of $\mathbf{(A_{s(i+1)},  B_{s(i+1)} , C_{s(i+1)}})$, we find the most similar column of $(\mathbf{A_{s(i)},  B_{s(i)} , C_{s(i)}})$, and store the correspondence. Finally, We can describe factor matrices as :

\begin{equation}
\mathbf{\widetilde{A}_s} = \begin{bmatrix}
\vdots\\
\end{bmatrix} *\mathbf{A} \Pi\lambda^{(1/N)} , \  \ \mathbf{\widetilde{B}_s} = \begin{bmatrix}
\vdots\\
\end{bmatrix} *\mathbf{B} \Pi\lambda^{(1/N)} , \  \ \mathbf{\widetilde{C}_s} = \begin{bmatrix}
\vdots\\
\end{bmatrix} *\mathbf{C} \Pi\lambda^{(1/N)}
\end{equation}

where $\widetilde{\mathbf{A}}_s,\widetilde{\mathbf{B}}_s,$ and $\widetilde{\mathbf{C}}_s$ are matrices of dimension $\widetilde{\mathbf{A}}_s \in \mathbb{R}^{pQ \times R}, \ \ \widetilde{\mathbf{B}}_s \in\mathbb{R}^{pQ \times R}$ and $\widetilde{\mathbf{C}}_s \in \mathbb{R}^{pQ \times R}$ respectively and N is number of dimensions of tensor. For 3-mode tensor , N is 3 and for 4-mode tensor, N is 4 and so on. Similarly, for 4-mode, we have $(\mathbf{\widetilde{A}_s , \widetilde{B}_s , \widetilde{C}_s,\widetilde{D}_s})$ and so on.   Even though for 3-mode tensor, $\mathbf{A}$ and $\mathbf{B}$ do not increase their number of rows over time, the incoming slices may contribute valuable new estimates to the already estimated factors.Thus, we update the all factor matrices in same way. This is able to partially achieve Goal G2.

Update results:  Final step is to remove the all singleton dimensions from the sets compression matrices {$\mathbf{U_p,V_p, W_p}$} and stack them together. A singleton dimension of tensor or matrix is any dimension for which size of matrix or tensor with given dimension becomes one. Consider the 5-by-1-by-5 array A. After removing its singleton dimension, the array A become 5-by-5. Now, we can write compression matrices as:

\begin{equation}
\widetilde{\mathbf{P}}_a= \begin{bmatrix}
\mathbf{U(:,:,1)^T}\\
\mathbf{	U(:,:,2)^T}\\
\vdots\\
\mathbf{	U(:,:,p)^T}
\end{bmatrix}  , \  \ \widetilde{\mathbf{P}}_b = \begin{bmatrix}
\mathbf{V(:,:,1)^T}\\
\mathbf{V(:,:,2)^T}\\
\vdots\\
\mathbf{V(:,:,p)^T}
\end{bmatrix} , \  \ \widetilde{\mathbf{P}}_c = \begin{bmatrix}
\mathbf{W(:,:,1)^T}\\
\mathbf{W(:,:,2)^T}\\
\vdots\\
\mathbf{	W(:,:,p)^T}
\end{bmatrix}
\end{equation}

where $\widetilde{\mathbf{P}}_a,\widetilde{\mathbf{P}}_b,$ and $\widetilde{\mathbf{P}}_c$ are matrices of dimension $\widetilde{\mathbf{P}}_a=\mathbb{R}^{pQ \times I},\ \ \widetilde{\mathbf{P}}_b=\mathbb{R}^{pQ \times J}$ and $\widetilde{\mathbf{P}}_c=\mathbb{R}^{pQ \times K}$ respectively. The updated factor matrices (A, B, and C) for 3-mode tensor $\tensor{X}$ (i.e. $\tensor{X}_{old} +  \tensor{X}_{new}$)  can be obtained by :

\begin{equation}
\mathbf{A}=\widetilde{\mathbf{P}}_a^{-1}*\widetilde{\mathbf{A}}_s,\ \
\mathbf{B}= \widetilde{\mathbf{P}}_b^{-1}*\widetilde{\mathbf{B}}_s , \ \ \mathbf{C}= \big[ \mathbf{C}_{old};\widetilde{\mathbf{P}}_c^{-1}*\widetilde{\mathbf{C}}_s\big]
\end{equation}

where $\mathbf{A},\mathbf{B}$ and $\mathbf{C}$ are matrices of dimension $\mathbb{R}^{I \times R}, \ \ \mathbb{R}^{J \times R}$ and $\mathbb{R}^{K_{1\dots n,(n+1) \dots m} \times R}$ respectively. Hence, we achieve Goal G2.

Complexity Analysis: As discussed previously, compression step's time and space complexity is $O(QSt_{new})$ and $O(Q^3)$ respectively.  Identifiability and update can be calculated in $O(pQI+pQR)$. Hence, time complexity is considered as $O(p^2QI+p^2QIR+QSt_{new})$. Overall, as S is larger than any other factors, the time complexity of can be written as $O(QSt_{new})$. In terms of space consumption, is quite efficient since only the compressed matrices, previous factor matrices and summaries need to be stored. Hence, the total cost of space is $pQ(pT+t_{new}+R)+(T+t_{old})R+Q^3$.

Dataset     Time Complexity     Space Complexity     Reference

$O(QSt_{new})$   $pQ(pT+t_{new}+R)+(T+t_{old})R+Q^3$

OnlineCP   $O(NSt_{new}R)$   $St_{new}+(2T+t_{old})R+(N-1)R^2$   <cit.>

SambaTen   $O(nnz(X)R(N-1)$    $nnz(X)+(T+t_{old}+2S)R+2R^2$   <cit.>

RLST   $O(R^2S)$   $St_{new}+(T+t_{old}+2S)R+2R^2$   <cit.>

ParaComp   $O(QSt_{new}R)$   $St_{new}+(T+t_{old})R+pQ^3$   <cit.>

Complexity comparison between and state-of-art methods.

§.§ Extending to Higher-Order Tensors

We now show how our approach is extended to higher-order cases. Consider N-mode tensor $\tensor{X}_{old} \in \mathbb{R}^{I_1 \times I_2 \times \dots \times I_{N-1} \times T_1}$. The factor matrices are $(\mathbf{A}^{(1)}_{old}, \mathbf{A}^{(2)}_{old},\dots, \mathbf{A}^{(N-1)}_{old}, \mathbf{A}^{(T_1)}_{old})$ for CP decomposition with $N^{th}$ mode as new incoming data. A new tensor $\tensor{X}_{new} \in \mathbb{R}^{I_1 \times I_2 \times \dots \times I_{N-1} \times T_2}$ is added to $\tensor{X}_{old}$ to form new tensor of $\mathbb{R}^{I_1 \times I_2 \times \dots  \times I_{N-1} \times T}$ where $T=T_1+T_2$. In addition, sets of compression matrices for old data are {$\mathbf{U}_p^{(1)},\mathbf{U}_p^{(2)},\dots, \mathbf{U}_p^{(N-1)}, \mathbf{U}_p^{(T)}$} and for new data it is {$\mathbf{U}_p^{'(1)},\mathbf{U}_p^{'(2)},\dots, \mathbf{U}_p^{'(N-1)}, U_p^{'(T)}$} for p number of summaries.

Each compression matrices are converted into compressed cubes i.e. for $\tensor{X}_{old}$ compressed cube is of dimension $\tensor{Y_{p}} \in \mathbb{R}^{Q^{(1)} \times Q^{(2)}  \dots \times Q^{(N-1)}  \times Q^{(N)} }$ and same follows for $\tensor{X}_{new}$. The updated summaries are computed using $\tensor{X_{p}}=\tensor{Y_{p}}+\tensor{Z_{p}}$ s.t. $\tensor{X_{p}}  \in \mathbb{R}^{Q^{(1)} \times Q^{(2)}  \dots \times Q^{(N-1)}  \times Q^{(N)} }$. After CP decomposition of  $\tensor{X_{p}}$, factor matrices and random compressed matrices are stacked as :

\begin{multline*}
(\widetilde{\mathbf{A}}_{s}^{(1)},\widetilde{\mathbf{A}}_{s}^{(2)},\dots, \widetilde{\mathbf{A}}_{s}^{(N-1)} , \widetilde{\mathbf{A}}_{s}^{(N)}) \leftarrow \Pi\big[(\widetilde{A}_{s(i)}^{(1)}, \widetilde{\mathbf{A}}_{s(i)}^{(2)} ,\dots, \widetilde{\mathbf{A}}_{s(i)}^{(N-1)},\widetilde{\mathbf{A}}_{s(i)}^{(N)} ) \ \  {;} \\
\ \ (\widetilde{\mathbf{A}}_{s(i+1)}^{(1)}, \widetilde{\mathbf{A}}_{s(i+1)}^{(2)} ,\dots, \widetilde{\mathbf{A}}_{s(i+1)}^{(N-1)},\widetilde{\mathbf{A}}_{s(i+1)}^{(N)} )\big] , i \in (1,p-1)
\end{multline*}

\begin{equation}
\& \  \ (\widetilde{\mathbf{P}}^{(1)},\widetilde{\mathbf{P}}^{(2)},\dots, \widetilde{\mathbf{P}}^{(N-1)} , \widetilde{\mathbf{\mathbf{P}}}^{(N)}) = \begin{bmatrix}
\widetilde{\mathbf{P}}_{(1)}^{(1)},\widetilde{\mathbf{P}}_{(1)}^{(2)}, \dots, \widetilde{\mathbf{P}}_{(1)}^{(N-1)} , \widetilde{\mathbf{P}}_{(1)}^{(N)}\\
\widetilde{\mathbf{P}}_{(2)}^{(1)},\widetilde{\mathbf{P}}_{(2)}^{(2)}, \dots, \widetilde{\mathbf{P}}_{(2)}^{(N-1)} , \widetilde{\mathbf{P}}_{(2)}^{(N)}\\
\vdots\\
\widetilde{\mathbf{P}}_{(p)}^{(1)},\widetilde{\mathbf{P}}_{(p)}^{(2)}, \dots, \widetilde{\mathbf{P}}_{(p)}^{(N-1)} , \widetilde{\mathbf{P}}_{(p)}^{(N)}
\end{bmatrix}
\end{equation}

Finally, the update rule of each non-temporal mode $\in (1,N-1)$ and temporal mode $\in (N)$ is :

\begin{multline*}
(\mathbf{A}^{(1)}, \mathbf{A}^{(2)} ,\dots ,\mathbf{A}^{(N-1)}, \mathbf{A}^{N}) \leftarrow
\widetilde{\mathbf{P}}^{(2)^{-1}}*\widetilde{\mathbf{A}}_{s}^{(2)}, \dots ,
\end{multline*}

\begin{equation}
\widetilde{\mathbf{P}}^{(N-1)^{-1}}*\widetilde{\mathbf{A}}_{s}^{(N-1)}, \big[ \mathbf{A}_{old}^{(N)}; \widetilde{\mathbf{P}}^{(N)^{-1}}*\widetilde{\mathbf{A}}_{s}^{(N)}) \big]
\end{equation}

\begin{bmatrix}
\vdots\\
\end{bmatrix}
\leftarrow
\begin{bmatrix}
\widetilde{P}^{(1)-1}*\widetilde{A}_{s}^{(1)}\\
\widetilde{P}^{(2)-1}*\widetilde{A}_{s}^{(2)}\\
\vdots\\
\widetilde{P}^{(N-1)-1}*\widetilde{A}_{s}^{(N-1)}\\
\widetilde{P}^{(N)-1}*\widetilde{A}_{s}^{(N)}
\end{bmatrix}

Finally, by putting everything together, we obtain the general version of our for 3-mode tensor, as presented in Algorithm <ref>.The matlab implementation for higher modes is also available at $link^{1}$.

§.§ Necessary characteristics for uniqueness

As we mention above, is able to identify the solution of the online CP decomposition, as long as the parallel CP decompositions on the compressed tensors are also identifiable. Empirically, we observe that if the decomposition was given data that have exact or near-trilinear structure (or multilinear in the general case), i.e., obeying the low-rank CP model with some additive noise, is able to successfully, accurately, and using much fewer resources than state-of-the-art, track the online decomposition. On the other hand, when given data that do not have a low trilinear rank, the results were of lower quality. This observation is of potential interest in exploratory analysis, where we do not know 1) the (low) rank of the data, and 2) whether the data have low rank to begin with (we note that this is an extremely hard problem, out of the scope of this chapter, but we refer the interested reader to previous studies <cit.> for an overview). If provides a good approximation, this indirectly signifies that the data have appropriate trilinear structure, thus CP is a fitting tool for analysis. If, however, the results are poor, this may indicate that we need to reconsider the particular rank used, or even analyzing the data using CP in the first place. We reserve further investigations of what this observation implies for future work.

§ EXPERIMENTS

We design experiments to answer the following questions:

Q1. Memory Efficient: How  much memory required for updating incoming data?

Q2. Fast and Accurate: How fast and accurately are updates in compared to incremental algorithms?

Q3. Scalability: How does the running time of increase as tensor data grow (in time mode)?

Q4. Sensitivity: What is the influence of parameters on ?

Q5. Effectiveness: How used in real-world scenarios?

For our all experiments, we used Intel(R) Xeon(R), CPU E5-2680 v3 @ 2.50GHz machine with 48 CPU cores and 378GB RAM.
In this section we extensively evaluate the performance of on multiple synthetic datasets, and compare its performance with state-of-the-art approaches.We experiment on the different parameters of , and how that affects performance

§.§ Evaluation Measures

We evaluate and the baselines using three criteria: fitness, processor memory used, and wall-clock time. These measures provide a quantitative way to compare the performance of our method. More specifically, Fitness measures the effectiveness of approximated tensor and defined as :

\begin{equation}
Fitness(\%)= 100* \Big(1- \frac{||\tensor{X}-\tensor{X}_{approx}||_F^2}{||\tensor{X} ||_F^2}\Big)
\end{equation}

higher the value, the better approximation.

CPU time (sec): indicates the average running time for processing for processing all slices for given tensor, measured in seconds, is used to validate the time efficiency of an algorithm.

Process Memory used(MB): indicates the average memory required to process each slices for given tensor, is used to validate the space efficiency of an algorithm.

§.§ Baselines

In this experiment, four baselines have been selected as the competitors to evaluate the performance.

OnlineCP<cit.>: it is online CP decomposition method, where the latent factors are updated when there are new data.

SambaTen<cit.>: Sampling-based Batch Incremental Tensor Decomposition algorithm is the most recent and state-of-the-art method in online computation of canonical parafac and perform all computations in the reduced summary space.

RLST<cit.>: Recursive Least Squares Tracking (RLST) is another online approach in which recursive updates are computed to minimize the Mean Squared Error (MSE) on incoming slice.

ParaComp<cit.>: an implementation of non-incremental parallel compression based tensor decomposition method. The model is based on parallel processing of randomly compressed and reduced size replicas of the data. Here, we simply re-compute decomposition after every update.

§.§ Experimental Setup

The specifications of each synthetic dataset are given in Table <ref>. We generate random tensors of dimension I = J = K with increasing I and other modes. Those tensors are created from a known set of randomly generated factors with known rank $R$, so that we have full control over the ground truth of the full decomposition. We dynamically calculate the size of incoming batch or incoming slice(s) for our all experiments to fit the data into 10% of memory of machine. Left over memory of machine  is used for computations for algorithm. We use 20 parallel workers for every experiment.

Dimension     $\#$number of non-zeros      Batch size     p   Q   shared

50 x 50 x 50    $125K$   5    20   30   5

100 x 100 x 100    $1M$   10    30   35   10

500 x 500 x 500    $125M$   50   40   30   6

1000 x 1000 x 1000     $1B$   20   50    40    10

5000 x 5000 x 5000    $7B$   10   90   70   25

10000 x 10000 x 10000    $1T$   10   110   100   20

50000 x 50000 x 50000    $6.25T$   4   140   150   30

Table of Datasets analyzed

Note that all comparisons were carried out over 10 iterations each, and each number reported is an average with a standard deviation attached to it. Here, we only care about the relative comparison among baseline algorithms and  it is not mandatory to have the best rank decomposition for every dataset. Although, there are few methods <cit.> available in literature to find rank of tensor. But for simplicity, we choose the rank R to 5 for all datasets. In case of  method-specific parameters, for the ParaComp algorithm, the settings are choose to give best performance. For OnlineCP, we choose the batchsize which gave best performance in terms of approximation fitness. For fairness, we also compare against the parameter configuration for that yielded the best performance for baselines. Also, during processing , for all methods we removed the unwanted variable from baselines to fairly compare with our methods.

§.§ Results

§.§.§ Memory efficient, Fast and Accurate

For all datasets we compute Fitness(%),CPU time (seconds) and Memory(MB) space required. For , OnlineCP, ParaComp,Sambaten and RLST we use 10% of the time-stamp data in each dataset as existing old tensor data. The results for qualitative measure for data is shown in Table <ref>-<ref> . For each of tensor data ,the best performance is shown in bold. All state-of-art methods address the issue very well. Compared with OnlineCP, ParaComp,Sambaten and RLST, give comparable fitness and reduce the mean CPU running time by up to 2x times for big tensor data. For all datasets, PARACOMP's accuracy (fitness) is better than all methods. But it is able to handle upto $\mathbb{R}^{10^4 \times 10^4 \times 250}$ size tensor only. For small size datasets, OnlineCP's efficiency is better than all methods. For large size dataset, outperforms the baseline methods w.r.t fitness, average running time (improved 2x-4x) and memory required to process the updates. It significantly saved 200-250% of memory as compared to Online CP, Sambaten and RLST as shown in Table <ref>. It saved 40-80% memory space compared to Paracomp. Hence, is comparable to state-of-art methods for small dataset and outperformed them for large dataset. These results answer Q1 & Q2 as the have comparable qualitative measures to other methods.

Size    5|c|Fitness(%)

I=J=K          OnlineCP    SambaTen   ParaComp   RLST

50    97.1$\pm$4.9   86.9$\pm$7.3   76.9$\pm$16.1   100.0$\pm$0    95.0$\pm$0.1

100    95.8$\pm$2.2   	93.4$\pm$6.9    72.1$\pm$3.7   	100.0$\pm$0 	   97.7$\pm$0.1

500   96.9$\pm$1.3   90.9$\pm$8.7   	84.1$\pm$0.6   100.0$\pm$0    	94.6$\pm$0.1

1000   96.1$\pm$0.1   71.6$\pm$10.7   	85.6$\pm$0.1	   100.0$\pm$0    	98.3$\pm$0.1

5000   90.8$\pm$14.7	   90.6 $\pm$42.7   	80.5$\pm$15.5   	100.0$\pm$0 	   98.8$\pm$0.1

10000   98.13$\pm$2.5   54.34$\pm$1.4   56.12$\pm$1.5   NA   97.3$\pm$0.2

50000   68.13$\pm$4.2   56.89$\pm$3.7   53.95$\pm$12.7   NA   NA

Dimension    5|c|CPU Time(sec)

I=J=K          OnlineCP    SambaTen   ParaComp   RLST

50    1.32$\pm$0.01   	0.95$\pm$0.1   2.53$\pm$0.1   1.09$\pm$0.1   1.35$\pm$0.1

100     1.92$\pm$0.3   1.35$\pm$0.3   5.19$\pm$0.3   	1.91$\pm$0.3   	4.97$\pm$0.3

500	   19.7$\pm$0.1   	20.73$\pm$4.4	   99.21$\pm$2.9	   26.3$\pm$0.1	   37.20$\pm$0.2

1000   	290.9$\pm$6.7	   455.7$\pm$24.5   	606.3$\pm$1.7   	511.37$\pm$1.6   	401.08$\pm$1.8

5000   4398.8$\pm$45.5   	5835.5$\pm$0.2K   	6779.7$\pm$0.3K   	5157.98$\pm$3.9   	6464.52$\pm$0.3K

10000   15406.5$\pm$156.9   21287.3$\pm$85.6   20589.81$\pm$69.4   NA   22974.76$\pm$0.1K

50000   78892.8$\pm$98.2   80159.6$\pm$23.5   79922.80$\pm$47.3   NA   NA

Experimental results for speed and accuracy of approximation of incoming slices. We see that gives comparable accuracy and speed to baseline.This answers our Q2.

Size     5|c|Process Memory (MB)

I=J=K         OnlineCP    SambaTen    ParaComp   RLST

50     	1.4$\pm$0.001   	1.6$\pm$0.002   	8.8$\pm$0.001   	5.5$\pm$0.003   	6.271$\pm$0.001

100   	13.2$\pm$0.04	   16.8$\pm$0.03	   30.1	$\pm$0.01   14.5$\pm$0.02	   19.255		$\pm$0.02

500   	 25.7$\pm$0.09   	2018.3$\pm$0.91   	2959.4$\pm$0.15   	26.3$\pm$0.01   	1070.138		$\pm$0.01

1000   	85.1$\pm$4.1   	16037.4$\pm$56.5   	13830.8$\pm$21.7   	110.9$\pm$5.1   	7789.841	$\pm$23.2