Robust Tensor Clustering with Non-Greedy Maximization

1 Introduction

Tensors can be taken as a natural representation of visual data. For example, an image can be represented by a second-order tensor [Ma et al., 2013] and a time-series video clip can be represented by a third-order tensor. Compared with the relaxation of tensor by vectorization an image into a long vector, tensor-based representation can preserve more spatial information of images. Therefore, tensor-based representation is more informative than the vectorization method. Additionally, the dimensionality of tensor-based representation is lower, which can effectively avoid high computational complexity, and large memory requirements generated by the vectorization of data [Guo et al., 2012].

Tensor clustering seeks to partition several m-order input tensors into groups by the minimization of certain criterion. It is a basic data analysis task and receives much attention recently in data mining [Zha et al., 2008], computer vision

[Guo et al., 2012], and so on. Like clustering on the vectors, tensor clustering is NP-hard to optimize, which needs additional optimization algorithms to approximatively solve it [Huang et al., 2008; Jegelka et al., 2009]. However, it is known that there exist outlying data points in the real-world datasets. When we use the above clustering algorithms to mine the structural information existing within these data, the outliers will reduce the performance of clustering algorithms, resulting in some misclustered sample data.

This urges us to develop a tensor clustering algorithm that is robust to the outliers. We call it the Robust Tensor Clustering (RTC). Our method firstly finds a lower rank approximation of the original tensor data using a L1 norm optimization function. Because the L1 norm doesn't exaggerate the effect of outliers compared with L2 norm, the minimization of the L1 norm approximation function makes RTC robust to outliers. Then we compute the HOSVD decomposition of this approximate tensor to obtain the final clustering results. We formulate the process of approximation into a framework of Tensor PCA with L1-norm (TPCA-L1) [Pang et al., 2010]. TPCA-L1 is a variant of Tensor PCA with Frobenius norm (TPCA-F) [De Lathauwer et al., 2000], where the Frobenius norm is replaced with the L1-norm. In this way, TPCAL1 will suppress the negative influence of outliers [Kwak, 2008]. However, [Pang et al., 2010] uses a greedy strategy like [Kwak, 2008] to maximize the function of TPCA-L1. Specifically, the projection directions are sequentially optimized one by one. This greedy method is easy to get stuck in a local solution. To overcome this drawback, we utilize a non-greedy maximization to directly solve the projection directions inspired by the work in [Nie et al., 2011], and thus much better solution can be obtained.

To sum up, the main contributions of this paper are as follows: 1) We propose a Robust Tensor Clustering (RTC) algorithm; 2) We utilize a non-greedy maximization method and then integrate it into the framework of RTC. Finally, we evaluate the performance of RTC on five public datasets and compare it with five state-of-the-art methods. Experimental results demonstrate the effectiveness of our method.

The rest of this paper is organized as follows. We present the framework of Robust Tensor Clustering with Non-Greedy Maximization in Section 2. The solution and algorithm are given in Section 3. Section 4 reports all experimental results. Finally, we summarize the conclusions in Section 5.

2 The Proposed Framework

Consider a sequence of input M-order (mode) tensor data Xj ∈ Rn1×n2×···×nM , j = 1, · · · , n, the goal of clustering is to find a K-way disjoint partitioning (S1, ..., SK) such that the following objective is minimized:

where Ck ∈ Rn1×n2×···×nM denotes the k-th cluster centroid. K is the number of clusters. ‖ · ‖ denotes the Euclidean distance between Xj and Ck. We introduce an indicator vector qk for the cluster Sk.

where qT k qk is the size of cluster Sk. Now the matrix Q is defined such that the k-th column of Q equal to qk/(qT k qk)1/2 . Q is an orthonormal matrix, QT Q = I, and Q ∈ Rn×K .

