Multi-Aspect Incremental Tensor Decomposition Based on Distributed In-Memory Big Data Systems

Before describing PARAFAC decomposition, we will present the relevant notations. The Kronecker product is denoted as ⊗. Given two matrices A and B, where the size of A is I × K and the size of B is J × R, the Kronecker product generates a matrix of size IJ × KR. The Khatri–Rao product is the column-wise Kronecker product of two matrices with the same number of columns. Given two matrices A and B, where the size of A is I × R and that of B is J × R, the Khatri–Rao product A ⊙ B generates a matrix of size IJ × R. The Hadamard product is the element-wise product, denoted as A * B, where the two matrices have the same size.

We will now describe tensor matricization, which unfolds a tensor into matrix format. Figure 1 shows a representative example, where a three-dimensional tensor X of size I × J × K is unfolded into three matrices. The unfolded matrices along the I, J, and K axes are denoted by X₍₁₎, X₍₂₎, and X₍₃₎, respectively. The sizes of the three unfolded matrices are X₍₁₎ ∈ R^{I × JK}, X₍₂₎ ∈ R^{J × IK}, and X₍₃₎ ∈ R^{K × IJ}.

Figure 1

However, when a new tensor is added (i.e., streaming or dynamic tensors), the total size of the tensor changes and the results of tensor matricization also change. Figure 2 shows the effect on matricization of the addition of a new tensor.

Figure 2

As tensors generally grow over time, in Figure 2, it is assumed that the K axis (time axis) is increasing and the other axes can also increase. The gray regions represent new values corresponding to the arrival of a new tensor. When the entire tensor is unfolded, the new tensor is located behind or under the existing tensor. Therefore, there is a need for a method in which only newly added tensors are decomposed, and the results of existing tensor decompositions are updated. The proposed method decomposes new values without decomposing the entire tensor; thereby, existing tensors are not unfolded, and new unfolded tensors are added to the existing matrix. The proposed method uses tensor indexing to concatenate new tensors and existing results.

PARAFAC decomposition, also known as CP decomposition, is one of the most widely used tensor decomposition algorithms. It decomposes a tensor into a sum of rank-one tensors (Gujral et al., 2018; Jeon et al., 2015; Kang et al., 2012; Yang & Yong, 2019a). A given three-dimensional tensor X is decomposed into R components to yield three factor matrices. This can be expressed as follows: (1)X≈∑r=1Rλrar ∘br ∘crX \approx \sum\limits_{r = 1}^R {{\lambda _r}{a_r}{\ ^ \circ }{b_r}{\ ^ \circ }{c_r}}

The alternating least squares (ALS) method is used for computing the factor matrices in PARAFAC decomposition. In the case of a three-dimensional tensor, this method fixes two factor matrices to obtain the third factor matrix (Yang & Yong, 2019a; 2019b). The PARAFAC-ALS method for computing the factor matrices is expressed as follows: (2)A=X(1)(C⊙B)(C⊤C*B⊤B)†B=X(2)(C⊙A)(C⊤C*A⊤A)†C=X(3)(B⊙A)(B⊤B*A⊤A)†\matrix{ {A = {X_{(1)}}(C \odot B){{({C^ \top}C*{B^ \top}B)}^\dagger }} \hfill \cr {B = {X_{(2)}}(C \odot A){{({C^ \top}C*{A^ \top}A)}^\dagger }} \hfill \cr {C = {X_{(3)}}(B \odot A){{({B^ \top}B*{A^ \top}A)}^\dagger }}}

