k-means法(k平均法)とは?アルゴリズムの仕組みと弱点を図解でわかりやすく解説

手元に1万件の顧客データがあるとします。正解ラベルは付いていません。それでも「似た者どうしをいくつかのグループに分けたい」— このとき最初に試す手法が k-means法(k平均法) です。

k-meansがやっていることは驚くほど単純です。「各点を一番近い代表点に割り当てる」「代表点を担当した点の平均に動かす」。この2つを交互に繰り返すだけ。それだけで、データが勝手にグループへ分かれていきます。

k-meansの概念模式図:割り当てステップと重心更新ステップを交互に繰り返す様子

この図がk-meansのすべてです。左は「いまの重心(星印)」、真ん中は各点を最も近い重心の色に塗り分けた状態、右は塗り分けた結果にもとづいて重心を平均位置へ動かした状態です。②と③をひたすら繰り返すと、重心はデータのかたまりの中心へ吸い寄せられていきます。重要なのは、初期の星印の位置が適当でも構わないという点です — なぜそれで正しい位置に落ち着くのかを、この記事で数式と実測の両面から確かめます。

なぜk-meansを学ぶのか

k-meansは「教科書に載っている基本手法」で終わる話ではありません。実務でも研究でも次のような場面で今も現役です。

  • 前処理としてのベクトル量子化 — 大量の特徴ベクトルを $K$ 個の代表ベクトルに置き換えて圧縮する。画像の減色処理、Bag-of-Visual-Wordsによる画像検索、音声コーデックのコードブック設計はすべてk-meansの応用です。
  • 探索的データ解析 — 主成分分析で次元を落としてからk-meansでグループ分けし、「そもそもこのデータには何種類のパターンがあるのか」を掴む。異常検知の前段としてもよく使われます。
  • 他手法の初期化 — 混合ガウスモデル(GMM)のEMアルゴリズムは初期値に敏感なので、まずk-meansを回して得た重心を初期値に使うのが定石です。

さらに、k-meansの構造そのものが学びの宝庫です。「隠れ変数(どのクラスタに属するか)とパラメータ(重心)を交互に最適化する」という骨格は、EMアルゴリズム・変分推論・行列分解など、機械学習の広い範囲で繰り返し現れます。k-meansを完全に理解しておくと、これらの手法が「あのパターンか」と見えるようになります。

この記事では、単にアルゴリズムを紹介するのではなく、なぜそのステップになるのかを導出し、なぜ必ず停止するのかを証明し、どこで失敗するのかを実測で示すところまで踏み込みます。

k-meansが解こうとしている問題

まず、k-meansが何を「良いクラスタリング」と考えているのかをはっきりさせましょう。

直感的には、良いクラスタリングとは「同じグループの点どうしが近く、違うグループの点は遠い」状態です。k-meansはこのうち前半だけを数式にします。すなわち、各グループに代表点(重心)を1つ置き、各点と自分の重心との距離の2乗の合計を最小にする、という基準です。

いま、$D$ 次元のデータが $N$ 個あるとします。

$$ \bm{X} = \{ \bm{x}_1, \bm{x}_2, \dots, \bm{x}_N \}, \quad \bm{x}_n \in \mathbb{R}^D $$

これを $K$ 個のクラスタに分けたい。クラスタ $k$ の重心を $\bm{\mu}_k \in \mathbb{R}^D$ と書きます。そして、点 $\bm{x}_n$ がクラスタ $k$ に属するかどうかを表す変数を導入します。

$$ r_{nk} \in \{0, 1\}, \quad \sum_{k=1}^{K} r_{nk} = 1 $$

$r_{nk}$ は、$\bm{x}_n$ がクラスタ $k$ に属するなら $1$、属さないなら $0$ を取ります。制約 $\sum_k r_{nk} = 1$ は「どの点もちょうど1つのクラスタに属する」ことを意味します。この表現は 1-of-K符号化(one-hot表現)と呼ばれ、$N \times K$ 個の $0/1$ 変数からなります。

この記法のもとで、k-meansの目的関数は次のように書けます。

$$ \begin{equation} J = \sum_{n=1}^{N} \sum_{k=1}^{K} r_{nk} \, \| \bm{x}_n-\bm{\mu}_k \|^2 \end{equation} $$

$r_{nk}$ が $0$ か $1$ しか取らないので、内側の和は「$\bm{x}_n$ が実際に属しているクラスタの重心との距離の2乗」だけを残します。つまり $J$ は、全データ点について「自分の重心までの距離の2乗」を足し上げたものです。統計の言葉では クラスタ内平方和(Within-Cluster Sum of Squares, WCSS) と呼ばれ、$J$ が小さいほど各クラスタがぎゅっとまとまっていることを意味します。

ここで気をつけたいのは、$J$ を小さくする変数が2種類あることです。

  • $r_{nk}$(どの点をどのクラスタに入れるか)— 離散変数
  • $\bm{\mu}_k$(重心をどこに置くか)— 連続変数

離散変数と連続変数が混ざった最小化問題で、しかも $r_{nk}$ の取り方は $K^N$ 通りあります。$N = 600$、$K = 3$ でも $3^{600}$ 通り — 全探索は絶対に不可能です。実際、k-meansの厳密解を求める問題はNP困難であることが知られています。

