Robust Manifold Nonnegative Tucker Factorization for Tensor Data Representation

In this paper, we aim to explore robust manifold NTF methods to address the above challenges. Our contributions are threefold.

• 1.

  Three robust manifold NTF algorithms are proposed. They first apply a half-quadratic optimization algorithm to transform the intractable problem into a weighted NTF model where the weights can handle the outliers. They can be implemented by incorporating different prior distribution for the measurement between the original data and reconstructed data. Unlike the Gaussian distribution of loss functions, other flexible distribution of metric can handle outliers more efficiently. The weights are adjusted adaptively with respect to the error. To our knowledge, this is the first work of NTF for robust weighted learning to handle outliers.

• 2.

  Three RMNTF algorithms reach the state-of-the-art performance are proposed. The three algorithms fall into the one major subclasses of NTF technologies, denoted as weighted NTF. Different distribution of loss function between original data and reconstructed data may influence the ability to handle outliers. Specifically, RMNTF with correntropy induced metric (RMNTF-CIM), RMNTF with Huber function (RMNTF-Huber), and RMNTF with Cauchy function (RMNTF-Cauchy) have been studied in this paper.

• 3.

  We demonstrate the convergence, robustness and invariance of RMNTF.

In this work, we first introduce some preliminaries and related work in the following two subsections, then present three RMNTF algorithms in Section 2. Section 3 presents the experimental results. Finally, Section 4 concludes our findings.

NMF: Given a non-negative data matrix $\mathbf{X}=[x_{1},\dots,x_{N}]\in\mathbb{R}_{\geq 0}^{N\times M}$, where $N$ refers to the number of data points and $M$ indicates the dimension of the feature. The objective of NMF is to find two nonnegative and low rank factor matrices $\mathbf{U}\in\mathbb{R}_{\geq 0}^{N\times r}$ and $\mathbf{V}\in\mathbb{R}_{\geq 0}^{M\times r}$ and the product of the two matrices approximates the original data matrix $\mathbf{X}$ by $\mathbf{X}\approx\mathbf{U}\mathbf{V}^{T}$, generally $r\ll\min\{M,N\}$. NMF incorporates the nonnegativity constraint and obtains the part-based representation as well as enhancing the interpretability [8] and these methods have a close relation with K-means [18]. However, the factorization of matrices is generally nonunique and many regularizors [19] or constraints [20] have been developed to alleviate nonuniqueness of decomposition. In particularly, NMF methods are proposed to utilize the priors to obtain a better representation [21] and achieve a robust clustering [22, 23].

However, the multiway structure information of high-dimensional data such as RGB images or videos cannot be represented by NMF. Vectorization of the high-dimensional data may lead to excess numbers of parameters and structure interruptions. Nonnegative tensor decomposition based methods have been proposed for the high-dimensional data.

Nonnegative CANDECOMP/PARAFAC (NCP) methods: The CP-based tensor model (Fig. 1) [24] decomposes $\mathcal{X}$ into a linear combination of $R$ rank-one tensors as follows:

where $\mathrm{Diag}^{N}(\boldsymbol{\lambda})$ is the $N$th-order diagonal tensor. The CP rank of $\mathcal{X}$ is given by $\mathrm{Rank}_{cp}(\mathcal{X})=R$ denotes the smallest number of the rank-one tensor decomposition [11]. The CP-based tensor model assumes each tensor element can be calculated by a summation of $R$ products and it is restrictive since it only considers $R$ possible interactions between latent factors.

Existing CP-based tensor models have a flexible subspace representation which may not consider the intrinsic manifold information of high-dimensional data and the performance of downstream tasks will be limited. A few regularization and prior knowledge strategies have been studied in these models [11]. For example, Zhao et al. [25] formulated CP factorization using a hierarchical probabilistic model and employed a fully Bayesian treatment by incorporating a sparsity-inducing prior over multiple latent factors and the appropriate hyperpriors over all hyperparameters to determine the rank of models. Zhao et al. [26] proposed a probabilistic model to recover the underlying low-rank tensor which modeled by multiplicative interactions among multiple groups of latent factors, and the additive sparse tensor modeling outliers. Zhou et al. [11] introduced concurrent regularizations which regularized the entire subspace in a concurrent and coherent way to avoid the strong scale restrictions of $L_{2}$ regularization. Chen et al. [27] proposed a generalized weighted low rank tensor factorization which represented the sparse component as a mixture of Gaussian, and unified the Tucker and CP factorization in a joint framework to handle complex noise and outliers. However, these methods neglect nonnegative constrains and may not learn the part-based and physical meaning representations.

