本記事では、2024年2月に arXiv へ投稿された次の論文の内容を概説する。

He, Yudong. "Imbalanced Data Clustering using Equilibrium K-Means." arXiv preprint arXiv:2402.14490 (2024).

この論文では、勾配法の観点から既存のクラスタリングアルゴリズム　K-Means, Fuzzy C-Means を一般化し、新しいクラスタリングアルゴリズム Equilibrium K-Means (EKM) を提案している。提案した EKM は、下図のような不均衡データに有効である。

クラスタリングの定式化と既存アルゴリズム

K-Means

K-means クラスタリングでは、$N$ 個のデータ ${\mathbf x}$ を $K$ 個のクラスタ $\mathbb S_1,\ldots,\mathbb S_K$ に分割する、以下の最適化問題を解く。

\newcommand{\minimize}{\mathop{\rm minimize}\limits}
\newcommand{\subjectto}{\mathop{\rm subject\ \rm to}\limits}
\begin{align}
\minimize_{\mathbb S_1,\ldots,\mathbb S_K}&\quad\sum_{k=1}^{K}\sum_{\mathbf x\in\mathbb S_k}\|\mathbf x-\mathbf c_k\|_2^2 \\
\subjectto&\quad \mathbf c_k = \frac{1}{|\mathbb S_k|}\sum_{\mathbf x \in \mathbb S_k}\mathbf x
\end{align}

この問題は NP 困難であることが知られており¹、以下のヒューリスティックなアルゴリズムが提案されている。

Lloyd's algorithm

Centroids $\mathbf c_1^{(1)}, \ldots, \mathbf c_K^{(1)}$ を初期化

以下を反復

$K$ 個のクラスタを更新する

\mathbb S_k^{(\tau)} = \left\{\mathbf x \mid \|\mathbf x-\mathbf c_k^{(\tau)}\|_2 \leq \|\mathbf x-\mathbf c_i^{(\tau)}\|_2, \forall i \in \{1,\ldots,K\}\right\}

Centroids を更新する

\mathbf c_k^{(\tau+1)} = \frac{1}{\left|\mathbb S_k^{(\tau)}\right|}\sum_{\mathbf x \in \mathbb S_k^{(\tau)}}\mathbf x

このアルゴリズムは、Python では、例えば scikit-learn の k-means のデフォルトアルゴリズムに指定されている。(algorithm 参照)

Fuzzy C-Means

K-means では各データを一つのクラスタに割り当てていたが、fuzzy c-means は、クラスタ $k$ ごとに各データ $\mathbf x_n$ の帰属度 $u_{kn}$ を計算し、クラスタ中心 $\mathbf c_k$ をデータの重み付き平均により更新する手法であり、以下の最小化問題を扱う。

\begin{align}
\minimize_{\mathbf c_1, \ldots, \mathbf c_K} &\quad \sum_{n=1}^N \sum_{k=1}^K \left(u_{kn}\right)^m \|\mathbf x_n - \mathbf c_k\|_2^2 \\
\subjectto &\quad u_{kn} \in [0, 1],\ \sum_{k=1}^K u_{kn} = 1,\ 0 < \sum_{n=1}^N u_{kn} < N
\end{align}

K-means と同様に、帰属度 $u_{kn}$ とクラスタ中心 $\mathbf c_k$ を交互に更新するアルゴリズムが知られている。

Centroids $\mathbf c_1^{(1)}, \ldots, \mathbf c_K^{(1)}$ を初期化

以下を反復

帰属度を更新する

u_{kn}^{(\tau)} = \left(\sum_{i=1}^K \left(\frac{\|\mathbf x_n - \mathbf c_k\|_2}{\|\mathbf x_n - \mathbf c_i\|_2}\right)^{\frac{2}{m-1}}\right)^{-1}

Centroids を更新する

\mathbf c_k^{(\tau+1)} = \frac{\sum_{n=1}^N\left(u_{kn}^{(\tau)}\right)^m\mathbf x_n}{\sum_{n=1}^N\left(u_{kn}^{(\tau)}\right)^m}

日本語の記事だと、この辺とか参考になるかも。

ファジィc-means法 - 機械学習の「朱鷺の杜Wiki」

Python では、scikit-fuzzy なるライブラリで実行できるらしい。