In Equation (2), X₍₁₎(C ⊙ B), X₍₂₎(C ⊙ A), and X₍₃₎(B ⊙ A) (i.e., the Khatri–Rao products) lead to intermediate data explosions. To alleviate this, several methods have been proposed for avoiding the calculation of Khatri–Rao products. For example, Kang et al. (2012), Jeon et al. (2015), and Park et al. (2016) calculated only non-zero values and avoided direct Khatri–Rao product calculations by multiplying the column vectors of the factor matrices. These methods frequently use map-and-reduce processes for each operation, and thus excessive shuffling may result in slow execution. Papalexakis, Faloutsos, and Sidiropoulos (2012) used random sampling for tensor decomposition. Jeon et al. (2016) proposed SCouT, which is a large and scalable coupled matrix–tensor factorization algorithm. However, all these methods are based on Hadoop. As Hadoop-based tensor decomposition methods are not suitable for iterative operations, several methods based on Spark have been proposed. For example, the algorithms in Yang and Yong (2018, 2019a) are both based on Apache Spark and reduce the required number of map-and-reduce operations. In the method in Yang and Yong (2018), memory overflow occurs when big data are processed in small distributed systems. This issue is resolved in Yang and Yong (2019a), where tensors are divided into sub-tensors, which are subsequently decomposed. In addition to tensor decomposition, methods for tensor completion are being extensively investigated. Ge et al. (2018) proposed a method called DisTenC; it is a distributed large-scale tensor completion algorithm based on Apache Spark. However, the main goals of tensor decomposition significantly differ from those of tensor completion. Tensor completion is based on low-rank approximations to obtain an accurate estimate of missing values in tensors, whereas tensor decomposition focuses on latent patterns and factors. Therefore, our study focuses on tensor decomposition rather than tensor completion.