NTF methods: The Tucker model (Fig. 2) assumes that an original tensor $\mathcal{X}$ can be well approximated as

where the tucker rank of $\mathcal{X}$ is denoted as $\mathrm{Rank}_{\mbox{tc}}(\mathcal{X})=(R_{1},\dots,R_{n},\dots,R_{N})$ with $R_{n}=\mathrm{Rank}(\mathbf{X}_{(n)})$.

We note that the CP decomposition is a special case of the Tucker decomposition. Although the tucker decomposition is invariant to rotations in the factor matrices, it shares parameters across latent factor matrices by a core tensor. In contrast, CP decomposition methods force each factor vector to capture potentially redundant information [28]. CP decomposition methods are more prone to overfitting than Tucker decompositions. The main interest in Tucker model is to find subspaces for tensor approximation.

Li et al. [29] introduced a manifold regularization into the core tensors of NTF which preserved the tensor geometry information. But the representation space of the core tensor will increase exponentially as tensor order increases, which results in high computational complexity. Jiang et al. [30] added a graph Laplacian regularization on a low-dimensional factor matrix to improve the robustness of tensor decomposition. Sun et al. [5] proposed a heterogeneous tensor decomposition for clustering by performing dimensionality reduction on the first $N-1$ order of the tensor, and incorporating some useful constraints on the last-mode factor matrix for clustering. However, these methods neglect nonnegative constraints and may loss physical meaning of low-rank representations. Yin et al. [31] incorporated Laplacian Eigenmaps and Locally Linear Embedding as the manifold regularization terms into the least square form of NTF model. Pan et al. [32] introduced the orthogonal constraint into the group of factor matrices of NTF, which not only helps to keep the inherent tensor structure but also well performs in data compression. Yin et al. [33] proposed Hypergraph Regularized NTF which preserved nonnegativity in tensor factorization and uncover the higher-order relationship among the nearest neighborhoods. However, these methods decompose tensor data by minimizing the Euclidean distance which fails on the not cleaned data.

Manifold learning: Manifold learning is a problem which encodes the geometric information of the data space. Its goal is to find a representation in which two objects are close to each other after dimension reduction if they are close in the intrinsic geometry of data manifold. Based on this idea, many types of manifold learning algorithms have been proposed, such as ISOMAP [34], LLE [35], Laplacian eigenmaps [36] and locality preserving projections [37]. The advantage of introducing manifold structures is that it can preserve the intrinsic geometric information of data points, and it has been shown to be useful in a wide-range of applications, such as face recognition [38], text mining [19], and multimedia interaction [39].

Recently, the idea of manifold learning has been employed to matrix and tensor analysis. For instance, Cai et al. [19] proposed graph regularized nonnegative matrix factorization (GNMF) which utilized the intrinsic geometric information. However, GNMF may not have optimal solutions due to the noises or outliers of data, consequently other extension methods based on the GNMF have been proposed. Moreover, NMF based on manifold learning methods may break the structure of mutliway data points. Hence, tensor based manifold learning methods have been proposed such as Graph regularized Nonnegative Tucker Decomposition [40], LLE based nonnegative tensor decomposition [31]. These manifold nonnegative tensor decomposition methods assume that the distribution of noise is Gaussian and may fail on grossly corrupted datasets.

Half-Quadratic minimization was pioneered by Geman and Reynolds [42], and it was used to alleviate the computational task in the context of image reconstruction with nonconvex regularization. Let $\mathcal{X}$ denote an original tensor and $\hat{\mathcal{X}}$ denote a reconstruction tensor of $\mathcal{X}$. Replacing the squared residual of data-fidelity terms in [43] on each entry with a generic function:

where $\mathcal{E}$ represents the residual error between original tensor $\mathcal{X}$ and reconstructed tensor $\hat{\mathcal{X}}$, $g(\cdot)$ is chosen to be robust to outliers or gross errors, and $\Phi(\cdot)$ denotes the regularization terms with respect to $\{\mathbf{A}^{(n)}\}_{n=1}^{N}$. The minimizer $\hat{\mathcal{X}}$ of cost function $\mathcal{J}(\mathcal{X},\hat{\mathcal{X}})$ involving the reconstruction error which is nonlinear with respect to $\mathcal{X}$ and the regularization term $\Phi(\{\mathbf{A}^{(n)}\}_{n=1}^{N})$. When factor matrices $\mathbf{A}^{(n)},n=1,\dots,N-1$ and core tensor $\mathcal{S}$ have many nonzero entries or ill-conditioned, the computation of factorization is costly. specifically, the loss function is possibly non-quadratic and non-convex, and it is difficult to optimize directly. Fortunately, the half-quadratic minimization [44] has been developed to solve the intractable optimization. According to the conjugate function and half quadratic theory [44], the reconstruction error term $g(\mathcal{E})$ can be performed as

