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

文章目录

谱聚类算法详解及代码实现
- 参考
- 关于谱聚类介绍
- - 谱聚类概述
  - 谱聚类前置知识
  - - 无向权重图
    - 邻接矩阵
    - 度矩阵
    - 拉普拉斯矩阵
    - 相似度矩阵
- 确定目标函数
- - 初始化目标函数（最小割MinCut方法）
  - - 我们先只考虑只有两个类的情况（c = 2)
    - 推广到多个类（c取任意值）
  - 目标函数改进（RatioCut 与 Ncut）
  - 目标函数求解
  - - RatioCut 目标函数近似求解
    - - 从聚两类开始（c = 2)
      - 多聚类(c 取任意值)
  - RatioCut 谱聚类算法流程
- Ncut 谱聚类算法流程
- Python-RatioCut 代码实现
- Python-Ncut 算法实现

参考

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=[w11w12w13w14w21w22w23w24w31w32w33w34w41w42w43w44](1)W = \left[ \begin{matrix} {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{matrix} \right] \tag{1} W=⎣⎢⎢⎡w11w21w31w41w12w22w32w42w13w23w33w43w14w24w34w44⎦⎥⎥⎤(1)

邻接矩阵

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

图一

W=(0001100111010001100011000)(2)W = \begin{pmatrix} {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{pmatrix} \tag{2} W=⎝⎜⎜⎜⎜⎛0001100111010001100011000⎠⎟⎟⎟⎟⎞(2)

度矩阵

度矩阵是对角矩阵，对角上元素为各个顶点的度。上述图一的度矩阵为式(3)所示。
D=(2000003000001000002000002)(3)D = \begin{pmatrix} {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{pmatrix} \tag{3} D=⎝⎜⎜⎜⎜⎛2000003000001000002000002⎠⎟⎟⎟⎟⎞(3)

拉普拉斯矩阵

拉普拉斯矩阵L= D - W，其中D为度矩阵， W为邻接矩阵。图一的拉普拉斯矩阵如式(4)所示。
L=(200−1−103−1−1−10−1100−1−1020−1−1002)(4)L = \begin{pmatrix} {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{pmatrix} \tag{4} L=⎝⎜⎜⎜⎜⎛200−1−103−1−1−10−1100−1−1020−1−1002⎠⎟⎟⎟⎟⎞(4)
用拉普拉斯矩阵来求解特征值，通过确定特征值来确定对应特征向量的个数，从而实现降维，然后再用K-means将降维的特征向量进行聚类。

相似度矩阵

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

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

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

方法一

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

wij={e∥xixj∥222σ2,if xi∈knn(xj)0,otherwise(方法1)wij = \begin{cases} e^\frac{\begin{Vmatrix}x_{i}&x_{j}\end{Vmatrix} _{2}^2}{2\sigma^2}, & \text{if $x_{i}$ $\in$ knn($x_{j}$)} \\[2ex] 0, & \text{otherwise} \end{cases} \tag{方法1} wij=⎩⎪⎨⎪⎧e2σ2∥xixj∥22,0,if xi ∈ knn(xj)otherwise(方法1)

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

方法二

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

wij=wji={e∥∥xixj∥222σ2,if xi∈knn(xj) or xj∈knn(xi) 0,otherwise(方法2)wij = wji = \begin{cases} e^\frac{\Vert\begin{Vmatrix}x_{i}&x_{j}\end{Vmatrix} _{2}^2}{2\sigma^2}, & \text{if $x_{i}$ $\in$ knn($x_{j}$) or $x_{j}$ $\in$ knn($x_{i}$) } \\[2ex] 0, & \text{otherwise} \end{cases} \tag{方法2} wij=wji=⎩⎪⎨⎪⎧e2σ2∥∥xixj∥22,0,if xi ∈ knn(xj) or xj ∈ knn(xi) otherwise(方法2)

方法三

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

wij=wji={e∥xixj∥222σ2,if xi∈knn(xj) and xj∈knn(xi) 0,otherwise(方法3)wij = wji = \begin{cases} e^\frac{\begin{Vmatrix}x_{i}&x_{j}\end{Vmatrix} _{2}^2}{2\sigma^2}, & \text{if $x_{i}$ $\in$ knn($x_{j}$) and $x_{j}$ $\in$ knn($x_{i}$) } \\[2ex] 0, & \text{otherwise} \end{cases} \tag{方法3} wij=wji=⎩⎪⎨⎪⎧e2σ2∥xixj∥22,0,if xi ∈ knn(xj) and xj ∈ knn(xi) otherwise(方法3)

最常用的方法一般是：

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

W=W+WT2W = \frac{W + W ^ T}{2} W=2W+WT

确定目标函数

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

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

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

对于无向图G(X,E) 进行切分的目标时将G划分为相互无连接的k个子图，每个子图包含点的集合A1A_{1}A1，A2A_{2}A2…AkA_{k}Ak,.(当然要保证每个点属于且只属于一个子图)。对于任意的两个子图点的集合A，B，可以定义A和B之间的切图权重为:
W(A,B)=∑i∈A,j∈Bwij(A-B 切图权重)W(A,B) = \sum_{i \in A, j \in B}{w_{ij}} \tag{A-B 切图权重} W(A,B)=i∈A,j∈B∑wij(A-B 切图权重)
那么我们可以把上述的切图权重公式扩展到k个子图的节点的集合种，即A1A_{1}A1,A2A_{2}A2…AkA_{k}Ak，定义切图cut为:
cut(A1,A2,A3...Ak)=12∑i∈A,j∈BkW(Ai,Ai‾)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}})} cut(A1,A2,A3...Ak)=21i∈A,j∈B∑kW(Ai,Ai)

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

其中Ai‾{\overline{A_{i}}}Ai 为Ai{A_{i}}Ai的补集（除了本身以外的所有子集的并集）。

因此我们可以进一步得到初步离散优化问题，最小切割目标函数:
minCut(V)=>min∑vi∈Ak,vj∈Ak‾wijmin Cut(V) => min \sum_{v_{i} \in A_{k}, v_{j} \in \overline{A_{k}} }{w_{ij}} minCut(V)=>minvi∈Ak,vj∈Ak∑wij

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

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

在这里我们定义只是向量f=(f1,f2,...fn)T∈Nn{f = (f_{1}, f{2},...f_{n})^T \in N^n}f=(f1,f2,...fn)T∈Nn
fi={1if vi∈A0if vi∈A‾f_{i} = \begin{cases} 1 & \text{if $v_{i}$ $\in$ A} \\\\ 0 & \text{if $v_{i}$ $\in$ $\overline{A}$} \end{cases} fi=⎩⎪⎨⎪⎧10if vi ∈ Aif vi ∈ A
因此目标函数转化为：
min∑vi∈Ak,vj∈Ak‾wij⟺min12∑i=1n∑j=1nwij(fi−fj)2min \sum_{v_{i} \in A_{k}, v_{j} \in \overline{A_{k}} }{w_{ij}} \Longleftrightarrow min \frac{1}{2} \sum_{i=1}^n\sum_{j=1}^n w_{ij}(f_{i} - f_{j})^2 minvi∈Ak,vj∈Ak∑wij⟺min21i=1∑nj=1∑nwij(fi−fj)2 {$}

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

我们将 12∑i=1sumj=1wij(fi−fj)2{\frac{1}{2} \sum_{i=1}sum_{j=1}w_{ij}(f_{i} - f_{j})^2}21∑i=1sumj=1wij(fi−fj)2 进行展开：

其中,D为度矩阵，W为邻接矩阵（相似矩阵），而L=D−W{L= D - W}L=D−W 为拉普拉斯矩阵，因此原目标函数等价于：
argminfi∈NnCut(V)=argminfi∈Nn12∑i=1n∑j=1nwij(fi−fj)2=argminfi∈NnfTLf\underset{f_{i} \in N^n}{arg \ min} Cut(V) = \underset{f_{i} \in N^n}{arg \ min}\frac{1}{2} \sum_{i = 1}^n\sum_{j=1}^n w_{ij}(f_{i} - f_{j})^2 = \underset{f_{i} \in N^n}{arg \ min} f ^ T L f fi∈Nnarg minCut(V)=fi∈Nnarg min21i=1∑nj=1∑nwij(fi−fj)2=fi∈Nnarg minfTLf
上面我们通过引入布尔指示向量将一个抽象的切割函数Cut(V){Cut(V)}Cut(V)转化为了就提可求解的fTLf{f^T L f}fTLf.但是这个是离散的优化函数，无法进行求导。

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

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

首先我们定义向量hk=(h(1,k),h2,k,...,hn,k)T∈Nnh_{k} = (h_{(1,k),h_{2,k}, ... ,h_{n,k} })^T \in N^{n}hk=(h(1,k),h2,k,...,hn,k)T∈Nn（其中：$k = 1,2, …,c $）并满足以下公式：
hik={1if vi∈Ak0otherwise(i=1,⋅⋅⋅,n;k=1,⋅⋅⋅,c)(指示向量)h_{ik} = \begin{cases} 1 & \text{if $v_{i} \in A_{k}$ } \\ \\ 0 & \text{otherwise} \end{cases} (i = 1, \cdot\cdot\cdot,n;k = 1,\cdot\cdot\cdot,c) \tag{指示向量} hik=⎩⎪⎨⎪⎧10if vi∈Ak otherwise(i=1,⋅⋅⋅,n;k=1,⋅⋅⋅,c)(指示向量)

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

我们将这c个指示向量hkh_{k}hk组合成一个指示矩阵H∈Nn∗cH \in N^{n * c}H∈Nn∗c, 矩阵HHH当中的每一列指示向量正交于其他任何一列非线性相关向量。
Tr(HTH)=n>0Tr(H^TH) = n > 0 Tr(HTH)=n>0
我们可以得到以下目标函数：
argminCut(A1,⋅⋅⋅,Ac)A1,A2,⋅⋅⋅,Ac⟺argminTr(HTLH)s.tTr(HTH)=n\underset{A_{1},A_{2},\cdot\cdot\cdot,A_{c}} {arg \ minCut(A_{1},\cdot\cdot\cdot,A_{c})} \Longleftrightarrow arg \ min \ Tr(H^TLH) \\ s.t \ \ Tr \ (H^TH) = n A1,A2,⋅⋅⋅,Acarg minCut(A1,⋅⋅⋅,Ac)⟺arg min Tr(HTLH)s.t Tr (HTH)=n
但是上面的目标函数存在一定局限，的哦到的分割结果往往倾向于将所连边最少且权值较低的孤立点分割出来。

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

目标函数改进（RatioCut 与 Ncut）

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

✨从顶点数出发，用子图的顶点数代表子图的规模。
∣A∣:=thenumberofverticesinA(子图A的势)|A| \ \ := \ the \ number \ of \ verticesin \ A \tag{子图A的势} ∣A∣ := the number of verticesin A(子图A的势)
✨从边的角度出发可以利用子图种边权重之和作为子图的规模。

vol(A):=∑vi∈Adi(子图A的体积)vol(A) \ \ := \ \sum_{v_{i} \in A}d_{i} \tag{子图A的体积} vol(A) := vi∈A∑di(子图A的体积)
对应的我们得到两种目标函数：

RatioCut(A1,A2,⋅⋅⋅,Ac):=∑k=1cCut(Ak,Ak‾)∣Ak∣=12∑k=1cW(Ak,Ak‾)∣Ak∣(比例切割)RatioCut(A_{1},A_{2},\cdot\cdot\cdot,A{c}) := \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}|} \tag{比例切割} RatioCut(A1,A2,⋅⋅⋅,Ac):=k=1∑c∣Ak∣Cut(Ak,Ak)=21k=1∑c∣Ak∣W(Ak,Ak)(比例切割)

Ncut(A1,A2,⋅⋅⋅,Ac):=∑k=1cCut(Ak,Ak‾)vol(Ak)=12∑k=1cW(Ak,Ak‾)vol(Ak)(归一化切割)Ncut(A_{1},A_{2},\cdot\cdot\cdot,A{c}) := \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})} \tag{归一化切割} Ncut(A1,A2,⋅⋅⋅,Ac):=k=1∑cvol(Ak)Cut(Ak,Ak)=21k=1∑cvol(Ak)W(Ak,Ak)(归一化切割)

目标函数求解

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

RatioCut 目标函数近似求解

从聚两类开始（c = 2)

在上面我门已经推导出了：
minCut(V)⟺minfTLfmin \ Cut(V) \Longleftrightarrow min \ f^TLf min Cut(V)⟺min fTLf
使用拉格朗日乘子将约束问题转化为无约束问题

我们将目标函数转化为LagrangianLagrangianLagrangian 函数L(f,λ)L(f,\lambda)L(f,λ)：
L(f,λ):=fTLf−λ(fT1)−λ(fTf−n)=fTLf−λ(fTf−n)L(f,\lambda) \ := \ f^TLf - \lambda(f^T1) - \lambda(f^Tf - n) \\ = f^TLf - \lambda(f^Tf - n) L(f,λ) := fTLf−λ(fT1)−λ(fTf−n)=fTLf−λ(fTf−n)
等价将目标约束问题转换：
minf∈RnL(f,λ)⟺minf∈RnfTLfs.tf⊥1(thatisfT1=0)fTf=n(fTf>0)\underset{f \in R^n}{min}L(f, \lambda) \Longleftrightarrow \underset{f \in R^n}{min}f^TLf \\ s.t f \bot 1 (that \ is \ f^T 1 = 0 ) \\ f^Tf \ = \ n \text{($f^Tf > 0$)} f∈RnminL(f,λ)⟺f∈RnminfTLfs.tf⊥1(that is fT1=0)fTf = n(fTf>0)
现在我们得到无约束且连续的目标函数，由于二次函数是天然的凸函数亦可导，可以快速找到全局最优解，只需要令器导数为0 即可得到极值。

求微分

求得导数

由常量微分dy=(dydx)Tdxdy = (\frac{dy}{dx})^Tdxdy=(dxdy)Tdx得：

dL(f,λ)df=[2fTLT−2λfT]T=2Lf−2λf=0⇒Lf=λf\frac{dL(f,\lambda)}{df} \ = \ [ 2f^TL^T - 2 \lambda f^T ]^T = 2Lf - 2 \lambda f = 0 \\ \Rightarrow Lf = \lambda f dfdL(f,λ) = [2fTLT−2λfT]T=2Lf−2λf=0⇒Lf=λf

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

求出极值

我们在等式两边，分别左乘fTf^TfT凑成目标函数形成后为：
fTLf=λfTf=λnf^TLf = \lambda f^Tf = \lambda n fTLf=λfTf=λn
由于n=∣V∣n = |V|n=∣V∣是常数，显然取决于目标函数值大小的元素时λ\lambdaλ，即：
minfTLffTf=minλmin \ \frac{f^TLf}{f^Tf} = min \ \lambda min fTffTLf=min λ
也就是说LLL的特征值λ\lambdaλ越大，目标函数值越大；特征值λ\lambdaλ越小，目标函数值越小。

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

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

{vi∈Aif fi≥0vi∈A‾if fi<0\begin{cases} v_{i} \in \ A & \text{if \ $f_{i} \geq 0 $} \\ \\ v_{i} \in \ \overline{A} & \text{if $f_{i} < 0 $} \end{cases} ⎩⎪⎨⎪⎧vi∈ Avi∈ Aif fi≥0if fi<0

多聚类(c 取任意值)

RatioCut(Ai,⋅⋅⋅,Ac)=∑k=1c(HTLH)ii=Tr(HTLH)RatioCut \ (A_{i},\cdot\cdot\cdot,A_{c}) \ = \ \sum_{k=1}^{c}(H^TLH)_{ii} \ = \ Tr(H^TLH) RatioCut (Ai,⋅⋅⋅,Ac) = k=1∑c(HTLH)ii = Tr(HTLH)

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

argminH∈Rn∗cTr(HTLH)s.tHTH=I[Tr(HTH)>0]\underset{H \in R^{n * c}}{arg \ min} \ Tr(H^TLH) \\ s.t \ H^T H = I \ \ \ \ \ \ \ \ \text{$[Tr(H^TH) > 0]$} H∈Rn∗carg min Tr(HTLH)s.t HTH=I [Tr(HTH)>0]

类似的，可以推出：

minTr(HTLH)⟺minλmin \ Tr(H^TLH) \Longleftrightarrow min \ \lambda min Tr(HTLH)⟺min λ
由于特征值λ\lambdaλ(拉格朗日乘子)在这里的物理意义代表切图的代价，因此我们可以将前kkk小特征值对应的kkk个特征向量拼接成一个矩阵U∈RN∗kU \ \in R^{N * k}U ∈RN∗k作为 H∈Rn∗cH \in R^{n*c}H∈Rn∗c的近似解。然后使用传统的k-means方法，将连续实数值离散化，将矩阵UUU的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_blobsplt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来正常显示负号class MyPen(object):"""## 用作绘图的装饰器"""def __init__(self, obj):obj.Pen = selfself.obj = objdef __call__(self, *args, **kwds):return self.obj(*args, **kwds)@staticmethoddef 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='+')@staticmethoddef viewResult(data,labels):"""## 绘制结果:pram ``data``: 数据集:pram ``clusterNum``: 聚类数量:pram ``labels``: 聚类结果:return None"""MyPen.draw(plt,data,labels)plt.show()@staticmethoddef viewCompare(data, trueLabels,predict_labels):"""## 绘制结果""",rows = 1cols = 2fig = 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 = dataself.cluster_num = cluster_numself.k = kdef 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 = distMatrixreturn distMatrixdef 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] = 1self.W = (W + W.T) / 2return (W + W.T) / 2def getD(self):"""## 获得度矩阵"""D = np.diag(np.sum(self.W,axis=1))self.D = Dreturn Ddef getL(self):"""## 获得拉普拉斯矩阵"""self.L = self.D - self.Wreturn self.Ldef 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_ = labelsreturn labelsdef loadData():"""# 导入数据集部分"""((nodesData, trueLabels)) = make_blobs()return nodesData, trueLabelsdef 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_blobsplt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来正常显示负号class MyPen(object):"""## 用作绘图的装饰器"""def __init__(self, obj):obj.Pen = selfself.obj = objdef __call__(self, *args, **kwds):return self.obj(*args, **kwds)@staticmethoddef 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='+')@staticmethoddef viewResult(data,labels):"""## 绘制结果:pram ``data``: 数据集:pram ``clusterNum``: 聚类数量:pram ``labels``: 聚类结果:return None"""MyPen.draw(plt,data,labels)plt.show()@staticmethoddef viewCompare(data, trueLabels,predict_labels):"""## 绘制结果""",rows = 1cols = 2fig = 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 = dataself.cluster_num = cluster_numself.k = kdef 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 = distMatrixreturn distMatrixdef 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] = 1self.W = (W + W.T) / 2return (W + W.T) / 2def getD(self):"""## 获得度矩阵"""D = np.diag(np.sum(self.W,axis=1))self.D = Dself.halfD = np.diag(np.sum(self.W,axis=1) ** -0.5)return Ddef getL(self):"""## 获得拉普拉斯矩阵"""self.L = self.D - self.Wself.L = self.halfD @ self.L @ self.halfD     ##两种算法的区别之处return self.Ldef 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_ = labelsreturn labelsdef loadData():"""# 导入数据集部分"""((nodesData, trueLabels)) = make_blobs()return nodesData, trueLabelsdef test():data,trueLabels = loadData()model = MySpectral(data,8)model.fit()model.Pen.viewCompare(data, trueLabels ,model.labels_)if __name__ == '__main__':test()