k はどう決めれば良い？

代表的なのはエルボー法（SSE が急激に下がりにくくなる k を選ぶ）とシルエット係数（クラスター内とクラスター間距離のバランスを評価）です。事前知識がある場合はそれを優先します。

k-means が収束しない・初期値で結果が変わる場合は？

k-means++ 初期化を使うことで初期重心の悪影響を抑えられます。scikit-learn の KMeans は既定で k-means++ を使い、n_init で複数回の初期化から最良を選びます。

k-means が苦手なデータは？

クラスターの形状が球状でない、密度差が大きい、クラスターサイズが極端に異なるデータでは性能が落ちます。そのような場合は DBSCAN、HDBSCAN、Gaussian Mixture などの代替を検討します。

k-means

入門

02.5.1

k-means

最終更新 2020-01-29 読了時間 4 分

まとめ

k-means は、正解ラベルのないデータを「距離が近い点同士」で $k$ 個のグループに分ける、もっとも基本的なクラスターリング手法です。やることは 2 つの繰り返しだけで、各点をいちばん近いセントロイド（重心）に割り当て、次に各セントロイドをそのグループの平均の位置へ動かします。この反復は、各点と所属セントロイドとの距離の 2 乗和（WCSS）を一度も増やさないので必ず止まりますが、止まった先が最良の分け方とは限らず、初期位置と $k$ の選び方で結果が変わります。このページでは make_blobs の 3 つの塊を KMeans(n_clusters=3, n_init=10) で分け、エルボー法とシルエット係数で $k$ を選ぶ手順、n_init が WCSS をどう安定させるかをコードで確かめます。

直感 #

学年全員の「勉強時間」と「睡眠時間」を散布図にすると、誰がどのタイプかは書いてなくても、点の集まりがいくつか見えます。k-means はこの集まりを機械的に見つけます。まず $k$ 個のセントロイドを適当な位置に置き、各点を「いちばん近いセントロイド」のグループに入れ、セントロイドをそのグループの平均の位置に動かす。この割り当てと移動を、セントロイドが動かなくなるまで繰り返すだけです。アルゴリズムは単純ですが、初期位置が悪いと別の答えに落ちることがあり、$k$ も人が決める必要があるので、n_init やシルエット係数のような補助と組み合わせて使います。

データ X k 個のセントロイドを初期化各点をいちばん近いセントロイドに割り当てセントロイドをクラスターの平均に更新割り当てが変わらない？最終クラスター

まだ変わる: 割り当てに戻る

何が起きているか（ステップで見る） #

コードを読む前に、アルゴリズムが図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと、k-means の 1 反復ごとの変化が見られます。

見るポイント：割り当てのステップでは重心（×印）は動かず、点の色だけが変わります。更新のステップはその逆です。色は変わらず、重心だけが「自分の色の点の平均」へ移動します。この 2 つを交互に繰り返しているだけです。右上の WCSS（各点と所属重心との距離の 2 乗和）は、どのステップでも増えません。これが「必ず止まる」理由です。最初の 3 つの重心はわざと偏った位置に置いています。それでも数回で安定しますが、初期位置が悪いと別の（局所的な）答えに落ちることがあり、k-means++ 初期化や n_init はこの弱点への対策です。

式で書くと何を最小にしているか #

目的関数 #

k-means が小さくしようとしている量は WCSS（Within-Cluster Sum of Squares）、つまり各サンプルと所属クラスターのセントロイドとの距離の 2 乗和です。ステップ図の右上に出ていた数値がこれです。

$$ J = \sum_{k=1}^{K} \sum_{\mathbf{x}_i \in C_k} \lVert \mathbf{x}_i - \boldsymbol{\mu}_k \rVert^2 $$

$C_k$ はクラスター $k$ に属するサンプルの集合、$\boldsymbol{\mu}_k$ はクラスター $k$ のセントロイドです。全ての分け方を調べて $J$ を最小にする問題は NP 困難なので、次の 2 ステップを交互に繰り返して局所最適解で妥協します。

Lloyd のアルゴリズム #

まず割り当てステップで、各サンプル $\mathbf{x}_i$ をもっとも近いセントロイドのクラスターに入れます。図では点の色だけが変わる手順です。

$$ c_i = \arg\min_{k} \lVert \mathbf{x}_i - \boldsymbol{\mu}_k \rVert^2 $$

次に更新ステップで、各クラスターのセントロイドを、そのクラスターに属するサンプルの平均に置き直します。図では × 印だけが動く手順です。