where $\phi(\mathcal{W})$ is the conjugate function of $g(\mathcal{E})$, $\mathcal{W}$ is the corresponded additional auxiliary variable, and $Q(\cdot,\cdot)$ is a quadratic term for $\mathcal{E}$ and $\mathcal{W}$. In this paper, we only consider the quadratic term of the multiplicative form [42]:

Substituting (4) and (5) into (3), we have the augmented cost function $\hat{\mathcal{J}}(\mathcal{X},\hat{\mathcal{X}};\mathcal{W})$:

The reconstruction error term involved in $\hat{\mathcal{J}}(\mathcal{X},\hat{\mathcal{X}};\mathcal{W})$ is half-quadratic. Hence, the minimizer $(\hat{\mathcal{X}},\mathcal{W})$ of $\min_{\hat{\mathcal{X}},\mathcal{W}}\left\{\hat{\mathcal{J}}(\mathcal{X},\hat{\mathcal{X}};\mathcal{W})\right\}$ is calculated by alternate minimization. At iteration $k$ we calculate

When $\hat{\mathcal{X}}$ is fixed, the minimization of the reconstruction error term $g(\cdot)$ is convex with respect to $\mathcal{W}$. The explicit optimum solution of $\mathcal{W}$ [45] can be determined as

where $\frac{[\cdot]}{[\cdot]}$ denotes element-wise division.

It is shown that the auxiliary variable $\mathcal{W}$ only depends on the loss function $g(\cdot)$. Since the outliers often cause large fitting errors, $\mathcal{W}$ is important for the objective functions to constrain overfitting. For the large outliers, the weights $\mathcal{W}$ should be small. On contrary, for the small errors, the weights $\mathcal{W}$ should be large. Therefore, the weights variable $\mathcal{W}$ can be seen as an outlier mask. The frequency used loss functions are shown in Fig. 3.

By using the nonnegative constraints and robust loss function for outliers, robust NTF can learn a part-based representation. Many robust NTF methods perform well in euclidean space. They fail to discover the intrinsic geometrical and discriminating structure of the data space. Here, we introduce a geometrically based regularization for our robust NTF framework.

First, supposing that the real data points $\mathcal{X}_{i}$ lies on a low-dimensional manifold $\mathcal{M}$ and $\mathbf{A}^{(N)}_{i}$ is the representation of $\mathcal{X}_{i}$ in the subspace. We make an assumption that if $\mathcal{X}_{i}$ and $\mathcal{X}_{j}$ are neighbors in data space, then their low-rank representations $\mathbf{A}^{(N)}_{i}$ and $\mathbf{A}^{(N)}_{j}$ are close enough to each other in $\mathcal{M}$. We build a regularization term as follows

where $f_{p}$ is the mapping function which project data point $\mathcal{X}_{i}$ to the low-rank representation $\mathbf{A}^{(N)}_{i}=f_{p}(\mathcal{X}_{i})$. $\mathcal{R}(\mathcal{M})$ is a measurement of the smoothness of $f_{p}$ along the geodesics in the intrinsic geometry of the data.

Based on [19], we use the similarity graph of $\mathcal{X}$. Suppose that $\mathbf{V}\in\mathbb{R}^{I_{N}\times I_{N}}$ defines the affinity of data points, then we use the Heat Kernel to describe the similarity between each pair of data points if nodes $i$ and $j$ are connected:

where $\tau$ is the width of the kernel used to control the similarity. Then, we calculate the diagonal matrix $\mathbf{D}\in\mathbb{R}^{I_{N}\times I_{N}}$, where $[\mathbf{D}]_{ii}=\sum_{j}[\mathbf{V}]_{ij}$ and the Laplacian matrix is $\mathbf{L}=\mathbf{D}-\mathbf{V}$. The graph regularization can be estimated as follows:

Liu et al. [46] proposed the concept of Cross correntropy which is a generalized similarity measure between two arbitrary scalar random variables $\mathbf{X}$ and $\mathbf{Y}$ defined by

where $\kappa_{\sigma}(\cdot)$ is the kernel function. In practice, the joint probabilistic density function is unknown and only a finite number of data points $\{(x_{i},y_{i})\}_{i=1}^{N}$ are available. The sample estimator of correntropy can be represented as

Based on the above definition of correntropy, Liu et al. [46] proposed correntropy induced metric (CIM) in the sample space which denoted as

