## 2 Basic Theory of Non-Negative Tensor Factorization
### Hierarchical Alternating Least Squares Algorithm for NTF  
The HALS algorithm [24] uses a set of local cost functions to obtain the simple update rules in vector-wise form, and is robust for both under-determined models and over-determined models. It is also suitable for solving problems relating to sparse representation and large-scale matrices [25; 26]. The extensions of HALS algorithms are also developed based on alpha divergence and beta divergence for rank-one tensor approximation [26].  
For a third order tensor \(\underline{\mathbf{Y}}\in\mathbb{R}^{I\times T\times Q}\), we consider the decomposition based on the parallel factor analysis (PARAFAC) model for efficiency, in which the tensor \(\underline{\mathbf{Y}}\) is estimated by the outer product of the vectors from three component factors \(\mathbf{A}=[\mathbf{a}_{1},\mathbf{a}_{2},\ldots,\mathbf{a}_{j}]\in\mathbb{R} ^{I\times J}\), \(\mathbf{B}=[\mathbf{b}_{1},\mathbf{b}_{2},\ldots,\mathbf{b}_{J}]\in\mathbb{R }^{T\times J}\), \(\mathbf{C}=[\mathbf{c}_{1},\mathbf{c}_{2},\ldots,\mathbf{c}_{J}]\in\mathbb{R }^{Q\times J}\); i.e.,  
\[\underline{\mathbf{Y}}=\sum_{j=1}^{J}\mathbf{a}_{j}\circ\mathbf{b}_{j}\circ \mathbf{c}_{j}+\underline{\mathbf{E}}\Leftrightarrow y_{itq}=\sum_{j=1}^{J}a _{ij}b_{tj}c_{qj}+e_{itq} \tag{9}\]  
The concept is to define residue matrices as:  
\[\begin{array}{ll}\underline{\mathbf{Y}}^{(j)}&=\underline{\mathbf{Y}}-\sum \limits_{p\
eq j}\mathbf{a}_{p}\circ\mathbf{b}_{p}\circ\mathbf{c}_{p}=\underline {\mathbf{Y}}-\sum\limits_{j=1}^{J}\mathbf{a}_{j}\circ\mathbf{b}_{j}\circ \mathbf{c}_{j}+\mathbf{a}_{j}\circ\mathbf{b}_{j}\circ\mathbf{c}_{j}\\ &=\underline{\mathbf{E}}+\mathbf{a}_{j}\circ\mathbf{b}_{j}\circ\mathbf{c}_{j}, (j=1,2,\ldots J),\end{array} \tag{10}\]  
and alternatively minimize the set of local cost functions according to the ALS algorithm:  
\[D_{F}^{(j)}(\underline{\mathbf{Y}}^{(j)}\|\mathbf{a}_{j}\circ\mathbf{b}_{j} \circ\mathbf{c}_{j})=\frac{1}{2}\|\underline{\mathbf{Y}}^{(j)}-\sum_{j=1}^{J} \mathbf{a}_{j}\circ\mathbf{b}_{j}\circ\mathbf{c}_{j}\|_{F}^{2}=\frac{1}{2}\| \underline{\mathbf{Y}}^{(j)}_{(1)}-\mathbf{a}_{j}\mathbf{c}_{j}^{T}\odot \mathbf{b}_{j}^{T}\|_{F}^{2}, \tag{11}\]  
where \(\underline{\mathbf{Y}}^{(j)}_{(1)}\) is the mode-1 unfolding of the residue \(\underline{\mathbf{Y}}^{(j)}\), and \(\odot\) denotes the Kahatri-Rao product.  
The stationary points of the above minimization are determined by the KKT optimality conditions given as  
\[\mathbf{a}_{j},\mathbf{b}_{j},\mathbf{c}_{j}>0,\]  
\[\
abla_{\mathbf{a}_{j}}D_{F}^{(j)}(\underline{\mathbf{Y}}^{(j)}\|\mathbf{a}_{ j}\circ\mathbf{b}_{j}\circ\mathbf{c}_{j})=0, \tag{12}\]  
\[\
abla_{\mathbf{b}_{j}}D_{F}^{(j)}(\underline{\mathbf{Y}}^{(j)}\|\mathbf{a}_{ j}\circ\mathbf{b}_{j}\circ\mathbf{c}_{j})=0, \tag{13}\]  
\[\
abla_{\mathbf{c}_{j}}D_{F}^{(j)}(\underline{\mathbf{Y}}^{(j)}\|\mathbf{a}_{ j}\circ\mathbf{b}_{j}\circ\mathbf{c}_{j})=0. \tag{14}\]Hence, the update rule for \(\mathbf{a}_{j}\) can be obtained as follows  
\[\mathbf{a}_{j}\leftarrow\left[\frac{\mathbf{Y}_{(1)}^{(j)}(\mathbf{b}_{j}\odot \mathbf{c}_{j})}{(\mathbf{b}_{j}\odot\mathbf{c}_{j})^{T}(\mathbf{b}_{j}\odot \mathbf{c}_{j})}\right]_{+}=\left[\frac{\mathbf{Y}_{(1)}^{(j)}(\mathbf{b}_{j} \odot\mathbf{c}_{j})}{(\mathbf{b}_{j}^{T}\mathbf{b}_{j}).*(\mathbf{c}_{j}^{T} \mathbf{c}_{j})}\right]_{+},(j=1,2,\ldots,J), \tag{15}\]  
where \(*\) represents the Hadamard product. Especially when the vectors are correspondingly normalized, i.e., \(\mathbf{b}_{j}/\|\mathbf{b}_{j}\|_{2}\) and \(\mathbf{c}_{j}/\|\mathbf{c}_{j}\|_{2}\), the above rule can be further simplified as:  
\[\mathbf{a}_{j}=\left[\mathbf{Y}_{(1)}^{(j)}(\mathbf{b}_{j}\odot\mathbf{c}_{j })\right]_{+}. \tag{16}\]  
Finally, the updaterules for \(\mathbf{b}_{j}\) and \(\mathbf{c}_{j}\) are analogously accessed.