The Sparse MinMax k-Means Algorithm for High-Dimensional Clustering

1 Introduction

Traditional clustering algorithms including k-means, kmedoids, and the hierarchical ones very often lose their effectiveness when the dataset contains significantly fewer data items compared to the dimensionality of the associated feature space [Witten and Tibshirani, 2010; Li et al., 2018; Chang et al., 2017; Jin and Wang, 2016]. Such curse of dimensionality for clustering, manifests in the following different forms [Pandove et al., 2018]: difficulty in global optimization of the clustering objective, distance concentration on Lp norms, the effect of irrelevant or noisy features, and extremely sparse data volume. In most scenarios, only a small proportion of the features can be assumed to be relevant for clustering [Witten and Tibshirani, 2010; Chang et al., 2017; Jin and Wang, 2016]. The goal of a clustering algorithm is to identify these features, avoid the negative influences of the noisy features, and thus, more accurately identify the underlying cluster structure [Witten and Tibshirani, 2010; Jin and Wang, 2016]. Conventional ways of handling the challenge of high-dimensional clustering include subspace clustering of various forms (CLIQUE, SUBCLUE, DBSCAN, DUSC etc.) [Kriegel et al., 2009], correlation clustering [Kriegel et al., 2008], bi-clustering [Pontes et al., 2015], and several dimensionality-reduction based approaches including Principal Component Analysis (PCA), independent component analysis, non-linear matrix factorization, isometric feature mapping, sparse feature weighting, and so on [Kriegel et al., 2009]. Among the competitive and recent approaches in this direction, the Influential Features PCA (IF-PCA) [Jin and Wang, 2016] is a spectral clustering method that selects features with a higher degree of relevance to the clustering task by using the largest Kolmogorov-Smirnov(KS) scores and then applies classical PCA to the post-selection data matrix. Sarkar and Ghosh [2019] recently suggested an approach to tackle high-dimensional clustering with a new data-driven measure of dissimilarity, referred by the authors as MADD (Mean of Absolute Differences of pairwise Distances) specifically tailored for the high-dimensional feature spaces.

∗Contact Author

The widely popular k-Means [Lloyd, 1982] algorithm suffers from a strong sensitivity to initialization. Its clustering solution is not at all robust against the initial seed points (candidate cluster centers) and thus, it often gets trapped in poor local minima [Pe˜na et al., 1999; Celebi et al., 2013]. To tackle this problem, Tzortzis and Likas [2014] suggested the MinMax k-means clustering algorithm, which begins with a random set of candidate cluster centers and attempts to minimize the maximum intra-cluster variance instead of the traditional sum. Minimizing the sum does not consider the relative differences among the cluster variances, but by minimizing the maximum intra-cluster variance, large variance clusters are avoided and high-quality solutions are produced.

The relaxed maximum variance objective of the MinMax k-Means clustering algorithm assigns weights to each cluster but it does not perform any weighting on the features. Sparse regularization can, therefore, be imposed to extend its clustering efficacy in the high-dimensional space. We first justify that MinMax k-Means objective can be reformulated into Witten and Tibshirani's [2010] sparse clustering framework. This framework offers a specific feature-weighting method, which optimizes the weighted cost objective function using a lassotype penalty (ℓ1-norm regularization), hence assigning exact zero weights to noisy features [Witten and Tibshirani, 2010]. In this work, we thus propose the novel Sparse MinMax kMeans algorithm. Our main contributions can be summarized in the following way:

• We justify that MinMax k-Means objective can be re-

formulated into Witten and Tibshirani's [2010] sparse clustering framework, which offers a specific featureweighting method and optimizes a weighted cost objective function by using a lasso type penalty.

• We propose the Sparse MinMax k-Means algorithm which maximizes a new weighted between-cluster sum of squares (BCSS) with the ℓ1-norm regularization to impose exact zero weights on the noisy features. We thus extend the advantages of the MinMax k-Means algorithm in the high-dimensional space.
• We compare the performance of our approach with other well-known high-dimensional clustering algorithms through extensive experiments on several real word datasets (especially, the high-dimensional gene microarray datasets).