勾配法による解釈

データ $\mathbf x_n$ と centroid $\mathbf c_k$ のユークリッド距離を $d_{kn}=\frac12 |\mathbf x_n - \mathbf c_k|_2^2$ とする。
また、K-means の最適化問題を $\min$ 関数を用いた以下の関数の最小化問題として定式化する (within-cluster sum of squares, WCSS)。

\mathrm{WCSS} = \sum_{n=1}^N \min(d_{1n},\ldots,d_{Kn})

ここで、$\min$ 関数をなめらかに近似した関数 $h$ を用いた以下の最適化問題を考えることで、クラスタリング問題を勾配法から見る。

\newcommand{\argmin}{\mathop{\rm arg\, \rm min}\limits}
\argmin_{\mathbf c_1,\ldots, \mathbf c_K} \sum_{n=1}^N h(d_{1n},\ldots,d_{Kn})

Smooth Minimum Functions

微分可能な単調増加関数 $f: [0, +\infty) \to [0, +\infty)$ が $\lim_{x\to+\infty}\frac{x}{f(x)}=0$ (equivalently $\lim_{x\to+\infty}\frac{1}{f'(x)}=0$) を満たすとする。$g(x)=1/f(x)$ とすることで、smooth minimum function $h_1: [0, +\infty)^K \to [0, +\infty)$ を以下で定義できる。

h_1(x_1,\ldots,x_K) = g^{-1}\left(g(x_1) + \cdots + g(x_K)\right)

例えば、$f(x)=\exp(\lambda x)$ とすると、$g(x)=\exp(-\lambda x)$, $g^{-1}(x) = -\frac1\lambda\log(x)$ であるので、

\begin{align}
h_1(x_1,\ldots,x_K) &= \mathrm{LSE}_\lambda(x_1,\ldots,x_K) \\
&= -\frac1\lambda \log\left(\sum_{k=1}^K\exp(-\lambda x_k)\right)
\end{align}

のように、LogSumExp 関数となる。$x_{\min}=\min(x_1,\ldots,x_K)$ とすると、$\lambda \to +\infty$ のとき、

\begin{align}
\mathrm{LSE}_\lambda(x_1,\ldots,x_K) &= -\frac1\lambda \log\left(\exp(-\lambda x_\min) \sum_{k=1}^K \exp\left(-\lambda(x_k - x_\min)\right)\right) \\
&= x_\min - \frac1\lambda \log\left(\sum_{k=1}^K \exp\left(-\lambda(x_k - x_\min)\right)\right) \\
&\to x_\min
\end{align}

となり、$\min$ に収束する。また、$f(x)=x^p\ (p>1)$ とすると、$g(x)=x^{-p}$, $g^{-1}(x)=x^{-\frac1p}$ であるので、

\begin{align}
h_1(x_1,\ldots,x_K) &= \mathrm{PN}_p(x_1,\ldots,x_K) \\
&= \left(\sum_{k=1}^K x_k^{-p}\right)^{-\frac1p}
\end{align}

のように、p-Norm 関数となる。$p \to +\infty$ のとき、$PN_p\to\min$ が確認できる。

\begin{align}
\mathrm{PN}_p(x_1,\ldots,x_K) &= \left(x_\min^{-p}\sum_{k=1}^K \left(\frac{x_k}{x_\min}\right)^{-p}\right)^{-\frac1p} \\
&= x_\min\left(\sum_{k=1}^K \left(\frac{x_k}{x_\min}\right)^{-p}\right)^{-\frac1p} \\
&\to x_\min
\end{align}

他にも、

h_2(x_1,\ldots,x_K) = \frac{x_1g(x_1)+\cdots+x_Kg(x_K)}{g(x_1)+\cdots+g(x_K)}

とすることでも smooth minimum function を構成できる。このとき、$f(x)=\exp(\alpha x)$ とすると、$g(x)=\exp(-\alpha x)$ であるので、

\begin{align}
h_2(x_1,\ldots,x_K) &= \mathrm{boltz}_\alpha(x_1,\ldots,x_K) \\
&= \frac{\sum_{k=1}^K x_k\exp(-\alpha x_k)}{\sum_{k=1}^K \exp(-\alpha x_k)}
\end{align}

のように、Boltzmann operator となる。$\alpha \to +\infty$ のとき、

\begin{align}
\mathrm{boltz}_\alpha(x_1,\ldots,x_K) &= \frac{\sum_{k=1}^K x_k\exp(-\alpha (x_k-x_\min))}{\sum_{k=1}^K \exp(-\alpha (x_k-x_\min))} \\
&\to x_\min
\end{align}

のとおり、$x_\min$ に収束する。

以下に $\min(x, 10)$, $\mathrm{LSE}_1(x, 10)$, $\mathrm{PN}_2(x, 10)$, $\mathrm{boltz}_1(x, 10)$ の例を示す。

Smoothed WCSS の最小化

WCSS 式中の $\min$ 関数を $\mathrm{LSE}$, $\mathrm{PN}$, $\mathrm{boltz}$ でそれぞれ近似した smoothed WCSS を勾配降下法で最小化することを考える。

LogSumExp
以下の最小化問題を考える。

\min_{\mathbf c_1,\ldots,\mathbf c_K} \quad J_{\mathrm{LSE}}(\mathbf c_1,\ldots,\mathbf c_K) = \sum_{n=1}^N-\frac1\lambda \log\left(\sum_{k=1}^K\exp(-\lambda d_{kn})\right)

この勾配は、

\begin{align}
\frac{\partial J_{\mathrm{LSE}}}{\partial \mathbf c_k} &= \sum_{n=1}^N \frac{\exp(-\lambda d_{kn})}{\sum_{i=1}^K\exp(-\lambda d_{in})}\frac{\partial d_{kn}}{\partial \mathbf c_k} \\
&= -\sum_{n=1}^N \frac{\exp(-\lambda d_{kn})}{\sum_{i=1}^K\exp(-\lambda d_{in})}(\mathbf x_n - \mathbf c_k)
\end{align}

となるので、勾配降下法による更新式は、

\mathbf c_k^{(\tau+1)} = \mathbf c_k^{(\tau)} - \gamma_k^{(\tau)}\frac{\partial J_{\mathrm{LSE}}}{\partial \mathbf c_k}\left(\mathbf c_k^{(\tau)}\right)

となる。ここで、$\mathbf c_{k}^{(\tau+1)}$ の更新式に $\mathbf c_k^{(\tau)}$ が現れないようにステップ幅 $\gamma_k^{(\tau)}$ を調整すると、

\gamma_k^{(\tau)}=\left(\sum_{n=1}^N\left(\frac{\exp\left(-\lambda d_{kn}^{(\tau)}\right)}{\sum_{i=1}^K\exp\left(-\lambda d_{in}^{(\tau)}\right)}\right)\right)^{-1}

である。$u_{kn}^{(\tau)}=\frac{\exp\left(-\lambda d_{kn}^{(\tau)}\right)}{\sum_{i=1}^K\exp\left(-\lambda d_{in}^{(\tau)}\right)}$ とすると、

\mathbf c_k^{(\tau+1)} = \frac{\sum_{n=1}^N u_{kn}^{(\tau)} \mathbf x_n}{\sum_{n=1}^N u_{kn}^{(\tau)}}

となる。さらに、$\lambda\to+\infty$ とすると、$d_{kn}\leq d_{in}\ (\forall i)$ のとき $u_{kn}^{(\tau)}\to1$, それ以外では $u_{kn}^{(\tau)}\to0$ となるので、更新式は K-means と一致する。

p-Norm
以下の最小化問題を考える。

\min_{\mathbf c_1,\ldots,\mathbf c_K} \quad J_{\mathrm{PN}}(\mathbf c_1,\ldots,\mathbf c_K) = \sum_{n=1}^N\left(\sum_{k=1}^K d_{kn}^{-p}\right)^{-\frac1p}

この勾配は、

\frac{\partial J_{\mathrm{PN}}}{\partial \mathbf c_k} = -\sum_{n=1}^N \frac{d_{kn}^{-p-1}}{\left(\sum_{i=1}^K d_{in}^{-p}\right)^{\frac{1+p}{p}}}(\mathbf x_n - \mathbf c_k)

となるので、勾配降下法による更新式は、

\mathbf c_k^{(\tau+1)} = \mathbf c_k^{(\tau)} - \gamma_k^{(\tau)}\frac{\partial J_{\mathrm{PN}}}{\partial \mathbf c_k}\left(\mathbf c_k^{(\tau)}\right)

となる。ここで、$\mathbf c_{k}^{(\tau+1)}$ の更新式に $\mathbf c_k^{(\tau)}$ が現れないようにステップ幅 $\gamma_k^{(\tau)}$ を調整すると、

\begin{align}
\gamma_k^{(\tau)} &= \left(\sum_{n=1}^N \frac{\left(d_{kn}^{(\tau)}\right)^{-p-1}}{\left(\sum_{i=1}^K \left(d_{in}^{(\tau)}\right)^{-p}\right)^{\frac{1+p}{p}}}\right)^{-1} \\
&= \left(\sum_{n=1}^N \frac{1}{\left(\sum_{i=1}^K \left(d_{kn}^{(\tau)}\right)^p\left(d_{in}^{(\tau)}\right)^{-p}\right)^{\frac{1+p}{p}}}\right)^{-1}
\end{align}

である。$u_{kn}^{(\tau)}=\left(\sum_{i=1}^K \left(d_{kn}^{(\tau)}\right)^p\left(d_{in}^{(\tau)}\right)^{-p}\right)^{-1}$ とし、$p=1/(m-1)$ とすると、

\mathbf c_k^{(\tau+1)} = \frac{\sum_{n=1}^N \left(u_{kn}^{(\tau)}\right)^m \mathbf x_n}{\sum_{n=1}^N \left(u_{kn}^{(\tau)}\right)^m}

となり、fuzzy c-means の更新式と一致する。

Boltzmann Operator
以下の最小化問題を考える。

\min_{\mathbf c_1,\ldots,\mathbf c_K} \quad J_{\mathrm{B}}(\mathbf c_1,\ldots,\mathbf c_K) = \sum_{n=1}^N\frac{\sum_{i=1}^K d_{in}\exp\left(-\alpha d_{in}\right)}{\sum_{i=1}^K \exp\left(-\alpha d_{in}\right)}

この勾配は、

\frac{\partial J_{\mathrm{PN}}}{\partial \mathbf c_k} = -\sum_{n=1}^N \frac{\exp\left(-\alpha d_{kn}\right)}{\sum_{i=1}^K\exp\left(-\alpha d_{in}\right)} \left[1 - \alpha\left(d_{kn} - \frac{\sum_{i=1}^Kd_{in}\exp\left(-\alpha d_{in}\right)}{\sum_{i=1}^K\exp\left(-\alpha d_{in}\right)}\right)\right](\mathbf x_n - \mathbf c_k)

となるので、勾配降下法による更新式は、

\mathbf c_k^{(\tau+1)} = \mathbf c_k^{(\tau)} - \gamma_k^{(\tau)}\frac{\partial J_{\mathrm{B}}}{\partial \mathbf c_k}\left(\mathbf c_k^{(\tau)}\right)

となる。ここで、$\mathbf c_{k}^{(\tau+1)}$ の更新式に $\mathbf c_k^{(\tau)}$ が現れないようにステップ幅 $\gamma_k^{(\tau)}$ を調整すると、

\gamma_k^{(\tau)} = \left(\sum_{n=1}^N \frac{\exp\left(-\alpha d_{kn}^{(\tau)}\right)}{\sum_{i=1}^K\exp\left(-\alpha d_{in}^{(\tau)}\right)} \left[1 - \alpha\left(d_{kn}^{(\tau)} - \frac{\sum_{i=1}^Kd_{in}^{(\tau)}\exp\left(-\alpha d_{in}^{(\tau)}\right)}{\sum_{i=1}^K\exp\left(-\alpha d_{in}^{(\tau)}\right)}\right)\right]\right)^{-1}

である。$w_{kn}^{(\tau)}=\frac{\exp\left(-\alpha d_{kn}^{(\tau)}\right)}{\sum_{i=1}^K\exp\left(-\alpha d_{in}^{(\tau)}\right)} \left[1 - \alpha\left(d_{kn}^{(\tau)} - \frac{\sum_{i=1}^Kd_{in}^{(\tau)}\exp\left(-\alpha d_{in}^{(\tau)}\right)}{\sum_{i=1}^K\exp\left(-\alpha d_{in}^{(\tau)}\right)}\right)\right]$ とすると、

\mathbf c_k^{(\tau+1)} = \frac{\sum_{n=1}^N w_{kn}^{(\tau)} \mathbf x_n}{\sum_{n=1}^N w_{kn}^{(\tau)}}

となる。この更新によるクラスタリングを equilibrium K-means (EKM) と呼ぶ。

数値実験

論文中の実験を K-means, fuzzy c-means, equilibrium K-means の3手法で再現実験する。

ライブラリのインポート

以下を使用する。

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial.distance import cdist
from sklearn.cluster import kmeans_plusplus

データ点の生成

正規分布を用いてランダムデータを生成する。可視化のため次元数は2とし、クラスタ数は3、各クラスタのデータ数は2000, 50, 50と不均衡になるように生成する。また、centroids の初期点は K-means++ アルゴリズムで生成する (sklearn.cluster.kmeans_plusplus を使用)。

seed = 777
np.random.seed(seed)

m1 = np.array([-2.0, 2.0], dtype=np.float32)
m2 = np.array([2.0, -2.0], dtype=np.float32)
m3 = np.array([4.0, 4.0], dtype=np.float32)
means = np.vstack((m1, m2, m3))
n1 = 2000
n2 = 50
n3 = 50
dim = 2

x1 = np.random.normal(m1, size=(n1, dim))
x2 = np.random.normal(m2, size=(n2, dim))
x3 = np.random.normal(m3, size=(n3, dim))
X = np.vstack((x1, x2, x3))
centers = kmeans_plusplus(X, 3, random_state=seed)[0]

colors = matplotlib.colormaps['tab10']
plt.scatter(x1[:, 0], x1[:, 1], alpha=0.5, color=colors(0))
plt.scatter(x2[:, 0], x2[:, 1], alpha=0.5, color=colors(2))
plt.scatter(x3[:, 0], x3[:, 1], alpha=0.5, color=colors(1))
plt.scatter(means[:, 0], means[:, 1],
            marker='+', color='k', label='True center')
plt.scatter(centers[:, 0], centers[:, 1],
            marker='x', color='k', label='Initial centroid')
plt.legend()
plt.show()

アルゴリズムの実装

勾配の計算と各反復の計算の関数を適当に実装する。

def lse_grad(X: np.ndarray, centers: np.ndarray,
             lambd: float = 1.0) -> np.ndarray:
    distances = 0.5 * cdist(centers, X, metric='euclidean') ** 2
    dmin = distances.min(axis=0, keepdims=True)
    expd = np.exp(-lambd * (distances - dmin))
    grads = expd / expd.sum(axis=0, keepdims=True)
    return grads


def pnorm_grad(X: np.ndarray, centers: np.ndarray,
               p: float = 1.0, eps: float = 1e-8) -> np.ndarray:
    distances = 0.5 * cdist(centers, X, metric='euclidean') ** 2
    n_clusters = centers.shape[0]
    m = (1.0 + p) / p
    p_distances = np.power(distances, p)
    # avoid zero division
    p_distances[p_distances == 0] = eps
    inv_p_distances = 1.0 / p_distances
    grads = 1.0 / (p_distances * inv_p_distances.sum(axis=0, keepdims=True))
    grads = np.power(grads, m)
    return grads


def boltzmann_grad(X: np.ndarray, centers: np.ndarray,
                   alpha: float = 1.0) -> np.ndarray:
    distances = 0.5 * cdist(centers, X, metric='euclidean') ** 2
    dmin = distances.min(axis=0, keepdims=True)
    expd = np.exp(-alpha * (distances - dmin))
    sumexpd = expd.sum(axis=0, keepdims=True)
    sumdexpd = (distances * expd).sum(axis=0, keepdims=True)
    grads = (expd / sumexpd) * (1 - alpha * (distances - (sumdexpd / sumexpd)))
    return grads


def lloyd_iter(X: np.ndarray, centers: np.ndarray) -> np.ndarray:
    distances = cdist(centers, X, metric='euclidean')
    index = np.argmin(distances, axis=0)
    n_clusters = centers.shape[0]
    new_centers = centers.copy()
    for i in range(n_clusters):
        new_centers[i] = np.mean(X[index == i], axis=0)
    return new_centers


def cmeans_iter(X: np.ndarray, centers: np.ndarray,
                m: float = 2.0, eps: float = 1e-8) -> np.ndarray:
    p = 1.0 / (m - 1.0)
    weights = pnorm_grad(X, centers, p, eps)
    new_centers = (weights @ X) / weights.sum(axis=1, keepdims=True)
    return new_centers


def equilibrium_iter(X: np.ndarray, centers: np.ndarray,
                     alpha: float = 1.0) -> np.ndarray:
    weights = boltzmann_grad(X, centers, alpha)
    new_centers = (weights @ X) / weights.sum(axis=1, keepdims=True)
    return new_centers

実行結果

各アルゴリズム 100 反復実行する。

lloyd_centers = centers.copy()
for _ in range(100):
    lloyd_centers = lloyd_iter(X, lloyd_centers)

cmeans_centers = centers.copy()
for _ in range(100):
    cmeans_centers = cmeans_iter(X, cmeans_centers)

equilibrium_centers = centers.copy()
for _ in range(100):
    equilibrium_centers = equilibrium_iter(X, equilibrium_centers)

plt.figure(figsize=(16, 3))

plt.subplot(141)
plt.scatter(x1[:, 0], x1[:, 1], alpha=0.5, color=colors(0))
plt.scatter(x2[:, 0], x2[:, 1], alpha=0.5, color=colors(2))
plt.scatter(x3[:, 0], x3[:, 1], alpha=0.5, color=colors(1))
plt.scatter(means[:, 0], means[:, 1], marker='+', color='k', label='True center')
plt.title('Ground Truth')

plt.subplot(142)
index = np.argmin(cdist(lloyd_centers, X, metric='euclidean'), axis=0)
for i in range(lloyd_centers.shape[0]):
    plt.scatter(X[index == i][:, 0], X[index==i][:, 1], alpha=0.5, color=colors(i))
plt.scatter(means[:, 0], means[:, 1], marker='+', color='k', label='True center')
plt.scatter(lloyd_centers[:, 0], lloyd_centers[:, 1], marker='x', color='k', label='Centroid')
plt.legend()
plt.title('K-Means')

plt.subplot(143)
index = np.argmin(cdist(cmeans_centers, X, metric='euclidean'), axis=0)
for i in range(cmeans_centers.shape[0]):
    plt.scatter(X[index == i][:, 0], X[index==i][:, 1], alpha=0.5, color=colors(i))
plt.scatter(means[:, 0], means[:, 1], marker='+', color='k', label='True center')
plt.scatter(cmeans_centers[:, 0], cmeans_centers[:, 1], marker='x', color='k', label='Centroid')
plt.legend()
plt.title('Fuzzy C-Means')

plt.subplot(144)
index = np.argmin(cdist(equilibrium_centers, X, metric='euclidean'), axis=0)
for i in range(equilibrium_centers.shape[0]):
    plt.scatter(X[index == i][:, 0], X[index==i][:, 1], alpha=0.5, color=colors(i))
plt.scatter(means[:, 0], means[:, 1], marker='+', color='k', label='True center')
plt.scatter(equilibrium_centers[:, 0], equilibrium_centers[:, 1], marker='x', color='k', label='Centroid')
plt.legend()
plt.title('Equilibrium K-Means')

plt.show()

確かに、不均衡データにおいて EKM が他手法よりも優れた性能を示している。

おまけ①

なぜ EKM が不均衡データに対して優れた性能を示すのかについて、私個人の見解を示す。

Fuzzy c-means では、帰属度 $u_{kn}^{(\tau)} = \left(\sum_{i=1}^K \left(\frac{|\mathbf x_n - \mathbf c_k|_2}{|\mathbf x_n - \mathbf c_i|_2}\right)^{\frac{2}{m-1}}\right)^{-1}$ は明らかに正の値である。そのため、centroid はメンバーが多いクラスタに近づきやすい。
以下の図は、上記の数値実験の初期点における右上の centroid についての勾配を示す。

青のクラスタ付近の勾配は小さいが、データが大量にあるため、トータルの勾配としては青のクラスタへと近づく方向になってしまう。

一方で、equilibrium K-means では、重み $w_{kn}^{(\tau)}=\frac{\exp\left(-\alpha d_{kn}^{(\tau)}\right)}{\sum_{i=1}^K\exp\left(-\alpha d_{in}^{(\tau)}\right)} \left[1 - \alpha\left(d_{kn}^{(\tau)} - \frac{\sum_{i=1}^Kd_{in}^{(\tau)}\exp\left(-\alpha d_{in}^{(\tau)}\right)}{\sum_{i=1}^K\exp\left(-\alpha d_{in}^{(\tau)}\right)}\right)\right]$ は $1 - \alpha\left(d_{kn}^{(\tau)} - \frac{\sum_{i=1}^Kd_{in}^{(\tau)}\exp\left(-\alpha d_{in}^{(\tau)}\right)}{\sum_{i=1}^K\exp\left(-\alpha d_{in}^{(\tau)}\right)}\right) < 0$ のときに負の値となる。これは、

1 - \alpha\left(d_{kn}^{(\tau)} - \frac{\sum_{i=1}^Kd_{in}^{(\tau)}\exp\left(-\alpha d_{in}^{(\tau)}\right)}{\sum_{i=1}^K\exp\left(-\alpha d_{in}^{(\tau)}\right)}\right) = 1 - \alpha\left(\frac{\sum_{i=1}^K \left(d_{kn}^{(\tau)} - d_{in}^{(\tau)}\right) \exp\left(-\alpha d_{in}^{(\tau)}\right)}{\sum_{i=1}^K\exp\left(-\alpha d_{in}^{(\tau)}\right)}\right)

より、データが centroid $\mathbf c_k$ から比較的遠い ($d_{kn}$ が $d_{in}$ より大きい) ときに負の値となる。Fuzzy c-means の例と同様に、実際に勾配を見てみると、確かにクラスタの境界付近で勾配の方向が反転していることがわかる。

その結果、大きいクラスタに属するデータには負の重みが割り当てられ、遠ざかるように更新されるため、不均衡なデータに強いアルゴリズムとなっている。

以下、勾配の可視化コード。

n_grids = 20
x_grid = np.linspace(-6.0, 6.0, n_grids)
y_grid = np.linspace(-5.0, 7.0, n_grids)
x_mesh, y_mesh = np.meshgrid(x_grid, y_grid)
X_mesh = np.dstack([x_mesh, y_mesh]).reshape(-1, 2)

# Fuzzy c-means
i = 1
grads = pnorm_grad(X_mesh, centers)
center = centers[i][None]
grad = grads[i][:, None]
grad = (grad * (X_mesh - center)).reshape(n_grids, n_grids, 2)

plt.figure(figsize=(16, 9))
for i in range(lloyd_centers.shape[0]):
    plt.scatter(X[index == i][:, 0], X[index==i][:, 1], alpha=0.2, color=colors(i))
plt.quiver(x_mesh, y_mesh, grad[..., 0], grad[..., 1])
plt.scatter(centers[:, 0], centers[:, 1], marker='x', color='r', label='Centroid')
plt.savefig('grad_cmeans.png')
plt.show()

# Equilibrium K-means
i = 1
grads = boltzmann_grad(X_mesh, centers)
center = centers[i][None]
grad = grads[i][:, None]
grad = (grad * (X_mesh - center)).reshape(n_grids, n_grids, 2)

plt.figure(figsize=(16, 9))
for i in range(lloyd_centers.shape[0]):
    plt.scatter(X[index == i][:, 0], X[index==i][:, 1], alpha=0.2, color=colors(i))
plt.quiver(x_mesh, y_mesh, grad[..., 0], grad[..., 1])
plt.scatter(centers[:, 0], centers[:, 1], marker='x', color='r', label='Centroid')
plt.savefig('grad_ekm.png')
plt.show()

おまけ②

活性化関数 ReLU と Swish (SiLU) の関係は、今回の $\min$ と $\mathrm{boltz}$ の関係と対応している。

SiLU の参考
SiLU — PyTorch 2.2 documentation

Aloise, Daniel, et al. "NP-hardness of Euclidean sum-of-squares clustering." Machine learning 75 (2009): 245-248. ↩