记录一些基本的基本聚类方法的入门笔记，以便后续查找。阅读材料中的工作仅截止到2007年。后续的一些新进展没有考虑。谱聚类本来也想多写一点，然而感觉对来龙去脉还没有对混合密度这么清楚，等实现了几个算法之后再说吧。

直觉上来说，聚类本身是一种对数据相似的衡量，同一类之间尽量相似，不同类之间尽量不一样。但由于没有标注，哪怕使用相同的数据集，对于“类”的定义根据使用场景也很可能不同。

1. 混合密度方法

在使用混合密度的想法中，我们简单假设，一个样本一旦属于某一个类，那么它就是从这个类自己的分布中产生的，每个类可以拥有不同的分布、不同的形式，并且不同类之间的分布相互独立，不会互相影响。而每个样本属于某个类的先验概率也是独立的.
这部分的详细内容可参见 [1]，这里采用与其相似的记号.

  • 样本数据集合 \(D = \{x_k \mid k = 1, 2, \dots, n\}\)，共 n 个
  • 类别集合 \(\{w_j \mid j = 1, 2, \dots, c\}\), 共 c 个
  • 每个类的先验概率 \(\alpha = \{ \alpha_j = P(w_j) \mid j = 1, 2, \dots, c \}\)
  • 每个类的概率密度 \(p(x \mid w_j, \theta)\)，其中 \(\theta = \{\theta_j \mid j = 1, 2, \dots, c\}\)

同时由于上述假设存在，不同类的分布相互独立，因此有

\[p(x \mid w_j, \theta) = p(x \mid w_j, \theta_j)\]

用简单好记的说法，就是当 \(w_j\) 出现在条件部分时， \(\theta\) 和 \(\theta_j\) 可以互换，否则就只能用 \(\theta\) 表示

混合密度

由于前面的假设可以得出整体的样本的密度函数：

\[p(x \mid \theta) = \sum_{j=1}^c P(w_j) p(x \mid w_j, \theta_j)\]

如果我们知道密度函数的形式（比如待会儿我们要假设它是高斯密度），那么对上面这给整体的分布中得出的一个独立同分布数据集，就可以使用参数估计的办法 [2] 来求得这个分布中每个未知参数. 假设使用简单的最大似然方法，容易得到它的对数似然函数：

\[\begin{align} l(\theta) &= p(D \mid \theta) \\ &= \prod_{k=1}^n p(x_k \mid \theta) \\ ll(\theta) &= \ln l(\theta) \\ &= \sum_{k=1}^n \ln p(x_k \mid \theta) \\ &= \sum_{k=1}^n \ln (\sum_{j=1}^c P(w_j) p(x \mid w_j, \theta_j)) \end{align}\]

为了让似然函数最大，令其对于各变量导数为零（这里就是关于各个参数 \(\theta_j\) 导数为零），同时由于个分布之间独立，对一个参数求导时其他项直接省略，那么有

\[\begin{align} 0 &= \frac{\partial}{\partial \theta_j} ll(\theta) \\ &= \sum_{k=1}^n \frac{1}{p(x_k \mid \theta)} \frac{\partial}{\partial \theta_j} \sum_{j=1}^c P(w_j) p(x \mid w_j, \theta_j) \\ &= \sum_{k=1}^n \frac{1}{p(x_k \mid \theta)} \frac{\partial}{\partial \theta_j} P(w_j) p(x \mid w_j, \theta_j) \\ &= \sum_{k=1}^n \frac{P(w_j)}{p(x_k \mid \theta)} \frac{\partial}{\partial \theta_j} p(x \mid w_j, \theta_j) \\ &= \sum_{k=1}^n \frac{P(w_j)}{p(x_k \mid \theta)} p(x \mid w_j, \theta_j) \frac{\partial}{\partial \theta_j} \ln p(x \mid w_j, \theta_j) \\ &= \sum_{k=1}^n \frac{p(x, w_j \mid \theta)}{p(x_k \mid \theta)} \frac{\partial}{\partial \theta_j} \ln p(x \mid w_j, \theta_j) \\ &= \sum_{k=1}^n p(w_j \mid x_k, \theta) \frac{\partial}{\partial \theta_j} \ln p(x \mid w_j, \theta_j) \\ \end{align}\]

最后的形式可以简单记为：对每个样本，其后验乘以对数似然的偏导，等于整体的似然函数的偏导.

此外，上面是假设我们知道了先验的概率 \(P(w_j)\) 的情形，如果连先验也是未知的，即我们一开始也不知道到底类别的概率分布如何，那么就可以扩充一下，把先验的概率也作为未知数之一，因此也需要求关于它的偏导并令之为0。但是，光求偏导是没意义的，因为它们自由度并不是 c，即不同类概率之和为1，所有的先验必须满足这一条件，即求：

\[\max ll(\theta, \alpha) \\ s.t. \sum_{j=1}^c \alpha_j = 1\]

