讲座笔记：k-sums聚类：一种可替代k-means的高效聚类算法

经典聚类方法

聚类问题是机器学习领域的一个重要的基础问题，已经研究了很多年，也提出了很多方法。

聚类目标

$X=[x_1,\cdots ,x_n]\in \mathbb{R}^{d\times n}$ 为待聚类数据，$x_i\mathbb{R}^{d}$ 为第 $i$ 个样本点。

聚类任务的目标：将 $X$ 划分为不同的类，并使

处于同一类之间的点距离较近（相似度较大）
处于不同类之间的点距离较远（相似度较小）

k-means

k-means算法的优化目标：最小化类内样本点与类均值的距离。
其优化问题可以表示为：

$\begin{equation} min_{A_1,\cdots ,A_c}\sum_{i=1}^c \sum_{x_j\in A_i}\Vert x_j-m_i\Vert_2^2 \end{equation}$

其中，$m_i$ 表示第 $i$ 个类的均值，$c$ 是类的个数。

K-means的优化：Lloyd算法

初始化$c$个类中心
While not converge do:

将每个样本点 $x_i$，分配到与其最接近的中心所属的类别

根据当前的划分情况，重新计算每类的中心（均值)

k-means的矩阵形式

标签指示矩阵：$Y\in \mathbb{R}^{d\times n}$，$n$ 代表数据个数，$c$ 代表类别个数，每行只有一个元素为1，其余为0。
借助标签矩阵 $Y$，k-means 可重新表述为：

$\begin{equation} \mathop{min}_{Y\in \Phi^{n\times c},M} \Vert X-MY^T \Vert_F^2 \end{equation}$

$M=[m_1,\cdots ,m_c]\in \mathbb{R}^{d\times c}$，$m_i$ 表示第 $i$ 类的样本均值
$\Phi^{n\times c}$ 表示所有的标签矩阵的集合

$M$ 和 $Y$ 都是优化变量。当一个公式中存在多个优化变量的时候，我们通常固定一个变量，然后去优化另一个变量，进行交替优化。
在这里固定 $Y$ 去优化 $M$ 就对应了 Lloyd 算法中的更新均值。更新了均值以后再去更新划分就相当于去求 $Y$ 的最优值了。
从这个角度来看，Lloyd 算法是必定收敛的，我们采用交替优化每次优化一个变量，使得目标值越来越小，而这个目标值存在一个下界。

k-means的优势与不足

优点

速度快，所需参数少（只需要聚类个数)。

缺点

适用范围较窄，只能处理线性可分的数据。进一步来说只适合每一类都是一个球型分布的数据。因此它能聚类的数据比较简单，不能太复杂。
这个优化问题是一个NP难问题，聚类结果受初始化影响较大。

虽然有这么多缺点，但是因为简单好用，在实际中k-means属于是用的最多的一个聚类算法。

经典图割聚类算法/谱聚类算法

我们先要将数据进行构图，使用图来进行分析可以直观的建模。数据是一个低维向量，每个节点是一个数据，节点和节点之间的边的权重是由节点向量的欧氏距离就算出来的一个相似度。之所以叫图割是因为在类和类之间的边做切割。

对于图G，寻求一个划分，以使类间节点之间连接权重最小化。

$\begin{equation} \textrm{Ratio Cut:}\qquad\sum_{i=1}^{c}\textrm{cut}(A_i,\bar{A}_i)/|A_i| \end{equation}$ $\begin{equation} \textrm{Normalized Cut:}\qquad\sum_{i=1}^{c}\textrm{cut}(A_i,\bar{A}_i)/\sum_{x_j\in A_i} d_{jj} \end{equation}$ $\begin{equation} \textrm{MinMax Cut:}\qquad\sum_{i=1}^{c}\textrm{cut}(A_i,\bar{A}_i)/ \sum_{x_i,x_j\in A_k} w_{ij} \end{equation}$

不同算法中定义归一化的值不同，有的用类别个数、有的用度数、有的用类内的相似度之和。