そこで使うのが「片方を固定してもう片方を最適化する」という戦略です。次のセクションで、この戦略がなぜ驚くほどうまく機能するのかを見ていきます。

2ステップの導出:なぜ「最近傍割り当て」と「平均」なのか

k-meansのアルゴリズムは、よく次の2ステップとして紹介されます。

  1. 割り当てステップ:各点を、最も近い重心のクラスタに割り当てる
  2. 更新ステップ:各クラスタの重心を、そのクラスタに属する点の平均に更新する

多くの入門記事はここで終わりますが、この2つは天下り的な手続きではなく、目的関数 $J$ の最小化から導かれるものです。順に導出しましょう。

ステップ1:$\bm{\mu}_k$ を固定して $r_{nk}$ を最適化する

重心 $\bm{\mu}_1, \dots, \bm{\mu}_K$ を固定したとき、$J$ を最小にする $r_{nk}$ を求めます。ここで $J$ の形をよく見ると、$n$ についての和が外側にあり、$n$ が違う項どうしは互いに影響しません。

$$ J = \sum_{n=1}^{N} \underbrace{\left( \sum_{k=1}^{K} r_{nk} \| \bm{x}_n-\bm{\mu}_k \|^2 \right)}_{\text{点 } n \text{ だけの寄与}} $$

つまり $J$ の最小化は $N$ 個の独立した小問題に分解できます。点 $n$ について、制約 $\sum_k r_{nk} = 1$ のもとで $\sum_k r_{nk} \|\bm{x}_n-\bm{\mu}_k\|^2$ を最小化すればよい。$r_{nk}$ はちょうど1つだけ $1$ にできるので、距離の2乗が最小になる $k$ を選ぶのが答えです。

$$ \begin{equation} r_{nk} = \begin{cases} 1 & \text{if } k = \displaystyle\arg\min_{j} \| \bm{x}_n-\bm{\mu}_j \|^2 \\ 0 & \text{otherwise} \end{cases} \end{equation} $$

「一番近い重心に割り当てる」という直感的なルールが、目的関数から厳密に導かれました。しかもこれはその時点での大域最適です($\bm{\mu}$ を固定した条件のもとで、これ以上良い $r$ は存在しない)。

ステップ2:$r_{nk}$ を固定して $\bm{\mu}_k$ を最適化する

次に割り当て $r_{nk}$ を固定し、$J$ を最小にする $\bm{\mu}_k$ を求めます。$J$ は $\bm{\mu}_k$ について2次関数(下に凸)なので、微分してゼロと置けば最小値が得られます。

まず距離の2乗を内積の形に開いておきます。これは後の微分を見通しよくするためです。

$$ \| \bm{x}_n-\bm{\mu}_k \|^2 = (\bm{x}_n-\bm{\mu}_k)^\top (\bm{x}_n-\bm{\mu}_k) $$

これを $\bm{\mu}_k$ で微分します。ベクトルの2次形式の微分公式 $\frac{\partial}{\partial \bm{a}} (\bm{a}^\top \bm{a}) = 2\bm{a}$ と連鎖律を使うと、

$$ \frac{\partial}{\partial \bm{\mu}_k} \| \bm{x}_n-\bm{\mu}_k \|^2 = -2 (\bm{x}_n-\bm{\mu}_k) $$

これを $J$ 全体に適用します。$J$ の $\bm{\mu}_k$ に関係する項は $k$ 番目のクラスタの分だけなので、

$$ \frac{\partial J}{\partial \bm{\mu}_k} = -2 \sum_{n=1}^{N} r_{nk} (\bm{x}_n-\bm{\mu}_k) $$

停留条件として、これをゼロベクトルと置きます。両辺を $-2$ で割ると、

$$ \sum_{n=1}^{N} r_{nk} (\bm{x}_n-\bm{\mu}_k) = \bm{0} $$

和を2つに分けて、$\bm{\mu}_k$ が $n$ に依存しない定数であることを使って外に出すと、

$$ \sum_{n=1}^{N} r_{nk} \bm{x}_n = \bm{\mu}_k \sum_{n=1}^{N} r_{nk} $$

ここで $\sum_n r_{nk}$ は「クラスタ $k$ に属する点の個数」$N_k$ にほかなりません。両辺を $N_k$ で割って $\bm{\mu}_k$ について解くと、

$$ \begin{equation} \bm{\mu}_k = \frac{\sum_{n=1}^{N} r_{nk} \bm{x}_n}{\sum_{n=1}^{N} r_{nk}} = \frac{1}{N_k} \sum_{n \in C_k} \bm{x}_n \end{equation} $$

右辺はまさに クラスタ $k$ に属する点の標本平均 です。「重心を平均に置く」というルールもまた、目的関数から導かれました。k-means(k平均法)という名前の “means”(平均)はここから来ています。

なお、この導出は距離の2乗を使っていることに本質的に依存しています。もし2乗ではなく距離そのものの和を最小化するなら、最適な代表点は平均ではなく中央値(多次元では幾何学的中央値)になり、これは k-medians と呼ばれる別の手法になります。「平均」が出てくるのは2乗誤差を選んだ結果なのです。

2つのステップがどちらも「その時点での最適解」を与えることがわかりました。では、これを交互に繰り返すと何が起きるのでしょうか。

Lloydアルゴリズム:交互最適化の実際の挙動