构造拉格朗日函数：

\[LL(\theta, \alpha, \lambda) = ll(\theta, \alpha) + \lambda(\sum_{j=1}^c \alpha_j - 1)\]

对这个式子求偏导，\(\theta\) 的部分和上面一样，就不考虑了，直接看 \(\alpha\). 注意到由于 \(\alpha_j = P(w_j)\) 其实决定了类分布 \(p(x \mid w_j, \theta_j)\)，后者根本不含 \(\alpha_j\) 项，因此只是一个系数而已，便有：

\[\begin{align} \frac{\partial}{\partial \alpha_j} LL &= \lambda + \sum_{k=1}^n \frac{1}{p(x_k \mid \theta)} \frac{\partial}{\partial \alpha_j} \sum_{j=1}^c \alpha_j p(x_k \mid w_j, \theta_j) \\ &= \lambda + \sum_{k=1}^n \frac{p(x_k \mid w_j, \theta_j)}{p(x_k \mid \theta)} \\ &= 0 \end{align}\]

这里用到一个小技巧，上面有 c 个导数等于0的等式，将两边同时乘以 \(\alpha_j\) 并分别加起来，就有：

\[\begin{align} 0 &= \lambda \sum_{j=1}^c \alpha_j + \sum_{k=1}^n \sum_{j=1}^c \frac{p(x_k \mid w_j, \theta_j)\alpha_j}{p(x_k \mid \theta)} \\ 0 &= \lambda + \sum_{k=1}^n \sum_{j=1}^c \frac{p(x_k, w_j \mid \theta)}{p(x_k \mid \theta)} \\ 0 &= \lambda + \sum_{k=1}^n \frac{p(x_k \mid \theta)}{p(x_k \mid \theta)} \\ 0 &= \lambda + n \\ \lambda &= -n \end{align}\]

再将 \(\lambda = -n\) 代入上面导数为0项中，并两边乘以 \(\alpha_j\)，便有

\[\begin{align} n &= \sum_{k=1}^n \frac{p(x_k \mid w_j, \theta_j)}{p(x_k \mid \theta)} \\ n \alpha_j &= \sum_{k=1}^n \frac{p(x_k \mid w_j, \theta_j)\alpha_j}{p(x_k \mid \theta)} \\ \alpha_j &= \frac{1}{n} \sum_{k=1}^n P(w_j \mid x_k, \theta) \end{align}\]

所以，结果简单来说就是，先验的最大似然估计等于样本后验的平均值.

可以发现，上面动不动两边要乘以 \(\alpha_j\) 是为了把分子上的似然函数凑出一个联合密度函数来，好和分母或者求和进行处理。

高斯混合密度

上面推导的是一般情况，而假设我们知道了混合密度中每一个都是高斯，那么所求参数就变成了高斯的参数. 回忆一下多元高斯的密度函数为：

\[p(x \mid \mu, \Sigma) = \frac{1}{\sqrt[d]{(2\pi)} \sqrt{\vert\Sigma^{-1}\vert}}\exp(-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu))\]

考虑一下只有 \(\mu\) 未知的情况，上面偏导等于0的项中，需要求密度函数的对数形式。

\[\begin{align} \frac{\partial}{\partial \mu_j} \ln p(x \mid \mu_j, \Sigma_j) &= -\frac{1}{2} \frac{\partial}{\partial \mu_j} (x - \mu_j)^T \Sigma_j^{-1}(x - \mu_j) \\ &= \Sigma^{-1}(x - \mu_j) \end{align}\]

上面涉及到对二次型复合函数求导，具体可以参考 [3] 里的矩阵微分操作记录，这里由于我们需要一个列向量，所以结果反了一下，另外协相关矩阵是对称矩阵，因此最终形式也比较简单。

而如果 \(\Sigma_j\) 也未知的话我就不知道怎么求了，总之周志华老师 [4] 直接给了个结果。
（20160611补充：搜到一个关于估计协相关矩阵的，将二次型视为一个矩阵的 trace 这招比较犀利：Wikipedia）

用这种方法最后实际是一个 NP 问题，因为需要对所有的聚类情形遍历一遍并找出能让似然最大者。实际都是用 EM 来求的。

注意到上面令梯度等于0后，只有后验、密度函数需要求，而后验可以通过贝叶斯公式换为密度函数和先验，因此通过初始化一个多元参数 \(\mu_j, \Sigma_j, j = 1, 2, \dots, c\)，E 步直接求出后验，M 步再迭代求密度函数的各项未知数.

k-means

k-means 目标是最小化所有点到各自聚类中心的距离之和，实际上就是最大化 \(\Sigma_j = I\) 时的似然函数，因此可以视为上面高斯混合密度的一个特例.

同时，据此可以发现 k-means 为什么对并不聚成一坨的数据聚类效果较差，因为 \(\Sigma_j = I\) 显然是个很强的条件.

