The SpectACl of Nonconvex Clustering: A Spectral Approach to Density-Based Clustering

摘要

1. 最小割和最大密度经常用来对非凸数据聚类。
2. 最小割对噪音敏感。
3. 基于最大密度的聚类方法很难聚类聚簇中心们的密度不同的情况。
4. 本文介绍的SPECTACL方法集两者之长，去两者之短。
5. 本方法有谱聚类的易于实现的特点，并且用理论证明给出了合理的聚簇中心密度标准。

一、简介

1. 谱聚类通过对数据相似性矩阵的频谱进行降维，解决了聚类任务公式中固有的维数诅咒。DBSCAN是基于密度的算法。两者都能聚类非线性可分的数据。 本文介绍的SPECTACL方法集两者之长，去两者之短。
2. SpectACL 发现聚簇中心的平均密度很大。还发现每个中心的密度由加权邻接矩阵的谱（矩阵特征值的集合）决定了。因此SpectACL不像DBSCAN一样对参数minPts（DBSCAN里影响密度的）敏感。并且还能解决各个聚簇中心密度不一致的问题。
3. SpectACL在embedding（投影到低维向量空间中）后也有个用k-mean的步骤，本文证明了应用K-mean的合理性：通过特征向量分解法我们得到了SpectACL的目标函数的上界。同时我们得出优化目标函数的上界和用k-mean对特征向量聚类是等价的。（谱聚类没有类似证明）。

二、聚类方法讨论

符号

1. 数据矩阵: $D\in R^{m\times n}$。m个样本，每个样本n个特征。
2. ϵ − n e i g h b o r h o o d : 对 每 个 点 j ， N ϵ ( j ) = { l ∣ ∥ D j ⋅ − D l ⋅ ∥ < ϵ } \epsilon-neighborhood :对每个点j，\mathit{N}_{\epsilon}(j)=\{l| \parallel D_{j\centerdot} - D_{l\centerdot} \parallel <\epsilon \} ϵ−neighborhood:对每个点j，Nϵ​(j)={l∣∥Dj⋅​−Dl⋅​∥<ϵ}
3. Forbenius norm : $\parallel M\parallel$.(各项元素的绝对值平方的总和的平方根)
4. d维低秩近似 ：$M\approx V^{(d)}\Lambda V^{(d)}$
5. $1^{m\ast n}$ : m*n的01矩阵。每行只有一个1
6. 图的权重矩阵W : 好多种表示方法，本文后面用的$ϵ-neighborhood$方法好像
7. 图的度数矩阵G：对角矩阵，每个元素为该点的权重和，也就是$G_{ii}={\sum }_{j=1}^m{W}_{ij}$
8. 拉普拉斯矩阵L : L=G-W

2.1 k-mean

1. k-mean的流行流行原因：(1). 简单有效（2）理论上建立了优化方法（循环迭代，交替最小）（3）由聚簇内的点的分散可以直观理解聚簇的意义

2. k-mean潜藏的假设是每个聚簇里的所有点彼此相似。后果就是只能处理凸形数据。

3. k-mean的目标方程:
   $$\underset{Y,X}{\min }\parallel D-Y{X}^T{\parallel }^{2}s.t.Y\in 1^{m\times r},X\in {\mathbb{R}}^{n\times r}.(1)$$
   $Y$为m*r的01矩阵，当第i个样本属于第j类时，$Y_{ij}=1$。$X$代表聚簇中心，每列为一个聚簇中心。$Y{X}^T$为$m\times n$的矩阵，每行为原始数据点所属于的聚簇中心。整个式子的意义就是整个数据集合聚簇中心方差。

   这个目标函数是非凸的（看上去也不像=_=）。但是$Y$或者$X$一旦确定了其中一个，上述目标函数就是凸的了（如果只看一个样本点，随着Y变化聚簇中心会变化，$X$中的中心又是排了序的(吧。。。)，大概是二次方的，所以大概是凸的(吧。。。)）。

   如果$Y$确定了，那么最优的$X$应该取$X=D^{T}Y(Y^{T}Y{)}^{-1}$。$Y^{T}Y$是一个r*r对角矩阵，对角线上是每个行的內积，也就是属于该聚簇中心的点的数量。$(Y^{T}Y{)}^{-1}$就是$Y^{T}Y$的每个元素取倒数，大概是个比率的意思。$Y(Y^{T}Y{)}^{-1}$就是$Y$矩阵每行的1变成对应的那个倒数，这样$Y(Y^{T}Y{)}^{-1}$每列的和为1。$D^{T}Y(Y^{T}Y{)}^{-1}$就是数据集中的每个数据点除该数据点属于的聚簇的点的个数，再加和。也就是相当于对D中的属于每个聚簇中心的数据点取平均，得到聚簇中心向量。