ステップ1とステップ2を交互に繰り返す手続きは、提案者の名前をとって Lloydアルゴリズム と呼ばれます(1957年にベル研究所のStuart Lloydがパルス符号変調の量子化のために考案し、1982年に公刊されました)。手順をまとめると次のとおりです。

  1. 重心 $\bm{\mu}_1, \dots, \bm{\mu}_K$ を何らかの方法で初期化する
  2. 割り当てステップ:式(2)により全点の $r_{nk}$ を更新する
  3. 更新ステップ:式(3)により全重心 $\bm{\mu}_k$ を更新する
  4. 重心が動かなくなるまで2と3を繰り返す

実際に動かしてみましょう。ここでは2次元の人工データを使います。3つの2次元ガウス分布から200点ずつ、合計600点をサンプリングし、$K = 3$ でクラスタリングします。初期重心はわざと悪い位置 — 3つとも左下の狭い範囲に固めて — 置いてみます。

Lloyd法の反復ごとの重心移動と割り当ての変化

黒い線が重心の移動軌跡です。初期状態(左上)では $J = 4537.6$ でしたが、1回目の反復で $J = 2037.1$ とほぼ半分に落ちています。ここで注目したいのは、反復2回目から6回目までほとんど絵が変わらないことです($J$ は $1665.2 \to 1648.9$ とわずかしか減らない)。3つのうち2つの重心が右側の大きなかたまりを分け合ってしまい、左上のかたまりが放置された「行き詰まり状態」に入っているのです。それでも重心はじりじりと動き続け、最終的に28回目の反復で $J = 473.6$ という良い解にたどり着きました。

この「長い停滞のあとに一気に改善する」挙動は、k-meansが単純な勾配降下ではなく離散的な割り当ての組み替えによって進むことに由来します。割り当てが変わらない限り $J$ はほとんど減らず、境界付近の点が1つ乗り換えた瞬間に状況が動き出すのです。

では、この反復は本当に必ず止まるのでしょうか。次はそれを保証する性質を見ます。

なぜ必ず止まるのか:単調減少と有限性

k-meansには、実用上きわめて重要な保証があります。目的関数 $J$ は反復ごとに決して増えず、有限回で必ず停止するというものです。証明は驚くほど短く、しかも前節の導出そのものが証明になっています。

(a) 割り当てステップで $J$ は増えない。 ステップ1は $\bm{\mu}$ を固定したうえで $J$ を最小にする $r$ を選んでいます。直前の $r$ も選択肢の1つに含まれているので、新しい $r$ での $J$ が古い $r$ での $J$ を上回ることはありません。

(b) 更新ステップで $J$ は増えない。 ステップ2は $r$ を固定したうえで $J$ を最小にする $\bm{\mu}$ を与えています。$J$ は $\bm{\mu}_k$ について下に凸な2次関数なので、式(3)の停留点は確かに最小点です。やはり直前の $\bm{\mu}$ も選択肢に含まれるので、$J$ は増えません。

(c) 有限回で停止する。 $J$ は下に有界(常に $J \geq 0$)で、かつ単調非増加です。ここで効いてくるのが、割り当て $r$ が決まれば重心 $\bm{\mu}$ が一意に決まるという事実です。割り当ての総数は高々 $K^N$ 通りなので、アルゴリズムが取りうる $J$ の値も有限個しかありません。有限個の値の上で狭義に減少し続けることは不可能なので、いつかは $J$ が変化しなくなり、そこで停止します。

この議論を数値で確かめてみましょう。図2と同じ「わざと悪い初期値」の場合について、割り当てステップ直後と更新ステップ直後の $J$ を半ステップ刻みで記録しました。

k-meansの目的関数が半ステップごとに単調減少し有限回で停止する様子

左のグラフは階段状に $J$ が落ちていく様子です。56個の半ステップすべてで $J$ は減少しており、増加した箇所は0回でした(実測)。右のグラフは各半ステップでの減少量です。ここで注意したいのは、減少量が単調に小さくなるわけではないことです。半ステップ番号10〜35あたりでは減少量が $0.3 \sim 1.0$ 程度の低空飛行を続け、40付近で $10^2$ オーダーへ急増しています。前節で見た「長い停滞のあとの急改善」が、減少量のグラフにもはっきり現れました。

ただし、この保証には重大な但し書きがあります。「$J$ が減らなくなった」ことは「$J$ が最小になった」ことを意味しません。保証されるのは局所解への収束だけです。次のセクションで、この弱点がどれほど深刻かを実測で見ていきます。

初期値依存:同じデータから違う答えが出る

k-meansを実務で使うときに最初にぶつかる壁が、初期値依存です。同じデータ、同じ $K$ でも、初期重心の置き方次第で別の解に落ちます。

これを見るために、局所解が現れやすいデータを用意しました。5つのかたまりが非対称に配置され、そのうち2つ(左下の2群)が互いに近い、という配置です。$K = 5$ として、初期重心をデータの範囲内の一様乱数で置き、seedだけを変えて3回走らせました。

初期値を変えると異なる局所解に落ちるk-meansの初期値依存性