$$ \boldsymbol{\mu}_k = \frac{1}{|C_k|} \sum_{\mathbf{x}_i \in C_k} \mathbf{x}_i $$

最後に、割り当てが変化しなくなるまでこの 2 ステップを繰り返します。scikit-learn の KMeans では max_iter（既定 300）と tol が打ち切りの条件です。

収束の保証 #

各ステップで目的関数 $J$ は単調に減少します。割り当てステップでは各サンプルをもっとも近いセントロイドに付け替えるだけなので、距離の 2 乗が増えることはありません。更新ステップも同じです。クラスター内の平均はユークリッド距離の 2 乗和を最小にする点なので、ここでも $J$ は減少または不変です。$J$ は 0 以上で単調非増加なので有限回のステップで必ず止まりますが、止まった場所が大域的最適解とは限りません。

k-means++ 初期化 #

ランダム初期化では 2 つのセントロイドが同じ塊の中に落ちて悪い局所解に陥りやすいので、k-means++ では初期セントロイドを互いに離れるように選びます。まず最初のセントロイド $\boldsymbol{\mu}_1$ をデータからランダムに 1 点選びます。次に各サンプル $\mathbf{x}_i$ について、もっとも近い既存セントロイドとの距離 $D(\mathbf{x}_i)$ を計算し、$D(\mathbf{x}_i)^2$ に比例する確率で次のセントロイドを選びます。これを $K$ 個そろうまで繰り返すと、遠くの点ほど選ばれやすいのでセントロイドがデータ全体に散らばります。この初期化には、期待値で最適解の $O(\log K)$ 倍以内という保証があります。

$$ \mathbb{E}[J_{\text{k-means++}}] \le 8(\ln K + 2) \cdot J_{\text{opt}} $$

k-means++ でセントロイド初期化割り当てステップ c_i = argmin ||x_i − μ_k||² 更新ステップ μ_k = mean(C_k) 割り当てに変化あり？収束: 局所最適解 n_init 回繰り返して最小の J を採用

変化あり: 割り当てステップに戻る

参考リンク

Lloyd, S. (1982). Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2), 129–137.

Arthur, D., & Vassilvitskii, S. (2007). k-means++: The Advantages of Careful Seeding. Proceedings of the 18th Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), 1027–1035.

Pythonで確かめる #

ライブラリと実験データ #

使うのは make_blobs（データ生成）、KMeans（本体）、silhouette_score（$k$ の評価）の 3 つだけです。

1
2
3
4
5
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

3 行目正解の塊の数が分かるデータを使う
make_blobs は指定した数の塊を作る関数で、正解が 3 つと分かっているデータに k-means をかけることで、エルボー法やシルエット係数が本当に 3 を指すかを検算できます。実データでは正解がないので、この検算を合成データで済ませてから使います。
4 行目KMeans は k-means++ と n_init を内蔵している
sklearn.cluster.KMeans は init="k-means++" が既定なので、上の数式で説明した初期化はこのクラスを使うだけで有効になります。n_init による複数回試行と inertia_（WCSS）の記録も内蔵しているので、自分で反復を書く必要はありません。
5 行目k を選ぶための指標を読み込む
silhouette_score は正解ラベルなしで分け方の良さを −1〜1 で評価する関数です。WCSS（inertia_）は $k$ を増やせば必ず下がるので $k$ の比較には使えず、この指標で補います。

3 つの塊を持つ 300 点の合成データを生成し、ラベルを使わずに散布図だけ描きます。

1
2
3
4
5
6
7
8
9
X, y_true = make_blobs(
    n_samples=300, centers=3,
    cluster_std=1.0, random_state=42
)

plt.figure(figsize=(6, 5))
plt.scatter(X[:, 0], X[:, 1], s=20, alpha=0.7)
plt.title("合成データ（3クラスター）")
plt.show()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1-4 行目散らばり 1.0 の塊を 3 つ作る
centers=3、cluster_std=1.0 で、塊同士が少し重なる程度のデータにしています。random_state=42 を固定しているのは、後のエルボー法やシルエット係数の図を本文の説明と一致させるためです。cluster_std を 2.5 に上げると塊が溶け合い、エルボーの位置が読みにくくなります。
6-9 行目正解ラベル y_true をあえて使わない
散布図を c= なしで単色に描いているのは、クラスターリングが「正解を知らない状態」で始まることを図で示すためです。y_true は make_blobs が返す正解ですが、k-means には渡しません。あとで labels と比較すれば、手法がどれだけ正解を復元したかを調べられます。
7 行目点を半透明にして重なりを見る
alpha=0.7 で点を半透明にしているので、塊が重なる場所は濃く見えます。この濃淡が、k-means が割り当てに迷う境界領域の目印になります。