4. 上述(1)式等价于：
   $$\underset{Y}{\max }tr(Y^{T}D{D}^{T}Y(Y^{T}Y{)}^{-1})s.t.Y\in 1^{m\times r}(2)$$
   这个式子中$D^{T}Y(Y^{T}Y{)}^{-1}$的意义已经知道了，也就是相当于$\underset{Y}{\max }tr(Y^{T}DX)$。$Y^{T}D$是一个$r\times n$的矩阵，每行为该聚簇的所有点的和。所以这个式子就是求对角线的最大和，该对角线每个元素为某一聚簇的所有点的累加和与聚簇中心的內积。

   注意到这个目标函数只和数据点的內积有关，也就是只和数据点的相似度有关，这个是kernel k-mean以及非凸聚簇的关键。

2.2 图的切割

1. 通过将数据点转换为核矩阵所反映的高维希尔伯特空间，可以克服凸簇的限制。由上述那个式子只和数据点的相似度有关的结论可以推出数据的图表示法：每个样本点对应一个节点，点之间的相似度对应于边的权值。这样，这个图中每个紧密连通的块就可以当做一个聚簇。分类的代价就是求割集的代价。也就是最小割问题:（最小割不是用网络流求的吗？？？）
   $$cut(Y;W)=\sum _{s=1}^r{Y}_{\cdot s}^{T}W(1-Y_{\cdot s})(3)$$
   上面这个式子效果一般不好，可能出现有的聚簇只有几个点得情况，于是可以加个约束：
   $$Rcut(Y;W)=\sum _{s=1}^r\frac{Y_{\cdot s}^{T}W(1-Y_{\cdot s})}{\mid Y_{\cdot s}\mid }(4)$$
   式子（4）又等价于：
   $$\underset{Y}{\max }tr(Z^T(-L)Z)s.t.Z=Y(Y^{T}Y{)}^{-1/2}(5)$$

2. 等价的推导：
   Z的意义 ：$m\times r$的矩阵，每列的意义是一个指示向量，当第$i$行数据点属于第$j$个聚簇时，那个位置的值为$\frac{1}{\sqrt{\mid A_j\mid }}$，$A_j$代表第j类聚簇的点的集合，$\mid A_j\mid$代表该集合点的数量。
   $$\begin{array}{rl}{Z}_{\cdot i}^{T}L{Z}_{\cdot i}& =Z_{\cdot i}^{T}G{Z}_{\cdot i}-Z_{\cdot i}^{T}W{Z}_{\cdot i}\\ & =\sum _{j=1}^m{G}_{jj}{Z}_{ji}^2-\sum _{j=1}^m\sum _{k=1}^m{W}_{jk}{Z}_{ji}{Z}_{ki}\\ & =\frac{1}{2}(\sum _{j=1}^m{G}_{jj}{Z}_{ji}^2-2\sum _{j=1}^m\sum _{k=1}^m{W}_{jk}{Z}_{ji}{Z}_{ki}+\sum _{k=1}^m{G}_{kk}{Z}_{ki}^2)\\ & =\frac{1}{2}(\sum _{j=1}^m\sum _{k=1}^m{W}_{jk}{Z}_{ji}^2-2\sum _{j=1}^m\sum _{k=1}^m{W}_{jk}{Z}_{ji}{Z}_{ki}+\sum _{k=1}^m\sum _{j=1}^m{W}_{kj}{Z}_{ki}^2)\\ & =\frac{1}{2}\sum _{j=1}^m\sum _{k=1}^m{W}_{jk}(Z_{ji}-Z_{ki}{)}^2\\ & =\frac{1}{2}(\sum _{j\in A_i,k\notin A_i}{W}_{jk}(\frac{1}{\sqrt{\mid A_i\mid }}{)}^2+\sum _{j\notin A_i,k\in A_i}{W}_{jk}(-\frac{1}{\sqrt{\mid A_i\mid }}{)}^2)\\ & =\frac{Y_{\cdot i}^{T}W(1-Y_{\cdot i})}{\mid A_i\mid }\end{array}$$

   所以（4）和（5）相当于多了个负号，最小（4）就是最大（5）。

