谱聚类

1. 信息增益的度量

由于数据集 $X=[{\mathbf{x}}_1,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_n]$ 的图结构不像图像的像素有像素值, 数据样本没有规整的结构, 也没有特定的函数 $f$, 因此, 我们需要为每个样本赋予一个函数值 { f i } i = 1 n \{f_i\}_{i=1}^n {fi​}i=1n​, 或表示成函数向量 $\mathbf{f}$. 为获得整个数据集的凸凹性(或差异性)

这种差异性可以通过拉普拉斯变换表示为
$$\mathbf{h}=L\mathbf{f}=(D-W)\mathbf{f}=D\mathbf{f}-W\mathbf{f}$$
则 $\mathbf{h}$ 的第 $i$ 个 分量 $h_i$ 可表示为
$$h_i=d_i{f}_i-w_{i,:}\mathbf{f}=\sum _{j\in N_i}{w}_{ij}{f}_i-\sum _{j\in N_i}{w}_{ij}{f}_j=\sum _{j\in N_i}{w}_{ij}(f_i-f_j)$$
即, $h_i$ 是第 $i$ 个顶点 $v_i$ 与其近邻点的差的和. 若构造一个量来表示所有顶点与其近邻点的差异的总和, 则可以由以下几种方式:

• 第一种(计算困难): ${\sum }_{i=1}^n{h}_i^2={\sum }_{i=1}^n({\sum }_{j\in N_i}{w}_{ij}(f_i-f_j){)}^2$
• 第二种(正负值会抵消): ${\sum }_{i=1}^n{h}_i={\sum }_{i=1}^n{\sum }_{j\in N_i}{w}_{ij}(f_i-f_j)$
• 第三种(非负求和): ${\sum }_{i=1}^n{f}_i\cdot h_i={\sum }_{i=1}^n{f}_i{\sum }_{j\in N_i}{w}_{ij}(f_i-f_j)$

综合上述三种情况, 第三种度量总的差异是最适合的.

$$\begin{array}{ll}{\sum }_{i=1}^n{f}_i\cdot h_i& ={\sum }_{i=1}^n{f}_i{\sum }_{j\in N_i}{w}_{ij}(f_i-f_j)\\ & \\ & ={\sum }_{i=1}^n{\sum }_{j\in N_i}{w}_{ij}{f}_i(f_i-f_j)\\ & \\ & =\frac{1}{2}{\sum }_{i=1}^n{\sum }_{j\in N_i}{w}_{ij}(2f_i{f}_i-2f_i{f}_j)\\ & \\ & =\frac{1}{2}{\sum }_{i=1}^n{\sum }_{j\in N_i}{w}_{ij}(f_i{f}_i-2f_i{f}_j+f_j{f}_j)\\ & \\ & =\frac{1}{2}{\sum }_{i=1}^n{\sum }_{j\in N_i}{w}_{ij}(f_i-f_j{)}^2\end{array}$$

因此, 数据集信息的总增益的度量可表示为

$$\begin{array}{l}{\sum }_{i=1}^n{f}_i\cdot h_i=\frac{1}{2}{\Sigma }_{i=1}^n{\Sigma }_{j\in N_i}{w}_{ij}(f_i-f_j{)}^2\\ \\ =\frac{1}{2}{\Sigma }_{i=1}^n{\Sigma }_{j\in N_i}{w}_{ij}(f_i^2+f_j^2-2f_i{f}_j)\\ \\ ={\Sigma }_{i=1}^n{\Sigma }_{j\in N_i}{w}_{ij}{f}_i^2-{\Sigma }_{i=1}^n{\Sigma }_{j\in N_i}{w}_{ij}{f}_i{f}_j\\ \\ ={\Sigma }_{i=1}^n(w_{i1}{f}_1^2+w_{i2}{f}_2^2+\cdots +w_{in}{f}_n^2)-{\mathbf{f}}^{\top }W\mathbf{f}\\ \\ =(f_1,f_2,\cdots ,f_n)\left(\begin{array}{cccc}{d}_1& & & \\ & d_2& & \\ & & \ddots & \\ & & & d_n\end{array}\right)\left(\begin{array}{c}{f}_1\\ f_2\\ ⋮\\ f_n\end{array}\right)-(f_1,f_2,\cdots ,f_n)\left(\begin{array}{cccc}{w}_{11}& w_{12}& \cdots & w_{1n}\\ w_{21}& w_{22}& \cdots & w_{2n}\\ \cdots & \cdots & \cdots & \cdots \\ w_{n1}& w_{n2}& \cdots & w_{nn}\end{array}\right)\left(\begin{array}{c}{f}_1\\ f_2\\ ⋮\\ f_n\end{array}\right)\\ \\ ={\mathbf{f}}^{\top }D\mathbf{f}-{\mathbf{f}}^{\top }W\mathbf{f}\\ \\ ={\mathbf{f}}^{\top }(D-W)\mathbf{f}\\ \\ ={\mathbf{f}}^{\top }L\mathbf{f}\end{array}$$

综上可知, 对于离散样本点的函数 $\mathbf{f}=(f_1,f_2,\cdots ,f_n{)}^{\top }$, $w_{ij}$ 为两个样本点之间的权重, 则数据集上的信息总增益为

$${\mathbf{f}}^{\top }L\mathbf{f}=\frac{1}{2}{\Sigma }_{i=1}^n{\Sigma }_{j\in N_i}{w}_{ij}(f_i-f_j{)}^2$$

2. 谱聚类: 寻找最优的函数向量 $\mathbf{f}$

2.1 : 寻找一个最优的函数向量 $\mathbf{f}$

为数据集的每个样本点寻找一个合适的函数 $f_i$, 完成函数的求值 $\{f_i({\mathbf{x}}_i){\}}_{i=1}^n$. 为使得样本点之间可以在同一尺度下进行大小的比较, 我们将函数值限定在0,1之间, 则可转化为约束 优化问题
$$\begin{gathered} \underset{\mathbf{f}}{\min }{\mathbf{f}}^{\top }L\mathbf{f} \\ s.t{\mathbf{f}}^{\top }\mathbf{f}=1 \end{gathered}$$

利用拉格朗日乘子法将等式约束优化问题转化成无约束优化问题
$$\underset{\mathbf{f}}{\min }Q(\mathbf{f})={\mathbf{f}}^{\top }L\mathbf{f}-\lambda {\mathbf{f}}^{\top }\mathbf{f}$$
求极值
$$\frac{\partial Q(\mathbf{f})}{\partial \mathbf{f}}=2L\mathbf{f}-2\lambda \mathbf{f}=2(L\mathbf{f}-\lambda \mathbf{f})=0$$

即
$$L\mathbf{f}=\lambda \mathbf{f}(转化为求拉普拉斯特征向量问题)$$

这是拉普拉斯矩阵 $L$ 的特征方程. 由于拉普拉斯矩阵的行和为 $0$, 上式有一个平凡解 $\lambda =0$, 其对应的特征向量为 $1$. 此解显然与数据集无关, 不是优化问题的最优解. 又因为拉普拉斯矩阵是半正定矩阵, 特征值非负, 因此, $\mathbf{f}$ 是拉普拉斯矩阵第二小的特征值对应的特征向量.

将上式两边左乘${\mathbf{f}}^{\top }$知
$$min{\mathbf{f}}^{\top }L\mathbf{f}\iff min\lambda {\mathbf{f}}^{\top }\mathbf{f}=min\lambda ={\lambda }_{min}$$
即
$$\underset{\mathbf{f}}{\min }\frac{{\mathbf{f}}^{\top }L\mathbf{f}}{{\mathbf{f}}^{\top }\mathbf{f}}={\lambda }_{min}$$
称为瑞利熵

2.2 寻找鲁棒性更强的多个函数向量

为获得更加鲁棒的结果, 可以寻找多个向量函数, 然后进行信息的融合. 则数据集上的信息总增益表示为

$$\underset{{\mathbf{f}}_1,\cdots ,{\mathbf{f}}_{\mathbf{k}}}{\min }{\Sigma }_{i=1}^k{{\mathbf{f}}_{\mathbf{i}}}^{\top }L{\mathbf{f}}_{\mathbf{i}}$$

$$\begin{gathered} {\Sigma }_{i=1}^k{{\mathbf{f}}_{\mathbf{i}}}^{\top }L{\mathbf{f}}_{\mathbf{i}}=\underset{{\mathbf{f}}_1,\cdots ,{\mathbf{f}}_{\mathbf{k}}}{\min }\left(\begin{array}{c}{\mathbf{f}}_1^{\top }\\ ⋮\\ {\mathbf{f}}_k^{\top }\end{array}\right)\left(\begin{array}{cccc}{L}_{11}& L_{12}& \cdots & L_{1n}\\ L_{21}& L_{22}& \cdots & L_{2n}\\ \cdots & \cdots & \cdots & \cdots \\ L_{n1}& L_{n2}& \cdots & L_{nn}\end{array}\right)\left(\begin{array}{c}{\mathbf{f}}_1\\ ⋮\\ {\mathbf{f}}_k\end{array}\right) \\ =\underset{\mathbf{F}}{\min }tr({\mathbf{F}}^{\top }L\mathbf{F}) \end{gathered}$$

结论：此问题等价于求多个较小的特征向量。

2.3 谱聚类(spectral clustering)算法

通过构造连接矩阵的方式获取拉普拉斯矩阵, 然后进行最优化求解

• 图拉普拉斯版本, 此版本效率有点儿低

class spectralClust_graph:
    
    def __init__(self, nClust=2, gamma=13.5, tau=0.1):
        
        self.nClust = nClust # 初始化类数
        self.gamma = gamma # 径向基核函数参数
        self.tau = tau # 近邻半径参数
    
    # 计算距离矩阵    
    def pairwise_distances(self, X):
        n = X.shape[1]
        G = X.T@X
        H = np.diag(G).reshape(-1,1)@np.ones((1,n))
        dist = H+H.T-2*G
        Dist = np.sqrt(dist)
        print(Dist.shape)
        return Dist
        return Dist
    
    # 计算权重矩阵
    def create_graph_weights(self, X, gamma, tau):
        
        # YOUR CODE HERE
        distance_matrix = self.pairwise_distances(X.T)
        n = distance_matrix.shape[0]
        graph_weights = []
        weights = np.exp(-gamma*distance_matrix)
        weights = weights*(weights>=tau)
        for i in range(n):
            for j in range(n):
                if i != j:
                    graph_weights.append(weights[i,j])
        graph_weights = np.array(graph_weights)
        return graph_weights
    
    # 计算连接矩阵
    def construct_incidence_matrix(self, X):
        
        weights = self.create_graph_weights(X, self.gamma, self.tau)
        
        no_of_samples = X.shape[0]
        # 为数据集连接边分配索引
        edges = []
        for i in range(no_of_samples):
            for j in range(no_of_samples):
                if i == j:
                    continue
                else:
                    edges.append([i,j])
        
        no_of_edges = len(edges)
        incidence_matrix = np.zeros(shape=(no_of_edges,no_of_samples))
        for index in range(no_of_edges):
            indices = edges[index]
            incidence_matrix[index, indices[0]] = -np.sqrt(weights[index])
            incidence_matrix[index, indices[1]] = np.sqrt(weights[index])
        id0 = (incidence_matrix==0).all(1)
        incidence_matrix = np.delete(incidence_matrix,id0,axis=0)
        return incidence_matrix
        
    
    def fit(self, X):
        incidence_matrix = self.construct_incidence_matrix(X)
        graph_laplacian = incidence_matrix.T@incidence_matrix
        eigenvalues, eigenvectors = linalg.eigs(graph_laplacian,k=self.nClust,which='SM')
        Labels = k_means(eigenvectors.real,self.nClust)
        
        return Labels[1]

import numpy as np
from scipy.sparse import linalg
from matplotlib.colors import ListedColormap
from sklearn.datasets import make_moons
from sklearn.cluster import k_means
import matplotlib.pyplot as plt