3つとも同じデータ・同じアルゴリズム・同じ $K$ です。それでも結果はまったく違います。左(seed=2)は5つのかたまりを正しく分離して $J = 336.4$。中央(seed=7)は左下の2群を1つにまとめてしまい、代わりに上のかたまりを2つに割って $J = 542.4$。右(seed=29)はさらに悪く、$J = 974.5$ で、しかも担当する点が0個の「空クラスタ」が発生しています(左上に浮いている星印がそれです)。最悪の解は最良の解の 2.90倍 の目的関数値になりました。

空クラスタは実装上も厄介です。式(3)の分母 $N_k$ がゼロになるため、素直に実装するとゼロ除算になります。実務では「最も遠い点をその重心に割り当て直す」「最大クラスタを分割する」といった対処を入れます。scikit-learnも内部で同様の再配置を行っています。

では、この初期値依存はどれくらいの頻度で起きるのでしょうか。同じデータで初期値を200通り変えて、収束後の $J$ の分布を調べました。

ランダム初期化とk-means++の目的関数値の分布比較

オレンジがランダム初期化の結果です。200回のうち最良解($J = 336.4$)に到達したのは42.5%にすぎません。平均は $J = 487.7$、最悪は $J = 1204.8$ で、標準偏差は $172.4$ にもなります。つまり素朴なランダム初期化でk-meansを1回だけ走らせると、6割近い確率で悪い解を掴まされるということです。これは「収束が保証されている」という安心感が、いかに当てにならないかを示しています。

同じ図の青いヒストグラムは、次に説明する k-means++ で初期化した場合です。まずはその仕組みを見ましょう。

k-means++:初期値の引き方を賢くする

初期値依存への最も効果的な対処が k-means++(Arthur & Vassilvitskii, 2007)です。アイデアは一言で言えば「初期重心どうしをできるだけ散らばらせる」。ただし、単に一番遠い点を選ぶ(最遠点法)と外れ値を掴んでしまうので、距離の2乗に比例した確率でランダムに選ぶという折衷を取ります。

手順は次のとおりです。

  1. 1個目の重心 $\bm{c}_1$ を、データ点から一様ランダムに選ぶ
  2. 各点 $\bm{x}$ について、すでに選んだ重心までの最短距離 $D(\bm{x}) = \min_j \|\bm{x}-\bm{c}_j\|$ を計算する
  3. 次の重心を、確率 $\dfrac{D(\bm{x})^2}{\sum_{\bm{x}’} D(\bm{x}’)^2}$ でデータ点から選ぶ
  4. $K$ 個そろうまで2〜3を繰り返す

この「$D^2$ 重み付きサンプリング」が実際にどう働くかを可視化しました。

k-means++が距離の2乗に比例した確率で次の重心を選ぶ仕組み

1個目の重心はたまたま3番目のかたまりの中から選ばれました。中央のパネルは各点の選択確率を色と大きさで表したもので、遠いかたまりほど濃く大きくなっています。実測では、最も確率の高い点は一様抽選($1/600 = 0.00167$)の2.8倍にあたる $0.0046$ でした。クラスタ単位で合計すると差はさらに明確で、1個目の重心があるかたまりが選ばれる合計確率はわずか 4.7% まで下がり、最も遠いかたまりが 58.3% を占めます(一様抽選ならどれも33.3%)。右のパネルで実際に60回抽選した結果を見ると、ほぼすべてが遠いかたまりから引かれています。

この効果が最終的な解の質にどう効くかは、先ほどの図5に戻ると読み取れます。

  • 最良解に到達した割合:ランダム 42.5% → k-means++ 69.5%(悪い局所解に落ちる確率が 57.5% → 30.5% とほぼ半減)
  • 平均の目的関数値:$487.7$ → $419.7$(13.9% の改善)
  • 平均反復回数:$8.77$ → $6.21$(約29%高速化)

初期化に少し手間をかけるだけで、解の質と収束速度が同時に改善しています。理論面でも、k-means++の期待目的関数値は最適解 $J_{\text{opt}}$ に対して $\mathbb{E}[J] \leq 8(\ln K + 2) J_{\text{opt}}$ という保証があることが証明されており(Arthur & Vassilvitskii, 2007)、素朴なランダム初期化には何の保証もないのと対照的です。

もう1つ、実務でほぼ必ず併用するのが 複数回試行(n_init) です。初期値を変えて何回か走らせ、$J$ が最小のものを採用します。今回のデータでの実測は次のとおりでした。

試行回数 ランダム初期化 k-means++
1回 42.5% 69.5%
5回 97.5% 100%
10回 100% 100%

5回まわすだけで、どちらの初期化でもほぼ確実に最良解が得られます。scikit-learnの KMeans がデフォルトで init="k-means++" かつ n_init を複数回に設定しているのは、この2つの効果を最初から享受させるためです。

初期値の問題に目処がついたところで、次はもう1つの大きな悩み — クラスタ数 $K$ をどう決めるか — に進みます。

クラスタ数 $K$ をどう決めるか

k-meansは $K$ を人間が与えなければ動きません。ところが実データでは「いくつに分かれているか」こそ知りたいことです。この鶏と卵の問題に対する定番の道具が2つあります。

エルボー法

素朴に「$J$ が小さくなる $K$ を選ぶ」とすると失敗します。なぜなら $J$ は $K$ を増やせば必ず減るからです。極端な話、$K = N$ とすれば各点が自分だけのクラスタになり $J = 0$ になります。そこで見るのは $J$ の絶対値ではなく 減り方の変化 — グラフが「肘(elbow)」のように折れ曲がる点を探します。