3. $L$是半正定矩阵，如果将（5）中的$-L$换成正定矩阵$(\mid {\lambda }_m\mid I-L)$不会改变优化目标。$\mid {\lambda }_m\mid$是$L$的最大的特征值的绝对值。

   由上面那个推导的第五行可以看出L是半正定的。
   $(\mid {\lambda }_m\mid I-L)$是正定的推导：
   $$\begin{array}{rl}{Z}_i^T(\mid {\lambda }_m\mid I-L)Z_i& =Z_i^T\mid {\lambda }_m\mid I{Z}_i+\frac{1}{2}\sum _{j=1}^m\sum _{k=1}^m{W}_{jk}(Z_{ij}-Z_{ik}{)}^2\\ & =\mid {\lambda }_m\mid Z_i^T{Z}_i+\frac{1}{2}\sum _{j=1}^m\sum _{k=1}^m{W}_{jk}(Z_{ij}-Z_{ik}{)}^2\\ & =\mid {\lambda }_m\mid +\frac{1}{2}\sum _{j=1}^m\sum _{k=1}^m{W}_{jk}(Z_{ij}-Z_{ik}{)}^2\\ & >0\end{array}$$
   至于为啥这个替换不影响目标函数的解是因为$Z_i^T{Z}_i=I$
   因此ratio-cut和用$\mid {\lambda }_m\mid I-L$作为核矩阵的kernel-kmeans就等价了。也就是用$\mid {\lambda }_m\mid I-L$代替（2）中的內积之后（2）和（4）（5）的优化目标相同。
   证明：
   $$\begin{array}{rl}& Z^T=(Y(Y^{T}Y{)}^{-1/2}{)}^T=(Y^T\ast Y{)}^{-1/2}{Y}^T\\ & 所以\\ & \underset{Y}{\max }tr(Z^T(-L)Z)=\underset{Y}{\max }tr((Y^T\ast Y{)}^{-1/2}{Y}^T(-L)\ast Y(Y^{T}Y{)}^{-1/2})\\ & 令A=(Y^T\ast Y{)}^{-1/2}，B=Y^T(-L)\ast Y:\\ & \underset{Y}{\max }tr(Z^T(-L)Z)=\underset{Y}{\max }tr(ABA)=\sum _{i=1}^m{A}_{ii}^2{B}_{ii}\\ & \\ & \mathbf{k}\mathbf{e}\mathbf{r}\mathbf{n}\mathbf{e}\mathbf{l}\mathbf{k}-\mathbf{m}\mathbf{e}\mathbf{a}\mathbf{n}:\\ & \underset{Y}{\max }tr(Y^T(\mid {\lambda }_m\mid I-L)Y(Y^{T}Y{)}^{-1})\\ & \underset{Y}{\max }tr(Y^T(-L)Y(Y^{T}Y{)}^{-1})\\ & 令C=(Y^{T}Y{)}^{-1}=A^2:\\ & \underset{Y}{\max }tr(BC)=\sum _{i=1}^m{B}_{ii}{C}_{ii}\end{array}$$
   这就证明出来了。

2.3 谱聚类

1. 可以证明目标函数（5）的解是$Z$取$(\mid {\lambda }_m\mid I-L)$的前r个特征向量。有了$Z$以后还要将它离散成$Y$, 通常用k-mean来实现。
2. 拉普拉斯矩阵特征值为零的特征向量表示连通着的各个聚簇（暂时令$1^{m\times 1}$表示元素全为1的向量。因为$L$的行和为0，所以$L\times 1^{m\times 1}=0\times 1^{m\times 1}$，所以0对应的特征向量应该是$1^{m\times 1}$，也就是所有元素都属于这个簇，表示各个聚簇都是连通的整体），这个特性给出了能embedding的理论证明（这个0特征值对应的全1的特征向量肯定是能降维降掉的，其它的为啥呢？？？，咋就提供理论证明了）。
3. 算法流程：
   • 选择一种合适的方法计算W，并计算L
   • 计算转换了的拉普拉斯矩阵的前r+1个特征值 : $\mid {\lambda }_m\mid I-L=V^{(r+1)}\Lambda V^{(r+1{)}^T}$
   • 去掉$V^{r+1}$第一个特征向量，用k-means返回r个聚簇。
