CS Sparse K-means: An Algorithm for Cluster-Specific Feature Selection in High-Dimensional Clustering
Xiangrui Zeng Hongyu Zheng


A Hierarchical Approach for Visual Storytelling Using Image Description
    Md Sultan Al Nahian Tasmia Tasrin Sagar Gandhi Ryan Gaines Brent Harrison
    
=============================================================================




Feature selection is an important and challenging task in high dimensional clustering. For example, in  genomics, there may only be a small number of genes that are differentially expressed, which are informative to the overall clustering structure.  Existing feature selection methods, such as Sparse K-means, rarely tackle the problem of accounting features that can only separate a subset of clusters. In genomics, it is highly likely that a gene can only define one subtype against all the other subtypes or distinguish a pair of subtypes but not others. In this paper, we propose a K-means based clustering algorithm that discovers informative features as well as which cluster pairs are separable by each selected features. The method is essentially an EM algorithm, in which we introduce lasso-type constraints on each cluster pair in the M step, and make the E step possible by maximizing the raw cross-cluster distance instead of minimizing the intra-cluster distance. The results were demonstrated on simulated data and a leukemia gene expression dataset.  









§ INTRODUCTION


Dimensionality reduction has been a perpetual
problem in machine learning and computational biology, especially in computational genomics when a typical set of gene expression data has a huge number of genes. Typical dimensionality reduction methods like PCA would work on such cases, resulting in a low dimensional space embedding, but are less biologically meaningful, since each feature vector is a weight map over all genes, so all genes contribute to each feature vector. 

In the field of computational genomics, for both classification and clustering problems, there usually exists a small set of so-called "signature genes" that are actually reverent to the problem. And that is why the more common approach for dimensionality reduction in the field is to perform feature selection, that is, to select a number of features (genes in our case), and perform classification or clustering only with the inferred set of genes (or in a more flexible manner, weight important genes higher and give zero weight to genes that are similarly expressed between clusters). For applications of feature selection in computational biology, see <cit.>.

In this study, we focus on clustering problems. Here, the cluster structure can represent subtypes of samples, and as described before, data points (samples) might only differ in a small subset of features (genes), which are important in the differentiation of the subtypes. Finding the important genes that distinguish subtypes also helps the interpretation of the disease states because selected genes can be used for further pathway and network analysis. 

Sparse K-means <cit.> is an algorithm that performs two tasks: it assigns a weight to each feature (gene) between 0 and 1, and clusters on high-dimensional data. These two steps are dependent on each other; the weight is calculated with clustering results, and the clustering is based on feature weights, and an Expectation-Maximization couples the tasks and results in both a good clustering that takes feature sparsity into account, and a good weight assignment that is based on the clustering. An illustration of Sparse K-means can be found in Figure <ref>.





There is one shortcoming of the aforementioned method. It assigns an importance score (weight) to a feature across all clusters, which is intuitive since in feature selection frameworks, the importance score is designed to be the importance of a feature to the whole cluster structure. Nevertheless, it is common for a gene to contribute to a particular subtype but not others, and the Sparse K-means does not take this into consideration. Ability to infer the importance for each feature between clusters can also be beneficial in practice, since knowing the subtype-specific gene importance helps the interpretation and diagnostics of different subtypes.

Therefore, We propose CS sparse K-means (CS stands for Cluster-Specific), which clusters the samples using not only the adaptively selected features but also their contribution to each pair of clusters. The idea of CS sparse K-means is illustrated in Figure <ref>. The weight assignment of CS sparse K-means is shown in Table <ref>.













§ RELATED WORK


Many feature selection methods have been developed <cit.>. In general, there are two ways to perform feature selection, namely filter approach and wrapper approach; The former approach calculates some relevance score on features and discards those with low scores. The latter approach uses a predictor as a black box and aims to obtain a subset of features that maximize the predictor score over the said subset. Wrapper methods are more common in practice due to more control over the selection process and direct association with final results, and in this paper we will mostly focus on such algorithms.

Our proposed method can be seen as an extension of classical Sparse K-Means algorithm, but it's worth noting there are other works that also focus on the same problem, but with a much different approach. An earlier attempt comes from the paper titled "Pairwise Variable Selection for High-Dimensional Model-Based Clustering <cit.>". They introduce the following regularization term to the traditional Gaussian Mixture models(k over all features, i < j over all cluster pairs, λ is a tunable hyperparameter):



    λ∑_k∑_i < j |μ_ik - μ_jk|


Since L1-norm enforces sparsity, if for a feature k and two clusters i < j we have μ_ik = μ_jk, then this feature is considered not useful in separating cluster i and j.

The authors develop an EM algorithm to iteratively update means and probabilistic assignments.

We note that the method is simple and intuitive, but lacks intuitive control over desired sparsity and only gives binary results(the feature is useful or is not useful instead of assigning a value for the situation).

The paper titled "CRAFT: ClusteR-specific Assorted Feature selecTion" <cit.> uses Bayesian method for this problem. The three main parts of their models are:

