数理科学メモ：クラスタリング（k-means法）

本記事の概要

教師なし学習におけるクラスタリングについて、その代表的なアルゴリズムの1つであるk-means法を試したいと思います。

k-means法（k平均法）

ある与えられたデータ群 $\mathbf{x}_k$に対して、各データがどのクラスタ$\omega_i$に所属するかを判定するのがクラスタリングです。

k-means法は、そのクラスタ数が既知の場合に有効なクラスタリングのアルゴリズムです。
クラスタが幾つあるのか不明な場合には使用できません。

アルゴリズムは以下の通りです。

1. データ群 $\mathbf{x}_k$ からランダムに各クラスタ$\omega_i$の重心を選ぶ
2. 全てのデータと重心の距離を計算する。各データについて、最も距離が短くなった重心（プロトタイプと書いてある本もある）と同じクラスタに属するとみなす
3. 各クラスタごとに、そのクラスに属するデータの平均ベクトル$\sum_{{\mathbf{x}_k \in \omega_i}} \mathbf{x}_k / n_i$を求めて、それを新しいクラスタの重心として更新する
4. データについて新しいクラスタに再割り当てされていれば、2番からやり直す。再割り当てがなければ、そこで終了する

理論的な詳細は教科書をご参照ください。
例えば、私は本記事を「続・わかりやすいパターン認識」の「第10章 クラスタリング」の10.2節を読んだ上で作成しています。

大域的最適解と局所的最適解

アルゴリズムを以下に実装します。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from scipy.spatial import Voronoi, voronoi_plot_2d

# k-meansの実装
class KMeans:
    def __init__(self, k=5, max_iters=100):
        self.k = k
        self.max_iters = max_iters

    def fit(self, X):
        # 初期の重心をランダムに選択
        idx = np.random.choice(len(X), self.k, replace=False)
        self.centroids = X[idx]
        
        for _ in range(self.max_iters):
            # 各データ点から重心までの距離を計算
            distances = np.linalg.norm(X[:, np.newaxis] - self.centroids, axis=2)
            # 最も近い重心のインデックスを取得
            self.labels = np.argmin(distances, axis=1)
            
            # 新しい重心を計算
            new_centroids = np.array([X[self.labels == i].mean(axis=0) if np.any(self.labels == i) 
                             else X[np.random.choice(len(X))] for i in range(self.k)])
            
            # 収束確認（重心が動かなくなったら終了）
            if np.all(self.centroids == new_centroids):
                break
            self.centroids = new_centroids

    def predict(self, X):
        distances = np.linalg.norm(X[:, np.newaxis] - self.centroids, axis=2)
        return np.argmin(distances, axis=1)

# 決定境界（ボロノイ図）の描画
def plot_geometric_boundaries(model, X):
    # 重心（Centroids）を取得
    centroids = model.centroids
    
    # ボロノイ図を計算
    vor = Voronoi(centroids)
    
    plt.figure(figsize=(10, 8))
    
    # 1. ボロノイ境界を描画
    # show_vertices=Falseで頂点を隠し、line_colorsで線の色を指定
    voronoi_plot_2d(vor, ax=plt.gca(), show_vertices=False, line_colors='orange', line_width=2, line_alpha=0.8)
    
    # 2. データ点の描画
    labels = model.predict(X)
    plt.scatter(X[:, 0], X[:, 1], c=labels, s=30, cmap='viridis', edgecolor='k', zorder=3)
    
    # 3. 重心の描画
    plt.scatter(centroids[:, 0], centroids[:, 1], s=200, marker='X', c='red', label='Centroids', zorder=4)
    
    # 4. 表示範囲の調整
    plt.xlim(X[:, 0].min() - 1, X[:, 0].max() + 1)
    plt.ylim(X[:, 1].min() - 1, X[:, 1].max() + 1)
    
    plt.title("K-means Decision Boundaries (Calculated Geometrically as Voronoi Diagram)")
    plt.legend()
    plt.grid(True, linestyle='--', alpha=0.5)
    plt.show()
    