The RTC aims to find a more accurate indicator matrix Q′ ∈ Rn×K via making the clustering algorithm robust to outliers. To this end, RTC maximizes the following problem to get the projection matrix, and then uses the projection matrix to solve approximation for the original tensor data:

where Ui ∈ Rni×ri (i = 1, · · · ,M) is called projection matrix (which is orthogonal, i.e., U T i Ui = I). Generally, ri < ni, so we can obtain a low-rank approximation Yj for Xj. Yj ≈ Xj ∏Mi=1 ×iU T i , where Xj ∏Mi=1 ×iU T i = Xj ×1 U T 1 · · · ×M U TM , and Xj ×i U T i denotes the i-mode product of the tensor Xj with the matrix Ui ∈ Rni×ri . Please refer to [Kolda and Bader, 2009] for the formal definition of i-mode product. ‖ · ‖1 is the tensor's L1-norm. The L1-norm of a M-order tensor X ∈ Rn1×n2×···×nM is defined as the sum of the absolute value of its elements.

The problem in (3) is called TPCA-L1 proposed in [Pang et al., 2010]. TPCA-L1 is a variant of TPCA-F [De Lathauwer et al., 2000], replacing the Frobenius norm with the L1-norm, and thus can suppress the negative influence of outliers [Kwak, 2008].

Once Ui ∈ Rni×ri , i = 1, · · · ,M are computed, the lowrank approximation Yj for Xj can be get. Thus we obtain a T-order tensor Z ∈ Rr1×···×rM×n comprising of n low-rank approximation Yj ∈ Rr1×r2×···×rM , where T = M+1. We can decompose the tensor Z to obtain the final clustering results for Xj, (j = 1, · · · , n) by solving the following optimization function:

where S ∈ Rd1×d2×···×K is called the core tensor. Pi ∈ Rri×di (i = 1, · · · ,T-1) and PT ∈ Rn×K are called factor matrices. The solution to the optimization problem of (5) can be found by HOSVD [De Lathauwer et al., 2000]. [Huang et al., 2008] demonstrates that PT ∈ Rn×K is the clustering results for n two-order tensors Xj. i.e., PT ∈ Rn×K is the indicator matrix of clustering n tensor Xj into K cluster centroids. In this way, we obtain Q′ = PT ∈ Rn×K .

3 Solution and Algorithm

In this section, we firstly discuss the non-greedy maximization algorithm for the vector case, and then extend this nongreedy strategy to handle tensor data.

Given a sequence of vectors (one-order tensor) xj ∈ Rd , j = 1, · · · , n, then problem (3) is formulated as follows:

where U ∈ Rd×m is the projection matrix, and d > m. Without loss of generality, xj, j = 1, · · · , n are assumed to be centralized, i.e., ∑nj=1 xj = 0.

The problem (6) should be replaced with Eq. (7):

where αj = sgn(U T xj). sgn(·) denotes the sign function defined as follows: sgn(x) = 1 if x > 0, sgn(x) = -1 if x < 0, and sgn(x) = 0 if x = 0. Denote M = ∑nj=1 xjαT j , [Nie et al., 2011] proves that U can be obtained by computing the SVD of M as M = WΛV T , and then U = WV T . This nongreedy strategy directly solve the projection direction U, and has been demonstrated to get much better solution than the non-greedy method [Kwak, 2008].

In the next paragraphs, we will discuss how to extend this non-greedy strategy to handle tensor data. Because the optimization algorithm is easy to be generalized to much higher order tensors from two-order tensors [Pang et al., 2010], we constrain to a sequence of second-order tensors from now on. Given a sequence of second-order tensors Xj ∈ Rn1×n2 , j = 1, · · · , n (i.e., a sequence of matrices, we use Xj to replace Xj in the following paragraphs). Problem (3)

is formulated as follows:

where U1 ∈ Rn1×r1 ,U2 ∈ Rn2×r2 are projection matrices and n1 > r1, n2 > r2. We use an alternative optimization technique to solve problem (8). Specially, we shall compute U1 while fixing U2 and compute U2 while fixing U1. The procedure is iteratively performed until the convergence criteria is meeted with.

