Incremental Nonnegative Tucker Decomposition with Block-Coordinate Descent and Recursive Approaches

Incremental versions of various tensor decomposition models have been extensively studied in the last decade, including online and incremental CP [52,53], online TT [54], dynamic approximation to HT [55], etc. There are many approaches and computational strategies for the Tucker decomposition.

The first versions of streaming Tucker decomposition were proposed by Sun et al. [56,57] in 2006. They were addressed for window-based tensor analysis [57] and dynamic tensor analysis [56]. Both can be regarded as extensions of the baseline Tucker decomposition model to incremental processing. However, orthogonal factor matrices in these models are computed with the standard singular value decomposition (SVD) applied to additively updated covariance matrices, which is their bottleneck. The above models have been next generalized and deeply studied in [58]. The Tucker decomposition with an incremental SVD was discussed in [59]. Gujral et al. [60] proposed the sampling-based batch incremental tensor decomposition (SamBaTen) that extends the CP model to an incremental version.

Another computational approach for processing large-scale streaming tensors with the Tucker decomposition was proposed by Malik and Becker [61]. It is the TensorSketch algorithm that is based on sketching Kronecker products in the Tucker decomposition using randomization techniques. This method is also addressed for updating factor matrices with orthogonality constraints. The Tucker decomposition was also extended for processing large-scale sparse tensors by Oh et al. [11] who proposed the P-TUCKER algorithm. It performs alternating least squares (ALS) updates with a row-wise rule, which is a memory-saving strategy for updating factor matrices and can be easily parallelized for multi-core processing. The updated factors with ALS are then orthogonalized using the alternating QR factorization. Another version of a scalable and incremental version of the Tucker decomposition was proposed by Xiao et al. [49]. In their approach, the input tensor may be dynamic and long in each mode. It is partitioned into timestamp-based small-scale subtensors along each mode, and then the factor matrices are updated with the ALS-based rule considering the factor matrices from previous timestamps and the auxiliary matrices from the current timestamp. The modified Gram–Schmidt orthogonalization is then used to impose the orthogonality constraints on the estimated factors. Other versions of the online Tucker decomposition were also studied in the unpublished work [62]. There are also statistical approaches to the online Tucker decomposition, e.g., Bayesian streaming sparse Tucker decomposition [63] and online stochastic decompositions [64]. Recently, Chachlakis et al. [65] proposed the dynamic $L_1$-Tucker decomposition algorithm for dynamic and outlier-resistant Tucker analysis of multi-way streaming data. Their approach is also not addressed for NTD.

To the best of our knowledge, there is no incremental version of NTD. In this study, we propose a new computational approach to NTD that is inspired by the online version of NMF proposed in [48].

Notations: tensors, matrices, vectors, and scalars are denoted by calligraphic uppercase letters (e.g., $\mathcal{Y}$), boldface uppercase letters (e.g., Y), lowercase boldface letters (e.g., y), and non-bold letters (e.g., y), respectively. For a matrix Y, ${\mathit{Y}}_{i,j}$ denotes the $(i,j)$-th element, and ${\mathit{y}}_j$ or ${\underline{\mathit{y}}}_j$ stand for the j-th column or row, respectively. The set of non-negative real numbers is denoted by ${\mathbb{R}}_{+}$. The subtensor obtained from $\mathcal{Y}$ by selecting its t-th batch subtensor across the N-th mode is denoted ${\mathcal{Y}}^{\left(t\right)}\in {\mathbb{R}}_{+}^{I_1\times \dots \times I_{N-1}\times I_N^{\left(t\right)}}$.

Let $\mathcal{Y}=\left[{\mathit{Y}}_{i_1,i_2,\dots ,i_N}\right]\in {\mathbb{R}}^{I_1\times I_2\times \dots \times I_N}$ be the N-way array, which will be referred to as the N-modal tensor.

Unfolding, also known as matricization, is the way of getting a matrix from a multi-way array by reshaping the data. The mode-n matricization of $\mathcal{Y}$, denoted by ${\mathit{Y}}_{\left(n\right)}\in {\mathbb{R}}^{I_n\times {\prod }_{\mathit{P}\ne n}{I}_{\mathit{P}}}$, rearranges mode-n fibers placing them as the columns of matrix ${\mathit{Y}}_{\left(n\right)}$. Tensor entries $(i_1,i_2,\dots ,i_N)$ are mapped to matrix’s entries $(i_n,j)$, where