4. $W$矩阵的确定方法：一种是节点度数，另一种是$ϵ-neighborhood$方法， 还有k近邻法以及全连接法。
   还建议对$W$归一化，$\tilde{W}=G^{-1/2}W{G}^{-1/2}$，对应的转换拉普拉斯矩阵是$L_{sym}=I-\tilde{W}$。这种叫N-cut切图，是目标函数（4）的变种，不再除以每个聚簇的点的个数，而是除以这个聚簇里每个点的度数和。具体的不想详细写了。
5. 谱聚类这一系列方法的一个缺点是优化目标不明确。上文介绍的那个ratio-cut，也就是目标函数（4）（5），和谱聚类的优化目标在某些样例上存在分歧。在分歧的时候就不能用k-mean了，也有论文提出了其他方法来离散化但是都没有理论的证明，

2.3 鲁棒性问题

谱聚类的噪声敏感度一直是一个问题，这可能归因于它与最小割范式的关系，在最小割范式中，很少的边就足以让两个集群作为一个出现。
光谱聚类研究的最新进展包括两种提高对噪声鲁棒性的方法。
一种方法是结合基于密度的集群的思想，例如要求集群中的每个节点具有最小的节点度。最小化受最小度约束的切割类似于寻找由DBSCAN执行的核心点的连接组件。该方法的缺点是引入了一个新的参数，影响结果的质量。
另外一种方法是聚类和图的表示同时搞，在每个迭代里更新聚类和矩阵表示。但是在每个迭代里需要计算特征向量，计算代价很大。
最后还有一种是每次选出一个最优聚簇，然后每次一次减少数据集大小。

三、Spectral Averagely-Dense Clustering

1. 这里提出一种基于 聚簇所组成的子图 的聚簇的定义，因为本算法关心节点之间度数的差，我们用$ϵ-neighborhood$方法定义矩阵$W$。

2. 我们的目标是最大化平均聚簇密度：
   $$\underset{Y\in 1^{m\times r}}{\max }tr(Y^{T}WY(Y^{T}Y{)}^{-1})=\sum _s\delta (Y_{\cdot s},W)(6)$$
   $Y^{T}WY$是度数和，$(Y^{T}Y{)}^{-1}$是点的数量和的倒数，$\delta (Y_{\cdot s},W)$就是子图的平均度数。
   $$\delta (Y_{\cdot s},W)=\frac{Y_{\cdot s}^{T}W{Y}_{\cdot s}}{\parallel Y_{\cdot s}{\parallel }^2}=\frac{1}{\mid Y_{\cdot s}\mid }\sum _{j:Y_{js=1}}{W}_{j\cdot }{Y}_{\cdot s}(7)$$
   目标函数（6）等价于$W$正则化的最小割：因为对应的拉普拉斯方程是$L=I-\tilde{W}$，又因为公式（6）可以让拉普拉斯方程去掉I，所以等价。
   对应的拉普拉斯方程为$L=I-\tilde{W}$的证明：
   $$\begin{array}{rl}& L_1表示原来的拉普拉斯方程\\ & L\\ & =G^{-1/2}{L}_1{G}^{-1/2}\\ & =G^{-1/2}(G-W)G^{-1/2}\\ & =G^{-1/2}G{G}^{-1/2}-G^{-1/2}W{G}^{-1/2}\\ & =I-\tilde{W}\end{array}$$
   $L=G^{-1/2}{L}_1{G}^{-1/2}$是让每个元素为$\frac{L_{1_{ij}}}{\sqrt{G_{ii}{G}_{jj}}}$，目的是为了让N-cut的指示向量正交，具体不详细写。

3. 下面我们将得到公式（6）与$W$的谱的关系，通过让k-mean和spectral embedding联系起来。最终表明本文的方法和谱聚类有相同的步骤，可以高效处理大规模数据。

3.1 密集簇与特征向量映射