To begin, let U2 be initialized and then fixed. Now the problem is to find the optimal U1 such that (8) is maximized. When U2 is fixed, problem (8) can be formulated as follows:

where (U2)i ∈ Rn2×1 denotes the i-th column in U2. Problem (9) can be formulated as follows:

where p = n × r2, zk = Xj(U2)i and zk ∈ Rn1×1 . Problem (10) is the same to (6), and thus can be easily solved.

After the optimal U1 is obtained, we view U1 as the known term, and the task becomes the problem of computing the optimal U2 that maximizes (8):

where (U1)i ∈ Rn2×1 denotes the i-th column in U1. Problem (11) can be formulated as follows:

where p = n×r1, zk = (U1)Ti Xj , and zk ∈ R1×n2 . Problem (12) can be formulated as follows:

which is the same to problem(6).

We use the relative reduction of the Root Mean Square Reconstruction Error with L1-norm (RMSREL1) to check the convergence. The RMSREL1 is defined as follows:

where Yj = U T 1 XjU2, and Yj ∈ Rr1×r2 . RMSREL1 is a modified version of RMSRE [Ye, 2005], replacing the Frobenius norm with the L1-norm. The convergence is determined by checking whether the following inequality holds:

for some small threshold η > 0. RMSREL1(t) denotes the RMSREL1 value in the t-th iteration step. In our experiments, we choose η = 10-3 . Results in Section 4 show that the algorithm converges within two to three iterations.

After computing Yj ∈ Rr1×r2 , j = 1, · · · , n, we obtain a three-order tensor Y = [Y1, · · · , Yn] ∈ Rr1×r2×n . We use formula (5) to decompose Y ∈ Rr1×r2×n , where T = 3. Note PT is a continuous solution, in order to obtain the discrete indicator matrix Q′ , we do K-means across all the rows in PT ∈ Rn×K . In this way, we get the final discrete indicator matrix Q′ ∈ Rn×K . The whole algorithm for Robust Tensor Clustering is summarized in Algorithm 1.

4 Experiments

4.1 Datasets

We have conducted experiments on five real-world large datasets, including three face datasets (Yale, ORL, PIE),

Algorithm 1 Robust Tensor Clustering (RTC) algorithm.

• 1: Input: X = [X1, · · · ,Xn] ∈ Rn1×n2×n , where X is centralized, r1, r2, cluster number K;
• 2: Initialize U 1 2 ∈ Rn2×r2 such that (U 1 2 )T U 1 2 = I, t = 1;
• 3: while not converge do
• 4: Compute U t 1 using (10) while fixing U t 2 ;
• 5: Compute U t+1 2 using (13) while fixing U t 1 ;
• 6: t = t + 1;
• 7: end while
• 8: Computing Yj = U T 1 XjU2, j = 1, · · · , n to obtain Y = [Y1, · · · , Yn] ∈ Rr1×r2×n ;
• 9: Computing HOSVD of Y using formula (5), where T=3;
• 10: Doing K-means across all the rows in PT .
• 11: Output:Indicator matrix Q′ ∈ Rn×K

a handwritten digit dataset (USPS) and an object dataset (COIL20). The brief description of the datasets is listed as follows:

Yale The Yale face database 1 contains 165 gray scale images of 15 individuals, each individual has 11 images. Each image is reshaped into 64 × 64 pixels in our experiments.

ORL The ORL face database 2 consists of a total of 400 face images, of a total of 40 people (10 samples per person). Similarly, each image is reshaped into 64 × 64 pixels.

PIE The PIE 3 is a database of 41,368 images of 68 people, each person under 13 different poses. We randomly choose one near frontal poses (C05 in our experiment), and use all the images under different illuminations and expressions, thus we get 49 images for each individual and totally 3332 images. Each image is reshaped into 64 × 64 pixels.