合成データの散布図

基本的な学習 #

KMeans(n_clusters=3, n_init=10) で 3 つに分け、各点を所属ごとに塗り分けて収束後のセントロイドを赤い × で重ねます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
kmeans = KMeans(n_clusters=3, random_state=42, n_init=10)
labels = kmeans.fit_predict(X)

plt.figure(figsize=(6, 5))
plt.scatter(X[:, 0], X[:, 1], c=labels, cmap="viridis", s=20, alpha=0.7)
plt.scatter(
    kmeans.cluster_centers_[:, 0],
    kmeans.cluster_centers_[:, 1],
    c="red", marker="x", s=200, linewidths=3, label="セントロイド",
)
plt.legend()
plt.title("k-means クラスターリング結果")
plt.show()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1 行目k と初期化のやり直し回数を決める
n_clusters=3 が「いくつに分けるか」で、上のステップ図で置いた重心の数に対応します。n_init=10 は初期重心を変えて 10 回やり直し、WCSS が最小の結果だけを採用する指定です。初期位置が悪いと別の答えに落ちる弱点への保険で、random_state は再現性のためだけに固定しています。
2 行目fit_predict の中で反復が走る
この 1 行で、ステップ図の「割り当て → 更新」の反復がすべて実行され、各点の所属番号（0〜2）が返ります。番号はただのラベルで大小に意味はなく、実行環境によって色の割り振りが入れ替わることもあります。
5 行目所属ごとに色を塗って結果を見る
c=labels で点を所属クラスターごとに塗り分けています。ここで色の境目が「重心からの距離が等しい線」になっているか、ひとつの塊が 2 色に割れていないかを見るのが目的です。割れていれば k が多すぎる兆候です。
6-10 行目重心を × で重ねて位置を確認する
cluster_centers_ が収束後の重心です。各色の点の重心（平均の位置）に × が来ているはずで、×が点のない場所にあるならクラスターが細長い・輪状など k-means が苦手な形をしています。

k-meansクラスターリング結果

実行結果の読み方 #

3 つの塊がそれぞれ 1 色に塗られ、赤い × が各塊の中央に乗っていれば成功です。色の境目は隣り合う × から等距離の直線になり、cluster_std=1.0 で塊が重なった場所の点は、どちらの × に近いかだけで機械的に振り分けられます。

エルボー法でクラスター数を選ぶ #

WCSS（inertia_）を $k$ = 1 から 9 まで計算し、減少が鈍くなる「肘」の位置を $k$ の候補にします。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
inertias = []
K_range = range(1, 10)
for k in K_range:
    km = KMeans(n_clusters=k, random_state=42, n_init=10)
    km.fit(X)
    inertias.append(km.inertia_)

plt.figure(figsize=(8, 4))
plt.plot(K_range, inertias, "o-")
plt.xlabel("k（クラスター数）")
plt.ylabel("WCSS（慣性）")
plt.title("エルボー法")
plt.grid(True)
plt.show()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

2-6 行目k を 1 から 9 まで変えて WCSS を記録
range(1, 10) の各 $k$ で KMeans を学習し、収束後の WCSS が入っている inertia_ を集めます。各回に n_init=10 を付けているのは、初期値の運で WCSS が跳ねると肘の位置が読めなくなるからです。$k=1$ は全点を 1 つにまとめた場合で、WCSS が最大の基準点になります。
6 行目inertia_ は k を増やせば必ず下がる
inertia_ は $k$ を増やすほど単調に減り、$k$ = 点の数なら 0 です。だから「最小の WCSS」ではなく「下がり方が鈍る点」を探します。この性質のため、エルボー法だけでは $k$ を決め切れないことがあります。
8-14 行目折れ線で肘を目で探す
横軸 $k$、縦軸 WCSS の折れ線を描くと、正解が 3 つのデータでは $k=1 \to 2 \to 3$ で急落し、4 以降は緩やかになります。この折れ曲がりが「肘」で、plt.grid(True) は $k$ の位置を読み取りやすくするためです。

エルボー法によるk選択

シルエット係数でクラスター数を評価 #

シルエット係数は、各点について「同じクラスターの点との平均距離」と「いちばん近い他クラスターとの平均距離」を比べた値で、1 に近いほどクラスターが明確に分かれています。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
sil_scores = []
K_range2 = range(2, 10)
for k in K_range2:
    km = KMeans(n_clusters=k, random_state=42, n_init=10)
    labels_k = km.fit_predict(X)
    sil_scores.append(silhouette_score(X, labels_k))

