Feature selection for k-means clustering stability: theoretical analysis and an algorithm

Abstract

Stability of a learning algorithm with respect to small input perturbations is an important property, as it implies that the derived models are robust with respect to the presence of noisy features and/or data sample fluctuations. The qualitative nature of the stability property enhardens the development of practical, stability optimizing, data mining algorithms as several issues naturally arise, such as: how “much” stability is enough, or how can stability be effectively associated with intrinsic data properties. In the context of this work we take into account these issues and explore the effect of stability maximization in the continuous (PCA-based) k-means clustering problem. Our analysis is based on both mathematical optimization and statistical arguments that complement each other and allow for the solid interpretation of the algorithm’s stability properties. Interestingly, we derive that stability maximization naturally introduces a tradeoff between cluster separation and variance, leading to the selection of features that have a high cluster separation index that is not artificially inflated by the features variance. The proposed algorithmic setup is based on a Sparse PCA approach, that selects the features that maximize stability in a greedy fashion. In our study, we also analyze several properties of Sparse PCA relevant to stability that promote Sparse PCA as a viable feature selection mechanism for clustering. The practical relevance of the proposed method is demonstrated in the context of cancer research, where we consider the problem of detecting potential tumor biomarkers using microarray gene expression data. The application of our method to a leukemia dataset shows that the tradeoff between cluster separation and variance leads to the selection of features corresponding to important biomarker genes. Some of them have relative low variance and are not detected without the direct optimization of stability in Sparse PCA based k-means. Apart from the qualitative evaluation, we have also verified our approach as a feature selection method for $k$-means clustering using four cancer research datasets. The quantitative empirical results illustrate the practical utility of our framework as a feature selection mechanism for clustering.

Appendix

Proof

(Proof of Lemma 1) Based on Ding and He (2004), the continuous solution for the instance clusters is derived by the $k-1$ dominant eigenvectors of matrix $X_{fc}^TX_{fc}$, where $X_{fc}$ is a feature-instance matrix with the rows (features) being centered. Since $X$ is double-centered the sum of rows and columns of $X$ will be equal to 0, i.e. $\sum _iX_{ij}=\sum _jX_{ij}=0$. Thus, the continuous solution of Spectral $k$-means (for instance clustering) will be derived by the $k-1$ dominant eigenvectors of matrix $X^TX$.

Analogously, the continuous solution for the feature clusters is derived by the $k-1$ dominant eigenvectors of matrix $X_{ic}X_{ic}^T$, where $X_{ic}$ is a feature-instance matrix with the columns (instances) being centered. Since $X$ is double-centered, we will have that the instance-centered matrix $X_{ic}$ will be equal to $X$, i.e. $X_{ic}=X$. Thus, the continuous cluster solution will be derived by the dominant eigenvectors of matrix $XX^T$.

Using basic linear algebra one can easily derive that the matrices $XX^T$ and $X^TX$ have exactly the same eigenvalues.

Thus $\lambda _{k-1}(XX^T)-\lambda _{k}(XX^T)=\lambda _{k-1}(X^TX) -\lambda _{k}(X^TX)$, and the stability of the relevant eigenspaces will be equivalent. $\square$

Proof

(Proof of Theorem 1) We will start by decomposing the components $\lambda _{1}(Cov)$ and $\mathbf{Trace}(Cov)$. For $\lambda _{1}(Cov)$ we have:

$$\begin{aligned} \lambda _{1}(Cov)&= \lambda _{1}(\mathbf{diag}(u)X_{fc}X_{fc}^T\mathbf{diag}(u))\\&= \lambda _{1}(X_{fc}^T\mathbf{diag}(u)X_{fc})\\&= \max \limits _{||v||=1}v^T(X_{fc}^T\mathbf{diag}(u)X_{fc})v\\&= \max \limits _{||v||=1}\sum \limits _{i=1}^mu_i(v^Tx_{fc}(i))^2 \end{aligned}$$

In the above derivations we have used the following, easily verifiable facts: $\lambda _1(AA^T)=\lambda _1(A^TA)$ and $\mathbf{diag}(u)=(\mathbf{diag}(u))^2$. For $\mathbf{Tr}(Cov)$ we have:

$$\begin{aligned} \mathbf{Tr}(Cov)&= \mathbf{Tr}(\mathbf{diag}(u)X_{fc}X_{fc}^T\mathbf{diag}(u))\\&= \sum \limits _{i=1}^m u_i x_{fc}(i)^Tx_{fc}(i) \end{aligned}$$

Based on the above we can write:

$$\begin{aligned} Obj_{(ows)}(I\cup \{m_s\}) =\max \limits _{||v||=1}\sum \limits _{i\in I\cup \{m_s\}}(v^Tx_{fc}(i))^2-\frac{1}{n}\sum \limits _{i\in I\cup \{m_s\}} x_{fc}(i)^Tx_{fc}(i) \end{aligned}$$

The vector $v$ that maximizes the expression $\max _{||v||=1}\sum _{i\in I\cup \{m_s\}}(v^Tx_{fc}(i))^2$ is the dominant eigenvector of matrix $Cov=\mathbf{diag}(u)X_{fc}X_{fc}^T\mathbf{diag}(u)$, where $u$ is defined by the feature-set $I\cup \{m_s\}$ (recall that this expression was derived from $\lambda _1(Cov)$). Thus, based on the fact that for all $v$ such that $||v||=1$ it holds that $\lambda _1(A)\ge v^TAv$, we can employ any other normalized vector and compute a lower bound to the expression. In order to derive a useful tight bound we employ here the eigenvector $v_{(I)}$ of matrix $Cov=\mathbf{diag}(u)X_{fc}X_{fc}^T\mathbf{diag}(u)$, where $u$ is defined by the feature-set $I$ (i.e. without feature $\{m_s\}$).

Based on the above we can write:

$$\begin{aligned} Obj_{(ows)}(I\cup \{m_s\})&\ge \sum \limits _{i\in I\cup \{m_s\}}(v_{(I)}^Tx_{fc}(i))^2-\frac{1}{n}\sum \limits _{i\in I\cup \{m_s\}} x_{fc}(i)^Tx_{fc}(i)\\&=\sum \limits _{i\in I}(v_{(I)}^Tx_{fc}(i))^2-\frac{1}{n}\sum \limits _{i\in I} x_{fc}(i)^Tx_{fc}(i)\\&+(v_{(I)}^Tx_{fc}(m_s))^2-\frac{1}{n} x_{fc}(m_s)^Tx_{fc}(m_s)\\&=Obj_{(ows)}(I)+(v_{(I)}^Tx_{fc}(m_s))^2-\frac{1}{n} x_{fc}(m_s)^Tx_{fc}(m_s) \end{aligned}$$

$\square$

Proof

(Proof of Theorem 2) We will start by decomposing the components $\lambda _{1}(Cov)$ and $\mathbf{Trace}(Cov)$.

For $\lambda _{1}(Cov)$ we have:

$$\begin{aligned} \lambda _{1}(Cov)&= \lambda _{1}(C_m^uX_{fc}X_{fc}^TC_m^u)\\&= \lambda _{1}(X_{fc}^TC_m^uX_{fc})\\&= \lambda _{1}\left( X_{fc}^T\mathbf{diag}(u)\left( I -\tfrac{1}{\mathbf{card}(u)}e_me^T_m\right) \mathbf{diag}(u)X_{fc}\right) \\&= \max \limits _{||v||=1}v^T\left( X_{fc}^T\mathbf{diag}(u) \left( I-\tfrac{1}{\mathbf{card}(u)}e_me^T_m\right) \mathbf{diag}(u)X_{fc}\right) v\\&= \max \limits _{||v||=1}v^T(X_{fc}^T\mathbf{diag}(u)X_{fc})v\\&-\tfrac{1}{\mathbf{card}(u)}v^T(X_{fc}^T\mathbf{diag}(u)e_me^T_m \mathbf{diag}(u)X_{fc})v)\\&= \max \limits _{||v||=1}\sum \limits _{i=1}^mu_i(v^Tx_{fc}(i))^2 -\tfrac{1}{\mathbf{card}(u)}\left( v^T\sum \limits _{i=1}^mu_ix_{fc}(i)\right) ^2 \end{aligned}$$

In the above derivations we have used the following, easily verifiable facts: $\lambda _1(AA^T)=\lambda _1(A^TA), \,C_m^u=(C_m^u)^T, \,C_m^u=(C_m^u)^2$ and $\mathbf{diag}(u)=(\mathbf{diag}(u))^2$.

For $\mathbf{Tr}(Cov)$ we have:

$$\begin{aligned} \mathbf{Tr}(Cov)&= \mathbf{Tr}(C_m^uX_{fc}X_{fc}^TC_m^u)\\&= \mathbf{Tr}(X_{fc}^TC_m^uX_{fc})\\&= \mathbf{Tr}\left( X_{fc}\mathbf{diag}(u) \left( I-\tfrac{1}{\mathbf{card}(u)}e_me^T_m\right) \mathbf{diag}(u)X_{fc}\right) \\&= \mathbf{Tr}(X_{fc}^T\mathbf{diag}(u)X_{fc})\\&-\tfrac{1}{\mathbf{card}(u)}\mathbf{Tr}(X_{fc}^T \mathbf{diag}(u)e_me^T_m\mathbf{diag}(u)X_{fc}) \end{aligned}$$

