Deterministic Feature Selection for $k$-means Clustering

Introduction

This paper is about feature selection for $k$-means clustering, a topic that received considerable attention from scientists and engineers. Arguably, $k$-means is the most widely used clustering algorithm in practice . Its simplicity and effectiveness are remarkable among all the available methods . On the negative side, using $k$-means to cluster high dimensional data with, for example, billions of features is not simple and straightforward ; the curse of dimensionality makes the algorithm very slow. On top of that, noisy features often lead to overfitting, another undesirable effect. Therefore, reducing the dimensionality of the data by selecting a subset of the features, i.e. feature selection, and optimizing the $k$-means objective on the low dimensional representation of the high dimensional data is an attractive approach that not only will make $k$-means faster, but also more robust .

The natural concern with throwing away potentially useful dimensions is that it could lead to a significantly higher clustering error. So, one has to select the features carefully to ensure that one can recover comparably good clusters just by using the dimension-reduced data. Practitioners have developed numerous feature selection methods that work well empirically . The main focus of this work is on algorithms for feature selection with provable guarantees. Recently, Boutsidis et al. described a feature selection algorithm that gives a theoretical guarantee on the quality of the clusters that are produced after reducing the dimension . Their algorithm, which employs a technique of Rudelson and Vershynin , selects the features randomly with probabilities that are computed via the right singular vectors of the matrix containing the data ( describes a similar randomized algorithm with the same bound but faster running time). Although Boutsidis et al. give a strong theoretical bound for the quality of the resulting clusters (we will discuss this bound in detail later), the bound fails with some non negligible probability, due to the randomness in how they sample the features. This means that every time the feature selection is performed, the algorithm could (a) fail, and (b) return a different answer each time. To better address the applicability of such feature selection algorithms for $k$-means, there is a need for deterministic, provably accurate feature selection algorithms. We present the first deterministic algorithms of this type. Our contributions can roughly be summarized as follows.

Deterministic Supervised Feature Selection (Theorem 2). Given any dataset and any input $k$-partition of this dataset, there is a small set of $O(k)$ feature dimensions in which any near optimal output $k$-partition of the data is no more than a constant factor worse in quality than the given input $k$-partition (the quality of the input and output $k$-partitions are compared in the original dimension). Moreover, this small set of feature dimensions can be computed in deterministic low-order polynomial time. This is the first deterministic algorithm of this type. Prior work gives a randomized algorithm that can only reduce to $\Omega(k\log k)$ dimensions and guarantee comparable clustering quality.

Existence of a small set of near optimal features (Corollary 3). We prove existence of a small set of near optimal features. That is, given any dataset and the number of clusters $k$, there is a set of $O(k)$ feature dimensions such that any optimal $k$-partition in the reduced dimension is no more than a constant factor worse in quality than the optimal $k$-partition of the dataset. The existence of such a small set of features was not known before. Prior work only implies the existence of $\Omega(k\log k)$ features with comparable performance.

Deterministic Unsupervised Feature Selection (Theorem 4). Given any dataset and the number of clusters $k$, it is possible, in deterministic low-order polynomial time, to select $r$ feature dimensions, for any $r>k$, such that the optimal $k$-partition of the dimension-reduced data is no more than $O(n/r)$ worse in quality than the optimal $k$-partition; here $n$ is the number of features of the dataset. This is the first deterministic algorithm of this type. Prior work offers a randomized algorithm with error $\left(3+O(k\log k/r)\right)$, for $r=\Omega(k\log k)$.

Unsupervised Feature Selection with a small subset of features (Theorem 5). Finally, given any dataset and the number of clusters $k$, in randomized low-order polynomial time it is possible to select a small number, $r$, of feature dimensions, with $k<r=o(k\log k)$, such that the optimal $k$-partition for the dimension-reduced data is no more than $O(k\log(k)/r)$ worse in quality than the optimal $k$-partition. In particular, this is the first (albeit randomized) algorithm of this type that can select a small subset of $O(k)$ feature dimensions and provide an $O(\log k)$-factor guarantee. Prior work is limited to selecting $r=\Omega(k\log k)$ features. The new algorithm combines ideas from this paper with the technique of .

In order to prove our results we prove a structural result in Lemmas 10 and 11. This structural result quantifies the tradeoffs when selecting feature dimensions in terms of how the feature selection process preserves the matrix norms of certain crucial matrices. This is a general structural result and may be of independent interest.