Incremental tensor decomposition algorithms are used to decompose dynamic tensors. Recalculation time can be reduced by such algorithms. Most incremental tensor decomposition algorithms assume that only one axis increases over time. The most representative example is online CP (Zhou et al., 2018), which is an incremental CP decomposition algorithm in which it is assumed that when new tensors are added to an existing tensor, only the time axis increases. Ma, Yang, and Wang (2018) proposed the randomized online CP decomposition algorithm, which is based on random sampling. This algorithm avoids calculating Khatri-Rao products and reduces memory usage. However, the algorithms mentioned above are based on MATLAB routines running on a single machine. Recently, distributed incremental tensor decomposition methods have been introduced. Gujral et al. (2018) proposed a sampling-based batch incremental CP tensor decomposition algorithm, which has been implemented in both MATLAB (Gujral et al., 2018) and Apache Spark (this can be download from https://github.com/lucasjliu/SamBaTen-Spark). The Spark version is a re-implementation of the MATLAB version. Yang and Yong (2019b) proposed InParTen, which is also an incremental PARAFAC decomposition algorithm based on Apache spark. Yang and Yong (2019b) and Gujral et al. (2018) both assume that a tensor grows only along one axis. However, real-world tensors can grow along all axes.

Recently, a few incremental tensor decomposition methods that consider multiaxial growth have been developed. Song et al. (2017) proposed a multi-aspect dynamic tensor completion framework called MAST which is based on CP decomposition. This method decomposes added tensors by partitioning them into sub-tensors. This affects the results of adjacent parts, implying that sub-tensor decomposition must be performed sequentially. Najafi, He, and Yu (2019) proposed an outlier-robust multi-aspect streaming tensor completion and factorization algorithm called OR-MSTC (Najafi et al., 2019). The incremental tensor decomposition method in OR-MSTC generates sub-tensors from incoming tensors. Nimishakavi et al. (2018) proposed a side-information-infused incremental tensor analysis method called SIITA which incorporates side information and can process general incremental tensors by Tucker decomposition. However, this method cannot handle scalable tensor datasets because it does not consider scalability issues (Nimishakavi et al., 2018). Additionally, most multi-aspect tensor decomposition algorithms either are based on MATLAB or focus on tensor completion.

In this study, we developed an incremental tensor decomposition algorithm called InParTen2. It considers both dynamic and static tensor datasets. Furthermore, the proposed method is implemented in Apache Spark and can therefore handle distributed in-memory big-data systems without dividing new tensors into subtensors. It can process new tensors based on their indexes and the results of previous tensor decompositions. Additionally, InParTen2 avoids the calculation of Khari–Rao products and minimizes shuffling by using the Apache Spark platform. Table 2 provides a summary of various incremental tensor decomposition algorithms.

InParTen2 decomposes new tensors without generating sub-tensors. We first perform initialization prior to executing incremental PARAFAC-ALS. In the decomposition of the added tensors, there are certain parts that affect the existing factor matrices. Therefore, we add the existing factor matrices when A and B are updated. For example, to update A, we add the existing factor matrix A₁ and the new factor matrix A. This process can be expressed as follows: (3.1)P=Xold(1)(C1⊙B1)Pnew=Xnew(1)(C2⊙B1)P=P+Pnew\eqalign{ & P = {X_{old(1)}}({C_1} \odot {B_1}) \cr & {P_{new}} = {X_{new(1)}}({C_2} \odot {B_1}) \cr & P = P + {P_{new}}}(3.2)Q=C1⊤C1*B1⊤B1Qnew=C2⊤C1*B1⊤B1\matrix{ {Q = C_1^ \top{C_1}*B_1^ \top{B_1}} \hfill \cr {{Q_{new}} = C_2^ \top{C_1}*B_1^ \top{B_1}}}(3.3)Q=Q+QnewA=PQ†\eqalign{& { Q = Q + {Q_{new}}} \cr & {A = P{Q^\dagger }}}

The initial values of P and Q are set to X_old(1)(C₁ ⊙ B₁) and C1⊤C1*B1⊤B1C_1^ \top{C_1}*B_1^ \top{B_1} , respectively. Here, we can calculate the Khatri–Rao products by using the existing tensor to obtain an initial value of P. We can avoid repeated calculations by using the existing factor matrices A₁, B₁, and C₁ because they were obtained from the existing tensor. This method is expressed by the following equation: (4)Xold(1)(C1⊙B1)=A1(C1⊙B1)⊤(C1⊙B1)=A1(C1⊤C1*B1⊤B1)=A1Q{X_{old(1)}}\left( {{C_1} \odot {B_1}} \right) = {A_1}{\left( {{C_1} \odot {B_1}} \right)^ \top}\left( {{C_1} \odot {B_1}} \right) = {A_1}\left( {C_1^ \top{C_1}*B_1^ \top{B_1}} \right) = {A_1}Q

To calculate X_old(1)(C₁ ⊙ B₁), we can solve the equation above by multiplying Q by the existing matrix A₁. Equation (4) can effectively reduce execution time and improve memory efficiency because the existing tensor is not unfolded, and intermediate data explosion is avoided by using the Hadamard product. Time complexity is thus reduced from O(RIJK) to O(R²(I + J + K)). Additional values, such as X_old(2)(C₁ ⊙ A₁) and X_old(3)(B₁ ⊙ C₁), can also be obtained by a similar method. Subsequently, we initialize the factor matrices C_updated and A based on the existing factor matrices and unfold the new tensors prior to decomposing them. In the initialization process, two cases are considered. In the first case, all axes increase over time, and in the second, the time axis does not change or exhibits relatively small changes as the other axes increase. We will now explain the initialization process for the first case. In Figure 4 above, all dimensions increase over time. It is evident that the factor matrices obtained from each sub-tensor have overlapping parts: C₂, A₁, A₂, and B₁. These parts affect each other. In particular, the factor matrix C₂ affects all decomposed sub-tensors in the new tensors. Therefore, the initial value of C₂ can be obtained from the unfolded tensor X_new(3) and the existing factor matrices A₁ and B₁. In this case, the size of K’ for the added tensor is greater than zero; thus, we initialize C_updated to the initial value of C₂, as shown in Equation (5.1). In the second case, where the time axis (K axis) does not change or changes little, we initialize C_updated using C₁, as shown in Equation (5.2). (5.1)(5.2)Cupdated={C2=Xnew(3)(B1⊙A1)(B1⊤B1*A1⊤A1)†if K'>0C1if K'≤0{C_{updated}} = \left\{ {\matrix{ {{C_2} = {X_{new(3)}}\left( {{B_1} \odot {A_1}} \right){{\left( {B_1^ \top{B_1}*A_1^ \top{A_1}} \right)}^\dagger }} \hfill & {if\,K' > 0} \hfill \cr {{C_1}} \hfill & {if\,K' \le 0} \hfill \cr } } \right.

Subsequently, we solve for A using the unfolded tensors X_new(1) and C_updated combined with the existing factor matrices B₁ and A₁ as follows: (6)A=A1+Xnew(1)(Cupdated⊙B1)(Cupdated⊤Cupdated*B1⊤B1)†=(P+Pnew)(Q+Qnew)†A = {A_1} + {X_{new(1)}}\left( {{C_{updated}} \odot {B_1}} \right){\left( {C_{updated}^ \top{C_{updated}}*B_1^ \top{B_1}} \right)^\dagger } = \left( {P + {P_{new}}} \right){\left( {Q + {Q_{new}}} \right)^\dagger }

As mentioned previously, we update the factor matrices A based on Equation (3). However, when the size of the I axis of Xnew(1) increases by I’, the size of X_new(1)(C_updated ⊙ B₁) is (I + I’) × R. Therefore, X_new(1)(C_updated ⊙ B₁) overlaps with the initial P, and the value of the initial P should be added to that of P_new. To this end, we use indexing to add the parts that should be updated in the existing P; subsequently, the part of size I’ × R is appended to the bottom rows of the initial P because P and X_new(1)(C_updated ⊙ B₁) have different size. Additionally, Cupdated⊤Cupdated*B1⊤B1C_{updated}^ \top{C_{updated}}*B_1^ \top{B_1} and the existing Q have the same size, and therefore simple matrix addition can be performed.

Following this initial process, we use the incremental PARAFAC-ALS process to update the factor matrices B, A, and C. This update method can be expressed as follows: (7)B=B1+Xnew(2)(Cupdated⊙A)(Cupdated⊤Cupdated*A⊤A)†A=A1+Xnew(1)(Cupdated⊙B)(Cupdated⊤Cupdated*B⊤B)†Cupdate=Xnew(3)(B⊙A)(B⊤B*A⊤A)†\matrix{ {B = {B_1} + {X_{new(2)}}\left( {{C_{updated}} \odot A} \right){{\left( {C_{updated}^ \top{C_{updated}}*{A^ \top}A} \right)}^\dagger }} \cr {A = {A_1} + {X_{new(1)}}\left( {{C_{updated}} \odot B} \right){{\left( {C_{updated}^ \top{C_{updated}}*{B^ \top}B} \right)}^\dagger }} \cr {{C_{update}} = {X_{new(3)}}\left( {B \odot A} \right){{\left( {{B^ \top}B*{A^ \top}A} \right)}^\dagger }}}

First, to update the factor matrix B, we use the factor matrices C_updated and A from the new tensors. Similarly, to update A, we use indexing to add only the parts that should be updated in the factor matrix B₁. At the end of the iteration, we finally update the factor matrix C when the value of the initialized C_updated is C₂. In this case, C_updated can be sequentially added to the rows of the existing matrix C₁, as shown in the following equation: (8)C=[C1Cupdated]=[C1Xnew(3)(B⊙A)(B⊤B*A⊤A)†]C = \left[ {\matrix{ {{C_1}} \cr {{C_{updated}}} \cr } } \right] = \left[ {\matrix{ {{C_1}} \cr {{X_{new(3)}}\left( {B \odot A} \right){{\left( {{B^ \top}B*{A^ \top}A} \right)}^\dagger }} \cr } } \right]

In this section, we describe the implementation of InParTen2 in Apache Spark to achieve fast and memory-efficient execution. First, we use the coordinate format to load a new tensor as an input because real-word tensors are generally sparse. The coordinate format is a resilient distributed data (RDD) type consisting of three indexes and values. Subsequently, this format is used to load existing three-factor matrices. These are constructed from row indexes and row vectors. Owing to the limited memory space, only non-zero values of the tensor and three-factor matrices are loaded. Prior to the execution of PARAFAC-ALS, we initialize the factor matrices C_updated and A using the existing tensor. This can be carried out based on Equations (5.1) and (6). However, as mentioned previously, computing Khatri–Rao products can lead to intermediate data explosions. We modified the HaTen2-DRI method for implementation in Apache Spark. In HaTen2-DRI, intermediate data explosions occur 2nnz(x)R times. However, InParTen2 can reduce this to 2nnz(x). Khatri–Rao products can be avoided by using Equation (9). Thereby, all steps involving Khatri–Rao product operations can be omitted. (9)X(1)(C⊙B)=∑(j,k)=(1,1)JKX(1)(i,j,k)B(j,::)C(k,::){X_{(1)}}\left( {C \odot B} \right) = \sum\limits_{(j,k) = (1,1)}^{JK} {{X_{(1)(i,j,k)}}B(j,::)C(k,::)}

Here, the non-zero tensor values are multiplied by the two-factor matrices corresponding to the non-zero indexes in the tensor. For example, if there is a non-zero value in X_(1)(2,1,2), then all column values in row 1 of B and row 2 of C are multiplied. In the parts where tensor decomposition involving Khatri–Rao products was initially required, we use Equation (7) instead. We use the RDD functions join, map, and reduceByKey to implement Equation (7). To update A and B, the existing factor matrix should be added. As mentioned previously, we do not generate a zero matrix; rather, we use indexing to add to the values of the existing factor matrix. If the existing factor matrix does not require updating, the addition operation is not performed. These processes use the join and map functions.

Overall, the time and space complexity of InParTen2 is O((6nnz(X_new)+R²((I+I’)+(J+J’)+K’)) and O(R(I+J+K’)+2nnz(X_new)), respectively. Where nnz(X_new) is the number of non-zero in new tensor, and R is the decomposition rank.

We used synthetic and real-world datasets to evaluate InParTen2. We used 10 types of synthetic datasets for each tensor size. In the synthetic datasets, we generated random tensors of dimension I = J = K, where I ∈ [100, 500, 800], and used different initialized densities: 20% and 60% in each initial tensor. Further, we generated incoming tensors with densities 5%, 10%, and 20% in the initial tensors, which are half of that in the total tensor. In addition, we used different initial size and density.

We used three real-world datasets containing sparse three-dimensional tensors: YELP (Jeon et al., 2016), MovieLens (Jeon et al., 2016), and Netflix (Bennett & Lanning, 2007). Table 3 lists these datasets. YELP contains users, locations, time indexes, and rating values from 1 to 5. MovieLens and Netflix contain users, movies, time indexes, and rating values from 1 to 5. In the incremental tensor decomposition by InParTen2, it is assumed that the three dimensions increase, and the initial dimension is fixed. In the SamBaTen (Gujral et al., 2018) algorithm, which considers the increase in a single dimension, only the time axis of the initial dimension is set to be equal to that of InParTen2.

We compared the execution time and accuracy of InParTen2 with those of existing tensor decomposition methods based on Apache Spark. Execution time is an important measure because one of our goals is to reduce time cost. We compared the execution time using synthetic and real-world datasets. In the former case, we performed various measurements. In the first measurement, the initial tensors were set to a density of 20% and 60%, and the incoming tensors had densities of 5%, 10%, and 20%. This experiment was conducted to study the effect of the density of the incoming data on performance. In the second measurement, the size of the initial tensor and incoming tensors were varying. Here, the density of the incoming tensor was set to 20%, and that of the initial tensor to 20% and 60%. In this experiment, the incoming data had the same density, but performance was evaluated in terms of the size of the incoming tensor data. In the real-world datasets, we set the initial size. We tested 10 iterations and 10 ranks in each dataset.

The second comparison is based on accuracy. We reconstructed tensors through decomposition and used the relative error to investigate their similarity to the original tensors. We compared these results to those obtained by existing tensor decomposition methods, particularly static tensor decomposition. Relative error can be defined as follows: (10)Relative Error=‖X^reconstruct−Xoriginal‖2‖Xorginal‖2{Relative\,Error} = {{{{\left\| {{{\hat X}_{reconstruct}} - {X_{original}}} \right\|}_2}} \over {{{\left\| {{X_{orginal}}} \right\|}_2}}} where X^reconstruct{\hat X_{reconstruct}} denotes the reconstructed tensor using the decomposition results. A small relative error implies better approximation.

We used Spark-based tensor decomposition methods. There are several Hadoop-based distributed tensor decomposition algorithms, such as HaTen2 (Jeon et al., 2015) and BigTensor (Park et al., 2016), as well as MATLAB-based decomposition tools. However, Hadoop-based decomposition tools require longer iteration time, and therefore their execution time cannot be compared with that of Spark-based algorithms. MATLAB-based algorithms run on a single machine, thus limiting their ability to handle large data. Therefore, we used IM-PARAFAC and SamBaTen (Gujral et al., 2018), which run in the Apache Spark framework, for the comparison. IM-PARAFAC (Yang & Yong, 2019a) can decompose static tensors and handle large datasets. SamBaTen can decompose both static and incremental tensors. The incremental tensor decomposition of SamBaTen considers the increase in the time-axis only; accordingly, the initial size of the tensor datasets was set only for the time axis. The incremental tensor decomposition of SamBaTen is based on sampling. However, in the case of the static CP-ALS provided by SamBaTen, the HaTen2-Naive algorithm was re-implemented instead of the sampling algorithm. As InParTen2 also performs both static and incremental decomposition, both methods were tested. In the tests, the static tensor decomposition of InParTen2 is denoted as InParTen2-Static.

Initially, we considered the case in which initial tensors were set to a density of 20% and 60%, and the incoming tensors to 5%, 10%, and 20%. This test demonstrated the effect of density on computational cost for fixed size. For this experiment, the size of the initial tensor was set to half of the total tensor size. Table 4 and Figure 5 show the comparison of execution time for different densities of the incoming tensors. The best results are highlighted. Overall, it was demonstrated that sparse datasets exhibit shorter execution time. Further, InParTen2 was faster in every respect. InParTen2 expediently handled both static and incremental tensor decomposition. The latter was at least twice as fast as that by SamBaTen.

Figure 5

Table 5 shows the comparison of the relative errors for different densities of the incoming tensor datasets. InParTen2 exhibits similar relative errors to those of the other methods. In general, datasets with an initial density of 60% exhibit better relative error than those with an initial density of 20%. If the original datasets are sparse, the reconstructed tensor is filled with values, thus increasing the relative error.

Subsequently, we compared the execution time and relative error when data size increased. We considered initial tensors of different sizes under the same tensor size. Here, the density of the initial tensor was set to 20% and 60%, and that of the incoming tensors to 20%. Table 6 and Figure 6 show the comparison of the execution time for different initial size. The best results are highlighted. When the initial density was 60%, the execution time of the static decomposition tools increased gradually, but the dataset with an initial density of 20% exhibited the best performance when the initial size was half. The incremental decomposition by SamBaTen exhibited increased execution time. When the initial density was 60%, its execution time increased with the initial dimension, but when the initial density was 20%, the execution time decreased as the initial dimension increased. However, the incremental tensor decomposition by InParTen2 was faster regardless of the initial density when the initial dimension was larger and the amount of added data was small.

Figure 6

Table 7 shows the comparison of the relative errors for different initial dimensions. It can be seen that the relative errors by InParTen2 and the other tensor decomposition tools are similar. Overall, the relative error increases with density. Datasets with an initial density of 20% were overall sparse, and therefore they had similar accuracy regardless of the initial size.

Finally, we compared the execution time and relative error using real-world datasets. Table 8 shows the results. It can be seen that InParTen2 is overall faster. InParTen2-static is approximately 1.03 to 3.6 times as fast as IM-PARAFAC, and 17 times as fast as CP-ALS of SamBaTen. InParTen2 exhibits smaller execution time than static tensor decomposition. In addition, it is approximately 1.5 to 18 times as fast as SamBaTen, and unlike SamBaTen, it can handle large tensors. Furthermore, InParTen2 exhibits similar relative errors to those of the other methods. SamBaTen-CPALS and SamBaTen exhibit a large relative error (larger than that of InParTen2) because they sample the additional data. However, InParTen2 has a difference of 0.04 to 0.05 from InParTen2-static and IM-PARAFAC in the MovieLens datasets. This apparently increases the relative error because the number of tensor values is increased to an excessive degree.