前言
如果你对这篇文章感兴趣,可以点击「【访客必读 - 指引页】一文囊括主页内所有高质量博客」,查看完整博客分类与对应链接。
k k k-Means 作为一种经典聚类算法,相信大家都比较熟悉,其将簇中所有的点的均值作为簇中心,整个过程采用欧式空间中的距离度量。不同于 k k k-Means, k k k-Medoids 将距簇中所有点距离之和最小的点作为簇中心,如下所示:
medoid ( C ) : = arg min x i ∈ C ∑ x j ∈ C d ( x i , x j ) , \operatorname{medoid}(C):=\underset{x_i \in C}{\arg \min } \sum_{x_j \in C} d\left(x_i, x_j\right), medoid(C):=xi∈Cargminxj∈C∑d(xi,xj),
其中 d ( ⋅ , ⋅ ) d(\cdot, \cdot) d(⋅,⋅) 为采用的度量。整个过程希望最小化:
Loss : = ∑ i = 1 k ∑ x c ∈ C i d ( x c , m i ) , \text{Loss}:=\sum_{i=1}^k\sum_{x_c\in C_i} d(x_c, m_i), Loss:=i=1∑kxc∈Ci∑d(xc,mi),
其中 k k k 表示 k k k 个簇, C i C_i Ci 表示第 i i i 个簇, m i m_i mi 为 C i C_i Ci 的簇中心。接下来介绍一些实现上述目标的算法。
PAM (Partitioning Around Medoids)
在最初版本的 PAM 中,整体流程分为两步:
- 第一步为 BUILD,即贪心选取 k k k 个点作为 Medoids;
- 第二步为 SWAP,需迭代多次,每一次选取一对点 ( m i , x o ) (m_i,x_o) (mi,xo),用 x o x_o xo 将中心点 m i m_i mi 替换掉。
具体来说,在得到预处理的距离矩阵后,第一步一共贪心地执行 k k k 次,每一次选择一个使 Loss \text{Loss} Loss 下降最多的点作为 Medoids,这一步总的复杂度为 O ( n 2 k ) O(n^2k) O(n2k)。
第二步需迭代多次,每一次遍历所有的 ( m i , x o ) (m_i,x_o) (mi,xo) 组合,并计算采用该组合后, Loss \text{Loss} Loss 下降的幅度,选取下降幅度最大的组合作为交换,每一次迭代的复杂度为 O ( k ( n − k ) 2 ) O(k(n-k)^2) O(k(n−k)2)。
上述流程为比较暴力的方式,如果经过合理优化,例如提前为每个点计算离它最近的点和次近的点,则 SWAP 步每一次的迭代复杂度可以降为 O ( n 2 ) O(n^2) O(n2) [FasterPAM]。
CLARA
CLARA (Clustering LARge Applications) 是一种通过采样方式来加速 PAM 的方法,具体如下:
- 从大小为 N N N 的数据集中采样 n n n 次,每次采样出 s s s 个点;
- 对这 s s s 个点使用 PAM 算法,得到 k k k 个 medoids candidates;
- 从所有的 medoids candidates(一共 s ∗ k s*k s∗k 个)中挑出 k k k 个作为最终的 medoids.
最后一步可以采用随机抽取,投票加权,以及对 s ∗ k s*k s∗k 个点再执行一遍 PAM 算法等方式,整体过程的伪代码如下所示:
CLARANS
上述 CLARA 有一个问题,即每次采样出一个子集后,该次采样最终选择的 medoids candidates 就被限制在了这个子集中。那有没有什么方法,使得 medoids 的挑选仍然在所有点中进行,而不是局限在一个固定的子集中。
基于上述想法,CLARANS (Clustering Large Applications based on RANdomized Search) 提出在 PAM 的 SWAP 步骤中加入随机化采样,使得整体复杂度下降,具体如下:
- 随机挑选 k k k 个点作为初始 medoids;
- 随机将 k k k 个点中某一个点换成其它 n − k n-k n−k 个点中任意一个,判断 Loss \text{Loss} Loss 有无下降,若下降则重新执行该步,若持续 m a x n e i g h b o r maxneighbor maxneighbor 次置换 Loss \text{Loss} Loss 均未下降,则认为当前这组 medoids 为局部最优,进入下一步;
- 将当前这组 medoids 记录下来,并重复执行上述两步 n u m l o c a l numlocal numlocal 次,并从得到的 n u m l o c a l numlocal numlocal 局部最优中选一组 Loss \text{Loss} Loss 最小的输出。
上述过程对应下述算法:
其中「an arbitrary node in G n , k G_{n,k} Gn,k」即「从 n n n 个点随机挑出 k k k 个点,并将 k k k 个点的集合视作一个 node」,「random neighbor of a node」即随机将 k k k 个点的集合中某一个点换成其它 n − k n-k n−k 个点中的任意一个,「calculate the cost differential of the two nodes」即计算任意置换一个点后 Loss \text{Loss} Loss 的变化情况(该步复杂度为 O ( n − k ) O(n-k) O(n−k))。
Asymmteric k-Medoids
上述算法均基于对称的距离度量 d ( ⋅ , ⋅ ) d(\cdot,\cdot) d(⋅,⋅),那么当 d ( ⋅ , ⋅ ) d(\cdot,\cdot) d(⋅,⋅) 非对称时( d ( x , y ) ≠ d ( y , x ) d(x,y)\not=d(y,x) d(x,y)=d(y,x)),上述聚类算法应该如何变化?
此处介绍一个在 [ISIS16 - Sadaaki Miyamoto] 中给出的方法,即在簇 G i G_i Gi 中分别设立两个簇中心 v i , w i v_i,w_i vi,wi,其满足:
v i = arg min z j ∈ G i ∑ x k ∈ G i d ( x k , z j ) , w i = arg min y j ∈ G i ∑ x k ∈ G i d ( y j , x k ) . \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} viwi=zj∈Giargminxk∈Gi∑d(xk,zj),=yj∈Giargminxk∈Gi∑d(yj,xk).
随后重新定义 x x x 距簇 G i G_i Gi 的距离: D ( x , G i ) = α d ( x , v i ) + ( 1 − α ) d ( w i , x ) D(x,G_i)=\alpha d(x,v_i)+(1-\alpha)d(w_i,x) D(x,Gi)=αd(x,vi)+(1−α)d(wi,x),其中超参数 α ∈ [ 0 , 1 ] \alpha\in [0,1] α∈[0,1]。
Trimed
此处介绍一种在 N 个点的集合中,快速确定一个 medoid 的方法 Trimed。该方法在 [ICML17 - James Newling] 中提出,其思想为「当采用的度量 dist ( ⋅ , ⋅ ) \text{dist}(\cdot, \cdot) dist(⋅,⋅) 为距离度量时(满足对称性、同一性、三角不等式),可以通过剪枝,忽略一些不可能成为 medoid 的点,进而实现加速」。
对于集合 S = { x ( 1 ) , . . . , x ( N ) } \mathcal{S}=\{x(1),...,x(N)\} S={
x(1),...,x(N)},定义 E ( i ) E(i) E(i) 为:
E ( i ) = 1 N ∑ k = 1 N dist ( x ( i ) , x ( k ) ) . E(i)=\frac{1}{N}\sum_{k=1}^N \text{dist}(x(i),x(k)). E(i)=N1k=1∑Ndist(x(i),x(k)).
我们的目的是找到使 E ( i ) E(i) E(i) 最小的 i i i,最直接的做法就是将所有的 E ( i ) E(i) E(i) 计算出来,得到最小值,其时间复杂度为 O ( N 2 ) O(N^2) O(N2),但明显复杂度过高。
由于我们只需要找到 E ( i ) E(i) E(i) 最小值,而不需要生成对 E ( i ) E(i) E(i) 的完整排序,因此存在可以剪枝的空间。具体来说,Trimed 主要利用了下述不等式进行剪枝:
E ( j ) ≥ ∣ E ( i ) − dist ( x ( i ) , x ( j ) ) ∣ . E(j)\geq |E(i)-\text{dist}(x(i),x(j))|. E(j)≥∣E(i)−dist(x(i),x(j))∣.
证明如下:
E ( j ) = 1 N ∑ k = 1 N dist ( x ( j ) , x ( k ) ) ≥ 1 N ∑ k = 1 N ( dist ( x ( i ) , x ( k ) ) − dist ( x ( i ) , x ( j ) ) ) = E ( i ) − dist ( x ( i ) , x ( j ) ) E ( j ) ≥ 1 N ∑ k = 1 N ( dist ( x ( i ) , x ( j ) ) − dist ( x ( i ) , x ( k ) ) ) = dist ( x ( i ) , x ( j ) ) − E ( i ) \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} E(j)E(j)=N1k=1∑Ndist(x(j),x(k))≥N1k=1∑N(dist(x(i),x(k))−dist(x(i),x(j)))=E(i)−dist(x(i),x(j))≥N1k=1∑N(dist(x(i),x(j))−dist(x(i),x(k)))=dist(x(i),x(j))−E(i)
通过上述剪枝思路,当数据维度较低时,Trimed 可以在一定假设下将 O ( N 2 ) O(N^2) O(N2) 降至 O ( N N ) O(N\sqrt N) O(NN)(时间复杂度随数据维度指数增长,理论证明见原论文),整体算法如下:
BanditPAM
该篇文章的主要目的是「提高 PAM 算法的效率」,具体的思路是「引入 Upper Confidence Bound (UCB) 算法,使用采样估计的方式,来替代精准计算,进而降低时间复杂度」,并基于上述思路提出了 BanditPAM [NIPS20 - Mo Tiwari] 以及对应的代码实现。
如前文所述,PAM 算法分为 BUILD 与 SWAP 两步,其具体数学形式如下:
BUILD: arg min x ∈ X \ M l 1 n ∑ j = 1 n [ ( d ( x , x j ) − min m ′ ∈ M l d ( m ′ , x j ) ) ∧ 0 ] , SWAP: arg min ( m , x ) ∈ M × ( X \ M ) 1 n ∑ j = 1 n [ ( d ( x , x j ) − min m ′ ∈ M \ { m } d ( m ′ , x j ) ) ∧ 0 ] , \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} BUILD:x∈X\Mlargminn1j=1∑n[(d(x,xj)−m′∈Mlmind(m′,xj))∧0],SWAP:(m,x)∈M×(X\M)argminn1j=1∑n[(d(x,xj)−m′∈M\{
m}mind(m′,xj))∧0],
其中 ∧ \wedge ∧ 表示 min \min min, M l \mathcal{M}_l Ml 表示当前已选择的 l l l 个点的集合,BUILD 步对应再选取一个点 x x x 加入 M l \mathcal{M}_l Ml 的数学准则,SWAP 步对应再选取一对点 ( m , x ) (m,x) (m,x),并将 x x x 与 m m m 进行交换的数学准则。
基于上述数学形式,该篇文章为其形式化了一个统一的形式:
Shared Problem: arg min x ∈ S t a r 1 ∣ S r e f ∣ ∑ x j ∈ S r e f g x ( x j ) , \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), Shared Problem: x∈Starargmin∣Sref∣1xj∈Sref∑gx(xj),
其与 BUILD 和 SWAP 的具体对应关系如下:
BUILD: S t a r = X \ M l , S r e f = X , g x ( x j ) = ( d ( x , x j ) − min m ′ ∈ M l d ( m ′ , x j ) ) ∧ 0 , SWAP: S t a r = M × ( X \ M ) , S ref = X , g x ( x j ) = ( d ( x , x j ) − min m ′ ∈ M \ { m } d ( m ′ , x j ) ) ∧ 0. \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} BUILD:Star=X\Ml,Sref=X,gx(xj)=(d(x,xj)−m′∈Mlmind(m′,xj))∧0,SWAP:Star=M×(X\M),Sref =X,gx(xj)=(d(x,xj)−m′∈M\{
m}mind(m′,xj))∧0.
得到统一形式之后,我们可以发现精确求解 arg min x ∈ S t a r \arg \min_{x \in \mathcal{S}_{\mathrm{tar}}} argminx∈Star 的时间复杂度为 O ( ∣ S t a r ∣ ∣ S r e f ∣ ) O(|\mathcal{S}_{\mathrm{tar}}||\mathcal{S}_{\mathrm{ref}}|) O(∣Star∣∣Sref∣),为降低其时间复杂度,一个直接的思路就是对 ∣ S r e f ∣ |\mathcal{S}_{\mathrm{ref}}| ∣Sref∣ 下手,即通过多次采样,每次使用一个子集 S r e f _ b a t c h \mathcal{S}_{\mathrm{ref\_batch}} Sref_batch 来替代 S r e f \mathcal{S}_{\mathrm{ref}} Sref。
确定好这个思路后,选择 Upper Confidence Bound 算法来平衡「探索」与「利用」,在此处「探索」对应「继续采样使估计更准」,「利用」对应「根据现有估计结果,删除一部分选项」。
具体来说,每轮采样结束后,每个点 x ∈ S t a r x \in \mathcal{S}_{\mathrm{tar}} x∈Star 对应一个估计的均值 μ ^ x \hat{\mu}_x μ^x 与置信区间 C x C_x Cx,即代表 x x x 真实值在区间 [ μ ^ x − C x , μ ^ x + C x ] [\hat{\mu}_x-C_x, \hat{\mu}_x+C_x] [μ^x−Cx,μ^x+Cx] 中,因此我们可以将「区间下界」大于「最小区间上界」的点丢弃,即丢弃满足 μ ^ x − C x ≤ min y μ ^ y + C y \hat{\mu}_x-C_x\leq \min_y \hat{\mu}_y+C_y μ^x−Cx≤minyμ^y+Cy 的点,因为这些点的真实值大概率不是最小值。
上述思路对应的具体算法如下:
其中 σ x \sigma_x σx 由标准差得到,即 σ x = STD y ∈ S ref batch g x ( y ) \sigma_x=\operatorname{STD}_{y \in \mathcal{S}_{\text {ref batch }}} g_x(y) σx=STDy∈Sref batch gx(y)。
最后,该篇文章通过引入假设「for a fixed target point x x x and a randomly sampled reference point x J x_J xJ, the random variable Y = g x ( x J ) Y = g_x(x_J) Y=gx(xJ) is σ x \sigma_x σx-sub-Gaussian for some known parameter σ x \sigma_x σx」,证明上述算法在 O ( n log n ) O(n\log n) O(nlogn) 的时间复杂度下可以高概率得到与 PAM 精确计算一样的结果:
Parallel k-Medoids
此处提供几篇通过使用并行化、分布式来加速 k-Medoids 的文章:
- [KDD17 - Hwanjun Song] PAMAE: Parallel k-Medoids Clustering with High Accuracy and Efficiency
- [INS21 - Anton V. Ushakov] Near-optimal large-scale k-medoids clustering
参考资料
- k-medoids - Wikipedia
- [arXiv21 - Erich Schubert] Fast and Eager k-Medoids Clustering: O(k) Runtime Improvement of the PAM, CLARA, and CLARANS Algorithms
- [Book - Leonard Kaufman] Finding Groups in Data An Introduction to Cluster Analysis
- Advanced Partitional clustering: medoids, PAM and CLARA and lite versions
- [TKDE02 - Raymond T. Ng] CLARANS: A Method for Clustering Objects for Spatial Data Mining
- [ISIS16 - Sadaaki Miyamoto] Hierarchical and Non-hierarchical Medoid Clustering Using Asymmetric Similarity Measures
- [ICML17 - James Newling] A Sub-Quadratic Exact Medoid Algorithm
- [NIPS20 - Mo Tiwari] BanditPAM: Almost Linear Time k-Medoids Clustering via Multi-Armed Bandits