优化算法（以N-cut为例）

与k-means类似，借助标签矩阵 $Y$，N-cut的优化目标可写为：

$\begin{equation} \sum_{i=1}^{c} \frac{cut(A_i,\bar{A}_i)}{\sum_{x_j\in A_i}d_{jj}}=\sum_{i=1}^{c} \frac{\sum_{x_i \in A_i,x_j \in \bar{A}_i}w_{ij}}{\sum_{x_j\in A_i}d_{jj}} \end{equation}$ $\begin{equation} =\mathop{min}_{Y\in \Phi^{n\times c}} Tr \left( (Y^T\Delta Y)^{-1/2} Y^T (\Delta-W)Y(Y^T \Delta Y)^{-1/2} \right) \end{equation}$

该问题难以求解，因为$Y$是一个离散的矩阵，每一行只有一个为1，其余都为0，因此传统方法采用先松弛后离散的两阶段求解方式。
令 $F=\Delta^{1/2}Y(Y^T\Delta Y)^{-1/2}$ ，扔掉离散的约束，只保留$F$的正交性质,将问题松弛为

$\begin{equation} \mathop{min}_{F^TF=I,F\in \mathbb{R}^{n\times c}}Tr \left( F^T\Delta^{-1/2}(\Delta-W)\Delta^{-1/2}F \right) \end{equation}$

(8)的解为 $\Delta^{-1/2}(\Delta-W)\Delta^{-1/2}$最小 $c$ 个特征值对应的特征向量。
该类问题都是将 $Y$ 矩阵松弛以后进行特征分解，因特征分解也称谱分解，因此该类算法也叫做谱聚类算法

离散化

将图割问题松弛后，得到的解并非聚类结果，而是数据$X$的嵌入，因此需要借助k-means 或谱旋转等手段，得到最终的聚类结果。

谱聚类算法的优势与不足

优点：

因为使用了一个近邻图来进行分析，因此可以处理非线性可分数据，有着更加出色的聚类性能。

缺点：

计算复杂度较高，一般为$O(n^2 c)$（既要构图，又要进行谱分解）
先嵌入再离散的求解过程，所得解与原问题存在偏差，性能上存在一定的影响
构图的过程中，如何利用距离去把相似度求出来，往往会引入一个参数。构图参数对聚类性能存在较大影响

因为复杂度高，同时参数敏感，所以在实际场景中非常不好用，其使用频率远远没有k-means多。

改进k-means算法

处理非线性可分数据

核化 k-means 算法：将输入数据非线性地映射至高维空间，然后执行k-means算法。

该类算法虽然可以提升k-means的性能，但时间复杂度较高。

复杂度达到了 $O(n^2)$，并不实用。

更好的初始化

k-means++:选取相距较远的点作为初始类中心。

它不但使初始化的结果相对固定，并在一定程度上提高了算法性能。

k-means对初始化比较敏感，k-means++的思想为寻找两两之间相距较远的点作为初始点。首先随机选择一个点，然后就可以计算其他所有点到这个点的距离，距离越近的点被选的概率越小，距离远的被选的概率越大。该方法理论和实际效果都非常好，用的也非常多。

改进谱聚类算法

NysT_rom

令$W$表示相似度矩阵，它首先将$W$表示为

$W=\left[ \begin{array}{cc} A & B \\ B^T & C \end{array} \right]$

其中$A$是可对角化矩阵 $A=U_a\wedge U_a^T$ 。
矩阵$W$的特征向量$U$即可近似计算如下

$U^T=[U_a^T \wedge^{-1}U_a^T B]$

基于NysT_rom的快速谱聚类算法

一些快速谱聚类方法使用NysT_rom方法来近似计算谱聚类中涉及的特征向量。

对数据先进行采样，然后对采样的数据进行特征分解。原矩阵的特征向量对子矩阵的特征向量是有关系的，是一个近似的计算。
该方法有两个缺点：

这种近似计算要满足一些严格的数学假设，这些假设在实际中是难以满足的。所以在很多实际场景中这个近似是误差很大的。