where we use the Gaussian kernel in this paper, i.e., $\kappa_{\sigma}(e)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-e_{i}^{2}/2\sigma^{2}\right)$, and $e_{i}$ is denoted as $e_{i}=x_{i}-y_{i}$.

Substituting the error on each entry in Tucker model with the CIM, we obtain the objective function of RMNTF-CIM:

which is equivalent to solving the following optimization problem

Then, we introduce the half-quadratic minimization, $\mathcal{J}(\mathcal{X},\hat{\mathcal{X}})=\frac{1}{2}Q(\mathcal{W}\circledast\mathcal{E}^{2})+\phi(\mathcal{W})$, where $\phi(\mathcal{W})$ is the conjugate function of $\mathcal{J}(\mathcal{X},\hat{\mathcal{X}})$.

Robust statistics work well on model reconstruction under some observation with noise or outliers. Some popular M-estimators [47] such as Huber loss function and Cauchy function have been proposed to solve noisy data mining.

In this section, we take the Huber function in reconstruction error term $g(\cdot)$ to measure the quality of approximation by considering the connection between $L_{1}$ norm and $L_{2}$ norm:

where $c$ is the cutoff parameter to tradeoff between the $L_{1}$-norm and $L_{2}$-norm.

Substituting the Huber function on each entry in (3), we have the RMNTF-Huber by minimizing the following objective function:

where $\mathcal{E}=\mathcal{X}-\hat{\mathcal{X}}=\mathcal{X}-\mathcal{S}\times_{n=1}^{N}\mathbf{A}^{(n)}$.

Following the equation (10), we obtain the optimization of weighted tensor $\mathcal{W}$:

The optimization of factor matrices $\mathbf{A}^{(n)}$ and the core tensor $\mathcal{S}$ are the same as RMNTF-CIM. Here, we note that the cutoff parameter $c$ is set to the median of reconstruction errors, i.e., $c=\mathrm{median}(|\mathcal{E}|)$.

For any tensor $\mathcal{X}$, we define the reconstruction error of RMNTF-Cauchy as follows:

where $g_{\rm{Cauchy}}(x)=\ln(1+x),x\geq 0$.

As equation (10), the Optimization of weighted tensor $\mathcal{W}$ can be represented as:

The optimization of factor matrices $\mathbf{A}^{(n)}$ and the core tensor $\mathcal{S}$ are the same as RMNTF-CIM.

The convergence of the algorithms are guaranteed by the following theorem:

The robustness of RMNTF is guaranteed by the following Theorem:

The uniqueness of RMNTF is guaranteed by the following Theorem:

We conducted experiments on the COIL100, USPS, FEI, ORL and FERET image datasets.

COIL100¹¹1https://www.kaggle.com/jessicali9530/coil100 is an object categorization image database including 100 classes of objective, each of which contains $72$ images with $128\times 128$ of difference observation angle. For reprocessing, we resize each image into $64\times 64$ with the RGB representation by nearest neighbor interpolation algorithm. As a result, each object is represented as a $64\times 64\times 3$ tensor. In total, we have 7200 tensor objects.

USPS Dataset²²2https://www.kaggle.com/bistaumanga/usps-dataset consists of $11000$ images of handwritten digits $0\sim 9$ with $16\times 16$. Each of handwritten digits contains $1100$ images.

FEI Part 1 Dataset³³3https://fei.edu.br/~cet/facedatabase.html is the subset of FEI database, which consists of 700 color images of size $480\times 640$ collected from 50 individuals. Each individual has 14 different images under different observed and facial expressions. In our experiment, the images are resized to $48\times 64$. These image finally construct a fourth-order $48\times 64\times 3\times 700$ tensor.

AT & T ORL Dataset⁴⁴4https://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html: is the subset of FEI database, which consists of 700 color images of size $480\times 640$ collected from 50 individuals. Each individual has 14 different images under different observed and facial expressions. In our experiment, the images are resized to $48\times 64$ pixels. These image finally construct a fourth-order $48\times 64\times 3\times 700$ tensor.

FERET Dataset⁵⁵5https://www.nist.gov/itl/products-and-services/color-feret-database: This dataset [48] collects $14125$ grayscale face images acquired from $1199$ subjects. It is widely used for evaluating face recognition and clustering problem. In our experiment, we use a subset of $1400$ images with a size of $80\times 80$ from $200$ subjects where each subject contains $7$ images with varying postures, genders, lighting condition, shooting directions and race. In total, we stack the images into a third-order tensor with size of $80\times 80\times 1400$.

For COIL100, we randomly selected $5,10,20,30,50$ categories from the whole COIL100 for evaluation. For USPS, we randomly selected $3,5,7$ categories and the whole $10$ categories for the evaluation. For FEI, ORL and FERET data sets, we randomly selected $5,10,15,20$ categories for the evaluation. For each comparison, we reported the average results over $50$ Monte-Carlo runs.

