k-Medoids 聚类系列算法：PAM, CLARA, CLARANS, Trimed, BanditPAM

$k$-Means 作为一种经典聚类算法，相信大家都比较熟悉，其将簇中所有的点的均值作为簇中心，整个过程采用欧式空间中的距离度量。不同于 $k$-Means，$k$-Medoids 将距簇中所有点距离之和最小的点作为簇中心，如下所示：

$$ \operatorname{medoid}(C):=\underset{x_i \in C}{\arg \min } \sum_{x_j \in C} d\left(x_i, x_j\right), $$

其中 $d(\cdot, \cdot)$ 为采用的度量。整个过程希望最小化：

$$ \text{Loss}:=\sum_{i=1}^k\sum_{x_c\in C_i} d(x_c, m_i), $$

其中 $k$ 表示 $k$ 个簇，$C_i$ 表示第 $i$ 个簇，$m_i$ 为 $C_i$ 的簇中心。接下来介绍一些实现上述目标的算法。

PAM (Partitioning Around Medoids)

在最初版本的 PAM 中，整体流程分为两步：

第一步为 BUILD，即贪心选取 $k$ 个点作为 Medoids；
第二步为 SWAP，需迭代多次，每一次选取一对点 $(m_i,x_o)$，用 $x_o$ 将中心点 $m_i$ 替换掉。

具体来说，在得到预处理的距离矩阵后，第一步一共贪心地执行 $k$ 次，每一次选择一个使 $\text{Loss}$ 下降最多的点作为 Medoids，这一步总的复杂度为 $O(n^2k)$。

第二步需迭代多次，每一次遍历所有的 $(m_i,x_o)$ 组合，并计算采用该组合后，$\text{Loss}$ 下降的幅度，选取下降幅度最大的组合作为交换，每一次迭代的复杂度为 $O(k(n-k)^2)$。

上述流程为比较暴力的方式，如果经过合理优化，例如提前为每个点计算离它最近的点和次近的点，则 SWAP 步每一次的迭代复杂度可以降为 $O(n^2)$ [FasterPAM]。

CLARA

CLARA (Clustering LARge Applications) 是一种通过采样方式来加速 PAM 的方法，具体如下：

从大小为 $N$ 的数据集中采样 $n$ 次，每次采样出 $s$ 个点；
对这 $s$ 个点使用 PAM 算法，得到 $k$ 个 medoids candidates；
从所有的 medoids candidates（一共 $s*k$ 个）中挑出 $k$ 个作为最终的 medoids.

最后一步可以采用随机抽取，投票加权，以及对 $s*k$ 个点再执行一遍 PAM 算法等方式，整体过程的伪代码如下所示：

CLARANS

上述 CLARA 有一个问题，即每次采样出一个子集后，该次采样最终选择的 medoids candidates 就被限制在了这个子集中。那有没有什么方法，使得 medoids 的挑选仍然在所有点中进行，而不是局限在一个固定的子集中。

基于上述想法，CLARANS (Clustering Large Applications based on RANdomized Search) 提出在 PAM 的 SWAP 步骤中加入随机化采样，使得整体复杂度下降，具体如下：

随机挑选 $k$ 个点作为初始 medoids；
随机将 $k$ 个点中某一个点换成其它 $n-k$ 个点中任意一个，判断 $\text{Loss}$ 有无下降，若下降则重新执行该步，若持续 $maxneighbor$ 次置换 $\text{Loss}$ 均未下降，则认为当前这组 medoids 为局部最优，进入下一步；
将当前这组 medoids 记录下来，并重复执行上述两步 $numlocal$ 次，并从得到的 $numlocal$ 局部最优中选一组 $\text{Loss}$ 最小的输出。

上述过程对应下述算法：

其中「an arbitrary node in $G_{n,k}$」即「从 $n$ 个点随机挑出 $k$ 个点，并将 $k$ 个点的集合视作一个 node」，「random neighbor of a node」即随机将 $k$ 个点的集合中某一个点换成其它 $n-k$ 个点中的任意一个，「calculate the cost differential of the two nodes」即计算任意置换一个点后 $\text{Loss}$ 的变化情况（该步复杂度为 $O(n-k)$）。

