アナタ ハ k-Means ヲ スキ デス カ？
ハイ ワタシ ハ k-Means ガ ダイ スキ デス
概要
- Subspace k-Means を実装したよ
- Subspace k-Means と普通の k-Means の違いを説明するよ
- Subspace（に限らず）k-Means での最適クラスタ数の見積もりをしたよ
実装したよ
ということで、既に実装をしております。
なるべく普通の k-Means (sklearn.cluster.KMeans) と同じように動くように作ってます。
（2021.04.21 時点で scikit-learn == 0.24.1 で動くことは確認してます）
違う点は以下。
- Sparse Matrix (scipy.sparse.csr_matrix) には対応してません
- scikit-learn が採用しているアルゴリズム（lloyd と elkan）には対応してないし、そもそも Cython を全く使ってなく、numpy だけで書いてます
- なので、分かりやすいけど、遅いです
Subspace k-Means とは
そもそも Subspace k-Means って何だって話ですが、元ネタは以下の論文です。
ぶっちゃけ Google Scholar でググれば PDF も出てきます。
Mautz, Dominik, et al. "Towards an Optimal Subspace for K-Means." Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. ACM, 2017. (link)
- データを、クラスタリングに意味のある「クラスタ空間」と、クラスタリングに意味のない雑音である「ノイズ空間」の２つの部分空間（= Subspace）に分ける
- その「クラスタ空間」上で k-Means を行う
単純です。単純だからこそ素晴らしい。
じゃあどうやって部分空間に分けて、どうやって k-Means するかですが、めっちゃ詳しく論文に載ってます。
ここまで載ってたら実装できないわけがない。
ということで、これを物凄く忠実に実装したのが、上の repository になります。
（論文の著者は "SubKmeans" って名付けてるけど、"Sub" じゃあ意味が分からないので、"Subspace" って付けてます。そこが一番論文に忠実じゃないところ（笑））
k-Means と Subspace k-Means の違い
私が説明するのも烏滸がましいですが、k-Means とは以下のようなアルゴリズム（Expectation-Maximization algorithm: EM-algorithm）です。
- （初期化）データから適当に中心を $k$ 個（クラスタ数）決める
- （E-step）各データと $k$ 個の中心の距離を計算し、最も近い中心のクラスタにそのデータが属するとして、すべてのデータを各クラスタに割り振る
- （M-step）各クラスタに属するデータから、そのクラスタの中心を再計算する
- 収束するまで E-step と M-step を繰り返す
単純。素敵。抱いて。
じゃあ Subspace k-Means は、これと比較して何をしているのかというと、
- （初期化）データから適当に中心を $k$ 個（クラスタ数）決める
- （初期化２）仮決めのクラスタ次元 $m\ (0<m\leq k)$ と、$k$ 次元から $m$ 次元だけ取り出す単位行列 $P_C$、そして $k$ 次元のデータを $k$ 次元に射影する適当な $k \times k$ 直行行列 $V$ を決める
- （E-step）データと中心を、$V$ かまして $P_C$ かますことによって、射影して $m$ 次元だけ取り出して、その空間上で 各データと $k$ 個の中心の距離を計算し、最も近い中心のクラスタにそのデータが属するとして、すべてのデータを各クラスタに割り振る
- （M-step）各クラスタに属するデータから、そのクラスタの中心を 元の空間で 再計算する
- （M-step２）$m$（とそれに応じた $P_C$）と $V$ を更新する
- 収束するまで E-step と M-step を繰り返す
以上。単純。濡れた。
Subspace k-Means は何が嬉しいのか
「データをクラスタリングに意味のあるものとノイズに分解する」 の一言に尽きる。
- k-Means はデータ全部を使っちゃうので、ノイズが混ざっているとそれだけクラスタリングもごちゃごちゃになる
- Subspace k-Means は、意味のある（都合の良い）データだけを取り出してクラスタリングするので、k-Means では綺麗にクラスタリング出来ないデータもクラスタリング出来る
- 見た目も綺麗！！
なぜそんな事が出来るのか
詳しくは論文を読んで欲しいのだが、キモは（当たり前だが）直交行列 $V$ である。
k-Means の損失関数は以下で与えられる。（ちなみにこの値は fit 後の KMeans.inertia_ に格納されている）
$$
\sum_{i=1}^{k}\sum_{\mathbf{x}\in C_i}||\mathbf{x}-\mathbf{\mu}_i||^2
$$
これを最小化する $V$ は以下を最小化する $V$ である。（行列のトレースのテクニックだけを使って導出できる！）
$$
\mathrm{Tr}\left(P_{C}P_{C}^{T}V^{T}\left(\left[\sum_{i=1}^{k}\sum_{\mathbf{x}\in C_i}(\mathbf{x}-\mu_i)(\mathbf{x}-\mu_i)^{T}\right]-\sum_{\mathbf{x}\in D}(\mathbf{x}-\mu)(\mathbf{x}-\mu)^{T}\right)V\right) \
= \mathrm{Tr}\left(P_{C}P_{C}^{T}V^{T}\left(\left[\sum_{i=1}^{k}S_{i}\right]-S_{D}\right)V\right) \
= \mathrm{Tr}\left(P_{C}P_{C}^{T}V^{T}\Sigma V\right)
$$
上の式で、$\mathbf{x}\in C_{i}$ はクラスタ $C_i$ に属するデータ、$D$ はデータ全体、$\mu_i$ はクラスタ $C_{i}$ の中心、$\mu$ はデータ全体の中心を表す。
この中で、$P_{C}P_{C}^{T}$ は最初の $m$ 次元だけ取り出す操作。とすると $\mathrm{Tr}(V^{T}\Sigma V)$ の最初の $m$ 個の和を最小化すれば良いことになる。
もうこの時点で明らかなのだが、$\Sigma=\left[\sum_{i=1}^{k}S_{i}\right]-S_{D}$ って、$S_i$ が $\mu_i$ を中心とした共分散行列、$S_D$ が全体の共分散行列なんだから、$\Sigma$ は対称行列に決まっている。$V^{T}\Sigma V$ は、$V$ が直交行列（$V^T=V^{-1}$）なんだから、$\Sigma$ の対角化に他ならない。
すると $\mathrm{Tr}(V^{T}\Sigma V)$ は $\Sigma$ の固有値の和、要するに $\mathrm{Tr}\left(P_{C}P_{C}^{T}V^{T}\Sigma V\right)$ を最小化する $V$ は、$\Sigma$ の固有ベクトルを固有値がちっちゃい順に並べたヤツだ！
データの共分散の対角化と言うと、まず思い浮かぶのがPCA（Principal Component Analysis:主成分分析）。あれは固有ベクトルを固有値が大きい順に並べることで、データをクルクル回転させて、分散が大きい順に、軸が直行するように変換する方法。分散が小さい軸を削除することで次元削減が出来る。
Subspace k-Means は $\Sigma$ の固有値が小さい軸、もっと言うと、$m$ として $\Sigma$ の固有値が $0$ より十分小さい軸（の数）だけ取ることによって、クラスタ中心に収束する（都合の良い）データに変換して（都合の良い数だけ）取り出し k-Means する、という方法になる。
（固有値が $0$ もしくは $0$ より大きい軸は、発散するデータなんで、ノイズ！）
k-Means の EM は、損失関数が最小となるように期待して最大効率化するアルゴリズムだけど、Supspace k-Means は同時にデータの収束が最も良くなるように逆PCAで都合よく変換して都合の良い所だけ取り出すという、一挙両得のアルゴリズムなのだ！
wine でやってみる
実際に wine のデータでやってみた例（ipynb）を先程の repository に入れてるので、見て欲しい。
Subspace k-Means との比較として、PCA してから k-Means するやり方をまずは載せている。
PCA で次元削減する場合、累積寄与率を計算して、例えば累積寄与率が 0.8 を超える最小の次元数に削減する。例えばこの wine の例だと 5 次元となる。
で、 この 5 次元データを k-Means するわけだが、あらかじめワインは３種類って分かっているので、クラスタ数 3 で k-Means する。
この結果を図示したいのだが、5 次元を平面に射影するために、t-SNE を使っている。その結果がコレ。
悪くはないが、１直線に並んでしまっているので、クラスタ数が分かんなかった場合にグチャグチャになりそうな予感がする。
これが Subspace k-Means だとこうなる。
３つの島が綺麗に現れるので、クラスタ数が分からなくても、きちんと分類できそう！
ちなみに上の図は Subspace のクラスタ空間の方だが、ノイズ空間はこうなる。
本当にホワイトノイズである。
この分類精度を評価するのに、「どのクラスタがどの Ground Truth と対応して」と考えるのはめんどいので、NMI (Normalized Mutual Info Score) と AMI (Adjusted Mutual Info score) を使う。
全く同じスコアである。Subspace k-Means は、k-Means に比べて、少なくともおかしなことはしてない（俺が実装したのに！）ということが分かる。
pendigit でやってみる
流石に wine は簡単すぎたので、もうちょっと難しい pendigit でやってみる。
「手書き数字の分類」なんでクラスタ数は 10 に決まっているのだが、それを分からない想定で、最適クラスタ数の推定もやってみることにする。
まず上記の UCI のサイトのデータは .tra とか .tes とか見たこと無い拡張子なんでびっくりするけど、中身は単なる CSV なんで、落ち着いて対応しよう。（俺はスゲーあせった）
で、クラスタ数推定のために、以下のパターンを試した。
- PCAした場合とPCAしてない場合で、それぞれ k-Means と Subspace k-Means の４通り
- クラスタ数を見積もる指標として、以下の３通り
- inertia (distortion) の elbow point（肘点）
- silhouette（シルエット） score の knee point（膝点）
- calinski-harabasz score のknee point
elbow point, knee point の計算は、実装がめんどくさいので、kneed を使っている。が、kneed は（本来不要な）MinMaxScaling をしていて気持ち悪かったので、KneeLocator を継承し正規化をやらない ElbowLocator を作って計算している。（ここはどうでも良い所）
ちなみに elbow point は凹な関数の一番凹な所、knee point は凸な関数の一番凸な所、という理解で十分。（適当すぎる説明）
クラスタ数の推定
まずは所要時間。
やっぱり Cython 化してないし計算が複雑な Subspace k-Means が遅い。まぁこれは致し方ない。
次に収束するまでの iteration 数。
一概に Subspace k-Means が多いというわけではなく、k-Means の方が多い場合もある。ただ単に育ってきた環境が違う（セロリ）（収束までのルートが異なる）だけだろう。
そしていよいよ、最もメジャーな inertia vs クラスタ数。
丸点が elbow point。クラスタ数が増えるほど基本的に inertia は減るので、漸減する曲線になる。
PCA 有りの方は、元の 16 次元を 5 次元に圧縮しているので、inertia は当然小さい。
我ながらビックリしたのが、k-Means と Subspace k-Means は収束後はまったく同じ inertia の値を示す。wine で全く同じ結果になったわけだ。
で、偶然だが、PCA 無しの方がぴったり正解を当てている。
（正規化があっても無くても） elbow/knee point の位置はクラスタ数の上限・下限の値でブレてしまうので、まさに偶然。
ただ、inertia の elbow point が、クラスタ数の良い指標となることは示せたのだろう。
最後に silhouette score と calinski-harabasz score vs クラスタ数。丸点が knee point。
シルエットスコア（フランス語由来の単語書くのがめんどくさくなった）は基本的にガタガタになるので、まぁこれは想定内。
calinski-harabasz スコアは、筋が良いデータならばピョンと良い所で飛び出す筈なのだが、あんまり上手くいかなかった。こういうこともあるだろう。（適当すぎる感想）
ということで、PCA で次元削減すると情報量が減るので、クラスタ数の推定には向かない。シルエットスコアや calinski-harabasz スコアだと微妙なケースがある。よっていつもの inertia の elbow point で見事正解の 10 が当たったね、ということで、クラスタリングの精度の評価に移る。
クラスタリングの精度
まずは、全部のデータを t-SNE したものと Subspace k-Means のクラスタ空間（を t-SNE したもの）を、Ground Truth で色分けしたデータで比較する。
両者とも、遜色なく島ごとに色分けされている。
一部ごちゃっとなっているが、「１と３」「１と６」が混ざっているところなんかは、「字が汚くて、どっちだか分かんないの、あるよねー」と逆に納得感がある。
（t-SNE で綺麗に分かれるんだったら、最初から t-SNE や、minCut を使う Spectral Clustering しろよ、という話もあるが、どちらも fit & predict は出来るけど、fit した以外のデータに対する predict は出来ないのよ！）
（Spectral Clustering に関しては、また別途記事を書きたいと思っている）
で、クラスタリング結果だが、今回は inertia が同じ値ということもあって、ただの k-Means と Subspace k-Means を比べる。
まぁ、両者とも k-Means の手法的限界で、島の途中で分断されちゃってる所もある。
NMI, AMI を比較してみる。
Subspace k-Means の方が良い！！！
Subspace k-Means はクラスタ空間で k-Means しているだけあって、同じ local minimum でも、より良い local minimum に落ち込むことが出来たのだろう。濡れた！抱いて！先っぽ（略
ちなみにノイズ空間は安定のホワイトノイズである。
最後に一言
という訳で、Subspace k-Means は本当に素晴らしいので、是非私のなんちゃって実装なんかじゃなくて、ちゃんと Cython とかで高速化した上で scikit-learn に取り込まれて欲しいし、scikit-learn はバージョンが上がる毎に内部が結構変わるので、キャッチアップするのがめんどい。（次のリリースで変更しなきゃいけないのも既に分かっている）
誰にお願いすれば良いんだろう。
（蛇足）
よく分かってないが、私はこの repository で実績を解除したらしい。
北極海の海底に眠るストレージに保管される？ やったー（ミーハーなので、とりあえず喜ぶ）