USPS USPS 4 is an image dataset consisting of 9298 handwritten digits of "0" through "9". The size of each image is 16 ×16 pixels in our experiment.

COIL20 COIL20 5 contains 20 objects. The images of each objects were taken 5 degrees apart as the object is rotated on a turntable and each object has 72 images. The size of each image is 32x 32 pixels, with 256 grey levels per pixel.

All the features (pixel values) are scaled into [0,1] in our experiments and no other preprocessing step is applied.

4.2 EvaluationMetric

We use two metrics to measure the clustering performance: Accuracy (AC) and Normalized Mutual Information (NMI) [Cai et al., 2005]. Given an image xi , let ri and si be the obtained cluster label and the ground-truth label, respectively. The AC is defined as follows:

where n is the total number of images, δ(x, y) is the delta function that equals one if x = y and equals zero otherwise. map(ri) is the permutation mapping function that maps each cluster label ri to the equivalent ground-truth label. The best mapping can be found by using the Kuhn-Munkres algorithm.

1 http://cvc.yale.edu/projects/yalefaces/yalefaces.html

2 http://www.uk.research.att.com/facedatabase.html

3 http://vasc.ri.cmu.edu/idb/html/face/

4 http://www.csie.ntu.edu.tw/∼ cjlin/libsvmtools/datasets/

5 http://www.cs.columbia.edu/CAVE/software/softlib/

Let C denote the set of clusters obtained from the ground truth and C′ denote the set of clusters obtained from our algorithm. The Normalized Mutual Information (NMI) is defined:

where H(C) and H(C′ ) are the entropies of C and C′ , respectively. MI(C, C′ ) denotes the mutual information of C and C′ . It is easy to check that MI(C, C′ ) ranges from 0 to 1. MI = 1 if the two sets of clusters are identical, and MI = 0 if the two sets are independent.

4.3 Results

In order to evaluate the performance of proposed Robust Tensor Clustering algorithm, we compare it with five state-ofthe-art approaches: two traditional clustering methods: Kmeans [Bishop, 2006] and PCA+K-means [Bishop, 2006]; two spectral clustering methods: Self-tune SC (SSC) [Chen et al., 2011] and LSC-K [Chen and Cai, 2011]. We also compare the method in [Huang et al., 2008], and call it HOSVD clustering in our experiments.

All the methods need to input the cluster number K, we let K equal to the ground-truth number of classes in different datasets. In addition, SSC needs to input ˆK, which is the same to that in KNN. We tune ˆK in the rang of 5, 10, 15, 20 in our experiments. LSC-K needs to input the number of landmarks p, which is related to the size of datasets. In our experiments, we let p equal to 0.2, 0.4, 0.6, 0.8 of the dataset's size, respectively. Because LSC-K is robust to another parameter r, we let r = 6 according to [Chen and Cai, 2011]. In RTC, there are two parameters: r1 and r2 . For simplicity, we let r1 = r2 = r, and then tune r versus different values. The best performance is selected as the results. For K-means, PCA+K-means, SSC, and LSC-K, we reshape each image into a vector according to the size of images in different datasets, and then input the vectors to the clustering algorithms. For HOSVD and RTC, we directly use a matrix to represent an image, and then cluster these matrices. The dimensionality reduction parameter for HOSVD is the same as that in RTC.