Asymmteric k-Medoids

上述算法均基于对称的距离度量 $d(\cdot,\cdot)$，那么当 $d(\cdot,\cdot)$ 非对称时（$d(x,y)\not=d(y,x)$），上述聚类算法应该如何变化？

此处介绍一个在 [ISIS16 - Sadaaki Miyamoto] 中给出的方法，即在簇 $G_i$ 中分别设立两个簇中心 $v_i,w_i$，其满足：

$$ \begin{aligned} v_i & =\underset{z_j\in G_i}{\arg \min} \sum_{x_k \in G_i} d\left(x_k, z_j\right), \\ w_i & =\underset{y_j\in G_i}{\arg \min} \sum_{x_k \in G_i} d\left(y_j, x_k\right) . \end{aligned} $$

随后重新定义 $x$ 距簇 $G_i$ 的距离：$D(x,G_i)=\alpha d(x,v_i)+(1-\alpha)d(w_i,x)$，其中超参数 $\alpha\in [0,1]$。

Trimed

此处介绍一种在 N 个点的集合中，快速确定一个 medoid 的方法 Trimed。该方法在 [ICML17 - James Newling] 中提出，其思想为「当采用的度量 $\text{dist}(\cdot, \cdot)$ 为距离度量时（满足对称性、同一性、三角不等式），可以通过剪枝，忽略一些不可能成为 medoid 的点，进而实现加速」。

对于集合 $\mathcal{S}={x(1),…,x(N)}$，定义 $E(i)$ 为：

$$ E(i)=\frac{1}{N}\sum_{k=1}^N \text{dist}(x(i),x(k)). $$

我们的目的是找到使 $E(i)$ 最小的 $i$，最直接的做法就是将所有的 $E(i)$ 计算出来，得到最小值，其时间复杂度为 $O(N^2)$，但明显复杂度过高。

由于我们只需要找到 $E(i)$ 最小值，而不需要生成对 $E(i)$ 的完整排序，因此存在可以剪枝的空间。具体来说，Trimed 主要利用了下述不等式进行剪枝：

$$ E(j)\geq |E(i)-\text{dist}(x(i),x(j))|. $$

证明如下：

$$ \begin{aligned} E(j)&= \frac{1}{N}\sum_{k=1}^N \text{dist}(x(j),x(k)) \\ &\geq \frac{1}{N}\sum_{k=1}^N (\text{dist}(x(i),x(k))-\text{dist}(x(i),x(j))) \\ &=E(i)-\text{dist}(x(i),x(j)) \\ \\ E(j)&\geq \frac{1}{N}\sum_{k=1}^N (\text{dist}(x(i),x(j))-\text{dist}(x(i),x(k))) \\ &=\text{dist}(x(i),x(j))-E(i) \end{aligned} $$

通过上述剪枝思路，当数据维度较低时，Trimed 可以在一定假设下将 $O(N^2)$ 降至 $O(N\sqrt N)$（时间复杂度随数据维度指数增长，理论证明见原论文），整体算法如下：

BanditPAM

该篇文章的主要目的是「提高 PAM 算法的效率」，具体的思路是「引入 Upper Confidence Bound (UCB) 算法，使用采样估计的方式，来替代精准计算，进而降低时间复杂度」，并基于上述思路提出了 BanditPAM [NIPS20 - Mo Tiwari] 以及对应的代码实现。

如前文所述，PAM 算法分为 BUILD 与 SWAP 两步，其具体数学形式如下：

$$ \begin{gathered} \text{BUILD:} \quad \underset{x \in \mathcal{X} \backslash \mathcal{M}_l}{\arg \min } \frac{1}{n} \sum_{j=1}^n\left[\left(d\left(x, x_j\right)-\min _{m^{\prime} \in \mathcal{M}_l} d\left(m^{\prime}, x_j\right)\right) \wedge 0\right], \\ \text{SWAP:}\quad \underset{(m, x) \in \mathcal{M} \times(\mathcal{X} \backslash \mathcal{M})}{\arg \min } \frac{1}{n} \sum_{j=1}^n\left[\left(d\left(x, x_j\right)-\min _{m^{\prime} \in \mathcal{M} \backslash\{m\}} d\left(m^{\prime}, x_j\right)\right) \wedge 0\right], \end{gathered} $$