The standard model of the Tucker decomposition of $\mathcal{Y}$ has the form: where $\mathcal{G}=\left[{\mathit{G}}_{j_1,j_2,\dots ,j_N}\right]\in {\mathbb{R}}^{J_1\times J_2\times \dots \times J_N}$, for $J_n\le I_n$ and $n=1,\dots ,N$, is the core tensor, and ${\mathit{U}}^{\left(n\right)}=[{\mathit{U}}_1^{\left(n\right)},\dots ,{\mathit{U}}_{J_n}^{\left(n\right)}]=\left[{\mathit{U}}_{i_n,j_n}\right]\in {\mathbb{R}}^{I_n\times J_n}$ is the n-th factor matrix containing the features across the n-th mode of $\mathcal{Y}$. The symbols ${\times }_n$ and ∘ stand for the tensor-matrix product along the n-th mode, and the outer product, respectively. Tensor $\mathcal{G}$ has the multi-linear rank $(J_1,J_2,\dots ,J_N)$, i.e., $\forall n:\mathrm{rank}\left({\mathit{G}}_{\left(n\right)}\right)=J_n$, where ${\mathit{G}}_{\left(n\right)}$ is the unfolded version of $\mathcal{G}$ along the n-th mode.

Nonnegative Tucker decomposition (NTD) is a special case of the Tucker decomposition where the nonnegativity constraints are imposed onto all the factor matrices and the core tensor, i.e., $\forall i_n,j_n:{\mathit{U}}_{i_n,j_n}\ge 0,{\mathit{G}}_{j_1,j_2,\dots ,j_N}\ge 0$ for $n=1,\dots ,N$.

To estimate the factors in NTD, the alternating optimization scheme is typically used, where model (3) is unfolded along the n-th mode in each step. Thus: where ${\mathit{U}}^{{\otimes }_{-n}}\in {\mathbb{R}}_{+}^{{\prod }_{\mathit{P}\ne n}{I}_{\mathit{P}}\times {\prod }_{\mathit{P}\ne n}{J}_{\mathit{P}}}$ is a nonnegative matrix, and ⊗ denotes the Kronecker product.

Let $D\left({\mathit{Y}}_{\left(n\right)}\left|\right|\left[{\mathit{U}}^{\left(n\right)}{\mathit{G}}_{\left(n\right)}{\left({\mathit{U}}^{{\otimes }_{-n}}\right)}^T\right]\right)$ be the objective function that expresses dissimilarity between ${\mathit{Y}}_{\left(n\right)}$ and the unfolded version of the Tucker model. The optimization problem for estimation of the n-th factor in NTD can be presented in the generative form as the nonnegatively constrained nonlinear optimization problem: where

Core tensor $\mathcal{G}$ can be easily estimated by vectorizing the model (3). Thus where ${\mathit{U}}^{\otimes }={\mathit{U}}^{\left(N\right)}\otimes \dots \otimes {\mathit{U}}^{\left(1\right)}\in {\mathbb{R}}_{+}^{{\prod }_{n=1}^N{I}_n\times {\prod }_{n=1}^N{J}_n}$, $\mathit{Y}=\mathrm{vec}\left(\mathcal{Y}\right)\in {\mathbb{R}}_{+}^{{\prod }_{n=1}^N{I}_n}$ and $\mathit{g}=\mathrm{vec}\left(\mathcal{G}\right)\in {\mathbb{R}}_{+}^{{\prod }_{n=1}^N{J}_n}$ is vectorized core tensor $\mathcal{G}$.

The problems (5) and (7) can be solved using various constrained optimization solvers. In particular, when $D(·||·)$ is expressed by the Euclidean distance, problems (5), and (7) become the least-squares problems with nonnegativity constraints. The choice of the Euclidean distance is also motivated by the assumption that the residual error has a normal distribution, which is usually justified in practice. Assuming a wider class of data distributions, the $\beta$-divergence seems to be the best choice since it combines many well-known dissimilarity measures, including the generalized Kullback–Leibler divergence and the Itakura–Saito distance [1]. However, the proposed methodology explicitly involves the formulation of normal equations, which limits the profits of using other statistical measures than the Euclidean distance.

Without loss of generality, let us assume observed tensor $\mathcal{Y}\in {\mathbb{R}}_{+}^{I_1\times I_2\times \dots \times I_N}$ can be regarded as a sequence of N-way arrays concatenated along the N-th mode. Thus:

Let $\tilde{\mathcal{Y}}=\mathrm{cat}\left\{{\mathcal{Y}}^{\left(1\right)},\dots ,{\mathcal{Y}}^{(t-1)},N\right\}\in {\mathbb{R}}_{+}^{I_1\times \dots \times I_{N-1}\times {\tilde{I}}_N}$ be the N-way array representing the historical data, i.e., the multi-way arrays concatenated since the first tensor ${\mathcal{Y}}^{\left(1\right)}$ until the $(t-1)$-th array. We assume that ${\mathcal{Y}}^{\left(t\right)}\in {\mathbb{R}}_{+}^{I_1\times \dots \times I_{N-1}\times I_N^{\left(t\right)}}$ is the latest tensor which arrives in the t-th time slot or the t-th subwindow. Following this assumption, model (8) takes the form:

For the data having the structure in (9), the INTD has the following form: where ${\tilde{\mathit{U}}}^{\left(N\right)}\in {\mathbb{R}}_{+}^{{\tilde{I}}_N\times J_N}$ is the factor containing the past features of $\mathcal{Y}$ along the N-th mode, and ${\mathit{U}}_t^{\left(N\right)}\in {\mathbb{R}}_{+}^{I_N^{\left(t\right)}\times J_N}$ contains the current features that are valid in the t-th time slot.

According to the unfolding rule in (1) and using (4), we have:

Formula (11) implies that the current sample ${\mathcal{Y}}^{\left(t\right)}$ after unfolding with respect to the n-th mode ($n<N$) is located as the last $I_N^{\left(t\right)}$ columns in ${\mathit{Y}}_{\left(n\right)}$, whereas this sample after unfolding with respect to the N-th mode is located as the last $I_N^{\left(t\right)}$ rows of ${\mathit{Y}}_{\left(n\right)}$ according to (12).

In NTD, factor matrices ${\mathit{U}}^{\left(n\right)}$ for $n<N$ can be easily estimated from the strongly over-determined system of linear equations in (4). Assuming that the objective function in (5) is expressed by the Euclidean distance, system (4) can be transformed to the normal equations as follows:

Let ${\mathit{P}}_n={\mathit{Y}}_{\left(n\right)}{\mathit{U}}^{{\otimes }_{-n}}{\mathit{G}}_{\left(n\right)}^T\in {\mathbb{R}}_{+}^{I_n\times J_n}$ and ${\mathit{Q}}_n={\mathit{G}}_{\left(n\right)}{\mathit{U}}^{{\otimes }_{-n}T}{\mathit{U}}^{{\otimes }_{-n}}{\mathit{G}}_{\left(n\right)}^T\in {\mathbb{R}}^{J_n\times J_n}$. Inserting (11) into (13), we have for the t-th sample: where and

The sign $<·,·{>}_{-n}$ stands for the inner product (Let $\mathcal{A}=\left[a_{i_1,\dots ,i_N}\right]\in {\mathbb{R}}^{I_1\times \dots \times I_N}$ and $\mathcal{B}=\left[b_{i_1,\dots ,i_N}\right]\in {\mathbb{R}}^{I_1\times \dots \times I_N}$. The inner product of $\mathcal{A}$ and $\mathcal{B}$ along all but the n-th mode is defined as $C=<\mathcal{A},\mathcal{B}{>}_{-n}={\sum }_{i_1=1}^{I_1}\cdots {\sum }_{i_{n-1}=1}^{I_{n-1}}{\sum }_{i_{n+1}=1}^{I_{n+1}}\cdots {\sum }_{i_N=1}^{I_N}{a}_{i_1,\dots ,i_N}{b}_{i_1,\dots ,i_N}\in {\mathbb{R}}^{I_n\times I_n}$.) of tensors along all, but the n-th mode. Note that Formula (14) can be regarded as an iterative update rule with respect to the step t.

Since ${\mathit{U}}^{\left(N\right)T}{\mathit{U}}^{\left(N\right)}=\left[{\tilde{\mathit{U}}}^{\left(N\right)T},{\mathit{U}}_t^{\left(N\right)T}\right]\left[\begin{array}{c}{\tilde{\mathit{U}}}^{\left(N\right)}\\ {\mathit{U}}_t^{\left(N\right)}\end{array}\right]={\tilde{\mathit{U}}}^{\left(N\right)T}{\tilde{\mathit{U}}}^{\left(N\right)}+{\mathit{U}}_t^{\left(N\right)T}{\mathit{U}}_t^{\left(N\right)}$, matrix ${\mathit{Q}}_n$ for the t-sample can be expressed in a similar manner: where ${\mathit{U}}^{{\otimes }_{-\{n,N\}}}={\mathit{U}}^{(N-1)}\otimes \dots \otimes {\mathit{U}}^{(n+1)}\otimes {\mathit{U}}^{(n-1)}\otimes \dots \otimes {\mathit{U}}^{\left(1\right)}$, and