Table 1 lists the comparison results on five datasets. We conduct 30 tests on each dataset and the best performance is reported. From the results in Table 1, we can draw the following three conclusions. First, PCA+K-means achieves better performance than K-means, which demonstrates the effectiveness of PCA reduction dimension. Second, RTC outperforms HOSVD on all five datasets. The improvement varies from 1.23% (USPS) to 12.73% (Yale), and averages at 5.744% for AC (the average is 4.6% for NMI). This demonstrates the effectiveness of our proposed robust tensor clustering algorithm, which takes HOSVD as its second step. Third, HOSVD achieves better performance than K-means on average. This demonstrates that the matrix-based method can exploit more spatial information than vector-based method, and thus obtain the better clustering performance. In addition, we note that spectral clustering methods (SSC and LSC-K) achieve better performance than RTC on USPS and COIL20 (even the PCA+K-means also shows slightly improvement than RTC on COIL20 dataset). This is because the size of images in both the two datasets is smaller compared with that on Yale and PIE. Moreover, the content of images in these two datasets is relatively simple. In this case, the spatial information within images isn't sufficient, and can't play a important role in clustering, leading to the worse performance of matrix-based method (RTC and HOSVD) than vector-based methods (PCA+K-means, SSC and LSC-K).

Now we discuss the influence of parameter r to the RTC algorithm. Figure 1 illustrates the parameter tuning results of RTC on five datasets. As expected, the performance of our method is refining when we increase the value of r, and reaches a relatively stable stage. This is as expected because the approximation can't retain sufficient information for the original tensor data when r is too small, with the result of poor clustering performance. Instead, the clustering performance of RTC doesn't improve with the increasing of r, which explains that only a little information within the original tensor data is useful, and thus we can use a low-rank approximation to replace the original tensor data to do clustering. In addition, we can see that the curves of parameter tuning on Yale and ORL are less smooth compared with that on USPS, PIE and COIL20. The reason for this may be that USPS, PIE and COIL20 have more images than that on Yale and ORL.

In order to evaluate the robustness of RTC to outliers, we add artificial occlusions on some randomly selected images of USPS dataset to generate simulative polluted images (outliers) like that in [Kwak, 2008] and [Pang et al., 2010]. Then we test the performance of different methods on this polluted dataset. 10 sample images with artificial occlusions is listed in Figure 3. We generate polluted images in terms of different ratios. The polluted ratio is defined as the number of polluted images divided by the number of all images in dataset. Figure 2 illustrates the AC and NMI comparison results for SSC, LSC-K and RTC in terms of different polluted ratios. To better show the robustness of RTC, we individually list the comparisons between RTC and HOSVD (see the right figure in Figure 2). From Figure 2, we can see that RTC gradually outperforms the SSC and LSC-K with the increasing of polluted ratio (0.5 for AC and 0.6 for NMI), which demonstrates that RTC can obtain better clustering performance when the dataset is met with a high polluted ratio. Further, we add the artificial occlusions on other four datasets like that on USPS (the polluted ratio is 0.1 for Yale, 0.5 for ORL, 0.1 for PIE, 0.6 for USPS, and 1 for COIL20, respectively). The results on these datasets in listed in Table 2. We can see that RTC obtains the best clustering performance on all datasets, which effectively demonstrates RTC is more robust to the outliers than other clustering algorithms.

We also evaluate the convergence of Algorithm 1, the experimental results are drawn in Figure 4. We can see that RTC converges within two to three iterations on all five datasets.

5 Conclusions

In this paper, we proposed an algorithm of Robust Tensor Clustering (RTC). Experimental results show that RTC usually gets better performance when the dataset meets a polluted ratio, which effectively demonstrates that RTC is more robust to the outliers than other clustering algorithms. Moreover, as a matrix-based method, RTC can exploit more spatial information of images, and thus has better performance when the images reserve sufficient spatial information. In addition, we utilize a non-greedy maximization algorithm to solve the RTC problem. Experiments show that our algorithm can effectively converges within two to three iterations.

Acknowledgments

We would like to thank prof. Zhi-Hua Zhou for his insightful comments on this work. This work is supported by National High-tech R&D Program of China (2013AA01A601), 100 Talents Programme of The Chinese Academy of Sciences, the NSFC (under Grant 61202166), and Doctoral Fund of Ministry of Education of China (under Grant 20120032120042).