2 Preliminaries

In this section, we briefly discuss the MinMax k-Means and Sparse k-Means objectives along with the notations leading to the formulation of our proposed method.

2.1 Notations

We consider X = (xij) ∈ Rn×p to be our data set in matrix format where xij represents the j th feature (column) of the ith observation (row). Here n and p denote the number of observations and the number of features respectively. We consider K clusters and the set of cluster centers c = (ckj) ∈ RK×p . ck represents the kth cluster center and Ck denotes the kth cluster. Vk denotes variance of the cluster k where the cluster variance is defined as the sum, and not the average, of the squared distances from the observations belonging to the cluster to its center. δik is a cluster indicator variable with δik = 1 if xi belongs to cluster Ck and 0 otherwise. εmax denotes maximum intra-cluster variance. εw represents the weighted formulation of sum of the intra-cluster variances and wk denotes the weight assigned to cluster Ck in the MinMax k-Means algorithm [Tzortzis and Likas, 2014].

2.2 The MinMax k-Means Objective

The MinMax k-Means algorithm minimizes the maximum intra-cluster variance εmax (1):

Minimizing εmax is a non-trivial optimization problem. Thus, Tzortzis and Likas [2014] came up with a relaxed maximum variance objective. A weighted formulation εw of the sum of intra-cluster variances was thus constructed as in (2), where a higher weight wk was given to clusters with high variance to follow the behavior induced by the maximum variance criterion.

The exponent α is a user defined constant. 1 To compensate for the formation of large clusters, a higher weight should be induced for a higher variance. Maximizing εw with respect to the weights provide a way to realize this. Thus, the min-max problem can be written as:

2.3 The Sparse Clustering Framework

Witten and Tibshirani [2010] reformulated k-Means and hierarchical clustering as an optimization problem in the following way:

where fj(Xj,Θ) is a function that involves only the jth feature of the data, and Θ is a model parameter that belongs to a set τ . They further defined a sparse clustering framework as a solution to the above problem which is as follows:

where s is a tuning parameter that determines the number of retained features for clustering and ‖ω‖1 ,‖ω‖2 are respectively the ℓ1 and ℓ2-norms of the weight vector ω. ωj determines the contribution of the jth feature to the objective function (5). The ℓ1-norm or Lasso penalty results in sparsity in different applications [Witten and Tibshirani, 2010; Li et al., 2018; Chang et al., 2017].

3 The Sparse MinMax k-Means Algorithm

3.1 Deriving the Formulation

We introduce a new variable zik, which is similar to the cluster indicator variable δik and can be defined as:

Here wk's are the cluster weights as defined in the MinMax k-Means algorithm [Tzortzis and Likas, 2014] in Section 2.2. Thus, the εw defined in (2) can be re-written as follows:

We now provide Lemma 1 to prove that the MinMax k-Means clustering model can also be reformulated in the framework (4). We use a method similar to the one used by Chang et al. [2017] in context to sparse fuzzy clustering.

1 The role of α is discussed in Section 3.3

Lemma 1. For the data matrix X,

where n′k = ∑ni=1 zα ik = nkwα k , nk is the number of observations in cluster k, and ck is the kth cluster center such that:

Since n′k = ∑ni′=1 zα i′ k , we have

Similarly, we have

Thus,

The left hand side of Lemma 1 is the objective function (2) that is optimized by the MinMax k-Means Algorithm [Tzortzis and Likas, 2014], while the right hand side evaluates the dissimilarity within a cluster, which can be referred to as the within-cluster sum of squares (WCSS) of MinMax k-Means. The right hand side of Lemma 1 can be further simplified as:

Note that the product zα ik zα i′ k is non-zero only when both xi and xi′ ∈ Ck and is zero otherwise. The WCSS measure for our Sparse MinMax k-Means model can be written as:

From the constraints of the objective function (2) of the MinMax k-Means algorithm [2014], we can deduce that wα k ≤ 1. Hence, we can write:

This implies that the WCSS of MinMax k-Means will always be less than or equal to the WCSS defined for classical kMeans. Therefore, we can model the new between-cluster sum of squares (BCSS) of MinMax k-Means as:

where BCSS(Θ) = (BCSS1, . . . ,BCSSp)⊤ with Θ = C = (C1, . . . , CK) and τ is a set of all possible partitions of the observations into K clusters.

Now we can rewrite the MinMax k-Means problem by using the Witten and Tibshirani's [2010] sparse clustering framework as shown below:

Note that in MinMax k-Means, we maximize εw with respect to the cluster weights (wk's), and hence we must minimize

∑pj=1 BCSSj with respect to cluster weights. BCSSj, j = 1, . . . , p is a function that only involves the jth feature. Thus, we may conclude that MinMax k-Means fits into the framework (4). According to Witten and Tibshirani's [2010] sparse clustering framework (5), the MinMax k-Means can be generalized to the following model:

We will call (12) as the Sparse MinMax k-Means model.

Denoting aj = BCSSj, a = (a1, . . . , ap)⊤, the objective function of (12) can be rewritten as ∑pj=1 ωjaj . We use the technique similar to the one used in [Witten and Tibshirani, 2010] to build an algorithm to solve the problem (12). We alternatively fix two of w, c, and ω and maximize the objective with respect to the other. The optimization problem that arises in the final step while maximizing the objective with respect to ω can be written as:

Following the Proposition stated in page 715 of [Witten and Tibshirani, 2010]), the convex problem bears a solution of the form (13) is ω = S(a+,△)/∥∥S(a+,△)∥∥2 , where x+ indicates the positive part of x and △ = 0, if that results in ‖ω‖1 < s; else, △ > 0 is taken to get‖ω‖1 = s. S is the softthresholding operator, defined as S(x, c) = sign(x)(|x|-c)+. The assumptions are the same as in [Witten and Tibshirani, 2010], i.e., there exists a unique maximal element of a, and 1 ≤ s ≤ √p.

3.2 The Algorithm

From our definitions of WCSS and BCSS in Section 3.1, we see that maximizing the BCSS is equivalent to minimizing the WCSS and vice-versa. This property can be used to simplify our algorithm (1). Steps 3 and 4 can be optimized together by performing MinMax k-means clustering on the data after scaling each feature j by √ωj, which is same as updating the data matrix with each (i, j) element = √ωjxij . In Step 3, the new cluster centres can be calculated in the following way.

In Step 4, the cluster weight constraints in (12) and the new cluster centres obtained in Step 3 are incorporated in the objective via a Lagrange multiplier and the derivatives with respect to wk are set to zero. Thus, the cluster weights are updated as:

Tzortzis and Likas [2014] observed that addition of a memory effect to the weights could be beneficial in terms of increasing the stability. Thus, for each iteration t,

The algorithm iterates through steps 3-5 until the stopping criterion

Algorithm 1 Sparse MinMax k-Means Algorithm

Input: Data matrix X and number of clusters k

Parameter: Tuning parameter s

Output: Clusters C1, C2, . . . , Ck

• 1: Initialize ω as ω1 = . . . = ωp = 1 √p .

• 2: while stopping criteria (17) is not satisfied do
• 3: Optimize (12) with respect to C1,C2, . . . ,CK, keeping w and ω fixed. That is,

Maximizing (10) is same as minimizing (9).

• 4: Optimize (12) with respect to w, keeping C1, C2, . . . , CK abd ω fixed. That is,

Minimizing (10) is same as maximizing (9).

• 5: Optimize (12) with respect to ω, keeping C1,C2, . . . ,CK and w1,w2, . . . ,wK fixed, which results in the optimization problem stated in (13) and can be solved using the Proposition stated in page 715 of [Witten and Tibshirani, 2010] to get ωnew .
• 6: end while
• 7: return the clusters given by C1,C2, . . . ,CK, the cluster weights by w1,w2, . . . , wk and the feature weights corresponding to this clustering given by ω1, ω2, . . . , ωp.