In order to get deterministic algorithms and be able to select $O(k)$ feature dimensions, we need to use techniques that are completely different from those used in . Our approach is inspired by a recent deterministic result for decompositions of the identity which was introduced in and subsequently extended in , while use the randomized technique of to extract the features. This general approach might also be of interest to the Machine Learning and Pattern Recognition communities, with potential applications to other problems involving subsampling, for example sparse PCA and matrix approximation .

We first provide the basic background on $k$-means clustering that is needed to describe our results in Section 2. We postpone the more technical background needed for describing our main algorithms and for proving our main results to Section 3. We begin with the definition of the $k$-means clustering problem.In Section 3, we provide an alternative definition using matrix notation, which will be useful in proving the main results of this work. Consider a set ${\cal P}$ of $m$ points in an $n$-dimensional Euclidian space,

and an integer $k$ denoting the desired number of clusters. The objective of $k$-means is to find a $k$-partition of $\mathcal{P}$ such that points that are “close” to each other belong to the same cluster and points that are “far” from each other belong to different clusters. A $k$-partition of $\mathcal{P}$ is a collection

of $k$ non-empty pairwise disjoint sets which covers ${\cal P}$. Let $s_{j}=|\mathcal{S}_{j}|,$ be the size of $\mathcal{S}_{j}$. For each set $S_{j}$, let $\bm{\mu}_{j}\in\R^{n}$ be its centroid (the mean point),

where $\bm{\mu}({\mathbf{p}}_{i})$ is the centroid of the cluster to which ${\mathbf{p}}_{i}$ belongs. The goal of $k$-means is to find a partition ${\cal S}$ which minimizes ${\cal F}$ for a given $\mathcal{P}$ and $k$. We will refer to any such optimal clustering as,

The goal of feature selection is to construct points

(for some $r<n$ specified in advance) by projecting each ${\mathbf{p}}_{i}$ onto $r$ of the coordinate dimensions. Consider the optimum $k$-means partition of the points in $\mathcal{\hat{P}}$,

The goal of feature selection is to construct a new set of points $\mathcal{\hat{P}}$ such that,

Here, $\alpha$ is the approximation factor and might depend on $m,n,k$ and $r$. In words, computing a partition $\hat{{\cal S}}_{opt}$ by using the low-dimensional data and plugging it back into the clustering metric ${\cal F}({\cal P},\cdot)$ for the high dimensional data, gives an $\alpha$-approximation to the optimal value of the $k$-means objective function. Notice that we measure the quality of $\hat{{\cal S}}_{opt}$ by evaluating the $k$-means objective function in the original space, an approach which is standard . Comparing $\hat{{\cal S}}_{opt}$ directly to ${\cal S}_{opt}$, i.e. the identity of the clusters, not just the clustering error, would be much more interesting but at the same time a much harder (combinatorial) problem. A feature selection algorithm is called unsupervised if it computes $\mathcal{\hat{P}}$ by only looking at $\mathcal{P}$ and $k$. Supervised algorithms construct $\mathcal{\hat{P}}$ with respect to a given partition ${\cal S}_{in}$ of the data. Finally, an algorithm will be a $\gamma$-approximation for $k$-means ($\gamma\geq 1$) if it finds a clustering ${\cal S}_{\gamma}$ with corresponding value ${\cal F}_{\gamma}\leq\gamma{\cal F}_{opt}.$ Such algorithms will be used to state our results in a more general way.

[k-means approximation algorithm] An algorithm is a “$\gamma$-approximation” for $k$-means clustering ($\gamma\geq 1$) if it takes as input the dataset $\mathcal{P}=\{{\mathbf{p}}_{1},{\mathbf{p}}_{2},...,{\mathbf{p}}_{m}\}\in\R^{m\times n}$ and the number of clusters $k$, and returns a clustering ${\cal S}_{\gamma}$ such that,

The simplest algorithm with $\gamma=1$, but exponential running time, would try all possible $k$-partitions and return the best. Another example of such an algorithm is in with $\gamma=1+\epsilon$ ($0<\epsilon<1$). The corresponding running time is $O(mn\cdot 2^{(k/\epsilon)^{O(1)}})$. Also, the work in describes a method with $\gamma=O(\log k)$ and running time $O(mnk)$. The latter two algorithms are randomized. For other $\gamma$-approximation algorithms, see as well as the discussion and the references in .

Statement of our main results

We first discuss guarantees for supervised feature selection when the user has a candidate input $k$-partition. We next present results for the unsupervised case. We end this section with a brief discussion of the results.