虽然计算特征向量加快了，但是构图还是很慢。首先要给定一个图以后才能进行采样。

Anchor graph

选取$m$个锚点，将相似度矩阵$W$近似表示为

$W=B\Delta^{-1} B^T$

其中$B\in \mathbb{R}^{n\times m}$ 是样本点与锚点之间的相似度矩阵。
矩阵W的特征分解可以用$Z=B\Delta^{-1/2}$的奇异值分解代替

基于Anchor图的快速谱聚类算法
使用锚点图代替原来的相似度矩阵，大大降低了传统谱聚类算法的复杂度。一般来说该类方法复杂度为$O(nm^2+m^3)$。

也是采样 m 个点，但是并不去计算数据之间的近似度，因为一旦计算时间复杂度就是 $O(n^2)$。他是计算 n 个数据数据到 m 个代表点之间的近似度，这样复杂度就从 $O(n^2)$ 变到 $O(n\times m)$ ，所以他的速度的以加快。
同时还有个好处同时 $W$ 这个相似矩阵还是一个低秩的矩阵，后续的特征分解就可以使用奇异值分解进行代替，奇异值分解的复杂度要比特征分解的复杂度小很多。这两步过程都可以做到加速。
所以Anchor graph 用以加速用的多一些。

k-means和谱聚类的本质联系

k-means的矩阵形式

在k-means 的优化问题中，$M$的最优解为$M= XY(Y^T Y)^{-1}$。
据此，k-means的优化问题可等价表示为

$\mathop{min}_{Y\in \Phi^{n\times c},M} \Vert X-MY^T \Vert_F^2 \\ \Leftrightarrow \mathop{min}_{Y\in \Phi^{n\times c}} \Vert X-XY(Y^TY)^{-1}Y^T \Vert_F^2 \\ \Leftrightarrow \mathop{max}_{Y\in \Phi^{n\times c}} Tr\left( (Y^TY)^{-1/2} Y^TX^TXY(Y^TY)^{-1/2} \right)$

Normalized-cut

若使用对称双随机矩阵，N-cut的优化问题可写为

$\begin{equation} \mathop{min}_{Y\in \Phi^{n\times c}} Tr\left( (Y^T\Delta Y)^{-1/2} Y^T(\Delta -W)Y(Y^T\Delta Y)^{-1/2} \right) \end{equation}$ $\begin{equation} \Leftrightarrow \mathop{max}_{Y\in \Phi^{n\times c}} Tr\left( (Y^TY)^{-1/2} Y^TWY(Y^TY)^{-1/2} \right) \end{equation}$

双随机矩阵
$W\in \mathbb{R}^{n\times n}$ 行和为1，列和也为1。即:

$\sum_{i=1}^n w_{ij}=1,\forall j=1,\cdots ,n$

$\sum_{j=1}^n w_{ij}=1,\forall i=1,\cdots ,n$

统一视角

对于k-means，我们有：

$\begin{equation} \mathop{max}_{Y\in \Phi^{n\times c}} Tr\left( (Y^TY)^{-1/2} Y^TX^TXY(Y^TY)^{-1/2} \right) \end{equation}$

对于N-cut，我们有：

$\begin{equation} \mathop{max}_{Y\in \Phi^{n\times c}} Tr\left( (Y^TY)^{-1/2} Y^TWY(Y^TY)^{-1/2} \right) \end{equation}$

不难发现，而这可以统一到以下形式

$\begin{equation} \mathop{max}_{Y\in \Phi^{n\times c}} Tr\left( (Y^TY)^{-1/2} Y^TGY(Y^TY)^{-1/2} \right) \end{equation}$

其中$G=XX^T$表示k-means，$G=W$表示图割模型。

如果是内积矩阵的话就是k-means，如果是相似度矩阵的话就是谱聚类。

k-means 和谱聚类的离散形式

k-means优化问题也可以表示为 $(d_{j1}=\Vert X_j-X_1 \Vert^2_2)$