In these derivations we have used the following properties of the matrix Trace: $\mathbf{Tr}(AA^T)=\mathbf{Tr}(A^TA), \,\mathbf{Tr}(A+B)=\mathbf{Tr}(A)+\mathbf{Tr}(B)$ and $\mathbf{Tr}(\beta A)=\beta \mathbf{Tr}(A)$.

Now for $\mathbf{Tr}(X_{fc}^T\mathbf{diag}(u)X_{fc})$ we have:

$$\begin{aligned} \mathbf{Tr}(X_{fc}^T\mathbf{diag}(u)X_{fc})&= \mathbf{Tr}(\mathbf{diag} (u)X_{fc}X_{fc}^T\mathbf{diag}(u))\\&= \sum \limits _{i=1}^m u_i x_{fc}(i)^Tx_{fc}(i) \end{aligned}$$

Finally, for $\frac{1}{\mathbf{card}(u)}\mathbf{Tr}(X_{fc}^T\mathbf{diag} (u)e_me^T_m\mathbf{diag}(u)X_{fc})$ we have:

$$\begin{aligned}&\tfrac{1}{\mathbf{card}(u)}\mathbf{Tr}(X_{fc}^T\mathbf{diag}(u)e_me^T_m \mathbf{diag}(u)X_{fc})\\&=\tfrac{1}{\mathbf{card}(u)}\mathbf{Tr}(e^T_m\mathbf{diag}(u)X_{fc}X_{fc}^T \mathbf{diag}(u)e_m)\\&=\tfrac{1}{\mathbf{card}(u)}\left( \sum \limits _{i=1}^n u_ix_{fc}(i)\right) ^T\left( \sum \limits _{i=1}^m u_ix_{fc}(i)\right) \end{aligned}$$

Based on the afore derivations we can write the objective function in Eq. 4 as:

$$\begin{aligned} \max \limits _{||v||=1}\max \limits _{u\in \{0,1\}^n}\left[ O_1-\frac{1}{\mathbf{card}(u)}O_2\right] \end{aligned}$$

(8)

Where

$$\begin{aligned} O_1&= \sum \limits _{i=1}^mu_i\left[ (v^Tx_{fc}(i))^2-\frac{1}{n} x_{fc}(i)^Tx_{fc}(i)\right] \end{aligned}$$

(9)

$$\begin{aligned} O_2&= \left( v^T\sum \limits _{i=1}^mu_ix_{fc}(i)\right) ^2-\frac{1}{n}\left( \sum \limits _{i=1}^m u_ix_{fc}(i)\right) ^T\left( \sum \limits _{i=1}^m u_ix_{fc}(i)\right) \end{aligned}$$

(10)

We can now employ the fact that for all $v$ such that $||v||=1$ it holds that $\lambda _1(A)\ge v^TAv$ and replace the $v$ vector in equations $O_1$ and $O_2$ with a fixed normalized vector. For deriving a tight bound we employ $v_{(I)}$, the dominant eigenvector of $C_m^uX_{fc}X_{fc}^TC_m^u$, where $u$ is defined based on the feature set $I$ (excluding feature $m_s$).

Based on the above (with a slight abuse of notation, since the global objective is lower bounded using $v_{(I)}$ and not the individual $O_1$ and $O_2$), we can write $O_1(I\cup \{m_s\})$ as:

$$\begin{aligned} O_1(I\cup \{m\})\ge O_1(I)+(v_{(I)}^Tx_{fc}(m_s))^2-\frac{1}{n} x_{fc}(m_s)^Tx_{fc}(m_s) \end{aligned}$$

Now $O_2(I\cup \{m_s\})$ can be written as (again with a slight abuse of notation, since the global objective is lower bounded using $v_{(I)}$ and not the individual $O_1$ and $O_2$):

$$\begin{aligned} O_2(I\cup \{m_s\})&\ge \left( v_{(I)}^T\sum \limits _{i\in I\cup \{m_s\}}x_{fc}(i)\right) ^2\!-\!\frac{1}{n}\left( \sum \limits _{i\in I\cup \{m_s\}} x_{fc}(i)\right) ^T\left( \sum \limits _{i\in I\cup \{m_s\}} x_{fc}(i)\right) \\&= O_2(I)+(v_{(I)}^Tx_{fc}(m_s))^2+2\left( \sum \limits _{i\in I}v_{(I)}^Tx_{fc}(i)\right) (v_{(I)}^Tx_{fc}(m_s))\\&-\frac{1}{n}x_{fc}(m_s)^Tx_{fc}(m_s)-\frac{2}{n} \left( \sum \limits _{i\in I}x_{fc}(i)\right) ^T x_{fc}(m_s) \end{aligned}$$