Since matrix ${\mathit{Q}}_n^{\left(t\right)}$ is square and symmetric, i.e., ${\mathit{Q}}_n^{\left(t\right)}={\mathit{Q}}_n^{\left(t\right)T}$, factor ${\mathit{U}}^{\left(n\right)}$ can be updated for $\forall n<N$ by solving the following equation:

To estimate factor matrix ${\mathit{U}}_t^{\left(N\right)}$, the system of linear equations in (12) is used. Note that where ${\mathit{Z}}_N=\left({\mathit{U}}^{(N-1)}\otimes \dots \otimes {\mathit{U}}^{\left(1\right)}\right){\mathit{G}}_{\left(N\right)}^T$. Hence: and

Transforming system (23) to the normal equations, we have: ${\mathit{Z}}_N^T{\mathit{Y}}_{\left(N\right)}^{\left(t\right)T}={\mathit{Z}}_N^T{\mathit{Z}}_N{\mathit{U}}_t^{\left(N\right)T}$. Following the similar calculations as in (14) and (17), we obtain:

Matrix ${\mathit{U}}_t^{\left(N\right)}$ can be readily estimated from the system: using any nonnegatively constrained linear least-squares solver.

Applying the concept of the Gauss–Seidel method for least-squares problems, factors $\left\{{\mathit{U}}^{\left(n\right)}\right\}$ can be updated with the following block-coordinate descent update rule: where $q_{j_n{j}_n}^{(n,t)}$ is the $j_n$-th diagonal entry of ${\mathit{Q}}_n^{\left(t\right)}$, and ${\mathit{u}}_{j_n}^{\left(n\right)}$, ${\mathit{p}}_{j_n}^{(n,t)}$, and ${\mathit{q}}_{j_n}^{(n,t)}$ are the $j_n$-th columns of ${\mathit{U}}^{\left(n\right)}$, ${\mathit{P}}_n^{\left(t\right)}$, and ${\mathit{Q}}_n^{\left(t\right)}$, respectively.

Since ${\mathit{Z}}_N$ is common in both systems (22) and (23), then ${\mathit{Q}}_N^{\left(t\right)}={\mathit{Q}}_N^{(t-1)}$, and

Following the similar strategy to estimate tensor $\mathcal{G}$, we have: which can be equivalently rewritten in the vectorized form: where $\mathit{g}=\mathrm{vec}\left\{\mathcal{G}\right\}\in {\mathbb{R}}_{+}^{{\prod }_{\mathit{P}}{J}_{\mathit{P}}}$. Vector g, estimated from (30) using any NNLS solver, is then tensorized to $\mathcal{G}$.

Due to the scaling ambiguity, the columns in factors must be scaled to the unit norm, i.e., $l_1$-norm. For $n<N$, the scaling was performed by mapping: $\forall n<N:{\mathit{U}}^{\left(n\right)}\leftarrow {\mathit{U}}^{\left(n\right)}\mathrm{diag}\left\{\left|\right|u_j^{\left(n\right)}{\left|\right|}_1^{-1}\right\}$, where $u_j^{\left(n\right)}$ is the j-th column of ${\mathit{U}}^{\left(n\right)}$. For $n=N$, the updated block is ${\mathit{U}}_t^{\left(N\right)}$, and we assume that block ${\tilde{\mathit{U}}}^{\left(N\right)}$ has been already normalized. Hence, we have the mapping: ${\mathit{U}}_t^{\left(N\right)}\leftarrow {\mathit{U}}_t^{\left(N\right)}\mathrm{diag}\left\{{\left(1+\left|\right|u_{j,t}^{\left(N\right)}{\left|\right|}_1\right)}^{-1}\right\}$, where $u_{j,t}^{\left(N\right)}$ is the j-th column of ${\mathit{U}}_t^{\left(N\right)}$.