1. 公式（7）的$\delta$函数就是通常所说的瑞利商。这个商或者说系数有几个性质：
   (1). 范围在$W$的最大和最小特征值之间。
   证明：
   $$\begin{array}{rl}& W=V\Lambda V^T,p=V^{T}x\\ & \delta (x,W)=\frac{x^{T}Wx}{x^{T}x}=\frac{x^{T}V\Lambda V^{T}x}{x^{T}x}=\frac{p^T\Lambda p}{x^{T}x}=\frac{{\sum }_{i=1}^m{\lambda }_i{p}_i^2}{{\sum }_{i=1}^m{x}_i^2}\\ & 因为{\lambda }_1\le {\lambda }_2...\le {\lambda }_m\\ & \frac{{\lambda }_1{\sum }_{i=1}^m{p}_i^2}{{\sum }_{i=1}^m{x}_i^2}\le \delta (x,W)\le \frac{{\lambda }_m{\sum }_{i=1}^m{p}_i^2}{{\sum }_{i=1}^m{x}_i^2}\\ & \sum _{i=1}^m{p}_i^2=\sum _{i=1}^m(V_{\cdot i}^{T}x{)}^2=\sum _{i=1}^m{x}^T{V}_{\cdot i}{V}_{\cdot i}^{T}x=x^{T}x=\sum _{i=1}^m{x}_i^2\\ & {\lambda }_1\le \delta (x,W)\le {\lambda }_m\end{array}$$
   (2). $\delta (x,W)=\delta (cx,W)$,$c$是非0常数。
   $$\delta (cx,W)=\frac{c{x}^{T}Wcx}{c{x}^{T}cx}=\delta (x,W)$$
   (3). 特征值们是瑞利商的极值或者驻点。 $\delta (x,W)$对$x$求导，懒得写证明了。

2. 前d个特征向量张起来的空间里所有向量的瑞利商都$\ge {\lambda }_{d+1}$。
   证明：
   $$\begin{array}{rl}& 设向量c\in R^{d\times 1}为各个正交基的系数，不为0向量，那么：\\ & \delta (c{V}^d,W)=\sum _{i=1}^m\delta (c_i{V}_i^d,W)=\sum _{i,c_i\ne 0}^m\delta (V_i^d,W)=...\ge {\lambda }_m\\ & 懒得写了\end{array}$$
   所以第一大的特征值对应的特征向量就能张出一个密度大于${\lambda }_2$的空间来投影$W$, 让我们来搜索最优解。这个发现引出一个初步的方法：用$W$的特征向量将$W$降维，也就是把（6）(7)改成：
   $$\underset{Y\in 1^{m\times r}}{\max }\sum _s\frac{Y_{\cdot s}^T{V}^d{\Lambda }^d{V}^d{Y}_{\cdot s}}{\mid Y_{\cdot s}\mid }(8)$$
   比较（8）和（2），发现当（2）中的$D=U=V^r({\Lambda }^r{)}^{-1/2}$时与（8）一样。

   然而但是对$U$用k-means得不到期望的聚簇，因为目标函数的特征向量都是正交的，所以是非凸的有很多局部最优解。

3. 特征向量的元素有正数也有负数。
   第一特征向量具有较高的密度函数值δ,但向量元素有正有负使得概念性的解释这个向量困难（如果全是正数就可以把特征向量的每个元素比如${\lambda }_{ij}$看成第i个元素在第j个聚类上的权重啥啥的，大概能有个关系,称作模糊聚类指标。）。将k-means应用于一组高度密集的模糊聚类指标，将进一步投影包含聚簇中心信息的向量到正象限空间。
   一种可能性是通过对称的非负矩阵分解近似矩阵W，但这用np难问题代替了多项式可解的特征分解
   我们进行了另一种方法，表明模糊聚类指标是由一个简单的投影特征值派生出来的：将特征向量取绝对值。。。。
   当然论文里不是这么写的，是这么写的：$W$是对称实数矩阵，$v,\lambda$是$W$的一对对应的特征向量和特征值。让$v=v^{+}-v^{-}，v^{+},v^{-}\in {\mathbb{R}}_{+}^{m\times 1}$。那个非负向量$u=v^{+}+v^{-}$有密度$\delta (u)\ge \mid \lambda \mid$
   证明（终于有个证明是给我的了，不用自己搞）：
   $$\begin{array}{rl}& \mid \lambda \mid u_j=\mid \lambda \mid \mid v_j\mid =\mid W_{j\cdot }v\mid =\mid \sum _l{W}_{jl}{v}_l\mid \le \sum _l{W}_{jl}\mid v_l\mid =\sum _l{W}_{jl}{u}_l=W_{j}u\\ & \mid \lambda \mid u_j\le W_{j}u\mid \lambda \mid u_{\le }Wu\\ & u^T\mid \lambda \mid u_{\le }{u}^{T}Wu\\ & \delta (u,W)=\frac{u^{T}Wu}{u^{T}u}\ge \mid \lambda \mid \end{array}$$
   称$u$这种特征向量取绝对值的向量叫做推测投影特征向量（projected eigenvectors）。

        \space\space\space\space     这个是原始数据的亚子。
   
   这个是投影特征向量的亚子，很神奇。于是可以认为这种向量是模糊聚类指标。
   然后观察图片可以认为发现投影特征向量是非正交的，因为各个图片有交集，这样可以增加聚类的鲁棒性。