此外，对一些非线性可分的数据（比如螺旋线）, k-means 或者混合高斯这样的假设就很不适用了，有人对 k-means 加上了 kernel [5]，使得映射到高维空间后线性可分并且可求内积.

2. spectral clustering

基本想法

谱聚类 [6] 的思想是把数据点都看成图的节点，图上的边则为两个点的某种相似度，例如可以取高斯函数 \(d(x_i, x_j) = \exp(\Vert x_i - x_j \Vert^2 / (2\sigma^2))\). 而边可以是全连接、取 k-nn 连接、某个 \(d(x_i, x_j) < \epsilon\) 邻域内连接等等. 一个类之间的相似度应该比较高，不同类间的比较低.

构建 Laplacian Matrix，标准的 (unnormalized) 的定义为 \(L = D - W\)，其中 D 是对角阵，每个元素为每个顶点的度，W 为邻接矩阵。

这个矩阵有一些重要特点：

  • 行和为0
  • 必有一个0特征值和全为1的特征向量
  • 是半正定矩阵
  • 0特征值的重数是图的连通子图的数目

前两个容易证，后两个可参考 [6].

如果这个图已经不连通了，就自然成类了，但是现实数据总会连起来的。对聚成两类的情形来说，相当于求一个最小分割，使得被割掉的几条边的权和最小，即能找到一个最小二分割.

但这种二分割算法往往会找出一个将某远处的点单独为一类，其他所有点为一类，我们希望聚出的类别差那么远，因此可以使用子图的 cardinality(顶点数) 或者 volume(度的和) 来归一化一下.

不同做法

此外，对于拉普拉斯矩阵还有不同的解释

  • 图的最大分割
  • 随机游走
  • 摄动分析 [8], Ng 除了提出了摄动分析的解释之外也做了比较好的归一化工作，可以说他对谱聚类的贡献是突破性的

在标准拉普拉斯矩阵基础上，对矩阵做不同归一化 (normalized) 可以有不同的形式：

  • \(L_{sym} = D^{-1/2}LD^{1/2}\)，相当于用顶点数归一化
  • \(L_{rw} = D^{-1/2}L\)，等价于随机游走的方法，相当于用 volume 归一化

与其他关系

由于谱聚类最后实际都是转化为离散优化的问题，最后与求矩阵的 trace 迹有关，有一个推导了 kernel k-means 和 谱聚类的关系的工作，证明了二者的等价 [9].

这意味着当我们求超大矩阵的特征向量特别困难时可以继续用 k-means 的搜索和改良方法.

3. 其他方法

  • LVQ: 见 [4] ，需要标注数据，通过标注数据来修正，不知道这还算不算聚类
  • DBSCAN: 见 [4], 并没有什么特别之处，只是不停地以密度范围聚类而已，只不过工程实现里需要考虑距离索引的问题，可以用 geohash 或者 kd-tree 等

更多工作

聚类的其他方法也太多了，例如 2014 年的 Science 上密度峰聚类 [10] 等等，其他方法也有一个 2005 年的 survey 可以查看 [11].

Bibliography

[1]R. O. Duda, P. E. Hart, and D. G. Stork, Pattern Classification. John Wiley & Sons, 2012.
[2]从参数估计到线性判别函数. http://libzx.so/chn/learning/2016/03/27/from-estimator-to-linear-discriminant.html 2016.
[3]Randal J. Barnes, Matrix Differentiation. http://www.atmos.washington.edu/~dennis/MatrixCalculus.pdf
[4]周志华, 机器学习. 清华大学出版社, 2016.
[5]B. Schölkopf, A. Smola, and K.-R. Müller, “Nonlinear component analysis as a kernel eigenvalue problem,” Neural computation, vol. 10, no. 5, pp. 1299–1319, 1998.
[6]U. Von Luxburg, “A tutorial on spectral clustering,” Statistics and computing, vol. 17, no. 4, pp. 395–416, 2007.
[7]J. Shi and J. Malik, “Normalized cuts and image segmentation,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 22, no. 8, pp. 888–905, 2000.
[8]A. Y. Ng, M. I. Jordan, Y. Weiss, and others, “On spectral clustering: Analysis and an algorithm,” Advances in neural information processing systems, vol. 2, pp. 849–856, 2002.
[9]I. S. Dhillon, Y. Guan, and B. Kulis, “Kernel k-means: spectral clustering and normalized cuts,” in Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining, 2004, pp. 551–556.
[10]A. Rodriguez and A. Laio, “Clustering by fast search and find of density peaks,” Science, vol. 344, no. 6191, pp. 1492–1496, Jun. 2014.
[11]R. Xu, D. Wunsch, and others, “Survey of clustering algorithms,” Neural Networks, IEEE Transactions on, vol. 16, no. 3, pp. 645–678, 2005.