Our first result is within the context of supervised feature selection. Suppose that we are given points $\mathcal{P}$ and some $k$-partition ${\cal S}_{in}$. The goal is to find the features of the points in $\mathcal{P}$ from which we can find a partition ${\cal S}_{out}$ that is not much worse than the given partition ${\cal S}_{in}$. Notice that if ${\cal S}_{in}$ is arbitrarily bad, then ${\cal S}_{out}$ might be much better than ${\cal S}_{in}$; however, our bound does not capture how much better the resulting partition would be. It only guarantees that it will not be much worse than the given partition.

There is an $O(mn\min\{m,n\}+rk^{2}n)$ time deterministic feature selection algorithm which takes as input any set of points $\mathcal{P}=\{{\mathbf{p}}_{1},{\mathbf{p}}_{2},...,{\mathbf{p}}_{m}\}\in\R^{m\times n}$, a number of clusters $k$, a number of features to be selected $k<r<n$, and a $k$-partition of the points ${\cal S}_{in}$. The algorithm constructs $\mathcal{\hat{P}}=\{\hat{{\mathbf{p}}}_{1},\hat{{\mathbf{p}}}_{2},...,\hat{{\mathbf{p}}}_{m}\}\in\R^{m\times r}$ such that ${\cal S}_{out}$, an arbitrary $\gamma$-approximation on $\mathcal{\hat{P}}$, satisfies

Asymptotically, as $r\rightarrow\infty$, the approximation factor is $\gamma\cdot O\left(1+\sqrt{k/r}\right)$. Essentially the clustering ${\cal S}_{out}$ is at most a constant factor worse than the original clustering ${\cal S}_{in}$; but, ${\cal S}_{out}$ was computed using the lower dimensional data. This means we can compress the number of the feature dimensions without destroying a given specific clustering in the data. This is useful, for example, in privacy preserving applications where one seeks to release minimal information of the data without destroying much of the encoded information . Notice that the feature selection part of the theorem (i.e. the construction of $\mathcal{\hat{P}}$) is deterministic. The $\gamma$-approximation algorithm, which can be randomized, is only used to describe the clustering that can be obtained with the features returned by our deterministic feature selection algorithm (same comment applies to Theorem 4). The corresponding algorithm is presented as Algorithm 4 along with the proof of the theorem in Section 5. Prior to this result, the best, and in fact the only method with theoretical guarantees for this supervised setting is randomizedWe should note that describes the result for the unsupervised setting but it’s easy to verify that the same algorithm and bound apply to the supervised setting as well. and gives

Further, requires $r=\Omega(k\log k)$, otherwise the analysis breaks. We improve this bound by a factor of $O(\log k)$; also, we allow the user to select as few as $r=O(k)$ features.

A surprising existential result is a direct corollary of the above theorem by setting ${\cal S}_{in}={\cal S}_{opt}$ and using an optimal clustering algorithm on the reduced dimension data ($\gamma=1$).

For any set of points $\mathcal{P}=\{{\mathbf{p}}_{1},{\mathbf{p}}_{2},...,{\mathbf{p}}_{m}\}\in\R^{m\times n}$, integer $k$, and $\epsilon>0$, there is a set of $r=O(k/\epsilon^{2})$ features $\mathcal{\hat{P}}=\{\hat{{\mathbf{p}}}_{1},\hat{{\mathbf{p}}}_{2},...,\hat{{\mathbf{p}}}_{m}\}\in\R^{m\times r}$ such that, if

In words, for any dataset there exist a small subset of $O(k/\epsilon^{2})$ features such that the optimal clustering on these features is at most a $\left(5+\epsilon\right)$-factor worse than the optimal clustering on the original features. Unfortunately, finding these features without the knowledge of ${\cal S}_{opt}$ is not obvious.

2 Unsupervised Feature Selection

Our second result is within the context of unsupervised feature selection. In unsupervised feature selection, the goal is to obtain a $k$-partition that is as close as possible to the optimal $k$-partition of the high dimensional data. The theorem below shows that it is possible to reduce the dimension of any high-dimensional dataset by using a deterministic method and obtain some theoretical guarantees on the quality of the clusters.

There is an $O(mn\min\{m,n\}+rk^{2}n)$ time deterministic algorithm which takes as input any set of points $\mathcal{P}=\{{\mathbf{p}}_{1},{\mathbf{p}}_{2},...,{\mathbf{p}}_{m}\}\in\R^{m\times n}$, a number of clusters $k$, and a number of features to be selected $k<r<n$. The algorithm constructs $\mathcal{\hat{P}}=\{\hat{{\mathbf{p}}}_{1},\hat{{\mathbf{p}}}_{2},...,\hat{{\mathbf{p}}}_{m}\}\in\R^{m\times r}$ such that ${\cal S}_{out}$, an arbitrary $\gamma$-approximation on $\mathcal{\hat{P}}$, satisfies

