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]。
对于集合
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)(时间复杂度随数据维度指数增长,理论证明见原论文),整体算法如下:
如前文所述,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。
具体来说,每轮采样结束后,每个点
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