シルエット係数

もう1つはクラスタの「まとまりの良さ」と「離れ具合」を直接測る指標です。点 $i$ について、

  • $a(i)$:同じクラスタ内の他の点との平均距離(小さいほど良い=まとまっている)
  • $b(i)$:他のクラスタのうち最も近いクラスタの点との平均距離(大きいほど良い=離れている)

と定義し、シルエット値を次で計算します。

$$ \begin{equation} s(i) = \frac{b(i)-a(i)}{\max\{a(i),\, b(i)\}} \in [-1, 1] \end{equation} $$

分母の $\max$ は値を $[-1, 1]$ に収めるための正規化です。$s(i)$ が $1$ に近ければ「自分のクラスタにしっかり収まっている」、$0$ 付近なら「境界上で曖昧」、負なら「隣のクラスタのほうが近い=割り当てが怪しい」ことを意味します。全点の平均をとった値でクラスタ数を比較します。

先ほどの3クラスタのデータで、$K$ を1から8まで変えて両方を計算しました。

エルボー法とシルエット係数によるクラスタ数Kの決め方

左のエルボー法では、$J$ は $K=1$ の $4086.5$ から単調に減り続けています。しかし減り方に注目すると、$K=2 \to 3$ で 72.7% も減ったのに、$K=3 \to 4$ では 19.5% しか減りません。この「折れ曲がり」が $K = 3$ に肘があることを示しています。右のシルエット係数はさらに明快で、$K=3$ で 0.714 と最大値を取り、$K=2$(0.576)や $K=4$(0.615)を明確に上回りました。真のクラスタ数が3であることを、両方の指標が正しく当てています。

平均値だけでなく、点ごとのシルエット値の分布を見るとさらに情報が得られます。

K=2,3,5のシルエット図の比較

各クラスタの点を、シルエット値の大きい順に横棒として並べた図です。$K=3$(中央)では3つのクラスタが同じような厚みで、どれも横棒が長く伸びています。一方 $K=2$(左)は片方のクラスタだけが極端に厚く、$K=5$(右)は細いクラスタが並んで平均も $0.510$ まで落ちています。クラスタの厚みが偏っていないか、短い棒(あるいは負の値)が多くないかを見るのがシルエット図の使い方です。

ただし、これらの指標には共通の限界があります。どちらも「球状のクラスタ」を暗黙に前提していることです。エルボー法の $J$ は距離の2乗和そのものですし、シルエット係数もユークリッド距離の平均で測っています。データが球状でなければ、両方とも人間の直感とはズレた $K$ を推す可能性があります。実データでは肘がはっきり見えないことも珍しくなく、最終的にはドメイン知識と組み合わせて決めることになります。

では、その「球状の前提」はどこから来ているのでしょうか。k-meansの決定境界の形を調べると、その正体がはっきり見えます。

決定境界はボロノイ図になる

k-meansが学習したモデルは、$K$ 個の重心だけです。新しい点が来たときの予測は「最も近い重心のクラスタ」— つまり空間全体が重心によって領域分割されます。この分割の形を数式で確かめましょう。

点 $\bm{x}$ が重心 $\bm{\mu}_1$ と $\bm{\mu}_2$ から等距離にある境界を考えます。

$$ \| \bm{x}-\bm{\mu}_1 \|^2 = \| \bm{x}-\bm{\mu}_2 \|^2 $$

両辺を内積で展開します。

$$ \| \bm{x} \|^2-2 \bm{x}^\top \bm{\mu}_1 + \| \bm{\mu}_1 \|^2 = \| \bm{x} \|^2-2 \bm{x}^\top \bm{\mu}_2 + \| \bm{\mu}_2 \|^2 $$

ここで両辺に共通の $\|\bm{x}\|^2$ が消えるのがポイントです。この項が消えることで、$\bm{x}$ の2次の項がなくなります。残りを整理すると、

$$ 2 \bm{x}^\top (\bm{\mu}_2-\bm{\mu}_1) = \| \bm{\mu}_2 \|^2-\| \bm{\mu}_1 \|^2 $$

これは $\bm{x}$ について 1次式 です。すなわち境界は直線(一般には超平面)になります。さらに右辺を変形すると、この超平面が2つの重心を結ぶ線分の垂直二等分線であることも確認できます。

$$ \begin{equation} \left( \bm{x}-\frac{\bm{\mu}_1 + \bm{\mu}_2}{2} \right)^\top (\bm{\mu}_2-\bm{\mu}_1) = 0 \end{equation} $$

左側の括弧は「中点から $\bm{x}$ へのベクトル」、右側は「重心どうしを結ぶ方向ベクトル」なので、この式は両者が直交することを述べています。こうして空間を最近傍の代表点で分割した図形を ボロノイ図(Voronoi diagram) と呼びます。

k-meansの決定境界がボロノイ分割(垂直二等分線)になること

左は実際に空間を色分けしたもので、境界が3本の直線分になっていることがわかります。右は2つの重心だけを取り出した拡大図で、境界(黒線)が重心を結ぶ灰色の線分と中点で直交している — まさに式(5)が示す垂直二等分線です。