Asymptotically, as $r\rightarrow\infty$ and $n/r\rightarrow\infty$, the approximation factor is $\gamma\cdot O\left(n/r\right)$. The corresponding algorithm is presented as Algorithm 5 along with the proof of the theorem in Section 5. Prior to this result, the best method for this task was given in , which replaces $O\left(n/r\right)$ with $\left(3+O\left(\sqrt{{k\log(k)}/{r}}\right)\right)$ and it is randomized, i.e. this bound is achieved only with a constant probability. Further, requires $r=\Omega(k\log k)$, so one cannot select $o\left(k\log k\right)$ features. Clearly, our bound is worse but we achieve it deterministically, and it applies to any $r>k$, even $r=k+1$.

Finally, it is possible to combine the algorithm of Theorem 4 with the randomized algorithm of and obtain the following result.

There is an $O\left(mnk+rk^{3}\log(k)+r\log r\right)$ time randomized algorithm which takes as input any set of points $\mathcal{P}=\{{\mathbf{p}}_{1},{\mathbf{p}}_{2},...,{\mathbf{p}}_{m}\}\in\R^{m\times n}$, a number of clusters $k$, and a number of features to be selected $k<r<4k\log k$. The algorithm constructs $\mathcal{\hat{P}}=\{\hat{{\mathbf{p}}}_{1},\hat{{\mathbf{p}}}_{2},...,\hat{{\mathbf{p}}}_{m}\}\in\R^{m\times r}$ such that ${\cal S}_{out}$, an arbitrary $\gamma$-approximation on $\mathcal{\hat{P}}$, satisfies with probability $0.4$,

Since $r<4k\log k$, the approximation factor is $\gamma\cdot O\left(k\log k/r\right)$; there is no result in the literature which provides an approximation guarantee for selecting $o(k\log k)$ features. Note that the theorem is stated with the assumption $r<4k\log k$, even though the corresponding algorithm (Algorithm 6) runs with any $r>k$. This is because the algorithm of Theorem 5 is essentially the main algorithm of when $r=\Omega(k\log k)$ (that is, we do not provide an improved result for $r=\Omega(k\log k)$). What we achieve is to break the barrier of having to select $r=\Omega\left(k\log k\right)$ features. So, Theorem 5, to the best of our knowledge, is the best available algorithm in the literature providing theoretical guarantees for unsupervised feature selection in $k$-means clustering using $O\left(k\right)$ features. The proof of Theorem 5 is given in Section 5. Comparing Theorem 5 with Theorem 4, we obtain a much better approximation bound at the cost of introducing randomization. We note here that the main algorithm of can be obtained as Algorithm 6 in Section 5 after ignoring the 4th and 5th steps and replacing the approximate SVD with the exact SVD in the first step (the algorithm in uses approximate SVD as we do in Algorithm 6).

3 Discussion

For unsupervised feature selection, the theoretical guarantee from the deterministic algorithm is weak, but we suspect that it may perform very well in practice on typical input data matrices. The empirical investigation of this algorithm would also be an interesting future direction.

All our algorithms require the top $k$ singular vectors of the data matrix (or an approximation to the top $k$ singular vectors). Then, the selection of the features is done by looking at the structure of these singular vectors and using the deterministic techniques of . We should note that takes a similar approach as far as computing the (approximate) singular vectors of the dataset; then employ the randomized technique of to extract the features.

The high level description of our results for unsupervised feature selection assumed that the output partition, obtained in the reduced-dimension space, is the optimal one: ${\cal S}_{out}=\hat{{\cal S}}_{opt}$. This is the case if we assume a $\gamma$-approximation algorithm with $\gamma=1$. Notice, though, that our theorems are actually more general and can accomodate any $\gamma\geq 1$.

Preliminaries

We now provide the necessary technical background that we will use in presenting our algorithms and proving our main theorems in Section 5.

2 Approximate Singular Value Decomposition

3 Linear Algebraic Definition of k𝑘k-means

4 Spectral Sparsification and Rudelson’s concentration Lemma

For a vector ${\mathbf{z}}$ and scalar $\delta>0$, define the function

At each iteration $\tau$, the algorithm selects $i_{\tau}$, $t_{\tau}>0$ for which

The running time of the algorithm is dominated by the search for an index $i_{\tau}$ satisfying

A structural bound for clustering error with feature selection