The hyperparameters of RMNTF-CIM, RMNTF-Huber, and RMNTF-Cauchy in all experiments were set as follows: $\lambda=10^{5},p=3$, unless otherwise stated. We compared RMNTF with three NMF methods, and six NTD methods. They are listed as follows:

• 1.

  Nonnegative Matrix Factorization (NMF) [8]: We unfold tensorial original data into a matrix to apply it to NMF algorithm and obtain a low rank activation matrix $\mathbf{V}$. Finally, we apply $\mathbf{V}$ to downstream tasks.

• 2.

  Graph Regularization Nonnegative Matrix Decomposition (GNMF) [19]: We unfold tensorial data into a matrix $\mathbf{X}$ and reduce the dimensionality of $\mathbf{X}$ into $K$ by GNMF, where the dictionary matrix $\mathbf{U}$ is used in graph regularization term. The coefficient of regularization $\lambda=10$, and the hyperparameter $p=3$.

• 3.

  Correntropy Induced Metric Nonnegative Matrix Factorization (CIMNMF) [49]:Similar to the above, we unfold the tensorial data into a matrix and apply it to CIMNMF. The loss function between the unfolding tensor and reconstruction data is based on the correntropy induced metric. It handles outlier rows by incorporating structural knowledge about outliers.

• 4.

  Nonnegative Tucker Decomposition (NTD) [1]: NTD works directly on the tensorial data. It is constructed by the tucker model with nonnegativity constraints and updated by multiplicative rules. This method can be considered as a special case of RMNTF when all the entries of weight tensor $\mathcal{W}$ are $1$, and $\lambda=0$.

• 5.

  Nonnegative Tucker Decomposition with alpha-Divergence ($\alpha$NTD) [50]: We directly use the tensorial data to $\alpha$NTD. It is considered $\alpha$-divergence as a discrepancy measure and derive multiplicative updating rules for NTD. The hyperparameter $\alpha=1$.

• 6.

  Sparse Nonnegative Tucker Decomposition (SparseNTD) [51]: The tensorial data is directly utilized in SparseNTD. It adopts the sparsity and nonnegativity constraints to a core tensor and several factor matrices. The hyperparameter $\lambda=0.5$.

• 7.

  Graph Regularized Nonnegative Tucker Decomposition (GNTD) [40]: GNTD use tensorial data and it constructs the nearest neighbor graph to maintain the intrinsic manifold structure of tensor and applies this constraint on the $N$th factor matrix. The hyperparameters $\lambda=10^{5}$, and $p=3$.

• 8.

  Locally Linear Embedding Regularized Nonnegative Tucker Decomposition (LLENTD) [31]: We directly use tensor in LLENTD. It incorporates Laplacian Eigenmaps as manifold regularization terms into the least square form of nonnegative tucker model. The hyperparameters $\lambda=10^{5}$, and $p=3$.

• 9.

  Manifold Regularization Nonnegative Tucker Decomposition (MRNTD) [29]: MRNTD utilizes tensor data directly. It employs the manifold regularization terms for the core tensors constructed in the NTD. The hyperparameters $\lambda=10^{5}$, and $p=3$.

• 10.

  Hypergraph regularized nonnegative tensor factorization (hyperNTD) [33]:We take tensorial data into hyperNTD directly. It incorporates a higher-order relationship among the nearest neighborhoods and the nonnegative tucker decomposition. The hyperparameters $\lambda=10^{5}$, and $p=3$.

To evaluate the clustering performance of these algorithms, we adopt two commonly used metrics: 1) accuracy(ACC) and 2) normalized mutual information (NMI). The ACC is defined by

where $n$ is the total number of samples in a dataset. We have the cluster labels ${\mathbf{\hat{l}}}=\{{\hat{l}_{1}},{\hat{l}_{2}},\cdots,{\hat{l}_{n}}\}$ and the ground-truth labels ${\mathbf{l}}=\{{l_{1}},{l_{2}},\cdots,{l_{n}}\}$. $\delta(x,y)$ is set to $1$ if and only if $x=y$, and $0$ otherwise. $\operatorname{map}(\cdot)$ is a displacement mapping function that maps each cluster label $\hat{l_{i}}$ to the equivalent label from the dataset.

Another metric NMI is defined by

where $\operatorname{H}(\mathbf{l})$ and $\operatorname{H}(\mathbf{\hat{l}})$ denote the entropy of $\mathbf{l}$ and $\mathbf{\hat{l}}$, respectively, and