1. Whether a feature should be selected for a cluster is indicated by binary variable v_kd which follows a Bernoulli distribution and the parameter for this Bernoulli is drawn from a Gaussian with predefined mean and variance.

2. The cluster assignments and the number of clusters are generated by a Dirichlet Process.


3. If a feature is selected for a cluster, its distribution within the cluster follows N(μ_kd, σ_kd), and if not, it follows N(μ_d, σ_d) (independent of clusters).

The authors then developed approximate inference algorithms for the parameters.

However, the process has a lot of parameters and requires several regularizations to work, while our methods are much simpler to interpret. Its upsides include being able to handle categorical features (not mentioned above) and inference over the number of clusters by virtue of Dirichlet Process.



§ METHODS




 §.§ Setup


Suppose we have a set of n samples and k features that take on real values (Denoted X∈ℝ^n× k), and a predefined number of clusters C. The classical clustering task asks for a partition z:[n]→ [C] that minimizes the inter-cluster distance:

    min∑_k∈[C]1/2|C_k|∑_i,j∈ C_k||X_i-X_j||^2


We are however inspired by the so-called un-normalized inter-cluster distance, which only differs by a scaling constant if the clusters are of similar size (C_i ≈ C_j):


    min1/2∑_k∈[C]∑_i,j∈ C_k||X_i-X_j||^2


For a fixed dataset, the quantity ∑_i,j ||X_i-X_j||^2 is the same regardless of the partitions. This means the above optimization problem is equivalent to the following, which we call un-normalized cross-cluster distance:


    max1/2∑_i,j, z(i)≠ z(j)||X_i-X_j||^2


We now introduce the parameters of our algorithm. w^(k)_i,j is a scalar that defines importance of feature k in separating cluster i and j. Given the above formulation, it's natural to scale the result based on the parameters we just introduced:


    max1/2∑_k∑_i,j, z(i)≠ z(j) w^(k)_z(i),z(j) (X_ik-X_jk)^2


And this will be the target function in our algorithm. Fix dataset X, we can write the function as F(w, z) = 1/2∑_k∑_i,j, z(i)≠ z(j) w^(k)_z(i),z(j) (X_ik-X_jk)^2; We next introduce constraints on w and design an E-M algorithm to perform clustering.



 §.§ Expectation-Maximization


In the E-step, we fix w and find an assignment z that will maximize F(w,z). As the title suggests, we are going to use K-Means at this step to get an assignment. However, since K-Means is not optimal, there is no formal guarantee(also as we mentioned above, K-Means was not designed for this objective, but they are similar if cluster sizes are close); Rather, we will later show that this strategy works well in practice.

First, note that we defined a distance metric between any two points a, b∈ℝ^k, only if their cluster assignment is given as z(a) and z(b):


    D(a, b) = ∑_k w^(k)_z(a), z(b) (a_k - b_k)^2


However, in K-Means, the only distance we will ever compute is the distance between a data point and current mean of a cluster(defined as μ_c∈ℝ^k for each cluster c), which (at the time of computation) are assumed to be in the same cluster, but we did not define D between two points in the same cluster. To circumvent this problem, we again take the idea of complementary: In original K-means, when assigning points, the rule is:


    z(i) = argmin_c ||X_i - μ_c||^2


This is actually equivalent to the following, since ∑_c ||X_i-μ_c||^2 is a constant:


    z(i) = argmax_c∑_c'≠ c ||X_i - μ_c'||^2


