A New Simplex Sparse Learning Model to Measure Data Similarity for Clustering

1 Background andMotivation

Graph theory is an important branch in mathematics and has many applications. There are many algorithms built on graphs: (1) Clustering algorithms [Ng et al., 2002; Hagen and Kahng, 1992; Shi and Malik, 2000; Nie et al., 2014], (2) Dimension reduction algorithms based on manifold, such as LLE [Roweis and Saul, 2000] and Isomap [Tenenbaum et al., 2000], (3) Semi-supervised learning algorithms, such as label propagation [Zhu et al., 2003; Wang et al., 2014], (4) Ranking algorithms, such as Page-Rank [Page et al., 1999] and Hyperlink-Induced Topic Search (HITS) [Kleinberg, 1999]. Many data sets also have a natural graph structure, web-pages and hyper-link structure, protein interaction networks and social networks, just to name a few. In this paper, we limit our discussion to the data sets that can be transformed to similarity graph by simple means.

∗Corresponding Author. This work was partially supported by NSF IIS-1117965, IIS-1302675, IIS-1344152, DBI-1356628.

The key to many applications mentioned lies in the appropriate construction of the similarity graphs. There are various ways to set up such graphs, which depend on both the application and data sets. However, there are still a few open issues: (i) Selecting the appropriate scale of analysis, (ii) Selecting the appropriate number of neighbors, (iii) Handling multi-scale data, and, (iv) Dealing with noise and outliers. There are already a few papers [Long et al., 2006; Li et al., 2007] solving one of these issues. In particular, the classic paper self-tuning spectral clustering [Zelnik-Manor and Perona, 2004] by L. Zelnik-Manor et al. solves issues (i) and (iii). However, there is no single method that solves all of these issues to the best of our knowledge.

In the remainder part of this section, we first provide a brief review about commonly used similarity graph construction methods. Next, we present a concise introduction to sparse representation so as to lay out the background of our method. In Section 2, we present our motivation and formulate our Simplex Sparse Representation (SSR) objective function Eq. (10). We propose a novel accelerated projected gradient algorithm to optimize the objective function in an efficient way. Our experiment part consists of two sections, one with a synthetic data set to highlight the difference between our method and previous method, the other with real data sets to demonstrate the impressive performance of SSR. We conclude this paper with the conclusion and future work.

1.1 Different Similarity Graphs

There are many ways to transform a given set x1,. . . ,xn of data points with pairwise similarities Sij or pairwise distance Dij into a graph. The following are a few popular ones:

1. The ε-neighbor graph: We connect all points whose pairwise distances are smaller than ε.
2. k-nearest neighbor graph: The vertex xi is connected with xj if xj is among the k-nearest neighbors of xi.
3. The fully connected graph: All points are connected with each other as long as they have positive similarities. This construction is useful only if the similarity function models local neighborhoods. An example of such function is the Gaussian similarity function S(xi, xj) = exp ( -‖xi-xj‖2 2σ2 ), where σ controls the width of neighborhoods.

Note that all the above three methods suffer from some of the open issues mentioned. The ε-neighbor graph method needs to find the appropriate ε to establish the edge connections, the k-nearest neighbor graph uses the specified k to connect neighbors for each node, the σ in Gaussian similarity method clearly determines the overall shape of the similarity graph. In short, they are not parameter free. Also, ε is clearly prone to the influence of data scale, k and σ would also be dominated by any feature scale inconsistence among these data vectors. Therefore, a parameter free and robust similarity graph construction method is desired and this is our motivation of this paper.

1.2 Introduction to Sparse Representation

Suppose we have n data vectors with size d and arrange them into columns of a training sample matrix

Given a new observation y, sparse representation method computes a representation vector

to satisfy

To seek a sparse solution, it is natural to solve

,where pseudo norm ℓ0 counts the number of non-zero elements in β and λ0 > 0 is the controlling parameter.