この事実は、k-meansの能力と限界を同時に説明します。予測が「重心との距離を測るだけ」なので、$N$ 個の学習データを保持する必要がなく、$K$ 個の重心さえあれば新規データを即座に分類できます(これが k-NN との決定的な違いです)。その一方で、境界が必ず直線である以上、曲がった形のクラスタは絶対に表現できません。次のセクションで、この限界が実際にどう表面化するかを見ます。

k-meansが失敗するデータ形状

k-meansは万能ではありません。むしろ「うまくいく条件」がかなり限定的です。典型的な失敗パターンを4つ、実測つきで見てみましょう。正解ラベルがわかっている人工データを使い、純度(各クラスタの多数決ラベルが当たっている割合、1.0で完全一致)で評価します。

同心円・三日月・大きさの違う群・細長い群でk-meansが失敗する様子

上段が「こう分けたい」正解、下段がk-meansの結果です。それぞれ何が起きているかを見ていきます。

同心円(純度 0.513) — 内側の円と外側の輪に分けたいのに、k-meansは縦にすっぱり切ってしまいました。純度0.513はコイン投げとほぼ同じで、完全な失敗です。内側と外側は重心が同じ位置になるため、重心ベースの手法では原理的に区別できません。

三日月(純度 0.737) — 2つの弧を分けたいのに、直線の境界では弧の一部が反対側に混ざります。純度0.737は「半分は当たっている」程度で、実用には耐えません。

大きさが極端に違う2群(純度 0.909) — 広く散らばった500点の群と、小さく密集した50点の群です。k-meansは大きい群を真っ二つに割り、小さい群を右半分に吸収してしまいました。これは目的関数の性質から必然です。$J$ は距離の2乗和なので、点数が多く広がった群を分割するほうが、点数が少ない群を独立させるより $J$ が小さくなるのです。k-meansにはクラスタのサイズを揃えようとするバイアスが組み込まれていると言えます。

細長い2群(純度 0.538) — 横に長く伸びた2つの帯を上下に分けたいのに、k-meansは縦に切りました。距離が全方向で等しく扱われる(等方的)ため、「横方向の広がりは無視して縦方向で分けたい」という意図が伝わらないのです。

これら4つの失敗はすべて、前節で見た「境界が直線」「距離が等方的」「サイズを揃えたがる」という3つの性質から説明できます。逆に言えば、同じくらいの大きさの球状のかたまりが分離しているというk-meansの前提が満たされるなら、k-meansは高速で頑健な良い選択です。

前提が満たされない場合の対処としては、密度ベースのDBSCAN(任意形状に対応)、スペクトラルクラスタリング(グラフの連結性で分ける)、そして後述する混合ガウスモデル(楕円形に対応)が候補になります。

ところで、これらの失敗のうちいくつかは「距離の測り方」を変えるだけで改善できます。次にその話をします。

距離とスケーリング:前処理で結果が変わる

k-meansはユークリッド距離の2乗しか見ません。ということは、各特徴量のスケールがそのまま結果に効いてしまいます。これは実データで最も頻繁に踏む地雷です。

例として、「身長のような特徴」と「体重のような特徴」を持つ2群のデータを考えます。本当は身長方向(横方向)で2つに分かれているデータです。ここで、身長の単位を cm から m に変えるとどうなるでしょうか。データの中身は何も変わらず、数字が100分の1になるだけです。

特徴量のスケールを揃えないとクラスタリング結果が単位で変わること

結果は劇的です。cm・kgのまま実行すると純度 0.870 で正しく横方向に分かれますが、身長だけを m に変えただけで純度は 0.597 まで崩壊しました。図の3枚目を見ると、横方向ではなく縦方向(体重)で分かれてしまっています。理由は標準偏差を見れば明らかで、cm表記なら特徴1が $9.22$、特徴2が $8.51$ とほぼ同等なのに対し、m表記では特徴1が $0.0922$ となり、体重のほうが92.2倍もばらつきが大きいことになるからです。距離の2乗和はばらつきの大きい軸に支配されます。

対処は標準化(z-score変換)です。各特徴から平均を引き、標準偏差で割ります。

$$ z_{nd} = \frac{x_{nd}-\bar{x}_d}{\sigma_d} $$

標準化してから実行すると、m表記でも純度 0.863 に回復しました。さらに重要なのは、cm表記を標準化した結果とm表記を標準化した結果のラベルが100%一致したことです(実測)。標準化は線形変換に対して不変なので、単位の選び方に結果が依存しなくなります。

ただし標準化が常に正解とは限りません。今回のデータでは、標準化前(0.870)と標準化後(0.863)で純度がわずかに下がっています。もともと2つの特徴のばらつきが同程度だった場合、標準化は情報を足すのではなく少し均してしまうことがあるのです。「単位が異なる特徴が混ざっている」「スケールが桁で違う」ときには必須、そうでなければ効果は限定的、と考えるとよいでしょう。

より根本的な対処として、距離の測り方自体を変える方法もあります。特徴間に相関がある場合はマハラノビス距離を使う、あるいは主成分分析で相関を落としてからk-meansをかける、といった前処理が定石です。

ここまでで、k-meansの弱点は「等方的な距離」と「ハードな割り当て」に集約されることが見えてきました。この2つを緩めた手法が、次に見る混合ガウスモデルです。

混合ガウスモデル(GMM)との関係

k-meansと混合ガウスモデル(GMM) は、しばしば「兄弟」として並べられます。両者の関係は感覚的な類似ではなく、k-meansはGMMの厳密な極限であるという数学的な関係です。