We can now naturally extend this idea for our case:


    z(i) = argmax_c∑_c'≠ c∑_k w^(k)_c, c' (X_ik-μ_c'k)^2


In the M-step, we fix assignment z and find a set of w that will maximize F(w,z). First of all, F is a linear function of w:


    F(w,z) = ∑_k ∑_a < b w^(k)_a, b (∑_i∈ C_a∑_j ∈ C_b (X_ik - X_jk)^2) = ∑_k, a<b w^(k)_a, b s^(k)_a, b


where s are independent of w and we know s≥ 0. Note that we can calculate s with sufficient statistics without iterating over all pairs in (C_a, C_b): s^(k)_a,b = |C_b|(∑_i∈ C_aX_ik^2) + |C_a|(∑_j∈ C_b X_jk^2) - 2(∑_i∈ C_aX_ik)(∑_j∈ C_b X_jk).

Without constraint on w we will just take it to infinity and it's not interesting. Inspired by Sparse K-Means, and Lasso in general, we introduce the constraint on w:


    ∑_k w^(k)_a, b≤ T,       ∀ a < b 
    ∑_k (w^(k)_a, b)^2 ≤ 1,       ∀ a < b 
    
    w ≥ 0


where T is a tunable parameter. Note that this can be broken down to c(c-1)/2 instances of following optimization problem:


    max        w^Tx 
    s.t.       ||w||_1 ≤ T 
           ||w||_2 ≤ 1 
           w ≥ 0


And the solution to this particular optimization problem is given by normalizing a parameterized soft-thresholding function:

    w     =  S(x, δ)/||S(x, δ)||_2
    
    S(x, δ)_i     = max (x_i - δ, 0)


The parameter δ is taken to be minimum non-negative value that satisfies ||w||_1 ≤ T(note that ||w||_2 = 1 is guaranteed). There is no compact representation for the optimal δ, but binary search would suffice.

As in general E-M algorithm, we alternative between E-step (updating z) and M-step (updating w) until the function F(w,z) converges, or maximum iteration reached. The next section gives a formal description of the algorithm.



 §.§ Implementation

We give a full description of our purposed Cluster-Specific Sparse K-Means on next page. We provide pseudo code with specifications of input and output, and both E and M steps necessary to perform the optimization. 

We don't explicitly define F(w,z) at line 22, since in practice, the M-step function also returns sum of s· f(δ) which is F(w,z).

The algorithm has been implemented in Python 3. In practice, we set the convergence threshold to 0.001 as the target value is rather huge in most cases (analogous to total cross-cluster distance, so the value scales in a quadratic way with n). In the binary search process, the starting interval is set to [0, max(s) - ϵ] to ensure at least one entry in f(δ) is positive. 






§ RESULTS



 §.§ Analysis of a simulated dataset

We simulated a dataset of 60 samples and 100 features. The background followed N(0,1). A small signal of 1.2 was added to feature 1 to 10 of sample 1 to 20, feature 11 to 20 of sample 21 to 40,  and feature 21 to 30 of sample 41 to 60. Therefore, we formed 3 clusters, each with 10 distinct features.



Sparse 3-means and CS sparse 3-means were performed with T = 5. Both algorithms obtained relatively high accuracy: Our algorithm misclassified one sample (accuracy 98.3%) while sparse K-means misclassified three (accuracy 95%). The assigned weight for each feature by sparse K-means is shown in heatmap (Figure <ref>).

However, our algorithm is able to generate cluster(pair) specific weights:



Therefore, CS sparse K-means achieved better accuracy as well as discovered informative features specific to each cluster pair.



 §.§ Analysis of a leukemia dataset with three subtypes

To demonstrate the efficacy of CS sparse K-means on a real-world dataset, we compared the performance of sparse 3-means and CS sparse 3-means on a leukemia gene expression dataset of 3 subtypes <cit.>.

The preprocessed dataset <cit.> consists 89 samples and 5135 genes. The 89 bone marrow aspirates or peripheral blood samples are of 3 subtypes: 33 cohorts with “inv(16)" abnormalities,  21 cohorts “t(15;17)" abnormalities, and 35 cohorts with “t(8;21)" abnormalities.

We applied sparse 3-means and CS sparse 3-means with T = 5 to this gene expression dataset. The sparse 3-means misclassified 4 samples (accuracy 95.5%): two “t(8;21)" abnormalities were misclassified as “inv(16)" abnormalities, one “t(8;21)" abnormality was misclassified as “t(15;17)" abnormality, and one `inv(16)" abnormality was misclassified as “t(8;21)" abnormality. The CS sparse 3-means misclassified 3 samples (accuracy 96.6%): two “t(8;21)" abnormalities were misclassified as “inv(16)" abnormalities, and one “t(8;21)" abnormality was misclassified as “t(15;17)" abnormality. Overall, CS sparse 3-means has higher accuracy.

Sparse 3-means selected 50 genes (with non-zero weights). CS sparse 3-means selected 99 genes (at least one cluster pair has non-zero weight). To assess the gene selection performance by both methods, we compared the assigned gene weights results with the signature genes reported in the original paper. The original paper reported one gene signature for “t(15;17)" abnormality: HGF, and ten gene signatures for “t(8;21)" abnormality: LAT2, C11orf21, C15orf39, CAPG, TRH, AIF1, HSPG2, PRAME, FAM698, and LCP. The 3 clusters were identified by their greatest overlapping true subtype labels.

The gene weight assignments by both methods were summarized in Table 2 and Table 3.





Out of the 11 signature genes, sparse 3-means discovered 4 of them out of its 50 selected genes (p < 2.5× 10^-6 by Fisher's exact test). CS sparse 3-means discovered 7 of them out if its 99 selected genes (p < 2.5× 10^-10 by Fisher's exact test). If we take into account of the correct cluster pair specific discovery by CS sparse 3-means, the p-value is even smaller (p < 9.0× 10^-13).

Therefore, on this leukemia gene expression dataset, CS sparse K-means achieved better clustering and gene selection accuracy than sparse K-means. Our algorithm also correctly discovered subtype specific informative genes which provide very important clues on the gene expression level alternation between specific disease subtypes. 















§ CONCLUSION


In this paper, we developed an algorithm for cluster-specific sparse clustering based on K-means. Our algorithm showed better clustering accuracy than sparse K-means on several simulated or real-world datasets. In genomics analysis, our algorithm can discover disease subtype specific informative genes by assigning weights on each gene with respect to each subtype pairs. 

In the future, we can improve the performance of our algorithm in several aspects: 1) Extend the process to Gaussian Mixtures; 2) Formally analyze convergence and correctness guarantees for this algorithm; 3) Introduce parallelism and early-stopping to speed up the process. We will also apply this algorithm to more real-world genomics dataset to assist disease subtype studies and to further develop our algorithm.


plain