$\begin{equation} \mathop{min}_{Y\in \Phi^{n\times c},M} \Vert X-MY^T \Vert^2_F \Leftrightarrow \mathop{min}_{A_1,\cdots ,A_c} \sum_{i=1}^{c} \frac{1}{\vert A_i \vert} \sum_{x_j,x_i\in A_i} d_{j1} \end{equation}$

N-cut的优化问题可以表示为

$\begin{equation} \mathop{max}_{Y} Tr\left( (Y^TY)^{-1} Y^TWY \right) \Leftrightarrow \mathop{max}_{A_1,\cdots ,A_c} \sum_{i=1}^{c} \frac{1}{\vert A_i \vert} \sum_{x_j,x_i\in A_i} d_{j1} \end{equation}$

二者模型完全一致，区别在于对相似性（不相似性）的度量方式。
谱聚类算法旨在最大化类内相似度，而k-means旨在最小化类内距离。
使用近邻图是谱聚类算法表现出良好性能的原因之一，但同时也导致了模型不易优化。

k-sums聚类方法

平衡数据（约束条件）

假设每类样本数相等，上述统一框架可等价为：

$\begin{equation} \mathop{max}_{Y\in \Phi^{n\times c}} Tr\left( (Y^TY)^{-1/2} Y^TGY (Y^TY)^{-1/2} \right) \qquad s.t.Y^TY=\bar{n}I \end{equation}$ $\begin{equation} \Leftrightarrow \mathop{max}_{Y\in \Phi^{n\times c}} Tr\left( Y^TGY \right) \qquad s.t.Y^TY=\bar{n}I \end{equation}$

因为 $-1/2$ 的存在导致问题比较难解，我们需要进行松弛和离散的操作。有一个想法是能否将这个 $-1/2$ 去掉，如果要去掉的话就要加入上面这个约束，其物理意义是每一类的数据个数都要一样。

Model

上述问题难以求解，因此我们舍弃了平衡约束，此时问题变为

$\begin{equation} \mathop{max}_{Y\in \Phi^{n\times c}} Tr\left( Y^TGY \right) \end{equation}$

模型一

$\begin{equation} \mathop{max}_{Y\in \Phi^{n\times c}} Tr\left( Y^TGY \right) \end{equation}$

不难发现，若$G=W$，则$Y$的任何一列全为1，其他列为零，都可以使问题(19)取得最大值。因此该问题存在$c$个平凡解。

模型二

有趣的是，若$G=D$，即采用不相似性度量，则可避免平凡解。
此时问题为:

$\begin{equation} \mathop{max}_{Y\in \Phi^{n\times c}} Tr\left( Y^TDY \right) \end{equation}$

这个问题又存在一个没有意义的解，他的最优解会使得Y的任何一列全为1，其他列都为0，即所有数据都聚成一类，会使得目标值达到最大，这就没有意义了。但将G换成距离就没有这个解了。

舍弃了平衡约束的模型是否有意义?
当每一类个数一样的时候是他的最优解，即使舍弃了平衡约束，模型依旧倾向于平衡划分。

平衡划分

在模型(20)中包含样本数过多的类会导致一个较大的目标函数值。具体来说：假设样本之间的距离近似相等，为$\alpha$。则有

$\begin{equation} \mathop{min}_{Y\in \Phi^{n\times c}} Tr\left( Y^TDY \right) \end{equation}$ $\begin{equation} \approx \mathop{min}_{Y\in \Phi^{n\times c}} Tr\left( Y^T\alpha(11^T-I)Y \right) \end{equation}$ $\begin{equation} \Leftrightarrow \mathop{min}_{Y\in \Phi^{n\times c}} \alpha (n_1^2+\cdots +n_c^2-n) \end{equation}$

$n_i=n/c,\forall i=1,\cdots,c$为问题(23)的最优解。

k-sums模型

受谱聚类算法启发，k-sums使用了一个近邻距离图$\tilde{D}$代替全距离矩阵$D$。
也就是说，我们的模型最终形式为