if __name__ == "__main__":
    # 构造数据集
    seed = 13
    np.random.seed(seed)
    no_of_samples = 1000
    X, y = make_moons(n_samples=no_of_samples, noise=0.1, random_state=seed)
    
    scg = spectralClust_graph()
    unsupervised_labels = scg.fit(X)
    
    # 画图
    colormap_bright = ListedColormap(['#FF0000', '#0000FF']) # 设置颜色
    plt.figure(figsize=(12,6))
    plt.subplot(121)
    plt.scatter(X[:, 0], X[:, 1])
    plt.xticks(fontsize=16)
    plt.yticks(fontsize=16)
    
    plt.subplot(122)
    plt.scatter(X[:, 0], X[:, 1], c=unsupervised_labels, cmap=colormap_bright, edgecolors='k')
    plt.xticks(fontsize=16)
    plt.yticks(fontsize=16)
    
    plt.show()    

• 通用的普聚类版本
  也可以通过度矩阵构造拉普拉斯矩阵 $L=D-W$, 然后进行优化.

class spectralCLUST:
    
    def __init__(self, nClust=2, gamma=13.5, tau=0.1, affinity=None):
        
        self.nClust = nClust # 初始化类数
        self.gamma = gamma # 径向基核函数参数
        self.tau = tau # 近邻半径参数
        self.affinity = affinity
    
    # 计算距离矩阵, 矩阵的每一列为一个样本点    
    def pairwise_distances(self, X):
        n = X.shape[1]
        G = X.T@X
        H = np.diag(G).reshape(-1,1)@np.ones((1,n))
        dist = H+H.T-2*G
        Dist = np.sqrt(dist)
        print(Dist.shape)
        return Dist
    
    # 计算拉普拉斯矩阵,矩阵的每一列为一个样本点 
    def create_Weights(self, X, gamma, tau):
        
        #
        distance_matrix = self.pairwise_distances(X)
        n = distance_matrix.shape[0]
        weights = np.exp(-gamma*distance_matrix)
        Weights = weights*(weights>=tau)
        return Weights
    
    def clustering(self,CKSym):
        N = CKSym.shape[1]
        n = self.nClust
        DN = np.diag(np.divide(1, np.sqrt(np.sum(CKSym, axis=0) + np.finfo(float).eps)))
        LapN = identity(N).toarray().astype(float) - np.matmul(np.matmul(DN, CKSym), DN)
        _, _, vN = np.linalg.svd(LapN)
        vN = vN.T
        kerN = vN[:, N - n:N]
        normN = np.sqrt(np.sum(np.square(kerN), axis=1))
        kerNS = np.divide(kerN, normN.reshape(len(normN), 1) + np.finfo(float).eps)
        km = KMeans(n_clusters=n).fit(kerNS)
        return km.labels_
        
    # 拟合函数,矩阵 X 的每一列为数据点
    def fit(self, X): 
        if self.affinity == None:
            gamma = self.gamma
            tau = self.tau
            W = self.create_Weights(X, gamma, tau)
            Labels = self.clustering(W)
        if self.affinity == 'precomputed':
            Labels = self.clustering(X)
        
        return Labels

import numpy as np
from scipy.sparse import linalg
from matplotlib.colors import ListedColormap
from sklearn.datasets import make_moons
from sklearn.cluster import KMeans
from scipy.sparse import identity
import matplotlib.pyplot as plt
    
if __name__ == "__main__":
    # 构造数据集
    seed = 13
    np.random.seed(seed)
    no_of_samples = 1000
    X, y = make_moons(n_samples=no_of_samples, noise=0.1, random_state=seed)
    
    clt = spectralCLUST()
    unsupervised_labels = clt.fit(X.T)    
    # 画图
    colormap_bright = ListedColormap(['#FF0000', '#0000FF']) # 设置颜色
    plt.figure(figsize=(12,6))
    plt.subplot(121)
    plt.scatter(X[:, 0], X[:, 1])
    plt.xticks(fontsize=16)
    plt.yticks(fontsize=16)
    
    plt.subplot(122)
    plt.scatter(X[:, 0], X[:, 1], c=unsupervised_labels, cmap=colormap_bright, edgecolors='k')
    plt.xticks(fontsize=16)
    plt.yticks(fontsize=16)
    
    plt.show()

小结

第二种编程的方法可以实现通用的功能, 既可以处理样本数据集, 同时也可以处理权重矩阵, 以备后续章节使用