Kernel Probabilistic K-Means Clustering

Abstract

Kernel fuzzy c-means (KFCM) is a significantly improved version of fuzzy c-means (FCM) for processing linearly inseparable datasets. However, for fuzzification parameter $m=1$, the problem of KFCM (kernel fuzzy c-means) cannot be solved by Lagrangian optimization. To solve this problem, an equivalent model, called kernel probabilistic k-means (KPKM), is proposed here. The novel model relates KFCM to kernel k-means (KKM) in a unified mathematic framework. Moreover, the proposed KPKM can be addressed by the active gradient projection (AGP) method, which is a nonlinear programming technique with constraints of linear equalities and linear inequalities. To accelerate the AGP method, a fast AGP (FAGP) algorithm was designed. The proposed FAGP uses a maximum-step strategy to estimate the step length, and uses an iterative method to update the projection matrix. Experiments demonstrated the effectiveness of the proposed method through a performance comparison of KPKM with KFCM, KKM, FCM and k-means. Experiments showed that the proposed KPKM is able to find nonlinearly separable structures in synthetic datasets. Ten real UCI datasets were used in this study, and KPKM had better clustering performance on at least six datsets. The proposed fast AGP requires less running time than the original AGP, and it reduced running time by 76–95% on real datasets.

1. Introduction

Clustering is an important unsupervised method, and the purpose of clustering is to divide a dataset into multiple clusters (or classes) with high intra-cluster similarity and low inter-cluster similarity. There have been many clustering algorithms, such as k-means (KM) and its variants [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]. Others are based on minimal spanning trees [17,18,19], density analysis [20,21,22,23,24,25], spectral analysis [26,27], subspace clustering [28,29], etc.

Generally, k-means minimizes the sum of squared Euclidean distances between each sample point and its nearest clustering center [1]. One variant of k-means is kernel k-means (KKM) [30,31,32,33], which is able to find nonlinearly separable structures by using the kernel function method. Another variant of k-means is fuzzy c-means (FCM) [2], which determines partitions by computing the membership degree of each data point to each cluster. The higher the membership degree, the greater the possibility of the data point belonging to the cluster. Although FCM is more flexible in applications [11,12,13,14,15,16], it is primarily suitable for linearly separable datasets. Kernel fuzzy c-means (KFCM) [34] is a significantly improved version of fuzzy c-means for clustering linearly inseparable datasets. However, the problem of KFCM with fuzzification parameter $m=1$ cannot be solved by existing methods.

To solve the special case of KFCM for $m=1$, a novel model called kernel probabilistic k-means (KPKM) is proposed. In fact, KPKM is a nonlinear programming model constrained on linear equalities and linear inequalities, and it is equivalent to KFCM at m = 1. In theory, the proposed KPKM can be solved by the active gradient projection (AGP) method [35,36]. Since the AGP method may take too much time on large datasets, we further propose a fast AGP (FAGP) algorithm to solve KPKM more efficiently. Moreover, we report experiments demonstrating its effectiveness compared with KFCM, KKM, FCM, and KM.

The paper is organized as follows: Section 2 reviews previous work. The KPKM algorithm is proposed in Section 3. Section 4 proposes a solution for KPKM. Section 5 presents descriptions and analyses of experiments. Conclusions and future work are mentioned in Section 6.

2. Background and Related Work

There has been a lot of work related to this paper, mainly including k-means, fuzzy c-means, kernel k-means and kernel fuzzy c-means.

K-means minimizes the sum of squared Euclidean distances between each sample point and its nearest cluster center. K-means first chooses initial clustering centers randomly or manually, and then partitions a dataset into several clusters (a data point belongs to the cluster whose clustering center is nearest to the data point), and computes the mean of a cluster as the clustering center. K-means repeatedly updates clustering centers and clusters until convergence. FCM has the same ideal as k-means. FCM introduces a membership degree $w_{ij}$ and a fuzzy coefficient m into the objective function. The higher $w_{ij}$ is, the greater possibility of the i-th data point belonging to the j-th cluster. K-means and FCM belong to partition-based clustering algorithms, and partition-based clustering algorithms usually are not able to cluster linearly inseparable datasets. Kernel method maps a linearly inseparable dataset into a linearly separable space, so kernel k-means (and FCM) using a kernel function can cluster linearly inseparable datasets.

2.1. K-Means and Fuzzy C-Means

Let $\mathit{X}=\left\{{\mathit{x}}_i|{\mathit{x}}_i\in {\mathbb{R}}^D,1\le i\le L\right\}$ represent a dataset. K-means divides it into k clusters. If ${\omega }_j$ denotes the j-th cluster, we have $\mathit{X}={\displaystyle \bigcup _{j=1}^K}{\omega }_j$ and $\forall 1\le i\ne j\le K,{\omega }_i\bigcap {\omega }_j=\varnothing$. Using $L_j=\left|{\omega }_j\right|$ to stand for the number of elements in ${\omega }_j$ with the center of ${\mathit{c}}_j$, k-means can be described as minimizing the following objective function:

$$J=\sum _{i=1}^K\sum _{{\mathit{x}}_i\in {\omega }_j}{∥{\mathit{x}}_i-{\mathit{c}}_j∥}^2$$

(1)

where

$${\mathit{c}}_j=\frac{1}{L_j}\sum _{{\mathit{x}}_i\in {\omega }_j}{\mathit{x}}_i.$$

(2)

Let $L={\displaystyle \sum _{j=1}^K}{L}_j$. By using c instead of k and assigning membership degree $w_{ij}$ to the i-th data point in the j-th cluster for $1\le i\le L$ and $1\le j\le C$, the fuzzy c-means clustering model can be formulated as follows:

$$\begin{array}{c}minJ={\displaystyle \sum _{j=1}^C}{\displaystyle \sum _{i=1}^L}{w}_{{}_{ij}}^m{∥{\mathit{x}}_i-{\mathit{c}}_j∥}^2\hfill \\ s.t.{\displaystyle \sum _{j=1}^C}{w}_{ij}=1,w_{ij}\ge 0\hfill \end{array}$$

(3)

where $m>1$, $w_{ij}$ and ${\mathit{c}}_j$ are computed alternately below [2].

$$\left\{\begin{array}{c}{w}_{ij}=\frac{{∥{\mathit{x}}_i-{\mathit{c}}_j∥}^{-\frac{2}{m-1}}}{{\displaystyle \sum _{k=1}^C}{∥{\mathit{x}}_i-{\mathit{c}}_k∥}^{-\frac{2}{m-1}}}\hfill \\ {\mathit{c}}_j=\frac{{\displaystyle \sum _{i=1}^L}{w}_{{}_{ij}}^m{\mathit{x}}_i}{{\displaystyle \sum _{i=1}^L}{w}_{{}_{ij}}^m}\hfill \end{array}\right..$$

(4)

2.2. Kernel K-Means and Kernel Fuzzy C-Means

To improve performance of k-means and fuzzy c-means in linearly inseparable datasets, we may develop their kernel versions by choosing a feature mapping $\phi (·):{\mathcal{R}}^D\to \mathcal{H}$ from data points to kernel Hilbert space [37]. Though $\phi$ is usually unknown, it must satisfy

$$K(\mathit{x},\mathit{y})={\left(\phi \left(\mathit{x}\right)\right)}^{\mathrm{T}}\phi \left(\mathit{y}\right)$$

(5)

where $K(\mathit{x},\mathit{y})$ is a kernel function. Commonly used kernel functions are presented in

Table 1. Commonly used kernel functions.

Name	Code
Linear kernel	$K(\mathit{x},\mathit{y})=〈\mathit{x},\mathit{y}〉$
Laplace Radial Basis Function kernel	$K_{\mathrm{lap}}(\mathit{x},\mathit{y})=exp\left(-\sigma ∥\mathit{x}-\mathit{y}∥\right)$
Gaussian Radial Basis Function kernel	$K_{\mathrm{gau}}(\mathit{x},\mathit{y})=exp\left(-\frac{{∥\mathit{x}-\mathit{y}∥}^2}{2{\sigma }^2}\right)$
Polynomial kernel	$K_{\mathrm{pol}}(\mathit{x},\mathit{y})={\left(\mathit{x}·\mathit{y}+\beta \right)}^{\alpha }$
Sigmoid kernel	$K_{tan}(\mathit{x},\mathit{y})=tanh\left(\alpha \mathit{x}·\mathit{y}+\beta \right)$

.

In Table 1, $〈\mathit{x},\mathit{y}〉=\mathit{x}·\mathit{y}$ denotes the inner product of $\mathit{x}$ and $\mathit{y}$, and $\sigma ,\alpha ,\beta$ are parameters of the kernel. The objective function of kernel k-means is defined as

$$J_k=\sum _{j=1}^K\sum _{{\mathit{x}}_i\in {\omega }_j}{∥\phi \left({\mathit{x}}_i\right)-{\mathit{c}}_j∥}^2$$

(6)

where

$${\mathit{c}}_j=\frac{1}{L_j}\sum _{{\mathit{x}}_i\in {\omega }_j}\phi \left({\mathit{x}}_i\right).$$

(7)

Using (5) and (7), we have

$$\begin{array}{c}{∥\phi \left({\mathit{x}}_i\right)-{\mathit{c}}_j∥}^2=K({\mathit{x}}_i,{\mathit{x}}_i)-\frac{2}{L_j}{\displaystyle \sum _{{\mathit{x}}_k\in {\omega }_j}^L}K({\mathit{x}}_k,{\mathit{x}}_i)+\frac{1}{{\left(L_j\right)}^2}{\displaystyle \sum _{{\mathit{x}}_l\in {\omega }_j}^L}{\displaystyle \sum _{{\mathit{x}}_h\in {\omega }_j}^L}K({\mathit{x}}_l,{\mathit{x}}_h)\hfill \end{array}.$$

(8)

Let $w_{ij}$ represent the membership degree of the i-th point belonging to the j-th class. Likewise, we can define the kernel FCM clustering model as follows.

$$\begin{array}{c}min{J}_f\left(\mathit{W}\right)={\displaystyle \sum _{j=1}^C}{\displaystyle \sum _{i=1}^L}{w}_{ij}^m{∥\phi \left({\mathit{x}}_i\right)-{\mathit{c}}_j∥}^2\hfill \\ s.t.{\displaystyle \sum _{j=1}^C}{w}_{ij}=1,w_{ij}\ge 0,m>1\hfill \end{array}$$

(9)

where

$${\mathit{c}}_j=\frac{{\displaystyle \sum _{l=1}^L}{w}_{lj}^m\phi ({\mathit{x}}_l)}{{\displaystyle \sum _{l=1}^L}{w}_{lj}^m},$$

(10)

$$\begin{array}{c}{∥\phi \left({\mathit{x}}_i\right)-{\mathit{c}}_j∥}^2=K({\mathit{x}}_i,{\mathit{x}}_i)-\frac{2}{{\displaystyle \sum _{l=1}^L}{w}_{lj}^m}{\displaystyle \sum _{k=1}^L}{w}_{kj}^{m}K({\mathit{x}}_k,{\mathit{x}}_i)+\frac{1}{{\left({\displaystyle \sum _{l=1}^L}{w}_{lj}^m\right)}^2}{\displaystyle \sum _{l=1}^L}{\displaystyle \sum _{h=1}^L}{w}_{lj}^m{w}_{hj}^{m}K({\mathit{x}}_l,{\mathit{x}}_h)\hfill \end{array}.$$

(11)

The membership degree is computed via

$$w_{ij}=\frac{{\left({∥\phi \left({\mathit{x}}_i\right)-{\mathit{c}}_j∥}^2\right)}^{-\frac{1}{m-1}}}{{\displaystyle \sum _{k=1}^C}{\left({∥\phi \left({\mathit{x}}_i\right)-{\mathit{c}}_k∥}^2\right)}^{-\frac{1}{m-1}}}.$$

(12)

3. Kernel Probabilistic K-Means

In this section, the kernel probabilistic k-means (KPKM) are proposed.

We first review the problem. When $m=1$, the KFCM model gets into a special case, namely,

$$\begin{array}{c}min{J}_f\left(\mathit{W}\right)={\displaystyle \sum _{j=1}^K}{\displaystyle \sum _{i=1}^L}{w}_{ij}{∥\phi \left({\mathit{x}}_i\right)-{\mathit{c}}_j∥}^2\hfill \\ s.t.{\displaystyle \sum _{j=1}^K}{w}_{ij}=1,w_{ij}\ge 0\hfill \end{array}$$

(13)

where

$${\mathit{c}}_j=\frac{{\displaystyle \sum _{l=1}^L}{w}_{lj}\phi \left({\mathit{x}}_l\right)}{{\displaystyle \sum _{l=1}^L}{w}_{lj}},$$

(14)

$$\begin{array}{c}{∥\phi \left({\mathit{x}}_i\right)-{\mathit{c}}_j∥}^2=K({\mathit{x}}_i,{\mathit{x}}_i)-\frac{2}{{\displaystyle \sum _{l=1}^L}{w}_{lj}}{\displaystyle \sum _{k=1}^L}{w}_{kj}K({\mathit{x}}_k,{\mathit{x}}_i)+\frac{1}{{\left({\displaystyle \sum _{l=1}^L}{w}_{lj}\right)}^2}{\displaystyle \sum _{l=1}^L}{\displaystyle \sum _{h=1}^L}{w}_{lj}{w}_{hj}K({\mathit{x}}_l,{\mathit{x}}_h)\hfill \end{array}.$$

(15)

This special case cannot be solved by Lagrangian optimization for $m>1$, because (12) cannot be computed when $m=1$ (when $m=1$, $\frac{1}{m-1}=\frac{1}{0}$ cannot be computed).

In this paper, we use the optimization methods to solve this problem, but the partial derivative of (13) with respect to $w_{ij}$ is ${∥\phi \left({\mathit{x}}_i\right)-{\mathit{c}}_j∥}^2$, and it does not contain $w_{ij}$.

In order to solve the problem, we introduce (14) into (13), and redefine the member degrees $w_{ij}$ as probability $p_{ij}$ for $1\le i\le L$ and $1\le j\le K$. Finally, we have

$$\begin{array}{c}J\left(\mathit{P}\right)={\displaystyle \sum _{j=1}^K}{\displaystyle \sum _{i=1}^L}{p}_{ij}{∥\phi \left({\mathit{x}}_i\right)-\frac{{\sum }_{l=1}^L{p}_{lj}\phi \left({\mathit{x}}_l\right)}{{\sum }_{l=1}^L{p}_{lj}}∥}^2\hfill \\ s.t.{\sum }_{j=1}^K{p}_{ij}=1,p_{ij}\ge 0\hfill \end{array}$$

(16)

where probability vector

$$\mathit{P}={\left\{p_{11},\dots ,p_{1K},p_{21},\dots p_{2K},\dots ,p_{LK}\right\}}^{\mathrm{T}}.$$

(16) is the proposed kernel probabilistic k-means.

The proposed KPKM is a soft clustering method. The $p_{ij}(0\le p_{ij}\le 1)$ is the probability of the i-th data point belonging to the j-th cluster. The higher $p_{ij}$ is, the greater possibility of the i-th data point belonging to the j-th cluster. In KKM, the membership degree has only two values (0 and 1). In the proposed KPKM, $w_{ij}\in \left[0,1\right]$ (although final $w_{ij}=0\mathtt{or}1$).

4. A Fast Solution to KPKM

KPKM is a nonlinear programming problem with linear equalities and linear inequalities constraints, and it is able to be solved by active gradient projection [35,36] theoretically. In this section, on the basis of AGP, we first calculate the gradient of of the objection function of KPKM, and then design a fast AGP algorithm to solve the KPKM model.

The fast AGP has two advantages: iteratively updating the projection matrix and estimating the maximum step length.

4.1. Gradient Calculation

For convenience, we define $F_{kj}$ as

$$F_{kj}={∥\phi \left({\mathit{x}}_k\right)-\frac{{\displaystyle \sum _{l=1}^L}{p}_{lj}\phi \left({\mathit{x}}_l\right)}{{\displaystyle \sum _{l=1}^L}{p}_{lj}}∥}^2.$$

(17)

Using the chain rule on (16), we obtain

$$\frac{\partial J}{\partial p_{ij}}={\displaystyle \sum _{k=1}^L}{p}_{kj}\frac{\partial F_{kj}}{\partial p_{ij}}+F_{ij}.$$

(18)

According to (17), we can further derive

$$\begin{array}{c}\frac{\partial F_{kj}}{\partial p_{ij}}=\frac{-2}{{\sum }_{i=1}^L{p}_{ij}}\left(K({\mathit{x}}_k,{\mathit{x}}_i)-\frac{{\sum }_{l=1}^L{p}_{lj}K({\mathit{x}}_k,{\mathit{x}}_l)}{{\sum }_{i=1}^L{p}_{ij}}\right.\left.-\frac{{\sum }_{l=1}^L{p}_{lj}K({\mathit{x}}_i,{\mathit{x}}_l)}{{\sum }_{i=1}^L{p}_{ij}}+\frac{{\sum }_{l=1}^L{\sum }_{h=1}^L{p}_{lj}{p}_{hj}K({\mathit{x}}_l,{\mathit{x}}_h)}{{\left({\sum }_{i=1}^L{p}_{ij}\right)}^2}\right)\hfill \end{array}.$$

(19)

By substituting (19) into (18), we finally get

$$\begin{array}{c}\frac{\partial J}{\partial p_{ij}}=K({\mathit{x}}_i,{\mathit{x}}_i)+\frac{1}{{\left({\sum }_{l=1}^L{p}_{lj}\right)}^2}{\displaystyle \sum _{l=1}^L}{\displaystyle \sum _{h=1}^L}{p}_{lj}{p}_{hj}K({\mathit{x}}_l,{\mathit{x}}_h)-\frac{2}{{\sum }_{l=1}^L{p}_{lj}}{\displaystyle \sum _{k=1}^L}{p}_{kj}K({\mathit{x}}_k,{\mathit{x}}_i)\hfill \end{array}.$$

(20)

and the gradient

$$\nabla J={\left[\begin{array}{ccccccc}\frac{\partial J}{\partial p_{11}}& \cdots & \frac{\partial J}{\partial p_{1K}}& \cdots & \frac{\partial J}{\partial p_{Lj}}& \cdots & \frac{\partial J}{\partial p_{LK}}\end{array}\right]}^{\mathrm{T}}.$$

(21)

4.2. Fast AGP

In the constraints of KPKM, there are L linear equalities and K×L linear inequalities,

$$\forall 1\le i\le L,{\displaystyle \sum _{j=1}^K}{p}_{ij}=1,$$

(22)

$$\forall 1\le i\le L,1\le j\le K,p_{ij}\ge 0.$$

(23)

Let $\varphi =K\times L$. Let ${\mathit{I}}_{\varphi \times \varphi }$ be the identity matrix of size $\varphi \times \varphi$. Define two matrices, inequality matrix $\mathit{A}$ and equality matrix $\mathit{E}$, where

$$\mathit{A}={\mathit{I}}_{\varphi \times \varphi },$$

(24)

$$\mathit{E}={\left[\begin{array}{ccc}\underset{K}{\underbrace{\begin{array}{ccc}1& \dots & 1\end{array}}}& \begin{array}{ccc}0& \dots & 0\end{array}& \begin{array}{ccc}0& \dots & 0\end{array}\\ \underset{K}{\underbrace{\begin{array}{ccc}0& \dots & 0\end{array}}}& \underset{K}{\underbrace{\begin{array}{ccc}1& \dots & 1\end{array}}}& \begin{array}{ccc}0& \dots & 0\end{array}\\ & \dots & \\ \begin{array}{ccc}0& \dots & 0\end{array}& \begin{array}{ccc}0& \dots & 0\end{array}& \underset{K}{\underbrace{\begin{array}{ccc}1& \dots & 1\end{array}}}\end{array}\right]}_{L\times \varphi }.$$

(25)

Note that each row of $\mathit{A}$ corresponds to one and only one inequality in (23), and each row of $\mathit{E}$ corresponds to one and only one equality in (22). Accordingly, the KPKM’s constraints can be simply expressed as

$$\mathit{AP}\ge \mathbf{0},\mathit{EP}=\mathbf{1}$$

(26)

where $\mathbf{1}={[1,\dots ,1]}_{L\times 1}^{\mathtt{T}}$. Let ${\mathit{P}}^{\left(0\right)}$ stand for a randomly initialized probability vector, and ${\mathit{P}}^{\left(k\right)}$ for the probability vector at iteration k. The rows of inequality matrix $\mathit{A}$ can be broken into two groups: one is active; the other is inactive. The active group is composed of all inequalities that must work exactly as an equality at ${\mathit{P}}^{\left(k\right)}$, whereas the inactive group is composed of the left inequalities. If ${\mathit{A}}_1^{\left(k\right)}$ and ${\mathit{A}}_2^{\left(k\right)}$ respectively denote the active group and the inactive group, we have

$${\mathit{A}}_1^{\left(k\right)}{\mathit{P}}^{\left(k\right)}=\mathbf{0},$$

(27)

$${\mathit{A}}_2^{\left(k\right)}{\mathit{P}}^{\left(k\right)}>\mathbf{0}.$$

(28)

At iteration k, the active matrix ${\mathit{N}}^{\left(k\right)}$ is defined as

$${\mathit{N}}^{\left(k\right)}=\left[\begin{array}{c}{\mathit{A}}_1^{\left(k\right)}\hfill \\ \mathit{E}\hfill \end{array}\right].$$

(29)

When $\mathit{N}={\mathit{N}}^{\left(k\right)}$ is not a square matrix, we can construct its projection matrix ${\mathit{G}}^{\left(k\right)}$ and the corresponding orthogonal projection matrix ${\mathit{Q}}^{\left(k\right)}$ as follows:

$${\mathit{G}}^{\left(k\right)}={\mathit{N}}^{\mathrm{T}}{\left(\mathit{N}{\mathit{N}}^T\right)}^{-1}\mathit{N},$$

(30)

$${\mathit{Q}}^{\left(k\right)}=\mathit{I}-{\mathit{G}}^{\left(k\right)}.$$

(31)

Suppose that $\mathit{n}$ from ${\mathit{A}}_2^{(k-1)}$ is the active row vector at iteration k, satisfying ${\mathit{n}}^{\mathrm{T}}{\mathit{P}}^{(k-1)}>\mathbf{0}$ and ${\mathit{n}}^{\mathrm{T}}{\mathit{P}}^{\left(k\right)}=\mathbf{0}$. According to matrix theory [38], we can more efficiently compute ${\mathit{G}}^{\left(k\right)}$ and ${\mathit{Q}}^{\left(k\right)}$ by

$$\begin{array}{c}{\mathit{G}}^{\left(k\right)}={\mathit{G}}^{(k-1)}+{\mathit{Q}}^{(k-1)}{\mathit{n}}^{\mathrm{T}}{〈{\mathit{Q}}^{(k-1)}{\mathit{n}}^{\mathrm{T}},{\mathit{Q}}^{(k-1)}{\mathit{n}}^{\mathrm{T}}〉}^{-1}\mathit{n}{\mathit{Q}}^{(k-1)}\hfill \end{array}.$$

(32)

Furthermore, we can compute the projected gradient by

$${\mathit{d}}^{\left(k\right)}=-{\mathit{Q}}^{\left(k\right)}\nabla J({\mathit{P}}^{\left(k\right)}).$$

(33)

Using ${\mathit{d}}^{\left(k\right)}$, we update the probability vector ${\mathit{P}}^{\left(k+1\right)}$ as

$${\mathit{P}}^{\left(k+1\right)}={\mathit{P}}^{\left(k\right)}+t^{\left(k\right)}{\mathbf{d}}^{\left(k\right)}$$

(34)

where $t^{\left(k\right)}$ is the step length. Usually, $t^{\left(k\right)}$ is chosen as a small number. For fast convergence, we estimate the maximum step length as follows.

$$t^{\left(k\right)}=t_{max}^{\left(k\right)}.$$

(35)

(1)

Let $p_{ij}^{(k+1)}=p_{ij}^{\left(k\right)}+t_{ij}{d}_{ij}^{\left(k\right)}=0$ for $p_{ij}^{\left(k\right)}>0$;

(2)

Compute $t_{ij}=-p_{ij}^{\left(k\right)}/\phantom{p_{ij}^{\left(k\right)}{d}_{ij}^{\left(k\right)}}{d}_{ij}^{\left(k\right)}$ for $p_{ij}^{\left(k\right)}>0$ and $d_{ij}^{\left(k\right)}<0$;

(3)

$t_{max}^{\left(k\right)}=min\left\{t_{ij}\right\}$.

As shown in Figure 1, maximum step length is compared with small step length.

When $\mathit{N}={\mathit{N}}^{\left(k\right)}$ is a square matrix, it must be invertible. In this case, we actually have ${\mathit{G}}^{\left(k\right)}={\mathit{N}}^{\mathrm{T}}{\left(\mathit{N}{\mathit{N}}^T\right)}^{-1}\mathit{N}={\mathit{I}}_{LK\times LK}$, and ${\mathit{Q}}^{\left(k\right)}=\mathbf{0}$. Thus, ${\mathit{d}}^{\left(k\right)}=-{\mathit{Q}}^{\left(k\right)}\nabla J\left({\mathit{P}}^{\left(k\right)}\right)$$=\mathbf{0}$ is not a feasible descent direction. A new descent direction can be computed as follows:

(1)

Compute a new vector,

$${\mathit{q}}^{\left(k\right)}={\left(\mathit{N}{\mathit{N}}^T\right)}^{-1}\mathit{N}\nabla J={\left({\mathit{N}}^{\mathrm{T}}\right)}^{-1}\nabla J.$$

(36)

(2)

Break ${\mathit{q}}^{\left(k\right)}$ into two parts ${\mathit{q}}_1^{\left(k\right)}$ and ${\mathit{q}}_2^{\left(k\right)}$, namely,

$${\mathit{q}}^{\left(k\right)}={\left[\begin{array}{cc}{\left({\mathit{q}}_1^{\left(k\right)}\right)}^{\mathrm{T}}& {\left({\mathit{q}}_2^{\left(k\right)}\right)}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}$$

(37)

where the size of ${\mathit{q}}_1^{\left(k\right)}$ is the number of rows of ${\mathit{A}}_1^{\left(k\right)}$, and that of ${\mathit{q}}_2^{\left(k\right)}$ is the number of rows of $\mathit{E}$.

(3)

If ${\mathit{q}}_1^{\left(k\right)}\ge \mathbf{0}$, stop. Otherwise, choose any element from ${\mathit{q}}_1^{\left(k\right)}$ that is less than 0 and delete the corresponding row of ${\mathit{A}}_1^{\left(k\right)}$; then use (29)–(31) and (33) to compute ${\mathit{d}}^{\left(k\right)}$.

The above fast AGP solution to KPKM is outlined in Algorithm 1. Compared with the original AGP, the fast AGP has two advantages: iteratively updating the projection matrix (shown in (32)) and estimating the maximum step length (shown in (35)).

4.3. Analysis of Complexity

In this section, the computational complexities for the traditional AGP algorithm and the proposed FAGP algorithm are analyzed. Let $T_A$ represent the iteration number of AGP. Let $T_F$ represent the iteration number of FAGP. In computations of all matrices, the computational complexity of the projection matrix $\mathit{G}$ is the highest. In the AGP algorithm, computing $\mathit{G}$ via (30) requires $O\left(K^3{L}^3\right)$, so the total computational complexity for AGP algorithm is $O\left(T_A{K}^3{L}^3\right)$. In the FAGP algorithm, $\mathit{G}$ is computed via (32), and it does not require to compute the inverse of matrices (i.e., ${\left({\mathit{NN}}^T\right)}^{-1}$), and computing $\mathit{G}$ via (32) requires $O\left(K^3{L}^2\right)$, so the total computational complexity for FAGP algorithm (i.e., Algorithm 1) is $O\left(T_F{K}^3{L}^2\right)$.

Algorithm 1 Fast active gradient projection (FAGP).

Input: X and K

Do:

(1) Let $k=0$, $\mathit{N}={\mathit{N}}^{\left(0\right)}=\mathit{E}$,${\mathit{G}}^{\left(0\right)}={\mathit{N}}^{\mathrm{T}}{\left(\mathit{N}{\mathit{N}}^T\right)}^{-1}\mathit{N}$,       initialize ${\mathit{P}}^{\left(0\right)}$, go to (4);

(2) Find $\mathit{n}\in {\mathit{A}}_2^{(k-1)}$ meeting ${\mathit{n}}^{\mathrm{T}}{\mathit{P}}^{(k-1)}>0$ and ${\mathit{n}}^{\mathrm{T}}{\mathit{P}}^{\left(k\right)}=0$;

(3) Compute ${\mathit{G}}^{\left(k\right)}$ by (32);

(4) Compute ${\mathit{Q}}^{\left(k\right)}$ by (31), and ${\mathit{d}}^{\left(k\right)}$ by (33);

(5) Compute $t^{\left(k\right)}$ by (35), and ${\mathit{P}}^{(k+1)}$ by (34);

(6) Let $k=k+1$, if ${\mathit{G}}^{\left(k\right)}\ne {\mathit{I}}_{LK\times LK}$, go to (2);

(7) Construct $\mathit{N}={\mathit{N}}^{\left(k\right)}$ by (29), and by (36);

(8) Break into ${\mathit{q}}_1^{\left(k\right)}$ and ${\mathit{q}}_2^{\left(k\right)}$ by (37);

(9) If all elements of ${\mathit{q}}_1^{\left(k\right)}$ are more than 0, stop;

(10) Choose any element less than 0 from ${\mathit{q}}_1^{\left(k\right)}$,       and delete the corresponding row of ${\mathit{A}}_1^{\left(k\right)}$;

(11) Reconstruct $\mathit{N}={\mathit{N}}^{\left(k\right)}$ by (29);

(12) Compute ${\mathit{G}}^{\left(k\right)}$ by (30), and ${\mathit{Q}}^{\left(k\right)}$ by (31);

(13) Compute ${\mathit{d}}^{\left(k\right)}$ by (33), and $t^{\left(k\right)}$ by (35);

(14) Compute ${\mathit{P}}^{(k+1)}$ by (34);

(15) Let $k=k+1$, and go to (7).

Output: the probability vector $\mathit{P}$

5. Experimental Results

In order to evaluate performance of the proposed KPKM model (equivalent to KFCM at m = 1) solved by the FAGP algorithm, we have conducted a lot of experiments on one synthetic dataset, ten UCI datasets (http://archive.ics.uci.edu/ml (accessed on 8 March 2021)) and the MNIST dataset (http://yann.lecun.com/exdb/mnist/ (accessed on 8 March 2021)). These datasets are detailed in Section 5.1, Section 5.2 and Section 5.3. In Section 5.1, we use a synthetic dataset to compare KPKM using a Gaussian kernel with KPKM using a linear kernel. In Section 5.2 and Section 5.3, we compare the proposed KPKM with KFCM, KKM, FCM, and KM to evaluate the performance of KPKM solved by the proposed fast AGP. Moreover, the Section 5.4 and Section 5.5 evaluate the descent stability and convergence speed of the proposed FAGP, respectively. In the experiments, we implemented our own MATLAB code for KPKM, KFCM, and KKM with two build-in functions, sparse and full, employed for matrix optimization. Moreover, we called MATLAB’s build-in functions fcm and kmeans for FCM and KM, respectively.

All experiments were carried out on a PC with an Intel(R) Core(TM) i7-4790 CPU at 3.60 GHz, 8.00 G RAM, running Windows 7 and MATLAB 2015a.

5.1. The Experiment on the Synthetic Dataset

In this experiment, we analyzed the influences of the Gaussian radial basis function kernel and the linear kernel on KPKM when clustering one synthetic dataset, which is shown in Figure 2. The synthetic dataset contained 300 points, with 100 and 200 being from two linearly inseparable classes, disc and ring, respectively. The result is displayed in Figure 3. Obviously, Gaussian radial basis function KPKM can cluster perfectly, whereas linear KPKM cannot.

5.2. Experiment on Ten UCI Datasets

In this experiment, we compare the clustering performance of KPKM with the performances of KFCM, KKM, FCM, and KM in terms of three measures: normalized mutual information (NMI), adjusted Rand index (ARI), and v-measure (VM) [39,40]. NMI is a normalization of the mutual information score to scale the results between 0 (no mutual information) and 1 (perfect correlation). The NMI is defined as

$$\mathrm{NMI}(U,V)=\frac{\mathrm{MI}(U,V)}{\mathrm{mean}\left(H\right(U),H(V\left)\right)}$$

(38)

where $H\left(U\right)$ is the entropy of U; MI is given by

$$\mathrm{MI}(U,V)=\sum _{i=1}^{\left|U\right|}\sum _{j=1}^{\left|V\right|}\frac{|U_i\cap V_j|}{N}log\frac{N|U_i\cap V_j|}{|U_i\left|\right|V_j|}$$

(39)

where $|U_i|$ is the number of the samples in cluster $U_i$. ARI is an adjusted similarity measure between two clusters. The ARI is defined as

$$\mathrm{ARI}=\frac{\mathrm{RI}-E\left(\mathrm{RI}\right)}{max\left(\mathrm{RI}\right)-E\left(\mathrm{RI}\right)}$$

(40)

where RI is the ratio of data points clustered correctly to all data points, and E(RI) is the expectation of RI. VM is a harmonic mean between homogeneity and completeness, where homogeneity means that each cluster is a subset of a single class, and completeness means that each class is a subset of a single cluster. The VM is defined as

$$\mathrm{VM}=\frac{(1+\beta )\times \mathrm{homogeneity}\times \mathrm{completeness}}{(\beta \times \mathrm{homogeneity}+\mathrm{completeness})}$$

(41)

The higher the NMI, ARI and VM, the better the performance.

We describe the ten UCI datasets represented in

Table 2. The ten UCI datasets used.

Name	Code	Instances	Classes	Dimensions
Iris	R1	150	3	4
Seeds	R2	210	3	7
Segmentation	R3	210	7	19
Glass	R4	214	6	9
Ionosphere	R5	351	2	33
Dermatology	R6	358	6	34
Breast-cancer	R7	683	2	9
Natural	R8	2000	9	294
Yeast	R9	2426	3	24
Waveform	R10	5000	3	21

by name, code, number of instances, number of classes and number of dimensions. A set of experiential parameters was selected. We set m = 2 for KFCM, and m = 1.3 for FCM. With the Gaussian radial basis function kernel, we also selected appropriate $\sigma$ (shown in

Table 3. Parameter $\sigma$ selected for kernel clustering methods.

Dataset	KPKM	KFCM	KKM
R1	1.08	1.22	0.9
R2	1.9	2	2
R3	510	530	540
R4	510	100	510
R5	1.5	1	1.3
R6	3.3	2	18
R7	12	10	15
R8	3	0.9	2.7
R9	15	10	25
R10	10	15	13

), and report the clustering results in

Table 4. Comparisons of KPKM with kernel fuzzy c-means (KFCM), kernel k-means (KKM), fuzzy c-means (FCM), and k-means (KM) on ten UCI datasets (DS) in terms of normalized mutual information (NMI), adjusted Rand index (ARI), and v-measure (VM).

DS	Method	KPKM	KFCM	KKM	FCM	KM
R1	NMI	0.8146	0.7900	0.7820	0.6723	0.6733
	ARI	0.8119	0.8015	0.7590	0.5763	0.5779
	VM	0.8146	0.7900	0.7820	0.7081	0.7149
R2	NMI	0.7073	0.6949	0.7038	0.6949	0.7025
	ARI	0.7141	0.7166	0.7231	0.7166	0.7135
	VM	0.7073	0.6949	0.7038	0.6949	0.6999
R3	NMI	0.5555	0.5503	0.5154	0.4678	0.5132
	ARI	0.3909	0.4141	0.3429	0.3172	0.3313
	VM	0.5553	0.5503	0.5139	0.5729	0.5252
R4	NMI	0.4436	0.3594	0.3943	0.3489	0.4178
	ARI	0.2796	0.2137	0.2542	0.2126	0.2551
	VM	0.4424	0.3593	0.3934	0.3807	0.3857
R5	NMI	0.2476	0.2390	0.2715	0.1299	0.1349
	ARI	0.1747	0.1098	0.1657	0.1727	0.1777
	VM	0.2476	0.2390	0.2715	0.1298	0.1348
R6	NMI	0.2919	0.2068	0.2778	0.1046	0.1032
	ARI	0.1795	0.1388	0.1698	0.0261	0.0266
	VM	0.2919	0.2065	0.2776	0.1095	0.1056
R7	NMI	0.7903	0.7825	0.7741	0.7478	0.7478
	ARI	0.8796	0.8741	0.8674	0.8464	0.8464
	VM	0.7903	0.7825	0.7741	0.7478	0.7478
R8	NMI	0.0531	0.0326	0.0556	0.0521	0.0536
	ARI	0.0253	0.0303	0.0283	0.0273	0.0261
	VM	0.0525	0.0312	0.0550	0.0495	0.0529
R9	NMI	0.0052	0.0031	0.0050	0.0043	0.0050
	ARI	0.0118	0.0084	0.0118	0.0109	0.0117
	VM	0.0052	0.0031	0.0050	0.0045	0.0045
R10	NMI	0.3654	0.3162	0.3637	0.3606	0.3622
	ARI	0.2546	0.2377	0.2541	0.2529	0.2536
	VM	0.3654	0.3162	0.3637	0.3559	0.3622

. The better results in each case is highlighted in bold. (The performance of many a method depends on parameters [41]. By experiments we found that for different algorithms with the same kernel function, the most appropriate parameters are usually different, so we used different settings to make the methods have the best clustering performances they could. We also tried to select the most appropriate parameter by using our experience.) We ran every algorithm 10 times, and average results were calculated.

From Table 4, we can see that

(1)

For NMI, KPKM had the best clustering results on nine datasets, so the clustering performance of KPKM was the best for NMI.

(2)

For ARI, KPKM had the best clustering results on six datasets. KFCM, KKM, FCM, and KM performed the best on 2, 1, 0, and 1 datasets, respectively, so KPKM performed the best for ARI.

(3)

For VM, KPKM, KFCM, KKM, FCM, and KM had the best clustering results on 7, 0, 2, 1, and 0 datasets, respectively. Thus, KPKM had the best clustering performance for VM.

Overall, KPKM is better than the other models.

5.3. Experiment on the MNIST Dataset

In this experiment, we used the MNIST dataset to compare KPKM with KFCM and KKM when clustering digital images. The MNIST dataset was composed of handwritten digits, with 60,000 examples for training and 10,000 examples for testing. All digits were size-normalized and centered in fixed-size images. To evaluate the clustering performances of KPKM, KFCM, and KKM, we randomly chose 1000 training examples; 100 are shown in Figure 4.

Moreover, we defined a CWSS kernel based on complex wavelet structural similarity (CWSS) [42]. CWSS may be regarded as a coefficient that does not change the structural content of an image. By using $CWSS(\mathit{x},\mathit{y})$ to denote the similarity between two images $\mathit{x}$ and $\mathit{y}$, we can express the CWSS kernel as

$$K_{\mathrm{cwss}}(\mathit{x},\mathit{y})=exp\left(-\frac{1-CWSS(\mathit{x},\mathit{y})}{2{\sigma }^2}\right)$$

(42)

where $\sigma =5$ is set for KFCM, and $\sigma =1$ for both KPKM and KKM. Additionally, m = 1.3 is set for KFCM. We present the results in

Table 5. Comparisons of KPKM with KFCM and KKM on the MNIST dataset in terms of NMI, ARI, and VM.

Measure	KPKM	KFCM	KKM
NMI	0.6830	0.5438	0.6613
ARI	0.5685	0.4047	0.5256
VM	0.6829	0.5411	0.6613

, with the examples of Figure 4 correspondingly being displayed in Figure 5. From Table 5, we can see that KPKM outperformed KFCM and KKM in terms of NMI, ARI, and VM. From Figure 5, we can observe that both KPKM and KKM found ten clusters, whereas KFCM found only seven clusters, although the number of clusters was set to ten.

5.4. Experiment for Descent Stability

In this experiment, we used 10 UCI datasets to evaluate descent stability of FAGP. AGP selects a small step length to converge. If AGP selects a large step length, the objective function value may descend with oscillation, and even does not converge. FAGP iteratively estimates a maximum step length at each iteration for speeding up its convergence. However, does it have any serious influence on the convergence? we ran the proposed FAGP on 10 UCI datasets to demonstrate the descent stability of FAGP. As shown in Figure 6, we can see that FAGP descends stably without oscillation at each iteration.

5.5. Convergence Speed Comparison between FAGP and AGP

In this experiment, we compared the convergence speeds of FAGP and AGP by using running time on 10 UCI datasets. $\eta$ is the ratio of FAGP’s running time to AGP’s. FAGP and AGP used the same initializations in each case. The results are presented in

Table 6. Running time (s) comparison of FAGP and AGP.

DS	FAGP	AGP	$\mathit{\eta }$
Ionosphere	0.719763	2.579700	27.90%
Iris	0.397728	9.410401	4.22%
Seeds	0.693613	8.650179	8.01%
Glass	3.482856	22.644511	15.38%
Segmentation	4.706316	19.125480	24.60%
Dermatology	10.659099	50.796225	20.98%
Breast	1.0153	6.1325	16.55 %
Natural	2309	–	–
Yeast	197.37	–	–
Waveform	213.64	–	–

and

Table 7. Comparison of the numbers of iterations used by FAGP and AGP.

DS	FAGP	AGP
Ionosphere	354	975
Iris	357	6506
Seeds	423	4061
Glass	1123	5616
Segmentation	1275	3986
Dermatology	1768	6784
Breast	682	694
Natural	16,378	–
Yeast	5609	–
Waveform	9939	–

(“–” means the running time was too long to obtain the final clustering results), and we can observe the proposed FAGP ran faster than AGP on all the 6 datasets. For Iris and Seeds datasets, FAGP used less than 10% of running time of AGP. FAGP also required fewer iterations than AGP. For large datasets, the running time of AGP was too long, but the proposed fast AGP could obtain the final clustering results.

6. Conclusions

In this paper, a novel clustering model (i.e., KPKM) was proposed. The proposed KPKM solves the problem of KFCM for m = 1, and this problem cannot be solved by existing method. The traditional AGP method can solve the proposed KPKM, but the efficiency of AGP is low. A fast AGP was proposed to speed up the AGP. The proposed fast AGP uses a maximum step length strategy to reduce the iteration number and uses an iterative method to update the projection matrix. The experimental results demonstrated that the fast AGP is able to solve the KPKM and the fast AGP requires less running time than AGP (the proposed FAGP requires 4.22–27.90% of the running time of AGP on real UCI datasets). The convergence of the proposed method was also analyzed by experiments. Additionally, in the experiments, the KPKM model could produce overall better clustering results than the other models, including KFCM, KKM, FCM, and KM. The proposed KPKM obtained the best clustering results on at least 6 real UCI datasets (a total of 10 real UCI datasets were used).

As future work, the proposed KPKM using other kernels will be evaluated in a variety of applications. For large datasets, the proposed method still has some disadvantages, so the next project will include speeding up the AGP on large datasets.