$\begin{equation} \mathop{min}_{Y\in \Phi^{n\times c}} Tr(Y^T\tilde{D}Y) \end{equation}$

近邻距离矩阵
对于样本点$x_i$，若$x_j$是样本点$x_i$的k近邻，或$x_i$是$x_j$的k近邻，则近邻距离矩阵$\tilde{D}$中，$\tilde{d}_{ij}=\Vert x_i-x_j\Vert^2_2$。否则，$\tilde{d}_{ij} = \gamma$。其中$\gamma$为所有近邻距离的最大值。

一开始这个D是一个距离矩阵，两两点之间的距离都要给定。使用传统全连接的距离图的话和k-means一样只能聚类线性可分的数据，但是现在是用那个近邻距离图的话则可以聚类更为复杂分布的数据。
使用近邻图的话模型的优化也会变得更为快速，甚至和类别个数无关，当类别个数特别大的时候该方法非常有优势。

模型优化

因为模型公式非常简单，所以优化也相对更为简单，可以采用标准坐标下降优化算法。
我们要优化 $Y$ ，而 $Y$ 是一个 $n$ 行的矩阵，我们可以一行一行的更新，比如只更新第一行，这个优化问题就是去数这个里面的最小值对应的索引。因为 $Y$ 一行只有一个为1其他都为0。

采用标准坐标下降优化算法

$\begin{equation} \mathop{min}_{Y_1} Tr(Y^T\tilde{D}Y) \end{equation}$ $\begin{equation} \Leftrightarrow \mathop{min}_{Y_1} \left( \left[ \begin{array}{cc} y_1 & Y_r \\ \end{array} \right] \left[ \begin{array}{cc} \tilde{d}_{1,1} & v^T \\ v & D_r \end{array} \right] \left[ \begin{array}{c} y_1^T \\ Y_r^T \end{array} \right] \right) \end{equation}$ $\begin{equation} \Leftrightarrow \mathop{min}_{Y_1} Tr( y_1 \tilde{d}_{1,1} y_1^T ) +2Tr(Y_r v y_1^T) \end{equation}$ $\begin{equation} \Leftrightarrow \mathop{min}_{Y_1} Tr(Y_r v y_1^T) \end{equation}$ $\begin{equation} \Leftrightarrow \mathop{min}_{Y_1} y_1^T \tilde{Y}^T d_1 \end{equation}$

$\tilde{d}_1\in \mathbb{R}^n$ 是样本 $x_1$ 与所有样本的平方欧氏距离。

问题(29)的解为

$\begin{equation} y_{ik}=\begin{cases} 1 & k=\textrm{argmin}_j(Y^T \tilde{d}_i)_j \\ 0 & \textrm{otherwise} \end{cases} \end{equation}$

其中$Y\in \mathbb{R}^{n\times c},\tilde{d}_i\in \mathbb{R}^n$

计算复杂度

计算 $Y^T\tilde{d}_i$ 需要 $O(nc)$ 时间，求其最小值需要 $O(c)$ 时间。
考虑 $Y$ 的稀疏性，则 $Y^T\tilde{d}_i$ 只需要 $O(n)$ 时间。
利用 $\tilde{d}_i$ 中大多值为常量 $\gamma$ 以及 $Y$ 的特点，可以进一步加速为 $O(k$)。

$\tilde{d}_i$ 是第 $i$ 个数据到其他所有数据距离组成的一个向量，因为这个向量是一个近邻图，所以只有 $k$ 个保留了他的距离，其他都为 $\gamma$ 常量。利用这个特点我们可以进一步加速，其复杂度从 $O(m)$ 变到 $O(k)$。

例子

假设需要聚类6个样本，当前聚类结果是：
1，2号为1类；3，4号各为1类；5，6号为1类。
假设1号样本与2号，3号为近邻；4，5，6号不是1号的近邻。
当更新 $y_1$ 时，令 $t=\tilde{d}_1^T Y$，则 $t_1,t_2$ 包含近邻值（红色方块），$t_3,t_4$只包含 $\gamma$ （蓝色圆）。