# サンプルデータの生成
n_samples = 300
n_clusters = 5
X, y_true = make_blobs(n_samples=n_samples, centers=n_clusters, cluster_std=0.60, random_state=0)

# モデルの学習
model = KMeans(k=5)
model.fit(X)

# 実行
plot_geometric_boundaries(model, X)

各点 $x$ とクラス$\omega_i$番目の重心 $c_i$ の距離は$|x - c_i|^2$で求まるので、np.linalg.normで実装します。

次に各点 $x$ に対して、最も距離が近い重心のインデックスを割り当てます。

self.labels = np.argmin(distances, axis=1)

重心の更新は、クラス $\omega_i$ に属する全てのデータ点の平均値$\mu_i = \frac{1}{|n_i|} \sum_{x \in \omega_i} x$を新しい重心とします。

クラスに属するデータが1つもない場合、平均値が0で割った不定値になってしまうので、そのような場合はランダムな値で重心を再設定しておきます。

new_centroids = np.array([X[self.labels == i].mean(axis=0) if np.any(self.labels == i) 
                else X[np.random.choice(len(X))] for i in range(self.k)])

最大でループを100回繰り返すように設定していますが、100回より前でも、重心の更新がなくなった時点で終了するようにしています。

if np.all(self.centroids == new_centroids):
    break

また、決定境界としてボロノイ境界を描画しておきます。
ボロノイ境界とは2つの重心 $c_i$ と $c_j$ について以下の条件を満たす点 $x$ の集合で垂直二等分線になります。
$$|x - \mu_i| = |x - \mu_j|$$

ここではライブラリscipy.spatial.Voronoi を用いてボロノイ境界を計算しています。
Voronoi(centroids)は入力された全ての重心座標から、境界線同士がぶつかる頂点（vor.vertices）と、どの頂点同士が繋がっているか（vor.ridge_vertices）を算出し、voronoi_plot_2dにより描画します。
繋がりあう頂点は実線で結ばれ、繋がる先のない頂点は無限遠へと伸びる破線で示されます。

以上を実行すると、以下のように5つのクラスが分類され、各クラスの中心に重心があります。

初期値を決めた時からループを回す途中経過の重心と決定境界も描画すると以下のようになります。
最初はランダムに決められた誤った重心が、ループを繰り返すことで修正されていく様子が分かります。

ですが、ランダムに取った重心の初期値によっては、このように大域的最適解が求まらないこともあります。
例えば、高い頻度で以下のように結果が出力されたりもします。
明らかに異なる２つのクラスタを1つのクラスタとみなし、1つのクラスタを２つに分割してしまいます。

このような場合の初期値を決めた時からループを回す途中経過の重心と決定境界も描画すると以下のようになります。
適切な分類がされる前に更新が止まり、誤った分類を局所最適解として学習しています。

局所最適解を避ける為に、複数回アルゴリズムを実行し、量子化誤差が最小の結果のみを採用するようにします。
量子化誤差とは以下のようにして求めます。
まず、各データ点に対し、すべての重心との距離を計算します。その中で最も短い距離を選びます。その最短距離を2乗して、全データについて和を取ったのが量子化誤差です。
数式で表現すると以下になります。

$$J = \sum_{i=1}^{n} \min_{k \in {1 \dots K}} |x_i - c_k|^2$$

最適な学習ができていれば、量子化誤差は小さくなっていきます。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from scipy.spatial import Voronoi, voronoi_plot_2d

