谱聚类

基本概念

定义无向图$G={V,E,W}$，其中，$V$为图的顶点，$E$为图的边，$W$为图中边的权值矩阵。

定义顶点的度为该顶点与其他顶点连接权值之和：

$$d_i = \sum_{j=1}^N w_{ij}$$

度矩阵$D$为：

$$D = \left [ \begin{matrix} d_1&&& \\ &d_2&& \\ &&\ddots& \\ &&&&d_N \end{matrix} \right ]$$

图的邻接矩阵（相似矩阵）$W$：

$$W=[w_{ij}]，1\leq i,j\leq N$$

对于$V$的一个子集$A\subset V$：

$$\begin{align} &V = \cup_{k=1}^K A_k\\ &A_i \cap A_j = \oslash，\forall i,j \in \{1,2,\cdots,N\} \\ &vol(A) = \sum_{i\in A} d_i \end{align}$$

相似矩阵

邻接矩阵$W$：由任意两点之间的权重值$w_{ij}$组成的矩阵。通常手动输入权重，但是在谱聚类中，只有数据点的定义，并没有直接给出这个邻接矩阵，那么如何构造邻接矩阵？

基本思想：距离较远的两个点之间的边权重值较低，而距离较近的两个点之间的边权重值较高，不过这仅仅是定性，需要定量的权重值。一般来说，我们可以通过样本点距离度量的相似矩阵$S$来获得邻接矩阵$W$。

构造方法：这里采用全连接法进行构造邻接矩阵。可以采用不同的核函数来定义边权重，这里采用高斯核函数RBF：

W_{ij} = S_{ij} = \left \{ \begin{array}{**lr**} exp(-\frac{||x_i-x_j||^2_2}{2\sigma^2})，& (i,j) \in E \\ 0，otherwise.& \end{array} \right.

拉普拉斯矩阵

拉普拉斯矩阵定义：

$$L = D - W$$

其中，$D$为度矩阵，$W$为邻接矩阵。

拉普拉斯矩阵的性质：

1. 拉普拉斯矩阵是对称矩阵，其所有特征值都是实数

2. 对于任意向量$f$，有：

$$f^TLf = \frac{1}{2} \sum_{i,j=1}^n w_{ij}(f_i-f_j)^2$$

​ 于是，可以得知：

$$\begin{align} f^TLf &= f^TDf - f^TWf \\ &= \sum_{i=1}^N d_i f_i^2 - \sum_{i,j=1}^N w_{ij}f_if_j \\ &= \frac{1}{2}(\sum_{i=1}^N d_i f_i^2 - 2\sum_{i,j=1}^N w_{ij}f_if_j+\sum_{j=1}^N d_j f_j^2) \\ &= \frac{1}{2}\sum_{i,j=1}^N w_{ij}(f_i-f_j)^2 \end{align}$$

1. 拉普拉斯矩阵是半正定的，且对应的n个实数特征值都大于等于0$(0=\lambda_1\leq \lambda_2 \leq \cdots \leq \lambda_N)$，且最小的特征值为0，对应的特征向量是单位向量。

由图分割问题到谱聚类问题

对于无向图$G$的切割图，目标是将图$G$切分成相互没有连接的$k$个子图，每个子图点的集合为：

$$\begin{align} &V = \cup_{k=1}^K A_k\\ &A_i \cap A_j = \oslash，\forall i,j \in \{1,2,\cdots,N\} \\ \end{align}$$

对于任意两个子图点的集合$A,B \subset V,A\cap B = \oslash$，定义$A$和$B$之间的切图之间的权重为：

$$W(A,B) = \sum_{i\in A,j\in B} w_{ij}$$

对于k个子图点的集合，定义切图$cut$为：

$$\begin{align} cut(V) &= cut(A_1,A_2,\cdots,A_k) \\ &= \sum_{k=1}^K W(A_k,\bar{A_k}) \\ &= \sum_{k=1}^K W(A_k,V) - W(A_k,A_k) \end{align}$$

其中，$\bar{A_i}$为$A_i$的补集（除去$A_i$子集以外其他的属于$V$的结点的集合）。

为了让切图可以让子图内的点权重和高(内部连接强)，子图间的点权重和低(外部连接弱)，需要使$cut$达到最小化。

优化目标：

$$\begin{align} \mathop{max}\limits_{\{A_k\}^K_{k=1}} cut(V) = \sum\limits_{i=1}^K W(A_i,\bar{A_i}) \end{align}$$

为了避免最小切图导致的切图效果不佳，需要对每个子图的规模做出限定。这里得到了一个新的优化目标：

$$Ncut(A_1,A_2,\cdots,A_k) = \frac{1}{2} \sum\limits_{i=1}^K \frac{W(A_i,\bar{A_i})}{vol(A_i)} = \frac{1}{2} \sum\limits_{i=1}^K \frac{W(A_i,\bar{A_i})}{\sum\limits_{i\in A_k}d_i}，d_i=\sum_{j=1}^Nw_{ij}$$

由上可知，算法的目标函数为：

$$\{\hat{A_k} \}_{k=1}^K = \mathop{argmin}_{\{A_k\}^K_{k=1}} \sum\limits_{k=1}^K\frac{W(A_k,v)-W(A_k,A_K)}{\sum\limits_{i\in A_k}\sum\limits_{j=1}^Nw_{ij}}$$

由于求解的时候，不方便计算，于是设法将目标函数转换为矩阵的形式表达。

引入指示向量

令

\left \{ \begin{array}{**lr**} y_i \in \{0,1\}^k,&\\ \sum\limits_{j=1}^K y_{ij} = 1，1 \leq i,j \leq N \end{array} \right.

其中，$y_i$是一个$k$维的向量，每一个维度的值为0或者1；$y_{ij}$表示第$i$个样本属于第$j$个类别，并且每个样本只能属于一个类别。

由上，可知目标函数为：

$$\begin{align} &Y = (y_1,y_2,\cdots,y_N)_{N\times K}^T \\ &\hat Y = \mathop{argmin}\limits_{Y} Ncut(V) = \mathop{argmin}\limits_Y \sum\limits_{k=1}^K \frac{W(A_k,\bar{A_k})}{\sum\limits_{i \in A_k}d_i} \end{align}$$

又，

$$\begin{align} \sum\limits_{k=1}^K \frac{W(A_k,\bar{A_k})}{\sum\limits_{i \in A_k}d_i} &= tr \left [ \begin{matrix} \frac{W(A_1,\bar{A_1})}{\sum\limits_{i \in A_1}d_i}&&& \\ &\frac{W(A_2,\bar{A_2})}{\sum\limits_{i \in A_2}d_i}&& \\ &&\ddots& \\ &&&&\frac{W(A_K,\bar{A_K})}{\sum\limits_{i \in A_K}d_i} \end{matrix} \right ] \\ &= tr \underbrace{\left [ \begin{matrix} W(A_1,\bar{A_1})&&& \\ &W(A_2,\bar{A_2})&& \\ &&\ddots& \\ &&&&W(A_K,\bar{A_K}) \end{matrix} \right ]_{K\times K}}_{O} \cdot \left[\underbrace{\left [ \begin{matrix} \sum\limits_{i \in A_1}d_i&&& \\ &\sum\limits_{i \in A_2}d_i&& \\ &&\ddots& \\ &&&&\sum\limits_{i \in A_K}d_i \end{matrix} \right ]_{K\times K}}_{P}\right]^{-1} \\ = tr(OP^{-1}) \end{align}$$

现在已知$W$，$Y$，求$O$，$P$，注意$Y$矩阵是$N\times K$维的，行代表结点标号，列代表属于的类别，其中$y_{ij}$表示第$i$个样本属于第$j$个类别。

有，

$$\begin{align} y_iy_i^T&= \left[\begin{matrix}0,0,\cdots,y_{ij}=1,\cdots,0\end{matrix}\right]\cdot \left[\begin{matrix}0\\0\\\cdots\\y_{ij}=1\\\cdots\\0\end{matrix}\right] \\ &= \left[ \begin{matrix} 0&\cdots&0&\cdots&0 \\ \vdots&\ddots&\vdots&\ddots&\vdots\\ 0&\cdots&y_{jj}=1&\cdots&0 \\ \vdots&\ddots&\vdots&\ddots&\vdots\\ 0&\cdots&0&\cdots&0 \\ \end{matrix} \right] \end{align}$$

于是，

$$\begin{align} Y^TY &= (y_1,y_2,\cdots,y_N)\cdot (y_1^T,y_2^T,\cdots,y_N^T) \\ &=\sum\limits_{i=1}^N y_iy_i^T \\ &= \left[ \begin{matrix} N_1&&&&\\ &N_2&&&\\ &&\ddots&\\ &&&&N_K\\ \end{matrix} \right]_{K\times K} \end{align}$$

其中，$N_K$表示在$N$个样本中，属于类别$k$的样本个数，且$\sum\limits_{k=1}^KN_k = N$，$N_k = |A_k|=\sum\limits_{i\in A_k} 1_N$，$1_N$表示元素全为$1$的$N\times 1$维向量。

有，

$$P = \sum\limits_{i=1}^Ny_i^Td_iy_i = Y^TDY = Y^T(W\cdot1_N)Y$$

使用拉普拉斯矩阵的性质，简化$O$矩阵，如下：

$$\begin{align} O &= \left [ \begin{matrix} W(A_1,\bar{A_1})&&& \\ &W(A_2,\bar{A_2})&& \\ &&\ddots& \\ &&&&W(A_K,\bar{A_K}) \end{matrix} \right ] \\ &= \left [ \begin{matrix} W(A_1,V)&&& \\ &W(A_2,V)&& \\ &&\ddots& \\ &&&&W(A_K,V) \end{matrix} \right ] - \left [ \begin{matrix} W(A_1,A_1)&&& \\ &W(A_2,A_2)&& \\ &&\ddots& \\ &&&&W(A_K,A_K) \end{matrix} \right ] \\ &= \left [ \begin{matrix} \sum\limits_{i \in A_1} d_i&&& \\ &\sum\limits_{i \in A_2} d_i&& \\ &&\ddots& \\ &&&&\sum\limits_{i \in A_K} d_i \end{matrix} \right ] - \left [ \begin{matrix} W(A_1,A_1)&&& \\ &W(A_2,A_2)&& \\ &&\ddots& \\ &&&&W(A_K,A_K) \end{matrix} \right ] \\ \end{align}$$

由于，

$$\begin{align} Y^TWT &= \left[ \begin{matrix} y_1 &y_2&\cdots&y_N \end{matrix} \right] \left[ \begin{matrix} w_{11}&\cdots&w_{1N} \\ \vdots&\ddots&\vdots \\ w_{N1}&\cdots&w_{NN} \\ \end{matrix} \right] \left[ \begin{matrix} y_1^T\\y_2^T\\\cdots\\y_N^T \end{matrix} \right] \\ &= \left[ \begin{matrix} \sum\limits_{i=1}^N y_1w_{i1} &\sum\limits_{i=1}^N y_2w_{i2}&\cdots&\sum\limits_{i=1}^N y_Nw_{iN} \end{matrix} \right] \left[ \begin{matrix} y_1^T\\y_2^T\\\cdots\\y_N^T \end{matrix} \right] \\ &=\sum\limits_{i=1}^N \sum\limits_{j=1}^N y_iw_{ij}y_j^T = \sum\limits_{i=1}^N \sum\limits_{j=1}^N y_iy_j^Tw_{ij} \\ &= \left[ \begin{matrix} \sum\limits_{i\in A_1,j\in A_1}w_{ij}&\sum\limits_{i\in A_1,j\in A_2}w_{ij}&\cdots& \sum\limits_{i\in A_1,j\in A_N}w_{ij} \\ \vdots&\vdots&\ddots&\vdots \\ \sum\limits_{i\in A_N,j\in A_1}w_{ij}&\sum\limits_{i\in A_N,j\in A_2}w_{ij}&\cdots& \sum\limits_{i\in A_N,j\in A_N}w_{ij} \end{matrix} \right] \end{align}$$

令，

$$O^{'} = Y^TDY - Y^WY$$

由$trace$性质可知，

$$tr(O^{'}P^{-1}) = tr(OP^{-1})$$

其中，$O^{'}$与$O$的对角元素一样。

于是，有

$$\hat Y = \mathop{argmin}_{Y} \ tr(OP^{-1}) = \mathop{argmin}_{Y} \ tr[Y^TLY\cdot (Y^TDY)^{-1}]$$

定义$H$为新的指示向量：

h_{ij}=\left \{ \begin{array}{**lr**} 0&v_i \not\in A_j&\\ \frac{1}{\sum\limits_{i\in A_j} d_i}&v_i \in A_j \end{array} \right.

有，

$$H^TDH = I$$

于是，

$$\hat H = \mathop{argmin}_{H}\ tr(H^TLH)，s.t.H^TDH=I$$

记，$H=D^{-\frac{1}{2}}F$，有$F^TF=I$，

$$\hat F = \mathop{argmin}_F \ tr(F^T{D^{-\frac{1}{2}}}L{D^{-\frac{1}{2}}}F)$$

注意$Y^TY$与$H^TH$表示的含义，两者的主对角线的值是相反的。

谱聚类实现

• 实现步骤

image-20210107164311389

image-20210107164257045

• 导入包

import warnings
import numpy as np
import matplotlib.pyplot as plt
from sklearn import preprocessing
from sklearn.cluster import KMeans
from sklearn.metrics import accuracy_score,recall_score
from sklearn.datasets.samples_generator import make_blobs
warnings.filterwarnings('ignore')

• 谱聚类

class SpectralClustering:
    def __init__(self,n_clusters=None,sigma=None,normalize=0):
        """
        
        """
        self.n_clusters = n_clusters
        self.W = None
        self.N = None # 样本数量/结点数量
        self.D = None # 度矩阵
        self.L = None # 拉普拉斯矩阵
        self.sigma = sigma # 高斯函数的sigma
        
        self.score = None
        
        self.normalize = normalize
        
        if self.n_clusters == None:
            self.n_clusters = 2
        
        if self.sigma == None:
            self.sigma = 3
    
    def init_params(self,X,y):
        self.X,self.y = X,y
        self.N = X.shape[0]
        self.W = self.similaryMatrix(X)
        self.D = self.diagMatrix()
        self.L = self.laplacianMatrix()
    
    def diagMatrix(self,):
        """
        对角矩阵
        """
        diag = np.sum(self.W,axis=1)
        diag = np.squeeze(diag)
        return np.diag(diag)
    
    def kernel(self,num,x,y):
        if num == 1: # 选择高斯核函数
            distance  = ((x - y)@(x - y))
            return np.exp(-distance/(2*self.sigma*self.sigma)) # 高斯函数
    
    def similaryMatrix(self,data):
        """
        相似矩阵/邻接矩阵/权值矩阵，权重值为距离的平方的倒数，表示越近的结点权重越大，越远的结点权重越小
        """
        S = np.zeros((self.N,self.N))
        for i in range(self.N):
            for j in range(self.N):
                S[i,j] = S[j,i] = self.kernel(1,data[i],data[j])
        S = S - np.diag(np.diag(S))
        return S
    
    def laplacianMatrix(self):
        L = self.D - self.W
        if self.normalize == 1:
            return np.sqrt(np.linalg.pinv(self.D)).dot(L).dot(np.sqrt(np.linalg.pinv(self.D)))
        return L
    
    def getEigenvectors(self,):
        eigenvalue,eigenvector = np.linalg.eig(self.L)
        index = np.argsort(eigenvalue)[:self.n_clusters]
        eigenvector = eigenvector[:,index]
        if self.normalize == 1:
#             eigenvector = preprocessing.StandardScaler().fit(
							eigenvector).transform(eigenvector) # 正则化
            eigenvector = preprocessing.normalize(eigenvector, norm='l2')
        return eigenvector
        
    def fit(self,X,y):
        self.init_params(X,y)
        vector = self.getEigenvectors()
        self.kmeans = KMeans(n_clusters=self.n_clusters)
        self.kmeans.fit(vector)
        self.y_pred = self.prediction()
        self.score = self.calc_score()
        self.plot()
        
    def prediction(self,):
        return self.kmeans.labels_
    
    def plot(self,):
        plt.scatter(X[:,0],X[:,1],c=self.y_pred)
        
    def calc_score(self,):
        a = accuracy_score(1-self.y_pred,self.y)
        b = accuracy_score(self.y_pred,self.y)
        return a if a>b else b

• 生成数据

def create_data():
    X,y = make_blobs(n_samples=500,centers=2,random_state=4)
    return X,y
X,y = create_data()
plt.scatter(X[:,0],X[:,1],c=y)

image-20210107164133043

• 运行

sc = SpectralClustering(n_clusters=2,normalize=1)
sc.fit(X,y)
sc.score

image-20210107164216858