plt.figure(figsize=(8, 4))
plt.plot(K_range2, sil_scores, "o-")
plt.xlabel("k（クラスター数）")
plt.ylabel("シルエット係数")
plt.title("シルエット係数によるk選択")
plt.grid(True)
plt.show()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

2 行目k=1 は評価できないので 2 から始める
シルエット係数は「他のクラスターとの距離」を使うので、クラスターが 1 つしかない $k=1$ では定義できません。range(2, 10) にしているのはそのためで、silhouette_score に $k=1$ のラベルを渡すとエラーになります。
5-6 行目ラベルと元データから係数を計算
silhouette_score(X, labels_k) は全点のシルエット係数の平均を返します。正解ラベルは使わず、X の距離と labels_k の分け方だけで計算する内部指標なので、実データでも使えます。計算量は点の数の 2 乗なので、数万点では sample_size 引数で間引きます。
8-14 行目最大になる k を読む
エルボー法と違い、シルエット係数は「最大の値を取る $k$」をそのまま候補にできます。正解が 3 つのこのデータでは $k=3$ が最大になり、エルボー法の肘と一致することを確認するのが図の目的です。両者がずれるときは、塊が重なっているか形が球状でない合図です。

シルエット係数によるk選択

n_initの役割 #

n_init を 1, 3, 10, 30 と変えて、初期値の試行回数が WCSS にどう影響するかを表示します。

1
2
3
4
5
6
results = []
for n in [1, 3, 10, 30]:
    km = KMeans(n_clusters=3, random_state=42, n_init=n)
    km.fit(X)
    results.append((n, km.inertia_))
    print(f"n_init={n:2d}  WCSS={km.inertia_:.2f}")

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

2-3 行目n_init だけを変えて比較する
n_clusters=3 と random_state=42 は固定し、n_init だけを 1 から 30 まで変えています。n_init=1 は初期化を 1 回しか試さないので初期位置の運がそのまま結果に出ます。30 回試せば WCSS が最小の結果を拾える確率が上がります。
4-5 行目inertia_ を集めて最良が選ばれたか見る
fit 後の inertia_ は、n_init 回の試行のうち WCSS が最小だったものです。n_init を増やしても値が変わらなければ既に最良の局所解に届いており、n_init=1 だけ値が大きければ初期値に当たり外れがあるデータだと分かります。
6 行目小数 2 桁で並べて差を読む
:.2f で WCSS を小数 2 桁に揃えて表示します。3 つの塊がはっきり分かれているこのデータでは 4 行とも同じ値になることが多く、差が出るのは塊が重なるデータや $k$ が大きい場合です。scikit-learn の既定値は n_init="auto" で、k-means++ なら 1 回になります。

k-meansはクラスターが球状で等分散であることを仮定しています。三日月型や密度が大きく異なるクラスターには適用できません。また、特徴量のスケールが異なる場合はStandardScalerで標準化してから適用してください。

エルボー法だけでは「肘」の位置が曖昧なことがあります。シルエット係数と併用し、複数の指標で $k$ を判断するのが実務的です。n_init=10以上に設定して初期値依存を減らすことも忘れずに行ってください。

使いどころと落とし穴 #

k-means が向いているのは、塊が球状で大きさがそろっており、点の数が数十万以上あって速さが要る場面です。計算は点の数 × $k$ × 反復回数に比例するだけなので、階層型クラスタリングや DBSCAN より大規模データに強い点が最大の利点です。落とし穴は 3 つあります。三日月や輪のような形は重心からの距離で割れないこと、特徴量のスケールが違うと値の大きい列だけで距離が決まること、そして外れ値 1 点が重心を引っぱることです。形が丸くないデータは DBSCAN、$k$ を決めたくないデータは X-means が次の選択肢です。

まとめ #

k-means は「各点をいちばん近いセントロイドに割り当てる」「セントロイドをその平均に動かす」の 2 手順を繰り返し、WCSS $J$ を減らしながらデータを $k$ 個に分けます。$J$ は増えないので必ず止まりますが局所解に落ちるため、init="k-means++" と n_init=10 以上で初期値の運を減らします。$k$ はエルボー法の肘とシルエット係数の最大値を併せて選び、正解が 3 つの make_blobs ではどちらも 3 を指しました。球状でないクラスターには DBSCAN など密度ベースの手法が適しています。

k-means++ — 初期セントロイドを D(x)² に比例する確率で離して選ぶ改良版
X-means — クラスター数を BIC で自動推定する発展手法
DBSCAN — 密度ベースのクラスターリング