谱聚类算法详解及代码实现

谱聚类算法详解及代码实现

参考

• https://www.163.com/dy/article/FQA23T6K0516EPQ9.html
• 文章1
• 图论绘图辅助工具

关于谱聚类介绍

谱聚类概述

谱聚类是从图论中演化出来的算法，在后来我们的聚类中得到了广泛的应用。闭气传统的K-Means算法，谱聚类对数据分布的适应性更强，聚类效果也很不错。其主要思想是将所有的数据看作是空间中的点，这些点之间可以用边连接起来。距离较远的两个点之间的边权重值较低，而距离较近的两个点之间的边权重值较高，通过对所有数据点组成的图进行切图，让切图后不同的子图间边的权重和尽可能低，而子图内的边权重和尽可能高，从而达到聚类的效果。
优点

• 谱聚类只需要数据之间的相似度矩阵，因此对于处理稀疏数据的聚类很有效
• 由于使用了降维，因此在处理高维数据聚类时的复杂度比传统聚类算法好
  缺点
• 聚类效果依赖于相似矩阵，不同的相似矩阵得到的最终聚类效果可能很不同。

谱聚类前置知识

无向权重图

对于一个图，我们一般用电的集合V和边的集合E来描述。即为G(V,E)。其中V即为我们数据集里面所有的点(v1,v2,……,vn)。对于V中的任意两个点vi和vj,我们定义权重wij为二者之间的权重。由于是无向图，所以wij=wji。
$$W=\left[\begin{array}{cccc}{w}_{11}& w_{12}& w_{13}& w_{14}\\ w_{21}& w_{22}& w_{23}& w_{24}\\ w_{31}& w_{32}& w_{33}& w_{34}\\ w_{41}& w_{42}& w_{43}& w_{44}\end{array}\right]\text{\hspace{1em}(1)}$$

邻接矩阵

邻接矩阵：是表示节点之间相邻关系的矩阵。如下图中的权重图中（假设各权重为1），其邻接矩阵可表示为事（2）所示。

图 一

$$W=\left(\begin{array}{ccccc}0& 0& 0& 1& 1\\ 0& 0& 1& 1& 1\\ 0& 1& 0& 0& 0\\ 1& 1& 0& 0& 0\\ 1& 1& 0& 0& 0\end{array}\right)\text{\hspace{1em}(2)}$$

度矩阵

度矩阵是对角矩阵，对角上元素为各个顶点的度。上述图一的度矩阵为式(3)所示。
$$D=\left(\begin{array}{ccccc}2& 0& 0& 0& 0\\ 0& 3& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 2& 0\\ 0& 0& 0& 0& 2\end{array}\right)\text{\hspace{1em}(3)}$$

拉普拉斯矩阵

拉普拉斯矩阵L= D - W，其中D为度矩阵， W为邻接矩阵。图一的拉普拉斯矩阵如式(4)所示。
$$L=\left(\begin{array}{ccccc}2& 0& 0& -1& -1\\ 0& 3& -1& -1& -1\\ 0& -1& 1& 0& 0\\ -1& -1& 0& 2& 0\\ -1& -1& 0& 0& 2\end{array}\right)\text{\hspace{1em}(4)}$$
用拉普拉斯矩阵来求解特征值，通过确定特征值来确定对应特征向量的个数，从而实现降维，然后再用K-means将降维的特征向量进行聚类。

相似度矩阵

一般通过K-近邻的方法构造相似度矩阵，一般是基于欧式距离进行构造的。其步骤一般分为以下两步：

1. ✔ 计算点对之间的欧式距离
2. 通过给定的参数k进行构造

其中第二步通常有以下三种方法实现

方法一

策略：通过给定参数k，选取距离当前点最近的k个点作为邻居（通常用高级核计算距离），令其余点到该点的距离为0。

$$wij=\{\begin{array}{ll}{e}^{\frac{{\Vert \begin{array}{cc}{x}_i& x_j\end{array}\Vert }_2^2}{2{\sigma }^2}},& \text{if xi ∈ knn(xj)}\\ 0,& \text{otherwise}\end{array}\text{\hspace{1em}(方法1)}$$

这样构造的问题：数据从有向图变成了无向图，即B是A的k近邻，但A不一定是B的k近邻。为了保证相似矩阵的对称性，有了下面两种构造方法。

方法二

策略：与方法一类似，不过在选取k近邻的时候，只要A是B的k近邻了我们也罢B当作A的k近邻（不管B是不是A的k近邻）,可以类似理解成或运算。

$$wij=wji=\{\begin{array}{ll}{e}^{\frac{\Vert {\Vert \begin{array}{cc}{x}_i& x_j\end{array}\Vert }_2^2}{2{\sigma }^2}},& \text{if xi ∈ knn(xj) or xj ∈ knn(xi) }\\ 0,& \text{otherwise}\end{array}\text{\hspace{1em}(方法2)}$$

方法三

策略：与方法一类似，不过在选取k近邻的时候，要保证A和B互为k近邻的时候这两个点的边才有权重。

$$wij=wji=\{\begin{array}{ll}{e}^{\frac{{\Vert \begin{array}{cc}{x}_i& x_j\end{array}\Vert }_2^2}{2{\sigma }^2}},& \text{if xi ∈ knn(xj) and xj ∈ knn(xi) }\\ 0,& \text{otherwise}\end{array}\text{\hspace{1em}(方法3)}$$

最常用的方法一般是：

将方法一算的的相似度矩阵与它的转置矩阵相加再除以2，这样就保证了求得的相似度矩阵也是对称矩阵。

$$W=\frac{W+W^T}{2}$$

确定目标函数

在确定好图的结构，问题就是将点聚成c个类的问题，可以将问题转化为切割无向图为c个子图。

根据上图，我们很容易想到一个方法，就是我们将距离较远的两个点切分到不同子图时，需要付出的代价很小。因此我们可以定义一个代价函数来作为我们初步的目标函数。

初始化目标函数（最小割MinCut方法）

对于无向图G(X,E) 进行切分的目标时将G划分为相互无连接的k个子图，每个子图包含点的集合$A_1$，$A_2$…$A_k$,.(当然要保证每个点属于且只属于一个子图)。对于任意的两个子图点的集合A，B，可以定义A和B之间的切图权重为:
$$W(A,B)=\sum _{i\in A,j\in B}{w}_{ij}\text{\hspace{1em}(A-B 切图权重)}$$
那么我们可以把上述的切图权重公式扩展到k个子图的节点的集合种，即$A_1$,$A_2$…$A_k$，定义切图cut为:
$$cut(A_1,A_2,A_3...A_k)=\frac{1}{2}\sum _{i\in A,j\in B}^{k}W(A_i,\overline{A_i})$$

这里除以2是因为按照求和公式我们会计算A和B之间的切图权重也会计算B和A之间的切图权重，而根据切图权重计算公式$W(A,B)=W(B,A)$,所以将最后结果除以2处理。

其中$\overline{A_i}$ 为$A_i$的补集（除了本身以外的所有子集的并集）。

因此我们可以进一步得到初步离散优化问题，最小切割目标函数:
$$minCut(V)=>min\sum _{v_i\in A_k,v_j\in \overline{A_k}}{w}_{ij}$$

我们先只考虑只有两个类的情况（c = 2)

通过从 c = 2 的情况理解引入指示向量的作用。

在这里我们定义只是向量$f=(f_1,f2,...f_n{)}^T\in N^n$
$$f_i=\{\begin{array}{ll}1& \text{if vi ∈ A}\\ & \\ 0& \text{if vi ∈ A}\end{array}$$
因此目标函数转化为：
$$min\sum _{v_i\in A_k,v_j\in \overline{A_k}}{w}_{ij}⟺min\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}(f_i-f_j{)}^2$$ {$}

可以重点理解引入指示向量前后式子两边的等价性。

我们将 $\frac{1}{2}{\sum }_{i=1}su{m}_{j=1}{w}_{ij}(f_i-f_j{)}^2$ 进行展开：

其中,D为度矩阵，W为邻接矩阵（相似矩阵），而$L=D-W$ 为拉普拉斯矩阵，因此原目标函数等价于：
$$\underset{f_i\in N^n}{argmin}Cut(V)=\underset{f_i\in N^n}{argmin}\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}(f_i-f_j{)}^2=\underset{f_i\in N^n}{argmin}{f}^{T}Lf$$
上面我们通过引入布尔指示向量将一个抽象的切割函数$Cut(V)$转化为了就提可求解的$f^{T}Lf$.但是这个是离散的优化函数，无法进行求导。

推广到多个类（c取任意值）

对于聚多个类的问题，我们要聚c个类$(k=1,2,....,c)$, 因此我们需要将原来的一个指示向量扩展为c个指示向量组成的指示矩阵$H$，其余步骤都一样。

首先我们定义向量$h_k=(h_{(1,k),h_{2,k},...,h_{n,k}}{)}^T\in N^n$（其中：$k = 1,2, …,c $）并满足以下公式：
$$h_{ik}=\{\begin{array}{ll}1& \text{if vi∈Ak }\\ & \\ 0& \text{otherwise}\end{array}(i=1,\cdot \cdot \cdot ,n;k=1,\cdot \cdot \cdot ,c)\text{\hspace{1em}(指示向量)}$$

上式其实就是给每个类分配了一个对应的布尔指示向量

我们将这c个指示向量$h_k$组合成一个指示矩阵$H\in N^{n\ast c}$, 矩阵$H$当中的每一列指示向量正交于其他任何一列非线性相关向量。
$$Tr(H^{T}H)=n>0$$
我们可以得到以下目标函数：
$$\begin{gathered} \underset{A_1,A_2,\cdot \cdot \cdot ,A_c}{argminCut(A_1,\cdot \cdot \cdot ,A_c)}⟺argminTr(H^{T}LH) \\ s.tTr(H^{T}H)=n \end{gathered}$$
但是上面的目标函数存在一定局限，的哦到的分割结果往往倾向于将所连边最少且权值较低的孤立点分割出来。

因此我们对原来的目标函数进行改进。

目标函数改进（RatioCut 与 Ncut）

其实原目标函数出现容易出现孤立点的原因是没有对子图规模进行限制，因此我们可以想到在原目标函数的基础上除以一个子图的规模。
子图规模的度量一般有以下两种方式：

1. ✨从顶点数出发，用子图的顶点数代表子图的规模。
   $$|A|:=thenumberofverticesinA\text{\hspace{1em}(子图A的势)}$$
2. ✨从边的角度出发可以利用子图种边权重之和作为子图的规模。

$$vol(A):=\sum _{v_i\in A}{d}_i\text{\hspace{1em}(子图A的体积)}$$
对应的我们得到两种目标函数：

$$RatioCut(A_1,A_2,\cdot \cdot \cdot ,Ac):=\sum _{k=1}^c\frac{Cut(A_k,\overline{A_k})}{|A_k|}=\frac{1}{2}\sum _{k=1}^c\frac{W(A_k,\overline{A_k})}{|A_k|}\text{\hspace{1em}(比例切割)}$$

$$Ncut(A_1,A_2,\cdot \cdot \cdot ,Ac):=\sum _{k=1}^c\frac{Cut(A_k,\overline{A_k})}{vol(A_k)}=\frac{1}{2}\sum _{k=1}^c\frac{W(A_k,\overline{A_k})}{vol(A_k)}\text{\hspace{1em}(归一化切割)}$$

目标函数求解

下面，我们先从RatioCut 的目标函数开始

RatioCut 目标函数近似求解

从聚两类开始（c = 2)

在上面我门已经推导出了：
$$minCut(V)⟺min{f}^{T}Lf$$
使用拉格朗日乘子将约束问题转化为无约束问题

我们将目标函数转化为$Lagrangian$ 函数$L(f,\lambda )$：
$$\begin{gathered} L(f,\lambda ):=f^{T}Lf-\lambda (f^{T}1)-\lambda (f^{T}f-n) \\ =f^{T}Lf-\lambda (f^{T}f-n) \end{gathered}$$
等价将目标约束问题转换：
$$\begin{gathered} \underset{f\in R^n}{\min }L(f,\lambda )⟺\underset{f\in R^n}{\min }{f}^{T}Lf \\ s.tf\perp 1(thatis{f}^{T}1=0) \\ {f}^{T}f=n\text{(fTf>0)} \end{gathered}$$
现在我们得到无约束且连续的目标函数，由于二次函数是天然的凸函数亦可导，可以快速找到全局最优解，只需要令器导数为0 即可得到极值。

求微分

求得导数

由常量微分$dy=(\frac{dy}{dx}{)}^{T}dx$得：

$$\begin{gathered} \frac{dL(f,\lambda )}{df}=[2f^T{L}^T-2\lambda f^T{]}^T=2Lf-2\lambda f=0 \\ \Rightarrow Lf=\lambda f \end{gathered}$$

我们发现，当Lagrange乘子$\lambda$是拉普拉斯矩阵$L$的特征值（eigenvalue）且指示向量$f$是$L$的特征向量(eigenvalue)时，函数有极值。

求出极值

我们在等式两边，分别左乘$f^T$凑成目标函数形成后为：
$$f^{T}Lf=\lambda f^{T}f=\lambda n$$
由于$n=|V|$是常数，显然取决于目标函数值大小的元素时$\lambda$，即：
$$min\frac{f^{T}Lf}{f^{T}f}=min\lambda$$
也就是说$L$的特征值$\lambda$越大，目标函数值越大；特征值$\lambda$越小，目标函数值越小。

那么$\lambda$的最小值时${\lambda }_{min}=0$, 令拉普拉斯矩阵的第二小特征值作为目标函数最优值。

对于二聚类($c=2$)问题，，我们最简单粗暴的定性方法是判断$f_i$是否大于等于0，即

$$\{\begin{array}{ll}{v}_i\in A& \text{if fi≥0}\\ & \\ v_i\in \overline{A}& \text{if fi<0}\end{array}$$

多聚类(c 取任意值)

$$RatioCut(A_i,\cdot \cdot \cdot ,A_c)=\sum _{k=1}^c(H^{T}LH{)}_{ii}=Tr(H^{T}LH)$$

由于指示矩阵$H$中向量之间的正交性，我们可以得到如下完整的目标函数：

$$\begin{gathered} \underset{H\in R^{n\ast c}}{argmin}Tr(H^{T}LH) \\ s.t{H}^{T}H=I\text{[Tr(HTH)>0]} \end{gathered}$$

类似的，可以推出：

$$minTr(H^{T}LH)⟺min\lambda$$
由于特征值$\lambda$(拉格朗日乘子)在这里的物理意义代表切图的代价，因此我们可以将前$k$小特征值对应的$k$个特征向量拼接成一个矩阵$U\in R^{N\ast k}$作为 $H\in R^{n\ast c}$的近似解。然后使用传统的k-means方法，将连续实数值离散化，将矩阵$U$的n行向量聚为c行向量，最终得到的label即为最终聚类结果。(推导这么多，看不懂的话可以直接看结果)

RatioCut 谱聚类算法流程

Ncut 谱聚类算法流程

Python-RatioCut 代码实现

import math
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来正常显示负号

class MyPen(object):
  """
    ## 用作绘图的装饰器
  """
  def __init__(self, obj):
    obj.Pen = self
    self.obj = obj

  def __call__(self, *args, **kwds):
    return self.obj(*args, **kwds)

  @staticmethod
  def draw(pen,data,labels):
    """
      ## 绘制
    """
    nodeNum = len(data)
    scatterColors = ['black', 'blue', 'green', 'yellow', 'red', 'purple', 'orange']
    for i in range(len(set(labels))):
      color = scatterColors[i % len(scatterColors)]
      x1= []; y1=[]
      for j in range(nodeNum):
        if labels[j] == i:
          x1.append(data[j, 0])
          y1.append(data[j, 1])
      pen.scatter(x1, y1, c=color, marker='+')

  @staticmethod
  def viewResult(data,labels):
    """
      ## 绘制结果

      :pram ``data``: 数据集
      :pram ``clusterNum``: 聚类数量
      :pram ``labels``: 聚类结果
      :return None
    """
    MyPen.draw(plt,data,labels)
    plt.show()
  
  @staticmethod
  def viewCompare(data, trueLabels,predict_labels):
    """
    ## 绘制结果
    """,
    rows = 1
    cols = 2
    fig = plt.figure(figsize=(10,10))
    ax1 = fig.add_subplot(rows, cols, 1)
    ax2 = fig.add_subplot(rows, cols, 2)
    ax1.set_title('正确')
    ax2.set_title('预测')
    MyPen.draw(ax1, data, trueLabels)
    MyPen.draw(ax2, data, predict_labels)
    plt.show()


@MyPen
class MySpectral(object):
  """
    ## 谱聚类算法实现

  """
  def __init__(self,data, cluster_num=8, k=8):
    self.data = data
    self.cluster_num = cluster_num
    self.k = k

  def distance(self,x1,x2):
    """
      # 计算两个点之间的距离
      
    """
    return np.sqrt(np.sum((x1 - x2) ** 2))
  
  def get_dist_matrix(self):
    """
      ## 获取距离矩阵
    """
    nodesNum = len(self.data)
    distMatrix = np.array([
      [self.distance(self.data[i], self.data[j]) for j in range(nodesNum)] 
      for i in range(nodesNum)
    ])
    self.distMatrix = distMatrix
    return distMatrix
  
  def getW(self):
    """
      ## 获取邻接矩阵 W
    """
    distMatrix = self.get_dist_matrix()
    W = np.zeros_like(distMatrix)
    for idx, item in enumerate(distMatrix):
      idxArr = np.argsort(item)[1:self.k + 1]
      W[idx,idxArr] = 1
    self.W = (W + W.T) / 2
    return (W + W.T) / 2
  
  def getD(self):
    """
      ## 获得度矩阵
    """
    D = np.diag(np.sum(self.W,axis=1))
    self.D = D
    return D
  
  def getL(self):
    """
      ## 获得拉普拉斯矩阵
    """
    self.L = self.D - self.W
    return self.L

  def getEigen(self):
    """
    ## 获得特征值和特征向量

    """
    eigenValues, eigenVectors = np.linalg.eig(self.L)
    idx = np.argsort(eigenValues)[:self.cluster_num]
    return eigenVectors[:,idx]
  
  def fit(self):
    """
      ## 谱聚类算法入口
    """
    W = self.getW()
    D = self.getD()
    L = self.getL()
    eigenVectors = self.getEigen()
    clf = KMeans(n_clusters= self.cluster_num)
    s = clf.fit(eigenVectors)
    labels = s.labels_
    self.labels_ = labels
    return labels


def loadData():
  """
  # 导入数据集部分
  """
  ((nodesData, trueLabels)) = make_blobs()
  return nodesData, trueLabels

def test():
  data,trueLabels = loadData()
  model = MySpectral(data,8)
  model.fit()
  model.Pen.viewCompare(data, trueLabels ,model.labels_)

if __name__ == '__main__':
  test()

代码实现思路

主要分成两部分把一部分是可视化的写了一些装饰器的类，另一个就是谱聚类算法的类（重点）。

算法效果
比如我们设置聚8类

Python-Ncut 算法实现

根据推导结论Ncut 只要在RatioCut 算法的基础上修改一点就ok了，而且根据自己的实验发现Ncut 能更好得保证子图不会过小

import math
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来正常显示负号

class MyPen(object):
  """
    ## 用作绘图的装饰器
  """
  def __init__(self, obj):
    obj.Pen = self
    self.obj = obj

  def __call__(self, *args, **kwds):
    return self.obj(*args, **kwds)

  @staticmethod
  def draw(pen,data,labels):
    """
      ## 绘制
    """
    nodeNum = len(data)
    scatterColors = ['black', 'blue', 'green', 'yellow', 'red', 'purple', 'orange']
    for i in range(len(set(labels))):
      color = scatterColors[i % len(scatterColors)]
      x1= []; y1=[]
      for j in range(nodeNum):
        if labels[j] == i:
          x1.append(data[j, 0])
          y1.append(data[j, 1])
      pen.scatter(x1, y1, c=color, marker='+')

  @staticmethod
  def viewResult(data,labels):
    """
      ## 绘制结果

      :pram ``data``: 数据集
      :pram ``clusterNum``: 聚类数量
      :pram ``labels``: 聚类结果
      :return None
    """
    MyPen.draw(plt,data,labels)
    plt.show()
  
  @staticmethod
  def viewCompare(data, trueLabels,predict_labels):
    """
    ## 绘制结果
    """,
    rows = 1
    cols = 2
    fig = plt.figure(figsize=(10,10))
    ax1 = fig.add_subplot(rows, cols, 1)
    ax2 = fig.add_subplot(rows, cols, 2)
    ax1.set_title('正确')
    ax2.set_title('预测')
    MyPen.draw(ax1, data, trueLabels)
    MyPen.draw(ax2, data, predict_labels)
    plt.show()


@MyPen
class MySpectral(object):
  """
    ## 谱聚类算法实现

  """
  def __init__(self,data, cluster_num=8, k=8):
    self.data = data
    self.cluster_num = cluster_num
    self.k = k

  def distance(self,x1,x2):
    """
      # 计算两个点之间的距离
      
    """
    return np.sqrt(np.sum((x1 - x2) ** 2))
  
  def get_dist_matrix(self):
    """
      ## 获取距离矩阵
    """
    nodesNum = len(self.data)
    distMatrix = np.array([
      [self.distance(self.data[i], self.data[j]) for j in range(nodesNum)] 
      for i in range(nodesNum)
    ])
    self.distMatrix = distMatrix
    return distMatrix
  
  def getW(self):
    """
      ## 获取邻接矩阵 W
    """
    distMatrix = self.get_dist_matrix()
    W = np.zeros_like(distMatrix)
    for idx, item in enumerate(distMatrix):
      idxArr = np.argsort(item)[1:self.k + 1]
      W[idx,idxArr] = 1
    self.W = (W + W.T) / 2
    return (W + W.T) / 2
  
  def getD(self):
    """
      ## 获得度矩阵
    """
    D = np.diag(np.sum(self.W,axis=1))
    self.D = D
    self.halfD = np.diag(np.sum(self.W,axis=1) ** -0.5)
    return D
  
  def getL(self):
    """
      ## 获得拉普拉斯矩阵
    """
    self.L = self.D - self.W
    self.L = self.halfD @ self.L @ self.halfD     ##两种算法的区别之处
    return self.L

  def getEigen(self):
    """
    ## 获得特征值和特征向量

    """
    eigenValues, eigenVectors = np.linalg.eig(self.L)
    idx = np.argsort(eigenValues)[:self.cluster_num]
    return eigenVectors[:,idx]
  
  def fit(self):
    """
      ## 谱聚类算法入口
    """
    W = self.getW()
    D = self.getD()
    L = self.getL()
    eigenVectors = self.getEigen()
    clf = KMeans(n_clusters= self.cluster_num)
    s = clf.fit(eigenVectors)
    labels = s.labels_
    self.labels_ = labels
    return labels


def loadData():
  """
  # 导入数据集部分
  """
  ((nodesData, trueLabels)) = make_blobs()
  return nodesData, trueLabels

def test():
  data,trueLabels = loadData()
  model = MySpectral(data,8)
  model.fit()
  model.Pen.viewCompare(data, trueLabels ,model.labels_)

if __name__ == '__main__':
  test()