The samples $\left\{{\mathcal{Y}}^{\left(t\right)}\right\}$ can be proceeded sequentially for $t=1,2,\dots$, however, this process must be initialized by performing a batch decomposition on a small initial sample $\tilde{\mathcal{Y}}$ using any NTF algorithm. Having the objects $\{\mathcal{G},{\mathit{U}}^{\left(1\right)},\dots ,{\tilde{\mathit{U}}}^{\left(n\right)}\}$ estimated from $\tilde{\mathcal{Y}}$, the matrices ${\mathit{P}}_N^{\left(0\right)}$, ${\mathit{Q}}_N^{\left(0\right)}$, ${\mathit{P}}_n^{\left(0\right)}$, and ${\mathit{Q}}_n^{\left(0\right)}$ are computed using (28), (25), (15), and (18) for $t=1$, respectively. The recursive incremental NTD (RI-NTD) algorithm for updating objects $\{\mathcal{G},{\mathit{U}}^{\left(1\right)},\dots ,{\mathit{U}}_t^{\left(n\right)}\}$ given $\left\{{\mathcal{Y}}^{\left(t\right)}\right\}$ for $t=1,2,\dots$ is expressed by Algorithm 1.

Algorithm 1:RI-NTD

Model (11) for $n<N$ can be equivalently expressed in the transposed form: where and

Matrix ${\mathit{Z}}^{\left(n\right)}=\left[\begin{array}{c}{\tilde{\mathit{Z}}}^{\left(n\right)}\\ {\mathit{Z}}_t^{\left(n\right)}\end{array}\right]$ has a block-structure. Assuming that $\mathrm{rank}\left({\tilde{\mathit{Z}}}^{\left(n\right)}\right)=J_n$ and $\mathrm{rank}\left({\mathit{Z}}_t^{\left(n\right)}\right)=J_n$, the k-th iterative update of matrix ${\mathit{U}}^{\left(n\right)T}$ for $k=0,1,\dots$ can be obtained using the block Kaczmarz method with nonnegativity constraints: where ${\left({\mathit{Z}}_t^{\left(n\right)T}\right)}^{†}={\mathit{Z}}_t^{\left(n\right)}{\left({\mathit{Z}}_t^{\left(n\right)T}{\mathit{Z}}_t^{\left(n\right)}\right)}^{-1}$ is the pseudo-inverse of ${\mathit{Z}}_t^{\left(n\right)T}$.

Matrix ${\mathit{U}}_t^{\left(N\right)}$ can be estimated from (23) using the same computational approach as presented in Section 2.2. Similarly, tensor $\mathcal{G}$ can be directly obtained from (30). The block Kaczmarz NTD (BK-NTD) for updating objects $\{\mathcal{G},{\mathit{U}}^{\left(1\right)},\dots ,{\mathit{U}}_t^{\left(n\right)}\}$ for $t=1,2,\dots$ is presented by Algorithm 2, which additionally has the inner iterations.

Algorithm 2:BK-NTD

Sweeping over the blocks across the N-th mode, the convergence rate of rule (36) can be analyzed using the concept of row paving in [50].

The proof of Lemma 1 can be easily obtained considering the result given in [50]. Lemma 1 shows that the convergence rate of rule (36) is linear.

The following benchmarks were used in the experiments:

The proposed algorithms: BK-NTD (Algorithm 2) and RI-NTD (Algorithm 1) were compared with the HALS-NTD [67], which is regarded in the literature as one of the most efficient algorithms for NTD. The maximum number of iterations in the HALS-NTD was set to 100 and the tolerance for early termination was equal to ${10}^{-5}$. In the incremental versions of NTD, we set ${\tilde{I}}_3=100$ (the number of samples in an initial block) and $I_n^{\left(t\right)}=100$ (the number of samples in a regular block). The number of inner iterations in BK-NTD was set to 5, i.e., $K_i=5$. The maximum number of iterations in the NNLS algorithm both for BK-NTD as well as for RI-NTD was set to 300.

The BK-NTD and RI-NTD were implemented in MATLAB 2016b, and the HALS-NTD was taken from the Matlab toolbox TENSORBOX (https://faculty.skoltech.ru/people/anhhuyphan (accessed on 1 August 2018). All of the algorithms were run on a machine supplied with a 4-core Intel Core i7 CPU, 64 GB RAM, and an SSD drive.

To analyze the susceptibility of the algorithms to noisy perturbations, the synthetic data were perturbed with an additive zero-mean Gaussian noise $\mathcal{N}(0,{\sigma }^2)$ with the variance ${\sigma }^2$ selected to satisfy a given level of signal-to-noise ratio (SNR).

The following tests were performed using the above mentioned benchmarks: