-
2021-09-11 10:51:47
学习心得
(1)laplacian matrix就是无向图中定义 L = D − A L=D-A L=D−A,其中D为邻接矩阵,A为度矩阵(是一个对角矩阵)。本文用的python计算拉普拉斯矩阵及其特征值、特征向量。
(2)numpy中文文档:https://www.numpy.org.cn/user/,可以配spyder看对应变量的矩阵元素。
(3)求图拉普拉斯矩阵的特征值:时间复杂度是 O ( n 3 ) O(n^3) O(n3)(n为节点数),所以当图很大时用该方法效率很低。文章目录
一、背景介绍
1.1 图论基础
定义一(图的邻接矩阵):
-
给定一个图 G = { V , E } \mathcal{G}=\{\mathcal{V}, \mathcal{E}\} G={V,E},其对应的邻接矩阵被记为 A ∈ { 0 , 1 } N × N \mathbf{A} \in\{0,1\}^{N \times N} A∈{0,1}N×N。 A i , j = 1 \mathbf{A}_{i, j}=1 Ai,j=1表示存在从结点 v i v_i vi到 v j v_j vj的边,反之表示不存在从结点 v i v_i vi到 v j v_j vj的边。
-
在无向图中,从结点 v i v_i vi到 v j v_j vj的边存在,意味着从结点 v j v_j vj到 v i v_i vi的边也存在。因而无向图的邻接矩阵是对称的。
-
在无权图中,各条边的权重被认为是等价的,即认为各条边的权重为 1 1 1。
-
对于有权图,其对应的邻接矩阵通常被记为 W ∈ { 0 , 1 } N × N \mathbf{W} \in\{0,1\}^{N \times N} W∈{0,1}N×N,其中 W i , j = w i j \mathbf{W}_{i, j}=w_{ij} Wi,j=wij表示从结点 v i v_i vi到 v j v_j vj的边的权重。若边不存在时,边的权重为 0 0 0。
一个无向无权图的例子:
其邻接矩阵为:
A = ( 0 1 0 1 1 1 0 1 0 0 0 1 0 0 1 1 0 0 0 1 1 0 1 1 0 ) \mathbf{A}=\left(\begin{array}{lllll} 0 & 1 & 0 & 1 & 1 \\ 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & 1 \\ 1 & 0 & 1 & 1 & 0 \end{array}\right) A=⎝⎜⎜⎜⎜⎛0101110100010011000110110⎠⎟⎟⎟⎟⎞
定义二(拉普拉斯矩阵,Laplacian Matrix):
- 给定一个图 G = { V , E } \mathcal{G}=\{\mathcal{V}, \mathcal{E}\} G={V,E},其邻接矩阵为 A A A,其拉普拉斯矩阵定义为 L = D − A \mathbf{L=D-A} L=D−A,其中度矩阵 D = d i a g ( d ( v 1 ) , ⋯ , d ( v N ) ) \mathbf{D=diag(d(v_1), \cdots, d(v_N))} D=diag(d(v1),⋯,d(vN))。
定义三(对称归一化的拉普拉斯矩阵,Symmetric normalized Laplacian):
- 给定一个图 G = { V , E } \mathcal{G}=\{\mathcal{V}, \mathcal{E}\} G={V,E},其邻接矩阵为 A A A,其规范化的拉普拉斯矩阵定义为
L = D − 1 2 ( D − A ) D − 1 2 = I − D − 1 2 A D − 1 2 \mathbf{L=D^{-\frac{1}{2}}(D-A)D^{-\frac{1}{2}}=I-D^{-\frac{1}{2}}AD^{-\frac{1}{2}}} L=D−21(D−A)D−21=I−D−21AD−21
1.2 拉普拉斯矩阵的变体
节点数为 n n n的简单图 G G G, A A A是 G G G的邻接矩阵,则如上面介绍的那样, G G G的拉普拉斯矩阵即 L = D − A L=D-A L=D−A。
(1)对称归一化后的拉普拉斯矩阵: L s y s = D − 1 / 2 L D − 1 / 2 = I − D − 1 / 2 A D − 1 / 2 L^{s y s}=D^{-1 / 2} L D^{-1 / 2}=I-D^{-1 / 2} A D^{-1 / 2} Lsys=D−1/2LD−1/2=I−D−1/2AD−1/2
(2)随机游走归一化的拉普拉斯矩阵: L r w = D − 1 L = I − D − 1 A L i , j r w : = { 1 i = j − 1 diag ( v i ) if i ≠ j and v i adjacent to v j 0 othewise \begin{aligned} &L^{r w}=D^{-1} L=I-D^{-1} A \\ &L_{i, j}^{r w}:= \begin{cases}1 & \mathrm{i}=\mathrm{j} \\ \frac{-1}{\sqrt{\operatorname{diag}\left(v_{i}\right)}} & \text { if } i \neq j \text { and } v_{i} \text { adjacent to } v_{j} \\ 0 & \text { othewise }\end{cases} \end{aligned} Lrw=D−1L=I−D−1ALi,jrw:=⎩⎪⎪⎨⎪⎪⎧1diag(vi)−10i=j if i=j and vi adjacent to vj othewise 1.3 拉普拉斯矩阵的优良性质:
- 拉普拉斯矩阵是半正定对称矩阵
- 对称矩阵有n个线性无关的特征向量,n是Graph中节点的个数 ⇒ \Rightarrow ⇒ 拉普拉斯矩阵可以特征分解
- 半正定矩阵的特征值非负
- 对称矩阵的特征向量构成的矩阵为正交阵 ⇒ U T U = E \Rightarrow U^{T} U=E ⇒UTU=E
1.4 GCN为啥要用拉普拉斯矩阵
- 拉普拉斯矩阵可以谱分解(特征分解)GCN是从谱域的角度提取拓扑图的空间特征的。
- 拉普拉斯矩阵只在中心元素和一阶相邻元素处有非零元素,其他位置皆为0.
- 传统傅里叶变换公式中的基函数是拉普拉斯算子,借助拉普拉斯矩阵,通过类比可以推导出Graph上的傅里叶变换公式。
二、Python代码实现
这是networkX库对稀疏矩阵A的处理方式,有少量改进(将所有内容保持稀疏):
n,m = A.shape diags = A.sum(axis=1) D = sps.spdiags(diags.flatten(), [0], m, n, format='csr') D - A
numpy
和scipy
两个库中模块中都提供了线性代数的库linalg
,但scipy
的linalg
的更全面一些。# laplacian矩阵 import numpy as np def unnormalized_laplacian(adj_matrix): # 先求度矩阵 R = np.sum(adj_matrix, axis=1) degreeMatrix = np.diag(R) return degreeMatrix - adj_matrix # 对称归一化的laplacian矩阵 def normalized_laplacian(adj_matrix): R = np.sum(adj_matrix, axis=1) R_sqrt = 1/np.sqrt(R) D_sqrt = np.diag(R_sqrt) I = np.eye(adj_matrix.shape[0]) return I - np.matmul(np.matmul(D_sqrt, adj_matrix), D_sqrt)
matlab的代码实现为:
L = diag(sum(A,2)) - A % or L=diag(sum(A))-A because A is symmetric
那么如何求矩阵的特征值呢——利用numpy的
linalg.eig
:# !/usr/bin/python ## -*-coding:utf-8 -*- import numpy as np A = np.array([[3,-1],[-1,3]]) print('打印A:\n{}'.format(A)) a, b = np.linalg.eig(A) print('打印特征值a:\n{}'.format(a)) print('打印特征向量b:\n{}'.format(b))
得到特征值和特征向量:
打印A: [[ 3 -1] [-1 3]] 打印特征值a: [4. 2.] 打印特征向量b: [[ 0.70710678 0.70710678] [-0.70710678 0.70710678]]
三阶矩阵试试,回归考研线代的题目:
import numpy as np A = np.array([[2,-3,1],[1,-2,1],[1,-3,2]]) print('打印A:\n{}'.format(A)) a, b = np.linalg.eig(A) print('打印特征值a:\n{}'.format(a)) print('打印特征向量b:\n{}'.format(b))
结果如下,看结果的特征向量矩阵,每一列代表一个一个特征向量,都是
打印A: [[ 2 -3 1] [ 1 -2 1] [ 1 -3 2]] 打印特征值a: [2.09037533e-15+0.00000000e+00j 1.00000000e+00+5.87474805e-16j 1.00000000e+00-5.87474805e-16j] 打印特征向量b: [[0.57735027+0.j 0.84946664+0.j 0.84946664-0.j ] [0.57735027+0.j 0.34188085-0.11423045j 0.34188085+0.11423045j] [0.57735027+0.j 0.17617591-0.34269135j 0.17617591+0.34269135j]]
三、关于图Fourier变换
根据卷积原理,卷积公式可以写成
f ∗ g = F − 1 { F { f } ⋅ F { g } } f * g=\mathcal{F}^{-1}\{\mathcal{F}\{f\} \cdot \mathcal{F}\{g\}\} f∗g=F−1{F{f}⋅F{g}}
正、逆Fourier变换
F ( v ) = ∫ R f ( x ) e − 2 π i x ⋅ v d x \mathcal{F}(v)=\int_{\mathrm{R}} f(x) e^{-2 \pi i x \cdot v} d x F(v)=∫Rf(x)e−2πix⋅vdxf ( x ) = ∫ R F ( v ) e 2 π i x ⋅ v d v f(x)=\int_{\mathbb{R}} \mathcal{F}(v) e^{2 \pi i x \cdot v} d v f(x)=∫RF(v)e2πix⋅vdv
一阶导数定义
f ′ ( x ) = lim h → 0 f ( x + h ) − f ( x ) h f^{\prime}(x)=\lim _{h \rightarrow 0} \frac{f(x+h)-f(x)}{h} f′(x)=h→0limhf(x+h)−f(x)
拉普拉斯相当于二阶导数
Δ f ( x ) = lim h → 0 f ( x + h ) − 2 f ( x ) + f ( x − h ) h 2 \Delta f(x)=\lim _{h \rightarrow 0} \frac{f(x+h)-2 f(x)+f(x-h)}{h^{2}} Δf(x)=h→0limh2f(x+h)−2f(x)+f(x−h)
在graph上,定义一阶导数为
f ∗ g ′ ( x ) = f ( x ) − f ( y ) f_{* g}^{\prime}(x)=f(x)-f(y) f∗g′(x)=f(x)−f(y)
对应的拉普拉斯算子定义为
Δ ∗ g f ′ ( x ) = Σ y ∼ x ( f ( x ) − f ( y ) ) \Delta_{*g} f^{\prime}(x)=\Sigma_{y \sim x} (f(x)-f(y)) Δ∗gf′(x)=Σy∼x(f(x)−f(y))
假设 D D D为 N × N N\times{N} N×N的度矩阵(degree matrix)
D ( i , j ) = { d i if i = j 0 otherwise D(i, j)=\left\{\begin{array}{ll}{d_{i}} & {\text { if } i=j} \\ {0} & {\text { otherwise }}\end{array}\right. D(i,j)={di0 if i=j otherwise
A A A为 N × N N\times{N} N×N的邻接矩阵(adjacency matrix)
A ( i , j ) = { 1 if x i ∼ x j 0 otherwise A(i, j)=\left\{\begin{array}{ll}{1} & {\text { if } x_{i} \sim x_{j}} \\ {0} & {\text { otherwise }}\end{array}\right. A(i,j)={10 if xi∼xj otherwise
那么图上的Laplacian算子可以写成
L = D − A L=D-A L=D−A
标准化后得到
L = I N − D − 1 2 A D − 1 2 L=I_{N}-D^{-\frac{1}{2}} A D^{-\frac{1}{2}} L=IN−D−21AD−21
定义Laplacian算子的目的是为了找到Fourier变换的基。传统Fourier变换的基就是Laplacian算子的一组特征向量
Δ e 2 π i x ⋅ v = λ e 2 π i x ⋅ v \Delta e^{2 \pi i x \cdot v}=\lambda e^{2 \pi i x \cdot v} Δe2πix⋅v=λe2πix⋅v
类似的,在graph上,有
Δ f = ( D − A ) f = L f \Delta{f}=(D-A)f=Lf Δf=(D−A)f=Lf
图拉普拉斯算子作用在由图节点信息构成的向量 f f f上得到的结果等于图拉普拉斯矩阵和向量 f f f的点积。那么graph上的Fourier基就是 L L L矩阵的 n n n个特征向量 U = [ u 1 … u n ] U=\left[u_{1} \dots u_{n}\right] U=[u1…un], L L L可以分解成 L = U Λ U T L=U \Lambda U^{T} L=UΛUT,其中, Λ \Lambda Λ是特征值组成的对角矩阵。
传统的Fourier变换与graph的Fourier变换区别
将 f ( i ) f(i) f(i)看成是第 i i i个点上的signal,用向量 x = ( f ( 1 ) … f ( n ) ) ∈ R n x=(f(1) \ldots f(n)) \in \mathbb{R}^{n} x=(f(1)…f(n))∈Rn来表示。矩阵形式的graph的Fourier变换为
G F { x } = U T x \mathcal{G F}\{x\}=U^{T} x GF{x}=UTx
类似的graph上的Fourier逆变换为
I G F { x } = U x \mathcal{I} \mathcal{G} \mathcal{F}\{x\}=U x IGF{x}=UxReference
(1)Python计算特征值与特征向量案例
(2)numpy的linalg
官方文档:https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.linalg.html
(3)用python求解特征向量和拉普拉斯矩阵
(4)图卷积网络(GCN)原理解析更多相关内容 -
-
拉普拉斯矩阵
2019-10-10 16:29:29拉普拉斯矩阵(Laplacian matrix) 也叫做导纳矩阵、基尔霍夫矩阵或离散拉普拉斯算子,主要应用在图论中,作为一个图的矩阵表示。 1.定义 给定一个有n个顶点的图G,它的拉普拉斯矩阵定义为: L=D-A 其中D为图的...拉普拉斯矩阵(Laplacian matrix) 也叫做导纳矩阵、基尔霍夫矩阵或离散拉普拉斯算子,主要应用在图论中,作为一个图的矩阵表示。
1.定义
给定一个有n个顶点的图G,它的拉普拉斯矩阵
定义为:
L=D-A
其中D为图的度矩阵,A为图的邻接矩阵。度矩阵在有向图中,只需要考虑出度或者入度中的一个。经过计算可以得
(1)若i =j,则
,
为顶点
的度;
(2)若i≠ j,但顶点
和顶点
相邻,则
;
(3)其它情况
。
也可以将这三种值通过除以
进行标准化。
2.举例
例图:不带权的无向图
度矩阵: (体现点与点之间关系)
邻接矩阵:(体现点与边之间关系)
拉普拉斯矩阵:(体现点与点 和 点与边之间关系)
3.性质
(1)拉普拉斯矩阵是半正定矩阵;
(2)特征值中0出现的次数就是图连通区域的个数;
(3)最小特征值是0,因为拉普拉斯矩阵每一行的和均为0;
(4)最小非零特征值是图的代数连通度。
-
解读拉普拉斯矩阵
2020-07-05 15:53:48拉普拉斯矩阵(Laplacian matrix) 也叫做导纳矩阵,这次笔记主要是记录下GCN学习时的注意点,在图论(Graph theory)中,对于图 G=(V,E): Laplacian 矩阵的定义为解读Laplace矩阵
1. 什么是Laplace矩阵?
拉普拉斯矩阵(Laplacian matrix) 也叫做导纳矩阵,这次笔记主要是记录下GCN学习时的注意点,在图论(Graph theory)中,对于图 G=(V,E):
Laplacian 矩阵的定义为 L = D - A (其中 L 是Laplacian 矩阵, D=diag(d)是对角矩阵,d=rowSum(A),对角线上元素依次为各个顶点的度, A 则是图的邻接矩阵)
若只考虑无向图,那么L就是对称矩阵。
对于无向图的Laplace矩阵,它有哪些性质?
- 半正定矩阵(特征值非负,且是对称矩阵);
- 对称矩阵(一定有n个线性无关的特征向量);
- 对称矩阵的不同特征值对应的特征向量相互正交,这些正交的特征向量构成的矩阵为正交矩阵;
- 由于是半正定矩阵,所以是对称阵,那么能特征值分解(EVD)。
由于 U 是正交矩阵,即UUT=I,所以特征值分解又可以写成:
2. 常见的Laplace矩阵
2.1 一般形式
2.2 对称归一化形式
2.3 随机游走归一化形式
Reference
-
带你玩转谱聚类及拉普拉斯矩阵
2020-02-21 10:29:06引言:在多变量统计和数据聚类中,谱聚类(Spectral Clustering)技术利用数据的相似矩阵的谱(特征值)进行降维。它将数据看成空间中的点,点对之间有边相连,距离越远的点对其边权值越小,距离越近的点对其边权值越大...引言:在多变量统计和数据聚类中,谱聚类(Spectral Clustering)技术利用数据的相似矩阵的谱(特征值)进行降维。它将数据看成空间中的点,点对之间有边相连,距离越远的点对其边权值越小,距离越近的点对其边权值越大。它将聚类问题转化为切图问题,使得切图后的总代价最小。即子图内点对之间边权值较大,子图间边权值较小。得到切图后子图的个数即为聚类的个数。最后本文力求用推理而非演绎的方式以加深大家对谱聚类的理解。
申明:本文严禁转载!
带着问题去阅读:
- 为什么叫谱聚类?谱的含义是什么?(见引言,矩阵的谱表示其特征值的集合)
- 谱聚类如何将聚类问题转化为图切割问题,常用的图构造方法是什么?
- 为什么要引入Laplacian矩阵?它有哪些性质?
- 为什么要取Laplacian矩阵第 2 \textbf{2} 2小或前 k k k小个特征值对应的特征向量?其特征值 λ \lambda λ有什么意义?它和瑞利商 (Rayleigh Quotient) \text{(Rayleigh Quotient)} (Rayleigh Quotient)有什么关系?前 k k k小中的 k k k一定和聚类个数 c c c相等吗?
- 为什么Laplacian矩阵0特征值的数量就是连通子图的个数?
- 谱聚类如何将一个难求解的离散优化问题松弛为易求解的连续优化问题?
- 为什么最后还要对特征向量进行k-means聚类?(需要将连续值离散化)
- 原目标函数聚类效果已经不错,为什么还要对拉普拉斯矩阵 L L L进行归一化,有几种归一化方法?分别又是什么?( L r w L_{rw} Lrw和 L s y m L_{sym} Lsym)
- 对于三种图拉普拉斯矩阵 L , L r w L,L_{rw} L,Lrw和 L s y m L_{sym} Lsym,它们的适用场景是什么?
- 我们发现在吴恩达( Ng \text{Ng} Ng)等人提出的算法中,在构造了 L s y m L_{sym} Lsym归一化矩阵后,又对前 k k k个特征向量构成的矩阵 U ∈ R n × k U\in R^{n\times k} U∈Rn×k按行再次进行了归一化,这样做是否多此一举? Ng \text{Ng} Ng等人的目的何在?
- 通过对谱聚类的学习,有哪些我们可以借鉴的思想可以用于我们的科研工作当中?
文章目录
- 一、问题抽象
- 二、目标函数的转化
- 三、目标函数的求解
- 四、谱聚类算法具体实现
- 五、谱聚类中值得借鉴的思想
- 参考文献
一、问题抽象
假设有 n n n个实数样本数据如下,每个样本有 d d d维,目标是要聚 c c c个类,并且数据分布并非云团,如何做?
X = { x 1 d , x 2 d , . . . , x n d } T , X ∈ R n × d X=\{x_1^d,x_2^d,...,x_n^d\}^T,X\in R^{n\times d} X={x1d,x2d,...,xnd}T,X∈Rn×d
将数据点抽象为无向带权图(图构造)
在图论中我们常用邻接矩阵W表示图(无向图是为了保证邻接矩阵是对称矩阵),因此我们只需按某个准则来计算数据点对之间的距离即可获得数据点的邻接矩阵。因此引出以下常用 3 3 3种准则来构造图。(对于3个数据点的邻接矩阵如下)
W = [ w 11 w 12 w 13 w 21 w 22 w 23 w 31 w 32 w 33 ] \begin{gathered} W=\begin{bmatrix} w_{11} & w_{12} & w_{13}\\ w_{21} & w_{22} & w_{23}\\ w_{31} & w_{32} & w_{33}\\ \end{bmatrix} \end{gathered} W=⎣⎡w11w21w31w12w22w32w13w23w33⎦⎤1.1 全连接法
所有点对之间的边权值均为正数(>0),采用核函数作为准则来计算点对之间的距离(也称为相似度或亲和度),进而直接得到实对称相似度矩阵 W = W T W=W^T W=WT,常用高斯核函数。令 s i j s_{ij} sij表示两点之间的边权重(也称两点相似度),则高斯核计算公式如下
w i j = s i j = s ( x i , x j ) = e − ∣ ∣ x i − x j ∣ ∣ 2 2 2 σ 2 (1.1) w_{ij} = s_{ij}= s(x_i,x_j)=e^{-\frac{||x_i-x_j||_2^2}{2\sigma^2}} \tag{1.1} wij=sij=s(xi,xj)=e−2σ2∣∣xi−xj∣∣22(1.1)
其中参数 σ \sigma σ控制着点对 x i x_{i} xi与 x j x_{j} xj之间的度量, σ \sigma σ越大,两点度量后的差距越小,聚类精度也越低; σ \sigma σ越小,两点度量后的差距越大。1.2 k k k-近邻法
计算点对之间的欧氏距离,通过给定参数 k k k,选取距离当前点最近的 k k k个点为邻居(常用高斯核计算距离),令其余点到该点距离为0。
w i j = { e − ∣ ∣ x i − x j ∣ ∣ 2 2 2 σ 2 i f x i ∈ k n n ( x j ) 0 o t h e r w i s e ( 可 能 存 在 w i j ≠ w j i ) (1.2) w_{ij} = \begin{cases} e^{-\frac{||x_i-x_j||_2^2}{2\sigma^2}}&if\; x_i\in knn(x_j)\\ 0& \ otherwise \end{cases} \quad (可能存在w_{ij} \not= w_{ji}) \tag{1.2} wij=⎩⎨⎧e−2σ2∣∣xi−xj∣∣220ifxi∈knn(xj) otherwise(可能存在wij=wji)(1.2)然而,这样做有一个显著的问题,数据从无向图变成了有向图。即你是我的邻居,但我不一定是你的邻居(即 w i j ≠ w j i w_{ij}\not= w_{ji} wij=wji)。为了保证相似度矩阵的对称性,论文[1]给出了两种解决方法将有向图变成无向图。
-
方法一:若两点间有两条有向边,则忽视一条,仅保留一条。具体做法如下:
w i j = w j i = { e − ∣ ∣ x i − x j ∣ ∣ 2 2 2 σ 2 i f x i ∈ k n n ( x j ) o r x j ∈ k n n ( x i ) 0 o t h e r w i s e (1.3) w_{ij} =w_{ji}= \begin{cases} e^{-\frac{||x_i-x_j||_2^2}{2\sigma^2}}&if\; x_i\in knn(x_j)\; or\; x_j\in knn(x_i)\\ 0& otherwise \end{cases} \tag{1.3} wij=wji=⎩⎨⎧e−2σ2∣∣xi−xj∣∣220ifxi∈knn(xj)orxj∈knn(xi)otherwise(1.3) -
方法二:当且仅当你是我的邻居,并且我也是你的邻居时,我俩之间才有边权值。具体做法如下
w i j = w j i = { e − ∣ ∣ x i − x j ∣ ∣ 2 2 2 σ 2 i f x i ∈ k n n ( x j ) a n d x j ∈ k n n ( x i ) 0 o t h e r w i s e (1.4) w_{ij} =w_{ji}= \begin{cases} e^{-\frac{||x_i-x_j||_2^2}{2\sigma^2}}&if\; x_i\in knn(x_j)\;and\; x_j\in knn(x_i)\\ 0& otherwise \end{cases} \tag{1.4} wij=wji=⎩⎨⎧e−2σ2∣∣xi−xj∣∣220ifxi∈knn(xj)andxj∈knn(xi)otherwise(1.4) -
方法三:采用平均化思想,先求每个点的 k k k近邻,得到不对称的有向图。然后采用如下做法变成对称的无向图,这也是实战常用的方法之一:
W = W + W T 2 (1.5) W = \frac{W + W^T}{2} \tag{1.5} W=2W+WT(1.5)
1.3 ϵ \epsilon ϵ-近邻法
给定一个阈值 ϵ \epsilon ϵ,如果两点距离大于这个阈值则认为这两点无边相连,否则这两点间的边权值为 ϵ \epsilon ϵ。具体做法如下:
w i j = w j i = { ϵ i f s ( x i , x j ) < ϵ 0 o t h e r w i s e (1.6) w_{ij} =w_{ji}= \begin{cases} \epsilon&if\; s(x_i,x_j)<\epsilon\\ 0& \ otherwise \end{cases} \tag{1.6} wij=wji={ϵ0ifs(xi,xj)<ϵ otherwise(1.6)
该方法得到的矩阵可视为无向无权图,因此精度不如其他方法。二、目标函数的转化
经过问题抽象后,我们将原始数据按k近邻法构造出如下拓扑图。现在我们可以将聚 c c c个类的问题转化为将无向图切割为 c c c个子图的问题。
由上图可知,切分的方法有很多,我们该如何切出 c c c个子图从而达到我们想要聚类的效果呢?很容易想到一个准则,即不同切图方法需要付出不同的代价,当我们将距离较近的看似为同一个类的两个点切分为不同类时,需要付出较大代价;当我们将距离较远的两个点切分到不同子图时,需要付出较小代价。因此我们可以定义一个代价函数来作为我们初步的目标函数。
2.1 初始目标函数的形成(最小割 Min Cut \text{Min Cut} Min Cut方法)
我们令无向图 G = < V , E > G=<V,E> G=<V,E>,其中 V ∈ R n × d V\in R^{n\times d} V∈Rn×d是 n n n个 d d d维顶点的集合, v i v_i vi表示第 i i i个顶点。 E E E是边的集合, e i j e_{ij} eij表示顶点 v i v_i vi到点 v j v_j vj的边。假设我们要将所有顶点 V V V切分成 c c c个子图,则令 A i A_i Ai表示第 i i i个子图所包含顶点的集合(称为图 G G G的一个分割), A ‾ i \overline A_i Ai表示第 i i i个子图的补图(其中 i = 1 , 2 , . . . , c i=1,2,...,c i=1,2,...,c,所以 V = A 1 ∪ A 2 ∪ . . . ∪ A c V=A_1\cup A_2 \cup...\cup A_c V=A1∪A2∪...∪Ac)。现在我们令 C u t ( V ) Cut(V) Cut(V)来表示切图的代价,则
C u t ( V ) = C u t ( A 1 , A 2 , . . . , A c ) (2.1) Cut(V)=Cut(A_1,A_2,...,A_c) \tag{2.1} Cut(V)=Cut(A1,A2,...,Ac)(2.1)
这是一个典型的 o n e vs a l l one\; \text{vs} \;all onevsall问题,由于我们不知道 A i A_i Ai中都有哪些顶点,大的分类问题我们不会做,我们可以先找最简单的 c a s e case case(即二分类)来做,然后再将多个小问题的解合并为一个大问题的完整解,显然该问题是可分且子问题的解也可合并,所以利用分治思想我们可将一个多分类问题转化为多个二分类问题。方便起见,我们用 W ( A i , A ‾ i ) W(A_i,\overline A_i) W(Ai,Ai)来表示子图 A i A_i Ai和其补图 A ‾ i \overline A_i Ai被切割后的代价,而二分类中两个子图的分割代价可用连接两个子图的所有边的权值之和来表示(前面提到边权值越大两点越近,反之越远)。因此进一步得到如下式子(左乘 1 / 2 1/2 1/2是因为对称性 w i j = w j i w_{ij}=w_{ji} wij=wji,相当于算了两遍):
C u t ( V ) = C u t ( A 1 , A 2 , . . . , A c ) = ∑ k = 1 c C u t ( A k , A ‾ k ) = 1 2 ∑ k = 1 c W ( A k , A ‾ k ) = 1 2 ∑ v i ∈ A k , v j ∈ A ‾ k , e i j ∈ E w i j ( i , j = 1 , 2 , . . . , n ) (2.2) Cut(V)=Cut(A_1,A_2,...,A_c)=\sum_{k=1}^cCut(A_k,\overline A_k)=\frac{1}{2} \sum_{k=1}^c W(A_k,\overline A_k)=\frac{1}{2} \sum_{v_i\in A_k,v_j\in \overline A_k,e_{ij}\in E}w_{ij} \quad (i,j=1,2,...,n) \tag{2.2} Cut(V)=Cut(A1,A2,...,Ac)=k=1∑cCut(Ak,Ak)=21k=1∑cW(Ak,Ak)=21vi∈Ak,vj∈Ak,eij∈E∑wij(i,j=1,2,...,n)(2.2)
因此,我们可以进一步得到初步离散优化问题,即最小割目标函数:
min C u t ( V ) ⇒ min ∑ v i ∈ A k , v j ∈ A ‾ k , e i j ∈ E w i j (2.3) \min \;\; Cut(V) \Rightarrow \min \;\;\sum_{v_i\in A_k,v_j\in \overline A_k,e_{ij}\in E}w_{ij} \tag{2.3} minCut(V)⇒minvi∈Ak,vj∈Ak,eij∈E∑wij(2.3)2.2 细化初始目标函数之引入指示向量( indicator vector \text{indicator vector} indicator vector)
得到最小割目标函数后,我们发现其约束 v i ∈ A k , v j ∈ A ‾ k , e i j ∈ E v_i\in A_k,v_j\in \overline A_k,e_{ij}\in E vi∈Ak,vj∈Ak,eij∈E过于笼统模糊,难以求解。因此,我们需要定性地引入指示向量来细化目标函数。在此我们先避开聚多类问题,从最简单的二聚类问题开始讨论,再逐步归纳到多聚类问题。
2.2.1 先讨论二聚类问题( c = 2 c=2 c=2)
对于聚二类问题,我们不再需要 A i A_i Ai,仅用 A A A足以表示 G G G的一个分割,另一个分割用 A ‾ \overline A A表示。因此原目标函数为
min C u t ( V ) ⇒ min C u t ( A , A ‾ ) ⇒ min ∑ v i ∈ A , v j ∈ A ‾ , e i j ∈ E w i j (2.4) \min \;\; Cut(V) \Rightarrow \min \;\; Cut(A,\overline A) \Rightarrow \min \;\;\sum_{v_i\in A,v_j\in \overline A,e_{ij}\in E}w_{ij} \tag{2.4} minCut(V)⇒minCut(A,A)⇒minvi∈A,vj∈A,eij∈E∑wij(2.4)① 目标函数的展开
对目标函数观察后我们发现,第三个约束 e i j ∈ E e_{ij}\in E eij∈E很好处理,如在构图方法中选择 k k k近邻法,使得点对之间有边相连时 w i j > 0 w_{ij}>0 wij>0,无边相连时 w i j = 0 w_{ij}=0 wij=0,因此选择 k k k近邻法可直接去除该约束。
对于剩下的两个约束 v i ∈ A , v j ∈ A ‾ v_i\in A,v_j\in \overline A vi∈A,vj∈A,我们定义指示向量 f = ( f 1 , f 2 , . . . , f n ) T ∈ N n f=(f_1,f_2,...,f_n)^T \in N^n f=(f1,f2,...,fn)T∈Nn(注:此处的 f i f_i fi是一个布尔值而不是向量)并满足如下条件:
f i = { 1 i f v i ∈ A 0 i f v i ∈ A ‾ (2.5) f_i = \begin{cases} 1 & if \; v_i \in A \\ 0 & if \; v_i \in \overline A \tag{2.5} \end{cases} fi={10ifvi∈Aifvi∈A(2.5)
这样的指示向量在论文[1]中被称为理想型指示向量(ideal indicator vector)。假设每个样本点以概率为 1 1 1属于某个类,不存在某个样本属于多个类的情况。我们的目的是,当满足约束 v i ∈ A , v j ∈ A ‾ v_i\in A,v_j\in \overline A vi∈A,vj∈A时 w i j w_{ij} wij有值,否则不对其进行求和。因此可以借助指示向量将目标函数细化如下:
min C u t ( A , A ‾ ) ⇔ min 1 2 ∑ v i ∈ A , v j ∈ A ‾ , e i j ∈ E w i j ⇔ min 1 2 ∑ i = 1 n ∑ j = 1 n w i j ( f i − f j ) 2 \min \;\; Cut(A,\overline A) \Leftrightarrow \min \;\; \frac{1}{2} \sum_{v_i\in A,v_j\in \overline A,e_{ij}\in E}w_{ij} \Leftrightarrow \min \;\; \frac{1}{2} \sum_{i=1}^n \sum_{j=1}^n w_{ij}(f_i-f_j)^2 minCut(A,A)⇔min21vi∈A,vj∈A,eij∈E∑wij⇔min21i=1∑nj=1∑nwij(fi−fj)2
现在先不着急求解目标函数,我们展开 w i j ( f i − f j ) 2 w_{ij}(f_i-f_j)^2 wij(fi−fj)2二次项如下:
1 2 ∑ i = 1 n ∑ j = 1 n w i j ( f i − f j ) 2 = 1 2 ∑ i = 1 n ∑ j = 1 n ( w i j f i 2 − 2 f i f j w i j + w i j f j 2 ) (2.6) \frac{1}{2} \sum_{i=1}^n \sum_{j=1}^n w_{ij}(f_i-f_j)^2 = \frac{1}{2} \sum_{i=1}^n \sum_{j=1}^n (w_{ij}f_i^2 -2f_if_jw_{ij}+w_{ij}f_j^2) \tag{2.6} 21i=1∑nj=1∑nwij(fi−fj)2=21i=1∑nj=1∑n(wijfi2−2fifjwij+wijfj2)(2.6)② 非标准化 Laplacian \text{Laplacian} Laplacian矩阵的引入
已知在图论中,一个无向图 G G G中第 i i i个顶点的(degree)度 d i d_i di表示如下:
d i = ∑ j = 1 n w i j (2.7) d_i = \sum_{j=1}^n w_{ij} \tag{2.7} di=j=1∑nwij(2.7)
因此上式可以继续展开如下,看到加权求和我们可以快速想到两向量的点积:
1 2 ∑ i = 1 n ∑ j = 1 n w i j ( f i − f j ) 2 = 1 2 ∑ i = 1 n ∑ j = 1 n ( w i j f i 2 − 2 f i f j w i j + w i j f j 2 ) = 1 2 ( ∑ i = 1 n d i f i 2 − ∑ i , j = 1 n 2 f i f j w i j + ∑ j = 1 n d j f j 2 ) = 1 2 ( 2 ∑ i = 1 n d i f i 2 − 2 ∑ i , j = 1 n f i f j w i j ) = ∑ i = 1 n d i f i 2 − ∑ i , j = 1 n f i f j w i j = ∑ i = 1 n f i d i f i − ∑ i , j = 1 n f i w i j f j = f T D f − f T W f = f T L f (2.8) \begin{aligned} \frac{1}{2} \sum_{i=1}^n \sum_{j=1}^n w_{ij}(f_i-f_j)^2 & =\frac{1}{2} \sum_{i=1}^n \sum_{j=1}^n (w_{ij}f_i^2 -2f_if_jw_{ij}+w_{ij}f_j^2) \\ & = \frac{1}{2} (\sum_{i=1}^n d_{i}f_i^2 - \sum_{i,j=1}^n 2f_if_jw_{ij}+\sum_{j=1}^n d_{j}f_j^2) \\ & =\frac{1}{2} (2\sum_{i=1}^n d_{i}f_i^2- 2\sum_{i,j=1}^n f_if_jw_{ij}) \\ & =\sum_{i=1}^n d_{i}f_i^2- \sum_{i,j=1}^n f_if_jw_{ij} = \sum_{i=1}^n f_id_{i}f_i- \sum_{i,j=1}^n f_iw_{ij}f_j \\ & =f^TDf-f^TWf \\ & =f^TLf \tag{2.8} \end{aligned} 21i=1∑nj=1∑nwij(fi−fj)2=21i=1∑nj=1∑n(wijfi2−2fifjwij+wijfj2)=21(i=1∑ndifi2−i,j=1∑n2fifjwij+j=1∑ndjfj2)=21(2i=1∑ndifi2−2i,j=1∑nfifjwij)=i=1∑ndifi2−i,j=1∑nfifjwij=i=1∑nfidifi−i,j=1∑nfiwijfj=fTDf−fTWf=fTLf(2.8)
其中, D D D为度矩阵, W W W为邻接矩阵(相似矩阵),而 L = D − W L=D-W L=D−W为Laplacian矩阵。没想到从一个二次项 w i j ( f i − f j ) 2 w_{ij}(f_i-f_j)^2 wij(fi−fj)2目标函数经过展开可以得到拉普拉斯矩阵。因此,原目标函数再次等价于:
arg min A ⊂ V C u t ( V ) ⇔ arg min f i ∈ N n 1 2 ∑ i = 1 n ∑ j = 1 n w i j ( f i − f j ) 2 ⇔ arg min f i ∈ N n f T L f s . t f i as defined in Eq. (2.5) \begin{aligned} \argmin_{A \subset V} \;\; Cut(V) \Leftrightarrow & \argmin_{f_i\in N^n} \;\; \frac{1}{2} \sum_{i=1}^n \sum_{j=1}^n w_{ij}(f_i-f_j)^2 \Leftrightarrow & \argmin_{f_i\in N^n} \;\; f^TLf \\ & s.t \;\;f_i \; \text{as defined in Eq. (2.5)} \\ \end{aligned} A⊂VargminCut(V)⇔fi∈Nnargmin21i=1∑nj=1∑nwij(fi−fj)2⇔s.tfias defined in Eq. (2.5)fi∈NnargminfTLf
现在,我们通过引入布尔指示向量将一个抽象的切割函数 C u t ( V ) Cut(V) Cut(V)转化为具体可求解的 f T L f f^TLf fTLf。但是我们通过约束 ( 2.5 ) (2.5) (2.5)发现这是一个离散优化函数,无法对其求导。不过,最小割问题已经有确定解法,只是这不是我们在此讨论的重点,最小割目标函数在谱聚类问题中的定义本身存在偏离,其分割结果往往更倾向于将所连边最少且边权值较低的孤立点分割出来。因此需要继续改进,具体如何改进我们将在2.3章节详细讨论。我们在此之所以对最小割目标函数的转换讨论这么详细是因为这个思想在后面会用到。
2.2.2 再回到多聚类问题( c c c取任意值)
讨论完二聚类的简单问题,我们再来归纳到多聚类问题就容易多了。对于多聚类问题,我们要聚 c c c个类( k = 1 , 2 , . . . , c k=1,2,...,c k=1,2,...,c),因此我们仅仅需要将原来一个指示向量扩展为 c c c个指示向量组合成的指示矩阵 H H H即可,其余步骤相似。
指示向量扩展到指示矩阵
对于多聚类问题,我们依然定义指示向量,只需要最后把 c c c个指示向量组合成一个指示矩阵即可。
我们定义指示向量 h k = ( h ( 1 , k ) , h ( 2 , k ) , . . . , h ( n , k ) ) T ∈ N n × 1 h_{k}=(h_{(1,k)},h_{(2,k)},...,h_{(n,k)})^T \in N^{n\times 1} hk=(h(1,k),h(2,k),...,h(n,k))T∈Nn×1(其中: k = 1 , 2 , . . . , c k=1,2,...,c k=1,2,...,c)并满足如下条件:
h i k = { 1 i f v i ∈ A k 0 o t h e r w i s e ( i = 1 , . . . , n ; k = 1 , . . . , c ) (2.9) h_{ik} = \begin{cases} 1 & if \; v_i \in A_k \\ 0 & otherwise \tag{2.9} \end{cases} \quad (i = 1,...,n;k=1,...,c) hik={10ifvi∈Akotherwise(i=1,...,n;k=1,...,c)(2.9)
现在,我们将这 c c c个指示向量 h k h_k hk组合成一个指示矩阵 H ∈ N n × c H\in N^{n\times c} H∈Nn×c,由于我们假设每个样本以概率为1属于某个类,因此矩阵 H H H当中的每一列指示向量正交(orthonormal)于其他任何一列向量。(注意:向量之间的正交不代表 H H H是正交阵)
【例如】我们要聚 c = 3 c=3 c=3个类,有 5 5 5个样本,聚类得到的情况是前两个样本属于一个类,接着两个样本属于一个类,最后一个样本独成一类。则
H T H = [ 1 1 0 0 0 0 0 1 1 0 0 0 0 0 1 ] × [ 1 0 0 1 0 0 0 1 0 0 1 0 0 0 1 ] = [ 2 0 0 0 2 0 0 0 1 ] \begin{gathered} H^T H= \begin{bmatrix} 1 & 1 & 0 & 0 & 0\\ 0 & 0 & 1& 1& 0\\ 0 & 0 & 0& 0& 1\\ \end{bmatrix} \times \begin{bmatrix} 1 & 0 & 0\\ 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ \end{bmatrix} = \begin{bmatrix} 2 & 0 & 0\\ 0 & 2 & 0\\ 0 & 0 & 1\\ \end{bmatrix} \end{gathered} HTH=⎣⎡100100010010001⎦⎤×⎣⎢⎢⎢⎢⎡110000011000001⎦⎥⎥⎥⎥⎤=⎣⎡200020001⎦⎤
用数学语言可表达为: Tr ( H T H ) = n > 0 \text{Tr}(H^TH)=n>0 Tr(HTH)=n>0。因此我们可以得到如下目标函数
arg min A 1 , . . . , A c C u t ( A 1 , . . . , A c ) ⇔ arg min H ∈ R n × c Tr ( H T L H ) s . t Tr ( H T H ) = n H as defined in Eq. ( 2.9 ) \begin{aligned} \argmin_{A_1,...,A_c} \; Cut(A_1,...,A_c) \Leftrightarrow \argmin_{H\in R^{n\times c}}& \;\; \text{Tr}(H^TLH) \\ s.t & \;\; \text{Tr}(H^TH) =n\\ & \;\; H \text{ as defined in Eq.} (2.9)\\ \end{aligned} A1,...,AcargminCut(A1,...,Ac)⇔H∈Rn×cargmins.tTr(HTLH)Tr(HTH)=nH as defined in Eq.(2.9)最终,我们按照二聚类的思路,结合扩展的指示矩阵可以得到相似的最小割目标函数。但是由于最小割目标函数的局限性。下面,我们将对最小割目标函数进行改进,并重新定义指示向量的值。
2.3 目标函数的改进( RatioCut \text{RatioCut} RatioCut与 Ncut \text{Ncut} Ncut)
由最小割目标函数可知,最小割得到的分割结果往往更倾向于将所连边最少且边权值较低的孤立点分割出来。如下图所示,我们更希望被分割的子图中顶点的数量相互均衡。因此我们需要加入其它限制以改进目标函数分割效果。
为了达到均衡的效果,我们很容易想到在原目标函数(切成每个子图的代价)的基础上除以切图后每个子图的规模。
因此,现在的问题就是如何度量切图后子图的规模? ①从顶点数的角度思考可以想到用子图的顶点数作为该子图的规模。②从边的角度思考可以想到用子图的边权和作为子图的规模。两种度量方式用数学语言有如下定义:
子 图 A 的 势 ∣ A ∣ : = the number of vertices in A 子 图 A 的 体 积 v o l ( A ) : = ∑ v i ∈ A d i (2.10) \begin{aligned} &子图A的势\;|A|:= \text{the number of vertices in } A \\ &子图A的体积\;vol(A):= \sum_{v_i\in A} d_i \\ \end{aligned} \tag{2.10} 子图A的势∣A∣:=the number of vertices in A子图A的体积vol(A):=vi∈A∑di(2.10)
可以看出第二种子图A的体积更能贴切描述子图的规模,因为它相比子图A的势含有边权值的信息,能更清楚的区别有相同顶点数的但边权和不同的规模。当原二聚类目标函数为
C u t ( V ) : = C u t ( A 1 , A 2 , . . . , A c ) = ∑ k = 1 c C u t ( A k , A ‾ k ) = 1 2 ∑ k = 1 c W ( A k , A ‾ k ) Cut(V):=Cut(A_1,A_2,...,A_c)=\sum_{k=1}^cCut(A_k,\overline A_k) =\frac{1}{2} \sum_{k=1}^c W(A_k,\overline A_k) Cut(V):=Cut(A1,A2,...,Ac)=k=1∑cCut(Ak,Ak)=21k=1∑cW(Ak,Ak)
根据前面规模的定义式 ( 2.10 ) (2.10) (2.10),我们可以得到如下两个改进后的目标函数 R a t i o C u t RatioCut RatioCut(比例切割)与 N c u t Ncut Ncut(归一化切割):
R a t i o C u t ( A 1 , A 2 , . . . , A c ) : = ∑ k = 1 c C u t ( A k , A ‾ k ) ∣ A k ∣ = 1 2 ∑ k = 1 c W ( A k , A ‾ k ) ∣ A k ∣ N c u t ( A 1 , A 2 , . . . , A c ) : = ∑ k = 1 c C u t ( A k , A ‾ k ) v o l ( A k ) = 1 2 ∑ k = 1 c W ( A k , A ‾ k ) v o l ( A k ) (2.11) \begin{aligned} &RatioCut(A_1,A_2,...,A_c) := \sum_{k=1}^c \frac{Cut(A_k,\overline A_k)}{|A_k|} = \frac12 \sum_{k=1}^c \frac{W(A_k,\overline A_k)}{|A_k|} \tag{2.11} \\ &Ncut(A_1,A_2,...,A_c) := \sum_{k=1}^c \frac{Cut(A_k,\overline A_k)}{vol(A_k)} = \frac12 \sum_{k=1}^c \frac{W(A_k,\overline A_k)}{vol(A_k)} \end{aligned} RatioCut(A1,A2,...,Ac):=k=1∑c∣Ak∣Cut(Ak,Ak)=21k=1∑c∣Ak∣W(Ak,Ak)Ncut(A1,A2,...,Ac):=k=1∑cvol(Ak)Cut(Ak,Ak)=21k=1∑cvol(Ak)W(Ak,Ak)(2.11)
得到改进后的目标函数后,我们将尝试对其进行求解。三、目标函数的求解
下面,我们先从非标准化 R a t i o C u t RatioCut RatioCut目标函数讲起。
3.1 RatioCut \text{RatioCut} RatioCut目标函数的近似解
3.1.1 先从二聚类说起( c = 2 c=2 c=2)
我们要求解的二聚类目标函数为
min A ⊂ V R a t i o C u t ( A , A ‾ ) \min_{A\subset V}\; RatioCut(A,\overline A) A⊂VminRatioCut(A,A)
现在,我们重新定义一下在 R a t i o C u t RatioCut RatioCut目标函数中的指示向量为 f = ( f 1 , f 2 , . . . , f n ) T ∈ R n f=(f_1,f_2,...,f_n)^T \in R^n f=(f1,f2,...,fn)T∈Rnf i = { ∣ A ‾ ∣ ∣ A ∣ i f v i ∈ A − ∣ A ∣ ∣ A ‾ ∣ i f v i ∈ A ‾ (3.1) f_i = \begin{cases} \sqrt{\frac{|\overline A|}{|A|}} & if \; v_i \in A \\ -\sqrt{\frac{|A|}{|\overline A|}} & if \; v_i \in \overline A \end{cases} \tag{3.1} fi=⎩⎨⎧∣A∣∣A∣−∣A∣∣A∣ifvi∈Aifvi∈A(3.1)
在之前,我们推导出
min C u t ( V ) ⇔ min f T L f \min \;Cut(V) \Leftrightarrow\min \; f^TLf minCut(V)⇔minfTLf
现在我们来反向推导一下看看能得出什么
f T L f = 1 2 ∑ i = 1 n ∑ j = 1 n w i j ( f i − f j ) 2 = 1 2 ∑ v i ∈ A , v j ∈ A ‾ w i j ( ∣ A ‾ ∣ ∣ A ∣ + ∣ A ∣ ∣ A ‾ ∣ ) 2 + 1 2 ∑ v i ∈ A ‾ , v j ∈ A w i j ( − ∣ A ∣ ∣ A ‾ ∣ − ∣ A ‾ ∣ ∣ A ∣ ) 2 = c u t ( A , A ‾ ) ( ∣ A ‾ ∣ ∣ A ∣ + ∣ A ∣ ∣ A ‾ ∣ + 2 ) = c u t ( A , A ‾ ) ( ∣ A ‾ ∣ + ∣ A ∣ ∣ A ∣ + ∣ A ∣ + ∣ A ‾ ∣ ∣ A ‾ ∣ ) = ∣ V ∣ ⋅ c u t ( A , A ‾ ) ( 1 ∣ A ∣ + 1 ∣ A ‾ ∣ ) = ∣ V ∣ ⋅ R a t i o C u t ( A , A ‾ ) (3.2) \begin{aligned} f^TLf & =\frac{1}{2} \sum_{i=1}^n \sum_{j=1}^n w_{ij}(f_i-f_j)^2 \\ & = \frac{1}{2} \sum_{v_i\in A,v_j\in \overline A} w_{ij}(\sqrt{\frac{|\overline A|}{|A|}}+\sqrt{\frac{|A|}{|\overline A|}})^2 + \frac{1}{2} \sum_{v_i\in \overline A,v_j\in A} w_{ij}(-\sqrt{\frac{|A|}{|\overline A|}}-\sqrt{\frac{|\overline A|}{|A|}})^2 \\ & =cut(A,\overline A)(\frac{|\overline A|}{|A|}+\frac{|A|}{|\overline A|}+2)\\ & =cut(A,\overline A)(\frac{|\overline A|+|A|}{|A|}+\frac{|A|+|\overline A|}{|\overline A|})\\ &=|V|\cdot cut(A,\overline A)(\frac{1}{|A|}+\frac{1}{|\overline A|})\\ &=|V|\cdot RatioCut(A,\overline A) \end{aligned} \tag{3.2} fTLf=21i=1∑nj=1∑nwij(fi−fj)2=21vi∈A,vj∈A∑wij(∣A∣∣A∣+∣A∣∣A∣)2+21vi∈A,vj∈A∑wij(−∣A∣∣A∣−∣A∣∣A∣)2=cut(A,A)(∣A∣∣A∣+∣A∣∣A∣+2)=cut(A,A)(∣A∣∣A∣+∣A∣+∣A∣∣A∣+∣A∣)=∣V∣⋅cut(A,A)(∣A∣1+∣A∣1)=∣V∣⋅RatioCut(A,A)(3.2)我们得到了 R a t i o C u t RatioCut RatioCut目标函数!!同时对于指示向量的正交性,我们有 f ⊥ 1 f \bot 1 f⊥1,证明如下
∑ i = 1 n f i ⋅ 1 = ∑ v i ∈ A ∣ A ‾ ∣ ∣ A ∣ − ∑ v i ∈ A ‾ ∣ A ∣ ∣ A ‾ ∣ = ∣ A ∣ ∣ A ‾ ∣ ∣ A ∣ − ∣ A ‾ ∣ ∣ A ∣ ∣ A ‾ ∣ = 0 \sum_{i=1}^n f_i\cdot 1 = \sum_{v_i\in A} \sqrt{\frac{|\overline A|}{|A|}} - \sum_{v_i\in \overline A} \sqrt{\frac{|A|}{|\overline A|}} = |A| \sqrt{\frac{|\overline A|}{|A|}} - |\overline A|\sqrt{\frac{|A|}{|\overline A|}} = 0 i=1∑nfi⋅1=vi∈A∑∣A∣∣A∣−vi∈A∑∣A∣∣A∣=∣A∣∣A∣∣A∣−∣A∣∣A∣∣A∣=0
同时我们还发现 f T f f^Tf fTf内积可得到一个常数。
f T f = ∑ i = 1 n f i ⋅ f i = ∑ v i ∈ A ∣ A ‾ ∣ ∣ A ∣ + ∑ v i ∈ A ‾ ∣ A ∣ ∣ A ‾ ∣ = ∣ A ∣ ∣ A ‾ ∣ ∣ A ∣ + ∣ A ‾ ∣ ∣ A ∣ ∣ A ‾ ∣ = ∣ V ∣ = n (3.3) f^Tf=\sum_{i=1}^n f_i\cdot f_i = \sum_{v_i\in A} \frac{|\overline A|}{|A|} + \sum_{v_i\in \overline A} \frac{|A|}{|\overline A|} = |A| \frac{|\overline A|}{|A|} + |\overline A|\frac{|A|}{|\overline A|} = |V|= n \tag{3.3} fTf=i=1∑nfi⋅fi=vi∈A∑∣A∣∣A∣+vi∈A∑∣A∣∣A∣=∣A∣∣A∣∣A∣+∣A∣∣A∣∣A∣=∣V∣=n(3.3)这是非常关键的条件!我们用拉格朗日乘子转化为无约束问题后可以得到一个有趣的结论。现在我们可以得到等价的完整目标函数
min A ⊂ V R a t i o C u t ( A , A ‾ ) ⇔ min f ∈ R n f T L f s . t f ⊥ 1 f i as defined in Eq.(3.1) f T f = n ( f T f > 0 ) \begin{aligned} \min_{A \subset V} \;RatioCut(A,\overline A) \Leftrightarrow & \min_{f \in R^n} \; f^TLf \\ & s.t \;\; f \bot 1 \\ & \quad \;\; f_i \text{ as defined in Eq.(3.1)} \\ & \quad \;\; f^Tf =n \quad(f^Tf>0) \\ \end{aligned} A⊂VminRatioCut(A,A)⇔f∈RnminfTLfs.tf⊥1fi as defined in Eq.(3.1)fTf=n(fTf>0)
在《A tutorial on spectral clustering》[1]这篇论文中(如下图),作者将第三个约束条件写为 ∣ ∣ f ∣ ∣ = n ||f||=\sqrt n ∣∣f∣∣=n,即 f T f f^Tf fTf得到的是一个常数。也是一个等式约束,即 f T f − n = 0 f^Tf-n=0 fTf−n=0。
由于约束 f i as defined in Eq.(3.1) f_i \text{ as defined in Eq.(3.1)} fi as defined in Eq.(3.1)的限制,即解向量 f i f_i fi的取值只有两个值,导致该问题是一个NP难的离散优化问题,目标函数不可导。如何解这个离散问题?
离散问题我们很难做或者不会做,但连续问题我们会做,所以我们想可不可以把离散问题转化成连续问题?在这里,我们可以放松一下约束条件,去掉 f i as defined in Eq.(3.1) f_i \text{ as defined in Eq.(3.1)} fi as defined in Eq.(3.1)这个令人头疼的限制,使得解向量 f i f_i fi可取实数域 R R R中的任意值。这样就将离散问题松弛为连续问题如下
min f ∈ R n f T L f s . t f ⊥ 1 ( that is f T 1 = 0 ) f T f = n ( f T f > 0 ) (3.4) \begin{aligned} & \min_{f \in R^n} \; f^TLf \\ & s.t \;\; f \bot 1 \;(\text{that is }f^T1 =0) \\ & \quad \;\; f^Tf =n \quad(f^Tf>0)\\ \end{aligned} \tag{3.4} f∈RnminfTLfs.tf⊥1(that is fT1=0)fTf=n(fTf>0)(3.4)使用拉格朗日乘子将约束问题转化为无约束问题
有约束的目标函数我们不会做,我们可以转化为无约束问题,因此目标函数 ( 3.4 ) (3.4) (3.4)可以转化为如下 L a g r a n g e Lagrange Lagrange函数 L ( f , λ ) L(f,\lambda) L(f,λ)
L ( f , λ ) : = f T L f − λ ( f T 1 ) − λ ( f T f − n ) = f T L f − λ ( f T f − n ) (3.5) \begin{aligned} L(f,\lambda):&= f^TLf -\lambda(f^T1)- \lambda(f^Tf-n) \\ & = f^TLf - \lambda(f^Tf-n)\\ \end{aligned}\tag{3.5} L(f,λ):=fTLf−λ(fT1)−λ(fTf−n)=fTLf−λ(fTf−n)(3.5)
因此
min f ∈ R n L ( f , λ ) ⇔ min f ∈ R n f T L f s . t f ⊥ 1 ( that is f T 1 = 0 ) f T f = n ( f T f > 0 ) \begin{aligned} \min_{f \in R^n} L(f,\lambda) \Leftrightarrow \;& \min_{f \in R^n} \; f^TLf \\ & s.t \;\; f \bot 1 \;(\text{that is }f^T1 =0) \\ & \quad \;\; f^Tf =n \quad(f^Tf>0)\\ \end{aligned} f∈RnminL(f,λ)⇔f∈RnminfTLfs.tf⊥1(that is fT1=0)fTf=n(fTf>0)
现在我们得到无约束且连续的目标函数,由于二次型函数是天然的凸函数亦可导,可以快速找到全局最优解,只需要令其导数为 0 0 0即可得到极值,对向量和矩阵求导详请见《机器学习常用矩阵求导方法》。①求微分
d L ( f , λ ) = d [ f T L f − λ ( f T f − n ) ] = d ( f T L f − λ f T f ) = d ( f T ) L f + f T L ( d f ) − λ d ( f T ) f − λ f T ( d f ) = ( d f ) T L f + f T L ( d f ) − λ ( d f ) T f − λ f T ( d f ) = t r [ ( d f ) T L f + f T L ( d f ) − λ ( d f ) T f − λ f T ( d f ) ] = t r [ ( d f ) T L f ] + t r [ f T L ( d f ) ] − λ ⋅ t r [ ( d f ) T f ] − λ ⋅ t r [ f T ( d f ) ] = t r [ f T L T ( d f ) ] + t r [ f T L ( d f ) ] − λ ⋅ t r [ f T ( d f ) ] − λ ⋅ t r [ f T ( d f ) ] = t r [ f T ( L T + L ) ( d f ) − 2 λ f T ( d f ) ] = [ f T ( L T + L ) − 2 λ f T ] d f = [ 2 f T L T − 2 λ f T ] d f (3.6) \begin{aligned} dL(f,\lambda) & = d [f^TLf -\lambda (f^Tf -n)] =d (f^TLf -\lambda f^Tf ) \\ & = d(f^T)Lf + f^TL(df) -\lambda d(f^T)f - \lambda f^T(df) \\ & = (df) ^TLf + f^TL(df) -\lambda (df) ^Tf - \lambda f^T(df) \\ & = tr[(df) ^TLf + f^TL(df) -\lambda (df) ^Tf - \lambda f^T(df)] \\ & = tr[(df) ^TLf] + tr[f^TL(df)] -\lambda \cdot tr[(df) ^Tf] - \lambda \cdot tr[ f^T(df)] \\ & = tr[f^TL^T(df) ] + tr[f^TL(df)] -\lambda \cdot tr[f^T(df) ] - \lambda \cdot tr[ f^T(df)] \\ &=tr[f^T(L^T+L)(df)- 2\lambda f^T(df)] \\ &=[f^T(L^T+L)- 2\lambda f^T]df \\ &=[2f^TL^T- 2\lambda f^T]df \\ \end{aligned} \tag{3.6} dL(f,λ)=d[fTLf−λ(fTf−n)]=d(fTLf−λfTf)=d(fT)Lf+fTL(df)−λd(fT)f−λfT(df)=(df)TLf+fTL(df)−λ(df)Tf−λfT(df)=tr[(df)TLf+fTL(df)−λ(df)Tf−λfT(df)]=tr[(df)TLf]+tr[fTL(df)]−λ⋅tr[(df)Tf]−λ⋅tr[fT(df)]=tr[fTLT(df)]+tr[fTL(df)]−λ⋅tr[fT(df)]−λ⋅tr[fT(df)]=tr[fT(LT+L)(df)−2λfT(df)]=[fT(LT+L)−2λfT]df=[2fTLT−2λfT]df(3.6)
最后一步 2 L T = L T + L 2L^T=L^T+L 2LT=LT+L,是因为 L = D − W L=D-W L=D−W,其中度矩阵 D D D和邻接矩阵 W W W均为对称矩阵,因此两个对称矩阵的差依然为对称矩阵,即 L T = L L^T=L LT=L②得到导数
由标量微分 d y = f ′ ( x ) d x dy = f'(x)dx dy=f′(x)dx推导出向量微分 d y = ( d y d x ) T d x dy = (\frac{dy}{dx})^Tdx dy=(dxdy)Tdx得d L ( f , λ ) d f = [ 2 f T L T − 2 λ f T ] T = 2 L f − 2 λ f = 0 ⇒ L f = λ f (3.7) \begin{aligned} \frac{d L(f,\lambda)}{d f} &= [2f^TL^T- 2\lambda f^T]^T =2Lf-2\lambda f =0\\ &\Rightarrow Lf=\lambda f \end{aligned} \tag{3.7} dfdL(f,λ)=[2fTLT−2λfT]T=2Lf−2λf=0⇒Lf=λf(3.7)
我们得到了一个惊人的结论!!当 L a g r a n g e Lagrange Lagrange乘子 λ \lambda λ是拉普拉斯矩阵 L L L的特征值(eigenvalue)且指示向量 f f f是 L L L的特征向量(eigenvector)时,函数有极值。那么,这个函数的极值到底是什么呢?我们继续分析,因为目标函数是 f T L f f^TLf fTLf,所以我们可以将等式 ( 3.7 ) (3.7) (3.7)两边分别左乘 f T f^T fT凑成目标函数形式后为
f T L f = λ f T f = λ n f^TLf=\lambda f ^Tf=\lambda n fTLf=λfTf=λn
因为 n = ∣ V ∣ n=|V| n=∣V∣是常数,显然取决于目标函数值大小的元素是 λ \lambda λ,即
min f T L f f T f = min λ (3.8) \min \; \frac{f^TLf}{f^Tf} = \min \; \lambda \tag{3.8} minfTffTLf=minλ(3.8)
也就是说 L L L的特征值 λ \lambda λ越大,目标函数值越大;特征值 λ \lambda λ越小,目标函数值越小!!随着入坑越深,我们已逐渐将连续优化问题转化成了特征分解问题。那么 λ \lambda λ的最小值是多少呢?我们从两个角度来分析①拉格朗日乘子角度
在高数中,我们之所以定义 L a g r a n g e Lagrange Lagrange乘子是为了方便求导,即将有约束的连续问题转化为无约束的连续问题,而定义后的 L a g r a n g e Lagrange Lagrange函数就是原函数的下界,即
min f T L f ≥ min L ( f , λ ) = f T L f − λ f T f ⇒ λ ≥ 0 \begin{aligned} \min \;f^TLf & \ge \min \;L(f,\lambda) \\ & = f^TLf -\lambda f^Tf \\ & \Rightarrow \lambda \ge 0 \end{aligned} minfTLf≥minL(f,λ)=fTLf−λfTf⇒λ≥0
由于 f T f = n > 0 f^Tf =n>0 fTf=n>0,则必然 λ ≥ 0 \lambda \ge 0 λ≥0,因此我们可满心欢喜地得到结论 λ min = 0 \lambda_{\min} = 0 λmin=0。②拉普拉斯矩阵性质角度
我们之前在式 ( 2.8 ) (2.8) (2.8)推导出,对于指示向量 f ∈ R n f\in R^n f∈Rn有
f T L f = 1 2 ∑ i = 1 n ∑ j = 1 n w i j ( f i − f j ) 2 f^TLf = \frac{1}{2} \sum_{i=1}^n \sum_{j=1}^n w_{ij}(f_i-f_j)^2 fTLf=21i=1∑nj=1∑nwij(fi−fj)2
由于 w i j ≥ 0 w_{ij}\ge0 wij≥0,二次项 ≥ 0 \ge0 ≥0,因此 f T L f ≥ 0 f^TLf \ge 0 fTLf≥0。又因为 min f T L f ⇔ min λ \min \; f^TLf \Leftrightarrow \min \; \lambda minfTLf⇔minλ,所以 λ min = 0 \lambda_{\min} = 0 λmin=0。令拉普拉斯矩阵的第二小特征值作为目标函数最优值
那我们是不是直接将 λ min = 0 \lambda_{\min} = 0 λmin=0这个解作为目标函数的解,求出对应特征向量 f f f(指示向量)就完事了呢?毫无疑问,它可以作为一个平凡解,但是否是我们想要的解?
显然不是!因为我们定义的目标函数值是有意义的,我们从最初的 C u t ( V ) Cut(V) Cut(V)到 R a t i o C u t ( V ) RatioCut(V) RatioCut(V)都代表着按某种方式切图后的代价。当这个代价为 0 0 0,则表示我们没有切图,所有的样本都属于一类,这种情况下得到的特征向量 f f f(指示向量)一定是常数向量 1 \boldsymbol 1 1。因此,我们一般不取最小而是取第二小特征值 λ \lambda λ作为最优值,其对应的特征向量 f f f即为我们所求的最优解。
实数解向量的定性
本节到这里就可以结束了吗?显然还不能TAT,别忘了,我们最初是将离散问题转化成连续问题求解的,也就是说当前目标函数所得到的最优解 f f f是任意实数解向量(real-valued solution vector),我们无法从中断言哪个样本确切的属于哪个类别。因此,我们还需要将实数解向量定性。
对于二聚类 ( c = 2 ) (c=2) (c=2)问题,这个很好做,毕竟我们的解向量 f ∈ R n f\in R^n f∈Rn,因此最简单粗暴的定性方法是判断 f i f_i f