其中 $\wedge$ 表示 $\min$，$\mathcal{M}_l$ 表示当前已选择的 $l$ 个点的集合，BUILD 步对应再选取一个点 $x$ 加入 $\mathcal{M}_l$ 的数学准则，SWAP 步对应再选取一对点 $(m,x)$，并将 $x$ 与 $m$ 进行交换的数学准则。

基于上述数学形式，该篇文章为其形式化了一个统一的形式：

$$ \text { Shared Problem: } \quad \underset{x \in \mathcal{S}_{\mathrm{tar}}}{\arg \min } \frac{1}{\left|\mathcal{S}_{\mathrm{ref}}\right|} \sum_{x_j \in \mathcal{S}_{\mathrm{ref}}} g_x\left(x_j\right), $$

其与 BUILD 和 SWAP 的具体对应关系如下：

$$ \begin{gathered} \text{BUILD:} \quad \mathcal{S}_{\mathrm{tar}}=\mathcal{X} \backslash \mathcal{M}_l, \mathcal{S}_{\mathrm{ref}}=\mathcal{X}, g_x\left(x_j\right)=\left(d\left(x, x_j\right)-\min _{m^{\prime} \in \mathcal{M}_l} d\left(m^{\prime}, x_j\right)\right) \wedge 0, \\ \text{SWAP:} \quad \mathcal{S}_{\mathrm{tar}}=\mathcal{M} \times(\mathcal{X} \backslash \mathcal{M}), \mathcal{S}_{\text {ref }}=\mathcal{X}, g_x\left(x_j\right)=\left(d\left(x, x_j\right)-\min _{m^{\prime} \in \mathcal{M} \backslash\{m\}} d\left(m^{\prime}, x_j\right)\right) \wedge 0. \end{gathered} $$

得到统一形式之后，我们可以发现精确求解 $\arg \min_{x \in \mathcal{S}_{\mathrm{tar}}}$ 的时间复杂度为 $O(|\mathcal{S}_{\mathrm{tar}}||\mathcal{S}_{\mathrm{ref}}|)$，为降低其时间复杂度，一个直接的思路就是对 $|\mathcal{S}_{\mathrm{ref}}|$ 下手，即通过多次采样，每次使用一个子集 $\mathcal{S}_{\mathrm{ref\_batch}}$ 来替代 $\mathcal{S}_{\mathrm{ref}}$。

确定好这个思路后，选择 Upper Confidence Bound 算法来平衡「探索」与「利用」，在此处「探索」对应「继续采样使估计更准」，「利用」对应「根据现有估计结果，删除一部分选项」。

具体来说，每轮采样结束后，每个点 $x \in \mathcal{S}_{\mathrm{tar}}$ 对应一个估计的均值 $\hat{\mu}_x$ 与置信区间 $C_x$，即代表 $x$ 真实值在区间 $[\hat{\mu}_x-C_x, \hat{\mu}_x+C_x]$ 中，因此我们可以将「区间下界」大于「最小区间上界」的点丢弃，即丢弃满足 $\hat{\mu}_x-C_x\leq \min_y \hat{\mu}_y+C_y$ 的点，因为这些点的真实值大概率不是最小值。

上述思路对应的具体算法如下：

其中 $\sigma_x$ 由标准差得到，即 $\sigma_x=\operatorname{STD}{y \in \mathcal{S}{\text {ref batch }}} g_x(y)$。

最后，该篇文章通过引入假设「for a fixed target point $x$ and a randomly sampled reference point $x_J$, the random variable $Y = g_x(x_J)$ is $\sigma_x$-sub-Gaussian for some known parameter $\sigma_x$」，证明上述算法在 $O(n\log n)$ 的时间复杂度下可以高概率得到与 PAM 精确计算一样的结果：

Parallel k-Medoids

此处提供几篇通过使用并行化、分布式来加速 k-Medoids 的文章：