# k-meansの実装（量子化誤差計算メソッドを追加）
class KMeans:
    def __init__(self, k=5, max_iters=100):
        self.k = k
        self.max_iters = max_iters
        self.centroids = None
        self.labels = None

    def fit(self, X):
        # 初期の重心をランダムに選択
        idx = np.random.choice(len(X), self.k, replace=False)
        self.centroids = X[idx]
        
        for _ in range(self.max_iters):
            # 各データ点から重心までの距離を計算
            distances = np.linalg.norm(X[:, np.newaxis] - self.centroids, axis=2)
            # 最も近い重心のインデックスを取得
            self.labels = np.argmin(distances, axis=1)
            
            # 新しい重心を計算（
            new_centroids = np.array([X[self.labels == i].mean(axis=0) if np.any(self.labels == i) 
                                      else X[np.random.choice(len(X))] for i in range(self.k)])
            
            # 収束確認（重心が動かなくなったら終了）
            if np.all(self.centroids == new_centroids):
                break
            self.centroids = new_centroids

    def compute_inertia(self, X):
        # 量子化誤差の計算
        distances = np.linalg.norm(X[:, np.newaxis] - self.centroids, axis=2)
        min_distances = np.min(distances, axis=1)
        return np.sum(min_distances**2)

    def predict(self, X):
        distances = np.linalg.norm(X[:, np.newaxis] - self.centroids, axis=2)
        return np.argmin(distances, axis=1)

# 決定境界（ボロノイ図）の描画
def plot_best_geometric_boundaries(model, X, inertia_val):
    centroids = model.centroids
    vor = Voronoi(centroids)
    
    plt.figure(figsize=(10, 8))
    
    # 1. ボロノイ境界を描画
    # show_vertices=Falseで頂点を隠し、line_colorsで線の色を指定
    voronoi_plot_2d(vor, ax=plt.gca(), show_vertices=False, 
                    line_colors='orange', line_width=2, line_alpha=0.8)
    
    # 2. データ点の描画
    labels = model.predict(X)
    plt.scatter(X[:, 0], X[:, 1], c=labels, s=30, cmap='viridis', edgecolor='k', zorder=3)
    
    # 3. 重心の描画
    plt.scatter(centroids[:, 0], centroids[:, 1], s=200, marker='X', 
                c='red', label='Centroids', zorder=4, edgecolor='white')
    
    # 4. 表示範囲の調整
    plt.xlim(X[:, 0].min() - 1, X[:, 0].max() + 1)
    plt.ylim(X[:, 1].min() - 1, X[:, 1].max() + 1)
    
    plt.title(f"Best K-means Geometric Boundaries (Min Inertia: {inertia_val:.2f})")
    plt.legend()
    plt.grid(True, linestyle='--', alpha=0.5)
    plt.show()

# サンプルデータの生成
n_samples = 300
n_clusters = 5
X, y_true = make_blobs(n_samples=n_samples, centers=n_clusters, cluster_std=0.60, random_state=0)

# 10回の試行を行い、量子化誤差が最小のモデルを選択
best_model = None
min_inertia = np.inf

print("Selecting the best model out of 10 runs...")
for i in range(10):
    # モデルの学習
    current_model = KMeans(k=5)
    current_model.fit(X)
    inertia = current_model.compute_inertia(X)
    
    if inertia < min_inertia:
        min_inertia = inertia
        best_model = current_model
        
# 実行
plot_best_geometric_boundaries(best_model, X, min_inertia)

上記の量子化誤差の最小を選ぶ実装ならば、運が悪くなければ、大域的最適解が求まります。

scikit-learnを用いた計算

scikit-learnを用いれば、k-means法は簡単に実装できます。

from sklearn.cluster import KMeans

model = KMeans(n_clusters = 3, random_state = 0)
model.fit(data)

dataはデータフレームや2次元numpy配列を用います。

scikit-learnのk-means法は、初期値のランダム性によりアルゴリズムのループ数が多くなってしまったり、局所最適解になってしまう問題を解決する為に、初期値をなるべく離れて選ぶk-means++法というものを用いています。

結果は

model.labels_

で取得できます。

詳細は機械学習の実務寄りの解説書をご覧ください。
例えば、「スッキリわかるPythonによる機械学習入門」の第15章「教師なし学習2：クラスタリング」では与えるデータの整形方法も含めてscikit-learnのk-means法について解説されています。