Recent discovery in [Donoho, 2004; Cand`es and Tao, 2006] found that the sparse solution in Eq. (3) could be approximated by solving the ℓ1 minimization problem:

or the penalty version:

This ℓ1 problem can be solved in polynomial time by standard linear programming methods [Roweis and Saul, 2000].

In the literature, sparse representation framework is well known for its robustness to noise and outliers in image classification [Wright et al., 2008]. It is noting that in such scheme, there is no restriction on the scale consistence for the data vectors. In other words, sparse representation has the potential to address the scale inconsistence and outlier issue presented in the introduction.

2 SimilarityMatrix Construction with Simplex Representation

Assuming we have the data matrix X = [x1, x2, . . . , xn] ∈ Rd×n , where each sample is (converted into) a vector of dimension d. The most popular method to compute the similarity matrix S is using Gaussian kernel function

One major disadvantage of this method is that the hyperparameter σ in the Gaussian kernel function is very sensitive and is difficult to tune in practice.

The sparse representation proposed in [Roweis and Saul, 2000] can be applied here to compute the similarity matrix S. Specifically, the similarities αi ∈ Rn-1 between the i-th feature and other features are calculated by a sparse representation as follows:

where X-i = [x1, . . . , xi-1, xi+1, . . . , xn] ∈ Rd×(n-1) , i.e, the data matrix without column i.

Since the similarity matrix is usually non-negative, we can impose the non-negative constraint on problem (6) to minimize the sum of sparse representation residue error:

The parameter λ is still needed to be well tuned here. It is easy to note that Eq. (7) is the sum of n independent variables, as a result, we will limit the discussion to a single term in the following context. More importantly, when the data are shifted by constants t = [t1, · · · , td]T , i.e., xk = xk + t for any k, the computed similarities will be changed. To get the shift invariant similarities, we need the following equation:

which indicates αTi 1 = 1. Imposing the constraint αTi 1 = 1 on the problem (7), we have

Interestingly, the constraints in problem (9) makes the second term constant. So problem (9) becomes

The constraints in problem (10) is simplex, so we call this kind of sparse representation as simplex representation.

Note that the above constraints (ℓ1 ball constraints) indeed introduce sparse solution αi and have empirical success in various applications. In order to solve the Eq. (10) for high dimensional data, it is more appropriate to apply first-order methods, i.e, use function values and their (sub)gradient at each iteration. There are quite a few classic first-order methods, including gradient descent, subgradient descent, and Nesterovs optimal method [Nesterov, 2003]. In this paper, we use the accelerated projected gradient method to optimize Eq. (10). We will present the details of the optimization and the algorithm in next section.

3 Optimization Details and Algorithm

Our accelerate projected gradient method introduces an auxiliary variable to convert the objective equation into an easier one, meanwhile the auxiliary variable approximates and converges to the solution αi in an iterative manner. The following part of this section provides the details. For the convenience of notations, we use f(αi) to represent the objective function in Eq. (10) and C to represent the associated constraints. In other words, we need to solve min αi∈C f (αi).

Let us assume that the auxiliary variable is zi and we will implement the alternative optimization method. We first solve α0i via solving Eq. (9) with λ = 1 without considering the constraints, then initialize z0i = α0i to start the iterative process, here zt-1 i denotes zi at iteration t - 1. c1 = 1 is the initial value of Newton acceleration coefficient, which will be updated at each iteration.

Next, we start the iterative approximation process. At iteration t, we approximate αi via Taylor expansion up to second order:

With extra terms independent of αi, we could write the previous objective function into a more compact form as follows:

The critical step of the our algorithm is to solve the following proximal problem:

Here we use v to represent zt-1 i - 1L f ′ (zt-1 i ) for notation convenience. We will introduce an efficient approach to solve Eq. (13) in a separate subsection below.

Then, we update the acceleration coefficient as follows:

Last the auxiliary variable is updated with the formula

Note that Eq. (14) and Eq. (15) together accelerate our algorithm.

The number of iteration is updated with t = t + 1.

We repeat the Eqs. (11)-(15) until the algorithm converges according to our criteria. Note that our algorithm converges at the order of O(t2 ), here t is number of iterations. We have to omit the proof due to space constrain, interested readers may refer to [Nesterov, 1983; 2003; 2005; 2007].

Our whole algorithm is summarized as follows:

Here the convergence criteria is that the relative change of ‖αi‖, the frobenius norm of αi, is less than 10-4 . The empirical experiments conducted in this paper show that our algorithm converges fast and always ends within 30 iteration.

One subtle but important point we need to point out here that the similarity matrix S = [ ˆα1, . . . , ˆαn] is not necessarily

Input: Data Matrix X

Output: α1, . . . , αn

1. Initialize αi and zi for i from 1 to n while Convergence Criteria Not Satisfied do

1. Optimize αi with the induced equation 13 using the method introduced in the corresponding subsection.

1. Update coefficient ct+1 = 1+√1+4c2t

Algorithm 1: Accelerated Projected Gradient Method

symmetric. Here, ˆαi is the result of inserting coefficient 0 for its i-th of αi, which is of dimension n. The coefficient αij, the similarity between xi and xj for X-i, does not have to be equal to αji. Our approach to deal with this is to get the final similarity matrix

i.e, get the average of S and its transpose. This is a conventional practice and the correctness is easy to see, so we omit the proof here. After we get the similarity matrix, we could then follow the standard spectral clustering procedure, i.e, calculate the Laplacian matrix, get the corresponding spectral vectors, then calculate the clustering based on K-means clustering. The key difference between our proposed method and standard spectral clustering method is how the similarity graph is constructed. Conventional similarity graph generally is constructed with one of the methods specified in the introduction part and clearly is not parameter-free. Also, these methods clearly are prone to outlier and data scale inconsistence influence.

We would like to summarize the main advantages of our method Simplex Sparse Representation (SSR).

First It is parameter free. There is no parameter tuning for either data vector scale or the number of neighbors.

Second The simplex representation adopts the merits of sparse representation. Our method is therefore robust to scale inconsistence and outlier noise issues.

Third Our optimization algorithm is easy to implement and converges fast.

3.1 Optimization Algorithm to Eq. (13)

We write the Lagrangian function of problem (13) as

where γ is a scalar and λ is a Lagrangian coefficient vector, both of which are to be determined. Suppose the optimal solution to the proximal problem (13) is α∗ , the associate Lagrangian coefficients are γ ∗ and λ∗ . Then according to the KKT condition [Boyd and Vandenberghe, 2004], we have the following equations:

here α∗ ij is the j-th scalar element of vector α∗i . Eq. (18) can be written as α∗ ij - vj - γ ∗1 - λ∗j = 0. According to the constraint 1T α∗i = 1, we have γ ∗ = 1-1T v-1T λ∗ n . So α∗ = (v - 11T n v+ 1n 1- 1T λ∗ n 1) +λ∗ .

Denote ¯λ∗ = 1T λ∗ n and u = v- 11T n v+ 1n 1, then we can write α∗ = u + λ∗ - ¯λ∗1. So ∀j we have

According to Eq. (19)-(22) we know uj + λ∗j - ¯λ∗ = (uj - ¯λ∗ )+, here x+ = max(x, 0). Then we have

So we can obtain the optimal solution α∗ if we know ¯λ∗ .

We write Eq.(22) as λ∗j = α∗ ij +¯λ∗ -uj. Similarly, according to Eqs. (19)-(21) we know λ∗j = (¯λ∗ -uj)+. As v is a n-

1-dimensional vector, we have ¯λ∗ = 1 n-1 n-1 ∑ j=1 (¯λ∗ - uj )+.

Defining a function as

so f(¯λ∗ ) = 0 and we can solve the root finding problem to obtain ¯λ∗ .

Note that λ∗ ≥ 0, f ′ (¯λ) ≤ 0 and f ′ (¯λ) is a piecewise linear and convex function, we can use Newton method to find the root of f(¯λ) = 0 efficiently, i.e,

4 Experiments on Synthetic Data Sets

In this part, we would like to compare our proposed method with Self-Tuning Spectral Clustering (STSC) proposed in [Zelnik-Manor and Perona, 2004]. STSC is also a parameter free spectral clustering method and relative new. We will first give a brief overview about STSC and then design a synthetic data experiment to highlight its difference between our method. STSC uses local scale parameters carefully and avoids using an uniform scale parameter, therefore the outlier influence is reduced. On the other hand, SSR focuses on representing every data vector using other vectors directly, therefore we don't worry about how to measure the similarity between vectors using appropriate metrics.

Next, we conduct a synthetic data experiment to demonstrate the difference in similarity graph construction between STSC and our algorithm SSR.

Input: Data Matrix X ∈ Rd×n

Output: α1, . . . , αn

1. Compute the local scale σi for each xi using the formula

where xK is xi's K-th neighbor

2. Construct the affinity matrix and normalized

Laplacian matrix L. 3. Find v1, . . . , vC the C largest eigenvectors of L and form the matrix

V = [v1, . . . , vC] ∈ Rn×C , where C is the largest possible group number

1. Recover the rotation R which best aligns V's columns with the canonical coordinate system using the incremental gradient descent scheme
2. Grade the cost of the alignment for each group

number, up to C, according to J = n∑ i=1 C∑ j=1 Z2 ij M2 i , where

Z = VR

6. Set the final group number Cbest to be the largest group number with minimal alignment cost.

7. Take the alignment result Z of the top Cbest eigenvectors and assign the original point xi to cluster c if and only if max j (Z 2 ij ) = Z 2 ic

Algorithm 2: Self-Tuning Spectral Clustering

We construct 3 groups of normal distributed vectors of dimension 10. First group consists of 10 vectors and with mean value 1, second, third group both contain 10 vectors but with mean 2 and 3 respectively. We try to do the clustering for these 3 groups vectors. To make the problem non-trivial, we also add noise with mean 0 and standard deviation 0.5 to the total 30 vectors. Such experiments are repeated 100 times. We now would like to introduce clustering measure metrics to evaluate clustering performance, which would be used to evaluate synthetic data set mentioned and real data sets in the next part.

4.1 EvaluationMetrics

To evaluate the clustering results, we adopt the three widely used clustering performance measures which are defined below.

Clustering Accuracy discovers the one-to-one relationship between clusters and classes and measures the extent to which each cluster contained data points from the corresponding class. Clustering Accuracy is defined as follows:

where ri denotes the cluster label of xi and li denotes the true class label, n is the total number of samples, δ(x, y) is the delta function that equals one if x = y and equals zero otherwise, and map(ri) is the permutation mapping function that maps each cluster label ri to the equivalent label from the data set.

NormalizedMutual Information(NMI) is used for determining the quality of clusters. Given a clustering result, the NMI is estimated by

where ni denotes the number of data contained in the cluster Ci(1 ≤ i ≤ c), ˆnj is the number of data belonging to the Lj(1 ≤ j ≤ c), and ni,j denotes the number of data that are in the intersection between cluster Ci and the class Lj .

Purity measures the extent to which each cluster contained data points from primarily one class. The purity of a clustering is observed by the weighted sum of individual cluster purity values, given as follows:

where Si is a particular cluster size of ni, nji is the number of the i-th input class that was assigned to the j-th cluster. K is the number of the clusters and n is the total number of the data points.

4.2 Synthetic Data Experiment Results

We first plot the similarity graph coefficients for the two methods in Fig. (1). Since the random vectors vary a lot during different times of experiments, it makes no sense plotting the average coefficient. As a result, we have to demonstrate the plot for a typical experiment, however, the plot does not vary noticeable according to our empirical experiments. It can be observed that our method yields a much sparse similarity graph.

Next, we summarize the average clustering performance of these two methods on this synthetic data set in Table 1. It can be observed that both methods achieve very good clustering results due to the clear structure of data matrix. However, SSR slightly outperforms STSC in terms of all three measure metrics.

5 Experiments on Real Data Sets

In this section, we will evaluate the performance of the proposed method on benchmark data sets. We compare the our method with K-means, NMF [Lee and Seung, 2001], Normalized Cut (NCut) [Shi and Malik, 2000] and SelfTuning Spectral Clustering (STSC) [Zelnik-Manor and Perona, 2004]. K-means and NMF are classic clustering methods and perfect for benchmark purpose. We also include NC and STSC here since their close connection with our proposed method.

Parameter setting is relatively easy for these methods. STSC and our method are both parameter free. For K-means and NMF, we set the number of clusters equal to the number of classes for each data set. For Ncut, we construct the similarity matrix via tuning the scale parameter σ from the list {10-2 , 10-1 , 1, 101 , 102 }. Under such setting of each method, we repeat clustering 30 times and report the average for each method.

There are in total 9 data sets used in our experiment section. Among those, 6 are image ones, AR 1 [Martinez and Kak, 2001], AT&T 2 [Samaria and Harter, 1994], JAFFE [Lyons et al., 1998], subset of MNIST [LeCun et al., ], subset of PIE 3 [Sim et al., 2002], subset of UMIST. The other 3 non-image ones are from UCI machine learning repository [Frank and Asuncion, 2010]: Abalone, Ecoli and Scale.

Table 1 summarizes the characteristics of the data sets used in the experiments.

5.1 Experiment Results and Discussions

Table 3 summarizes the results for all methods on the benchmark data sets specified above. We bold the corresponding result if it is statistically significant better than results from other methods via T-test. It can be observed that our proposed method outperforms other methods on all these data sets with the specified metrics. In particular, our method significantly gets a better result than Ncut (standard spectral clustering) and STSC (parameter-free spectral clustering) on some image data sets. AR and PIE, which both contain a large number of classes and noisy samples, are widely used benchmark data sets for sparse representation related classification methods. Note that Ncut could get better result if we tune the scale parameter in a refined way. However, as we mentioned in the introduction, choosing an appropriate scale parameter is difficult when the ground truth is unknown. Choosing an uniform scale parameter is also prone to outlier influence. The results demonstrate that our method provides the potential solution to the open issue in spectral clustering.

1 http://www2.ece.ohio-state.edu/ aleix/ARdatabase.html, other downloads 1.

2 http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase. html

3 http://www.zjucadcg.cn/dengcai/Data/data.html,we use first 10 images in each class

5.2 Computational Complexity and Scalability Discussion

Our method provides a parameter-free way to construct the similarity graph for subsequent clustering task. Its computational complexity is comparable to conventional spectral clustering methods, indeed significantly faster. For each data vector, the sparse coefficient vector will converge in the quadratic order in each iteration due to integrated Newton acceleration idea. The algorithm usually converges after a limited number of iterations. Our method runs efficiently on individual machines.

It is noting that our method does not have scalability issue on distributed computing system either. In addition to its impressive performance on individual machine, our method can be easily extended to distributed computing platform. With Apache Spark, we can hash all data vectors in the memory, compute the individual coefficient vectors on work nodes, and transmit the αs back to name node. This is due to the fact the sparse vector α for each individual vector can be easily paralleled.

6 Conclusion and FutureWork

In this paper, we proposed a parameter free spectral clustering method that is robust to data noise and scale inconsistence. Our framework addresses the potential solution to many open issues of spectral clustering. Our simplex representation objective function is derived in a natural way and solved via a novel algorithm. The projected gradient method was accelerated via a combination of auxiliary variable and Newton root finding algorithm. Empirical experiments on both synthetic and real data sets demonstrate the effectiveness of our method.

In the future, we have two further research directions. First, we will try to encode label information in our graph

(b) NMI

(c) Purity

construction, and therefore extend our framework to semisupervised case. In the literature, there have been extensive work in semi-supervised graph learning [Kulis et al., 2005; Huang et al., 2013b]. In terms of our framework, we may consider using ℓ2,1 norm in Eq. (7) to take advantage of structural sparsity instead of flat sparsity. Second, we may further look into how to improve the clustering performance after we get the similarity graph, prior work include [Huang et al., 2013a].