GMMは、データが $K$ 個のガウス分布の混合から生成されたと考えます。EMアルゴリズムのEステップでは、点 $\bm{x}_n$ がクラスタ $k$ に属する確率(負担率)を計算します。

$$ \gamma_{nk} = \frac{\pi_k \, \mathcal{N}(\bm{x}_n \mid \bm{\mu}_k, \bm{\Sigma}_k)}{\sum_{j=1}^{K} \pi_j \, \mathcal{N}(\bm{x}_n \mid \bm{\mu}_j, \bm{\Sigma}_j)} $$

k-meansの $r_{nk}$ が $0$ か $1$ しか取らない(ハード割り当て)のに対し、$\gamma_{nk}$ は $0$ と $1$ の間の連続値を取ります(ソフト割り当て)。この違いが何をもたらすかを見てみましょう。広がりの大きさが極端に違う2つの群を用意しました。

k-meansのハード割り当てと混合ガウスモデルのソフト割り当ての比較

左のk-meansは、中央の密集した群を左右に割ってしまい純度 0.850 にとどまりました。中央のGMMは各成分の広がり(共分散)も学習するため、純度 0.967 と大きく改善しています。実際、GMMが推定した2成分の標準偏差は $0.37$ と $1.44$ で、真の値($0.35$ と $1.5$)をほぼ正確に捉えていました。右のヒストグラムは各点の「所属確率の最大値」の分布です。多くの点は確率1に近い一方、6.7%の点は確率0.9未満で、最小値は $0.5025$ — つまり「どちらとも言えない」点にGMMは中間的な確率を与えています。k-meansはこれらをすべて強制的に0か1に潰していることになります。

なぜk-meansがGMMの極限なのか

ここからが本題です。GMMの共分散をすべて $\bm{\Sigma}_k = \varepsilon \bm{I}$(等方的で、全クラスタ共通、分散は固定値 $\varepsilon$)に制限してみます。ガウス分布の密度は

$$ \mathcal{N}(\bm{x}_n \mid \bm{\mu}_k, \varepsilon \bm{I}) = \frac{1}{(2\pi\varepsilon)^{D/2}} \exp\left( -\frac{\| \bm{x}_n-\bm{\mu}_k \|^2}{2\varepsilon} \right) $$

となります。これを負担率の式に代入すると、分子・分母に共通する正規化定数 $(2\pi\varepsilon)^{-D/2}$ が約分されて消えます。混合比も等しい($\pi_k = 1/K$)とすれば、それも約分されて、

$$ \gamma_{nk} = \frac{\exp\left( -\| \bm{x}_n-\bm{\mu}_k \|^2 / 2\varepsilon \right)}{\sum_{j=1}^{K} \exp\left( -\| \bm{x}_n-\bm{\mu}_j \|^2 / 2\varepsilon \right)} $$

という、距離の2乗を温度 $2\varepsilon$ のソフトマックスに通した形になります。ここで $\varepsilon \to 0$ とすると何が起きるでしょうか。分子・分母を「最も近い重心の項」で割って整理すると、最近傍でない項はすべて $\exp(-(\text{正の値})/2\varepsilon) \to 0$ となって消えます。したがって、

$$ \gamma_{nk} \to \begin{cases} 1 & \text{if } k = \arg\min_j \| \bm{x}_n-\bm{\mu}_j \|^2 \\ 0 & \text{otherwise} \end{cases} $$

これは式(2)の $r_{nk}$ そのものです。分散を0に近づけた極限で、GMMのEステップはk-meansの割り当てステップに一致します。Mステップも同様に、$\gamma_{nk}$ が $0/1$ になれば重み付き平均は単純平均に退化し、k-meansの更新ステップになります。

この極限を実際に数値で確かめました。共分散を $\varepsilon \bm{I}$ に固定したGMMを、$\varepsilon$ を変えながら同じ初期値から回し、k-meansの結果と比較したものが次の図です。

共分散を等方固定した混合ガウスモデルが分散ゼロの極限でk-meansに一致すること

左のグラフは負担率の最大値(平均)です。$\varepsilon = 4.0$ では $0.499$ — つまり「どちらのクラスタか五分五分」というぼやけた状態ですが、$\varepsilon$ を小さくすると $0.930$($\varepsilon=2$)、$0.992$($\varepsilon=1$)と急速に1へ近づき、$\varepsilon = 0.25$ では数値上ぴったり $1.0000$ になります。右のグラフは推定された重心とk-meansの重心との距離で、$\varepsilon = 4.0$ で $1.69$ だったずれが、$\varepsilon = 1.0$ で $1.4 \times 10^{-2}$、$\varepsilon = 0.1$ では $1.4 \times 10^{-8}$ — 実質的に完全一致です。理論が予言したとおり、$\varepsilon \to 0$ でGMMはk-meansに収束しました。

この関係を理解しておくと、両者の使い分けが自然に決まります。

k-means GMM
割り当て ハード(0/1) ソフト(確率)
クラスタの形 球状のみ 楕円(共分散を学習)
クラスタの大きさ 揃えたがる 混合比 $\pi_k$ で調整可能
パラメータ数 $KD$ $KD + KD(D+1)/2 + K$
計算コスト 軽い 重い($D$ が大きいと特に)
所属の曖昧さ 表現できない 確率として表現できる