4. SPECTACL算法流程总结：

   1. 计算邻接矩阵$W$
   2. 计算截断特征成分$W\approx V^d{\Lambda }^d{V^d}^T$
   3. 计算投影特征向量$U_{jk}=\mid V_{jk}^d\mid \mid {\lambda }_k{\mid }^{-1/2}$
   4. 用k-means对$U$聚类出$r$个聚簇。

5. 当$d=m$算法取复杂度的上界。
   绝大部分情况$d\ll m$的时候，我们优化的优化目标的近似函数。
   根据式子（8），可以讨论本算法的分块计算编码方法。
   谱聚类会在$d\to m$的时候得不到好结果，而本方法的实验结果证明平均下来本算法不受影响。
   $W$的计算方法不受限制（但是要实对称吧）。谱聚类的k近邻法就可以用。
   本方法只有参数$k$和$ϵ(有时候是k近邻的k)$。
   不能决定$k$挺伤的。
   下面有一种启发式决定$ϵ$的方法。还有一种计算$k和ϵ$敏感度的方法。

四、实验

1. 在很多综合生成的数据集上实验。真实数据上比较Normalized Mutual Information。并利用手写数据集对所得到的聚类进行了定性分析。

2. Mutual Information，互信息，MI， 表示两个变量X与Y是否有关系，以及关系的强弱。定义 为：
   $$I(X;Y)=\sum _{x\in X}\sum _{y\in Y}p(x,y)log(\frac{p(x,y)}{p(x)p(y)})$$
   $p(x,y)、p(x)、p(y)$分别是联合密度函数和边缘密度函数。
   互信息量$I(xi;yj)$在联合概率空间$P(XY)$中的统计平均值。 平均互信息$I(X;Y)$克服了互信 息量$I(xi;yj)$的随机性,成为一个确定的量。

   直观上，互信息度量 X 和 Y 共享的信息：它度量知道这两个变量其中一个，对另一个不确定 度减少的程度。例如，如果 X 和 Y 相互独立，则知道 X 不对 Y 提供任何信息，反之亦然，所以它 们的互信息为零。在另一个极端，如果 X 是 Y 的一个确定性函数，且 Y 也是 X 的一个确定性函数，那么传递的所有信息被 X 和 Y 共享：知道 X 决定 Y 的值，反之亦然。因此，在此情形互信息与 Y（或 X）单独包含的不确定度相同，称作 Y（或 X）的熵。而且，这个互信息与 X 的熵和 Y 的熵相同。

   python有库，要是因为数据量太大库用不了再自己写。

3. 实验比较五种算法:
   $ϵ-neighborhood$的SPECTACL。$ϵ$要让90%的点都有是个邻居以上。
   正则化的基于$KNN$的SPECTACL（$G_{knn}^{-1/2}{W}_{knn}{G}_{knn}^{-1/2}$）。$k$取10.
   DBSCAN，$minPts=10$
   normalized Spectral Clustering ，就是n-cut的谱聚类好像。
   Robust Spectral Clustering，RSC，$d=10$

4. 用算法得到的聚类结果和Ground Truth算$F-mesure$作为算法的评估标准

5. 对噪音的敏感度：。。。

6. 一种启发式的确定$ϵ$的方法