In this section, we develop two lemmas which are general structural bounds for the clustering error when using feature selection. These lemmas are the basis for all of our algorithms, and hence are crucial to proving our results, specifically Theorems 2, 4, and 5 (recall that there is no algorithm for Corollary 3 since that result is non-constructive).

(See Section 3.1 for the definition of the pseudo-inverse operator.) Then,

In the first 3 steps, we have used spectral submultiplicativity, and in the last step we have used the definition of the spectral norm of the pseudo-inverse. Combining all these bounds together (and using $\gamma\geq 1$):

Proofs of Main Theorems

This rank requirement follows from Lemma 7, because

2 Proof of Theorem 4

To bound the second term on the right hand side, we use spectral submultiplicativity to obtain

The result now follows by using the bounds from Lemma 8,

3 Proof of Theorem 5

We state the corresponding algorithm as Algorithm 6. To prove Theorem 5, we start with the general bound of Lemma 10; to apply the lemma, we need to satisfy the rank assumption, which will become clear shortly, during the course of the proof.

We can now apply the Markov’s inequality, to obtain that with probability at least $0.9$,

To bound the denominator of the second term, observe that

Thus, using (7) and (8) in (10), we get that with probability at least $0.9$,

By a simple application of Markov’s inequality and the union bound, both of the equations below hold with probability at least $0.6$,

Using (11),(12),(13) and (14) in (15), together with a union bound, we obtain that with probability at least $0.5$

The running time of Algorithm 6 is $O\left(mnk+rk^{3}\log(k)+r\log r\right)$, since it employs the approximate SVD of Lemma 6, the randomized technique of Lemma 9, and the deterministic technique of Lemma 7.

Related work

Feature selection has received considerable attention in the machine learning and pattern recognition communities. A large number of different techniques appeared in prior work, addressing feature selection within the context of both clustering and classification. Surveys include , as well as , which reports the results of the NIPS 2003 challenge in feature selection. Popular feature selection techniques include the Laplacian scores , the Fisher scores , or the constraint scores . None of these feature selection algorithms have theoretical guarantees on the performance of the clusters obtained using the dimension-reduced features.

The result most closely related to ours is the work in . This work provides a randomized algorithm which offers a theoretical guarantee. Specifically, for $r=\Omega(k\epsilon^{-2}\log k)$, it is possible to select $r$ features such that the optimal clustering in the reduced-dimension space is a $(3+\epsilon)$-approximation to the optimal clustering. Our result improves upon this in two ways. First, our algorithms are deterministic; second, by using our deterministic algorithms in combination with this randomized algorithm, we can select $r=O(k)$ features and obtain a competitive theoretical guarantee.

Next, we should mention that if one allows linear combinations of the features (feature extraction rather than feature selection), then there are algorithms that offer theoretical guarantees. First there is the SVD itself, which constructs $k$ (mixed) features for which the optimal clustering in this feature space is a $2$-approximation to the optimal clustering . It is possible to improve the efficiency of this SVD algorithm considerably by using the approximate SVD (as in Lemma 6) instead of the exact SVD to get nearly the same approximation guarantee with $k$ features. The exact statement of this improvement can be found in . Boutsidis et al. show how to select $O(k\epsilon^{-2})$ (mixed) features with random projections and also obtaining a $(2+\epsilon)$-guarantee. While these algorithms are interesting, they do not produce features that preserve the integrity of the original features. The focus of this work is on what one can achieve while preserving the original features.

A complementary line of research approaches the $k$-means problem by sub-sampling the points of the dataset; such a subset of points is called coreset. The idea here is to select a small subset of the points and by using only this subset obtain a partition for all the points that is as good as the partition that would have been obtained by using all the points. offer algorithms for $(1+\epsilon)$ approximate partitions. Note that we were able to give only constant factor approximations. For example, shows the existence of an $(1+\epsilon)$-approximate coreset of size $r=O(k^{3}/\epsilon^{n+1})$ ($n$ is the number of features). provides a coreset of size $r=poly(k,\epsilon^{-1})$. The techniques used for all these coresets are different from the techniques we used for feature selection; further, there is no clear way on how to use a coreset selection algorithm to select features. Moreover, the authors of use the coreset-based algorithm from to design a PTAS for k-means and other clustering problems. It would be interesting to understand whether the techniques from are useful for feature selection as well. In particular, it appears that there is potential to obtain relative error feature selection $k$-means algorithms by using such approaches.

Acknowledgements

Christos Boutsidis acknowledges the support from XDATA program of the Defense Advanced Research Projects Agency (DARPA), administered through Air Force Research Laboratory contract FA8750-12-C-0323.

References