まずk-meansで大まかな構造を掴み、必要ならGMMへ進むというのが定石です。実際、GMMのEMアルゴリズムはk-meansの結果を初期値として使うのが一般的で、両者は競合ではなく補完関係にあります。

計算量と大規模データへの適用

最後に実務上の性能特性に触れておきます。1回の反復のコストは次のとおりです。

  • 割り当てステップ:$N$ 個の点それぞれについて $K$ 個の重心との距離を計算 → $O(NKD)$
  • 更新ステップ:各点を1回ずつ足し込むだけ → $O(ND)$

したがって全体で $O(NKD)$ per iteration、反復回数を $T$ とすれば $O(NKDT)$ です。$N$ に対して線形なので、階層的クラスタリングの $O(N^2)$ や $O(N^3)$ と比べて桁違いにスケールします。今回の実験でも、$N=600$、$K=5$ で平均6〜9回の反復で収束しました。

さらに大規模なデータには次の工夫があります。

  • Mini-batch k-means — 毎回全点を使わず、ランダムな小バッチで重心を更新する。数百万点でも実用的な時間で動きます。精度はわずかに落ちますが、多くの用途では許容範囲です。
  • 三角不等式による枝刈り(Elkanの方法) — 前回の距離と重心の移動量から「計算しなくても最近傍でないとわかる」重心をスキップする。厳密解を保ったまま高速化できます。
  • 次元削減との併用 — $D$ が大きいと距離が意味を失う「次元の呪い」が効いてきます。主成分分析などで数十次元まで落としてからk-meansをかけるのが定石です。

これらはいずれもk-meansの本質(交互最適化)は変えず、計算の仕方だけを工夫したものです。アルゴリズムの骨格が単純だからこそ、こうした高速化が素直に効きます。

まとめ

k-means法について、アルゴリズムの導出から失敗パターンまでを見てきました。要点を整理します。

アルゴリズムの本質

  • 目的関数はクラスタ内平方和 $J = \sum_n \sum_k r_{nk} \|\bm{x}_n-\bm{\mu}_k\|^2$
  • 「最近傍に割り当てる」も「平均に更新する」も、$J$ の最小化から導かれる。天下りの手続きではない
  • 各ステップが $J$ を減らし、$J$ が取りうる値は有限個なので、必ず有限回で停止する(実測でも56半ステップすべてで単調減少、増加0回)

弱点と対処

  • 初期値依存:ランダム初期化では最良解に届いたのは200試行中 42.5% だけ。k-means++で 69.5% に改善し、平均 $J$ も 13.9% 向上。さらに5回試行すれば実質100%
  • $K$ の決定:エルボー法とシルエット係数。今回は $K=3 \to 4$ で $J$ の減少率が72.7%から19.5%へ落ち込み、シルエットも $K=3$ で最大 0.714 と、両者が一致して正解を指した
  • 形状の制約:境界が必ず直線(ボロノイ分割)なので、同心円は純度0.513、細長い群は0.538と完全に失敗する
  • スケール依存:単位を cm から m に変えただけで純度が 0.870 → 0.597 に崩壊。標準化すれば単位に依存しなくなる(実測でラベル100%一致)

位置づけ

k-meansはGMMの共分散を $\varepsilon \bm{I}$ に固定して $\varepsilon \to 0$ とした極限です。実測でも $\varepsilon = 0.1$ で重心のずれが $1.4 \times 10^{-8}$ まで縮み、理論どおり一致しました。この見方に立つと、k-meansの制約(球状・等方・ハード割り当て)がどこから来ているのかが一目瞭然になり、GMM・スペクトラルクラスタリング・DBSCANといった他手法との関係も整理できます。

単純だからこそ速く、単純だからこそ限界がはっきりしている — それがk-meansです。まずk-meansを回して構造を掴み、その結果が前提(球状・同程度の大きさ・分離している)に合っているかを図で確認する。合っていなければ別の手法へ進む。この使い方ができれば、k-meansは今でも最初の一手として十分に強力です。

関連記事

画像なし
混合ガウスモデルとEMアルゴリズムを理解して実装する
k-meansをソフト割り当てへ一般化した混合ガウスモデル。EMアルゴリズムの導出と実装を解説します。
【機械学習】EMアルゴリズムをゼロから理解する
k-meansの交互最適化を一般化した枠組み。潜在変数を含むモデルの最尤推定を扱います。
画像なし
DBSCAN密度ベースクラスタリングの理論と実装
同心円や三日月のような任意形状のクラスタを扱える密度ベース手法。k-meansが失敗する場面での代替です。
画像なし
グラフラプラシアンとスペクトラルクラスタリングの理論と導出と実装
データをグラフとみなして分割する手法。非凸なクラスタ形状に対応できます。
【scikit-learn】データの標準化と正規化する方法とスケールを元に戻す方法
k-meansの前処理に必須のスケーリング。標準化と正規化の使い分けと実装を解説します。
画像なし
マハラノビス距離の幾何学的な意味を徹底的に理解する
特徴間の相関を考慮した距離。等方的なユークリッド距離の限界を超えるための道具です。
主成分分析(PCA)完全ガイド — 固有値分解からSVD・カーネルPCA・実装まで
k-meansの前処理として定番の次元削減。相関を落として次元の呪いを緩和します。