In order to complete the proof we must express $\frac{1}{\mathbf{card}(I)+1}O_2(I)$ in the form of $\frac{1}{\mathbf{card}(I)}O_2(I)+C$. In order to achieve this goal we take into account the fact that $\frac{1}{\mathbf{card}(I)+1}=\frac{1}{\mathbf{card}(I)} -\frac{1}{\mathbf{card}(I)(\mathbf{card}(I)+1)}$ and write $Obj(I\cup \{m\})$ as:

With some simple algebraic derivations the final bound can be derived. $\square$

Proof

(Proof of Theorem 3) Recall that in Schur complement deflation, the deflation step is performed as follows:

$$\begin{aligned} A_t=A_{t-1}-\frac{A_{t-1}x_tx_t^TA_{t-1}}{x_t^TA_{t-1}x_t} \end{aligned}$$

Now if we consider that $A_t=X_{fc}^{(t)}(X_{fc}^{(t)})^T$ and also that $x_t$ is the dominant eigenvector of matrix $Cov^{(t-1)}$, we can write:

$$\begin{aligned}&A_t=A_{t-1}-\frac{A_{t-1}x_tx_t^TA_{t-1}}{x_t^TA_{t-1}x_t} \Rightarrow X_{fc}^{(t)}(X_{fc}^{(t)})^T=X_{fc}^{(t-1)}(X_{fc}^{(t-1)})^T\nonumber \\&\quad -\frac{X_{fc}^{(t-1)}(X_{fc}^{(t-1)})^Tx_tx_t^TX_{fc}^{(t-1)} (X_{fc}^{(t-1)})^T}{x_t^TX_{fc}^{(t-1)}(X_{fc}^{(t-1)})^Tx_t} \end{aligned}$$

(11)

Recall that $x_t$ is the dominant eigenvector of $Cov^{(t-1)}$ that can be written as $Cov^{(t-1)}=C_m^{u(t-1)}X_{fc}^{(t-1)}(X_{fc}^{(t-1)})^TC_m^{u(t-1)}$ (i.e. it is based on the selected feature subset $u(t-1)$).

Since $Cov^{(t-1)}$ is a double-centered matrix (its rows and columns are centered through the multiplication with $C_m^{u(t-1)}$), its dominant eigenvector will also be centered, i.e. $(C_m^{u(t-1)})x_t=x_t$. Based on this property, we can write:

$$\begin{aligned} x_t^TX_{fc}^{(t-1)}(X_{fc}^{(t-1)})^Tx_t&= x_t^TC_m^{u(t-1)} X_{fc}^{(t-1)}(X_{fc}^{(t-1)})^TC_m^{u(t-1)}x_t\nonumber \\&= x_t^TCov^{(t-1)}x_t=\lambda _{\max }^{(t-1)} \end{aligned}$$

(12)

It should be noted that if we are analyzing one-way stable Sparse PCA, we can directly derive that $x_t^TX_{fc}^{(t-1)}(X_{fc}^{(t-1)})^Tx_t=\lambda _{\max }^{(t-1)}$. Now, $(X_{fc}^{(t-1)})^Tx_t$ can be written as:

$$\begin{aligned} (X_{fc}^{(t-1)})^Tx_t=(X_{ft-1)})^TC_m^{u(t-1)}x_t= \sqrt{\lambda _{\max }^{(t-1)}}v_t, \end{aligned}$$

(13)

where $v_t$ is the dominant eigenvector of $(X_{fc}^{(t-1)})^TC_m^{u(t-1)}X_{fc}^{(t-1)}$ and $\lambda _{\max }^{(t-1)}$ is the dominant eigenvector of $Cov^{(t-1)}$.

We should again note that in the one-way stable case, we can directly derive that $(X_{fc}^{(t-1)})^Tx_t=\sqrt{\lambda _{\max }^{(t-1)}}v_t$.

Based on Eqs. 11–13 we can write

$$\begin{aligned} A_t&= A_{t-1}-\frac{A_{t-1}x_tx_t^TA_{t-1}}{x_t^TA_{t-1}x_t}\\ \Rightarrow X_{fc}^{(t)}(X_{fc}^{(t)})^T&= X_{fc}^{(t-1)} (X_{fc}^{(t-1)})^T-X_{fc}^{(t-1)}v_tv_t^T(X_{fc}^{(t-1)})^T\\ \Rightarrow X_{fc}^{(t)}&= X_{fc}^{(t-1)}(I-v_tv_t^T) \end{aligned}$$

$\square$