is satisfied, where the precision level ϵ was chosen as 10-4 as suggested in [Witten and Tibshirani, 2010].

3.3 Selection of the Tuning Parameters

With the increase (decrease) of the exponent α, the similarity among the weight values also reduces (enhances) pertaining to the enhancement (suppression) of the relative difference of the variances among the clusters. α is chosen by using a data driven approach to automatically adapt to the data set as in [Tzortzis and Likas, 2014]. It initiates with a small α (αinit), which is increased by αstep after each iteration, till the attainment of a maximum value α (αmax). Following a procedure similar to [Witten and Tibshirani, 2010], we employ a permutation technique and determine the gap statistic [Tibshirani et al., 2001] to select the value of s.

3.4 On Convergence of the Proposed Algorithm

An exact theoretical proof of the local convergence of the minmax objective function in our proposed algorithm is quite difficult and outside the scope of this paper. Nevertheless, we provide some empirical demonstrations of the convergence behavior by recording the variation of the weighted BCSS objective with iterations of the alternating optimization procedure. In all our experiments, the value of the weighted BCSS objective (12) is seen to converge within 10 iterations. Figure 1 empirically demonstrates the gradual convergence of the proposed algorithm (1) to a local stationary point on 5 selected row can be treated as statistically significant. Using the mean CERs we obtain a similar performance as the one obtained in Table 2.

To further analyze the performance of our algorithm, we use two more criteria (Retained Features (RF) and Dunn Index (DI)). DI [Dunn, 1974] denotes the ratio of the minimum distance between observations from two different clusters to the maximum intra-cluster dissimilarity. It can be mathematically defined as DI ≜ min1≤k≤l≤K δ(Gk, Gl)/ max1≤m≤K △m, where δ(Gk, Gl) is the inter-cluster distance between clusters Gk and Gl, and △m calculates the maximum distance between all items within cluster Gm. RF is the number of retained features in each of the algorithms or basically the number of non-zero feature weights (for sparse algorithms). We show the comparative results against IF-HCT-PCA and Sparse k-Means in Table 5. Here we do not compare with MADD as it is not one of the feature selection oriented clustering methods. Higher DI values would indicate better clustering and our method obtains the highest DI values among the three methods in 8 out of the 10 cases. The lowest RF values were, however, obtained in only 4 cases by our method although it has the best CER values in most of the cases (indicating a better recovery of the ground truth cluster structure).

In addition, we provide the average run-time (in seconds) over 20 independent runs for 4 main algorithms used in this study on the gene microarray datasets in Table 6. The IFHCT-PCA method is seen to have the highest run-time, while Sparse k-Means and our algorithm have a run-time of similar order (though, experimentally, the later takes a slightly higher execution time due to the extra cluster-weight update step). The MADD k-Means approach has the least run-time.

5 Conclusion

We proposed a Sparse MinMax k-Means approach to detect meaningful clusters in higher dimensional feature spaces. Our approach attempts to extend the reach of traditional k-Means clustering to the high-dimension and low sample-size (i.e. p ≫ n) situations. The experimental results obtained in Section 4 supports the improved performance of our approach over other approaches in general. Our algorithm obtained the best average rank in all the experiments among the other algorithms. In the cases where any other algorithm has produced a better CER, our scheme remains only marginally behind.

The proposed clustering scheme can be readily extended to sparse multi-view settings by using multiple kernels for handling more complicated vision data. It may also be useful to employ other alternative dissimilarity measures, especially the divergence-based ones (see for example [Chakraborty and Das, 2017]) in the proposed clustering framework.

Acknowledgements

This work was supported by the Institute of Information & communications Technology Planning & Evaluation (IITP) grant funded by the Korea government (MSIT) (2016-0-00564, Development of Intelligent Interaction Technology Based on Context Awareness and Human Intention Understanding).