kmeans, breathing kmeans
今日はkmeansとその派生のbreathing kmeansについて解説します。クラスタリング手法の最も基本的なアルゴリズムの一つです。breathing kmeansはarXivで偶然投稿されているのを見つけましたが、実用性についてはかなり薄いのかなとおもいます。面白そうだとは思ったので、解説してみます。
kmeasn
k平均法とも呼ばれます。CARTなどのように分類・回帰を行うアルゴリズムではなく、何らかの基準を用いてデータをクラスタリングし、その結果を人間が解釈しやすくすることが目的だと思います。あるユーザーグループをクラスタリングし、その特徴（性年代、居住地など）を集計し、どのような施策を打つことが有効かなどの検討に用いることができます。ただ、実データでやってみると結構な頻度でなんだかあんまり特徴づけられないクラスターが出てくることがあり、お客様への説明にどうしようかなと悩むことも多いです。
クラスタ数を学習前に設定しなければなりませんが、これと決まった方法はありません。データの特徴がわからないからkmeansをやるわけなので当然です。そうはいってもいくつかの決定のための基準は開発されています。後ほど紹介します。
アルゴリズムとしては、下記を最小化する問題です。
L = \mathop{\rm argmin}\limits_{C_1,...,C_K} \sum_{i=1}^{n} \mathop{\rm min}\limits_{j} || x_i - C_j||^2
各サンプルをクラスタにアサインし、各サンプルと属するクラスタとの距離の自乗の合計が最も小さくなるようにすると言い換えることができます。
アルゴリズム詳細
最も重要なハイパーパラメーターはクラスタ数Kです。ほかはIteration回数などの主に学習が収束するまでの制御パラメーターです。
以下の手順で学習が行われます。
- クラスタの初期値選択
- 各サンプルの最近接クラスタを決定し、クラスタ番号を更新
- 各クラスタに属するサンプルの重心にクラスタ中心を更新
- ２と３を収束条件に抵触するまで繰り返し
クラスタの初期値選択
最も単純な方法はランダムに選ぶことです。64bit浮動小数点の表現可能な範囲で選んでしまうととんでもないことになるので、実際には各特徴量の最大値最小値を上下限として選ぶことになると思います。
import numpy as np
"""
N_CLUSTERS # クラスタ数
X # データ点。np.array
C # クラスタの座標
"""
N, D = 1000, 10
N_CLUSTERS = 5
X = np.random.uniform(low=-10.0, high=10.0, size=(N, D))
N, D = X.shape
C = np.random.uniform(low=0.0, high=1.0, size=(N_CLUSTERS, D))
X_MAX = X.max(axis=0)
X_MIN = X.min(axis=0)
C = (X_MAX - X_MIN) * C + X_MIN
しかし、これには問題があります。それは完全にランダムに選ばれるため、複数のクラスタ中心が近いところに選ばれる可能性があります。収束性が悪いだけでなく、本来は異なるクラスタに属すべきであるサンプルが同じクラスタに属してしまう可能性があることです。
%matplotlib inline
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
import numpy as np
# ダミーデータの作成
N = 1000
N_C = 5
C, L = make_blobs(n_samples=N,
centers=N_C,
cluster_std=1,
random_state=42)
plt.figure(figsize=(6,6))
plt.scatter(C[:, 0], C[:, 1], c=L)
plt.show()
上のスクリプトの実行結果は以下です。5個のクラスタに分けられるようにデータを生成しています。
このデータに対して下記のスクリプトを実行してみましょう。
kmeans = KMeans(n_clusters=5, init="random", n_init=1, random_state=2)
kmeans.fit(C)
plt.figure(figsize=(6,6))
plt.scatter(C[:, 0], C[:, 1], c=kmeans.labels_)
plt.show()
すると異なるクラスタに属してほしい左上の部分のデータが同じクラスタになってしまっています。さらに、左下は一つのクラスタに属してほしいデータが二つのクラスタに分かれてしまっています。かなり無理やりですが、このようなことが発生する懸念が残ってしまいます。
これを解決する手法としてkmeans++が開発されました。アイデアとしては非常に単純なもので、
- 各クラスタの初期値をなるべく遠くなるように選択する
だけです。最初の1個だけは完全ランダムですが、以降は既存クラスタ中心から離れるように選んでいきます。実装も非常に単純です。先ほどのmake_blobのデータで試してみます。初期値の段階でこの状態に近いと嬉しいわけです。
def roulette_selection(probas):
"""
ルーレット選択アルゴリズム。選択確率の配列（和が1に正規化されていなくともよい）から選択されたインデックスを返す。
Parameters
----------
probas : 1-dim numpy.array or list
それぞれのサンプルの選択確率が格納された配列。和が1に正規化されていなくともよい。
"""
pick = np.random.uniform(0, sum(probas))
current = 0
for i, prob in enumerate(probas):
current += prob
if current > pick:
return i
def compute_euclid_distance(X, Cluster):
"""
データ（N個のサンプル）とクラスタ一つとの間のユークリッド距離をすべて計算する
Parameters
----------
X : numpy.array, NxD matrix
データ点
Cluster : 1-dim numpy.array
クラスタ中心
"""
N, D = X.shape
dist = np.zeros(N)
for n in range(N):
dist[n] = np.linalg.norm(X[n,:]-Cluster)
return dist
import sys
proba = np.ones(N) / N # 選ばれる確率。最初はすべてが同じ確率で選ばれるように同じ値で初期化する。
Clusters = [] # クラスタ中心を格納するリスト。
np.random.seed(seed=141)
for i in range(N_C):
# 選んで格納
idx = roulette_selection(proba)
Clusters.append(C[idx,:])
# probaを更新する。最も近いクラスタ中心との距離を格納する。
proba = np.ones(N) * sys.maxsize # 初期化
for c in Clusters:
distance_to_c = compute_euclid_distance(C, c)
proba = np.minimum(proba, distance_to_c)
Clusters = np.array(Clusters)
plt.figure(figsize=(6,6))
plt.scatter(C[:, 0], C[:, 1], c=L)
plt.scatter(Clusters[:, 0], Clusters[:, 1], c="red",s=200)
plt.show()
上記のスクリプトではseedを固定してうまくいくような結果を出していますが、この方法でも毎回うまくいくわけではありません。下の図の赤い点がkmeans++で選択された初期値です。
また最適な解近いものが得られることが理論的に保障されているため、これ以外を初期値選択に選ぶ理由は実用的な範囲ではないのではないでしょうか。
各サンプルの最近接クラスタを決定し、クラスタ番号を更新
次の手順をみてみます。これは非常に単純で、まず全サンプルと全クラスタ間の距離行列を計算します。その後、各サンプルに最も近い位置にあるクラスタのラベルを割り振ります。とはいえ、これを愚直に実装してしまうと、非常に低速です。Kmeansの処理時間はほぼここにかかっているので、高速化することには意味があります。ここでは3つ紹介したいと思います。（Elkanの手法のようなもの以外で）ほかにもありましたらご教授ください。ユークリッド距離以外には触れません。
1. sqrtを外す
文字通りの変更です。Kmeansで（ほかのクラスタリングアルゴリズムでも？）重要なことは、どのサンプルがどのクラスタにより近いかを判定することです。順序関係さえわかればよいので、距離を正確に定義通り求める意味はありません。sqrtのあるなしで距離の順序に変更は起きません。
2. ユークリッド距離計算に展開した公式を用いる
D次元空間上の点 $x$ と $c$ のユークリッド距離 $d(p, q)$は定義から展開すると以下で与えられます。
\begin{align}
d^2(x,c)
&= \sum_{i=1}^{D}(x_i - c_i)^2 \\
&= \sum_{i=1}^{D} x_i^2 + \sum_{i=1}^{D} c_i^2- 2 \sum_{i=1}^{D} (x_i \times c_i)
\end{align}
ここで $x$ を全入力データ中のあるサンプル、$c$ をある一つのクラスタ中心とします。したがって、距離行列を作成するためには上記の $d^2$ の計算を合計で $N\times C$ 回行わなければなりません。1行目は、差分をD回、自乗をD回、和をD回の合計で3D回の計算を行います。したがっておおよそ $N \times C \times 3D$ 回の処理が行われます。これを展開を行うことで少し計算回数を削減することができます。
まず、サンプルデータは学習中に変化することがありません。したがって2行目第1項は一度だけ計算しておけばよいわけです（D回の自乗とD回の和をN回繰り返すので、$2ND$回の計算です）。そして第2項はクラスタ中心に対して同様の計算を行いますが、一般的にサンプル数はクラスタ数に比べて非常に多いためほとんど無視できる数です。第3項はD回の積とD回の和をN回繰り返すので、$2ND$回の計算です。そしてこれらを収束するまで何度も計算をすることになるので、結局第3項の計算が支配的になります。演算の回数が3ND回から2ND回に削減することができました。
3. 行列積サブルーチンを用いる
１も２もサンプルの一つとクラスとの一つをそれぞれペアにして計算するかのように記述をしました。しかし、それを解決する方法があります。行列積を利用する方法です。展開したユークリッド距離の第3項はつまり、サンプルとクラスタのベクトルの内積をとっています。サンプルの行列をクラスタの行列の転置したものの行列積をとればよいだけになります。
import numpy as np
from scipy.linalg.blas import sgemm
N = 1000
D = 10
K = 5
X = np.random.uniform(-10, 10, (N, D)).astype(np.float32)
C = np.random.uniform(-10, 10, (K, D)).astype(np.float32)
X_sq = np.sum(X*X, axis=1).reshape((N,1)) # 前処理で1回だけ算出しておくだけでよい
C_sq = np.sum(C*C, axis=1).reshape((1,K)) # イテレーションごとに計算する
Dist = np.zeros((N, K)) + X_sq + C_sq
Dist = Dist - sgemm(2.0, X, C, trans_b=True)
単純なユークリッド距離の計算を自分で実装してしまうとほぼ最適化されていないものができてしまうでしょう。もしかしたら上記の実装より良いものがあるのかもしれませんが、そこに時間をかけるとすさまじい労力となってしまいます。それであればOpenBlasやMKLなど非常によく最適化された既存のライブラリを用いたほうがよいでしょう（こんなライブラリを自分で実装していおいてなんですが）。
各クラスタに属するサンプルの重心にクラスタ中心を更新
次にクラスタ重心を更新します。単純に考えると、クラスターごとに属するサンプルの座標の和とカウントを取得、最後に平均をとることで更新できます。しかしこれも高速化の余地があります。各イテレーションでその前のイテレーションにおける各サンプルが属したクラスターを比較してみるとよくわかります。最初のイテレーションでは多くのサンプルが頻繁にクラスターを移動しますが、イテレーションが進むに従いほぼ更新がなくなります。更新されるのは主に二つのクラスターと同じような位置にあるサンプルのみです。したがって次のように改良できます。
- １回目のイテレーションは単純に和とカウントから重心を計算する。クラスターラベルとカウントは古いものとして配列に格納
- ２回目以降は、まず現在の重心にカウントを掛け戻す。以下はサンプルごとにループを行う
- 新しいクラスターラベルと古いものを比較してラベルが異なる場合は、古いクラスターから引き新しいへ足す。カウントも同様
breathing kmeans
次はbreathing kmeansです(github)。精度と速度が改善するという部分だけをみて面白そうだなと思って読んだのですが、かなり利用範囲が制限されそうです。以下のような点が適用に適しているようです。
- 各クラスターは同じような形で、クラスタ内で分散が少ない
- 各クラスターは同じようなサンプル数
- クラスターは分割性が良い
- クラスター数は既知
特に最後の部分が実用的でないかと思います。
アルゴリズム
通常のkmeansと異なり、breathingの名前は学習中にクラスター数が変化することから来ています。以下のような手順で進めます。
- 通常のkmeansを1回実施
- breathing in （k個クラスター追加）
- breathing out （k個クラスター削除）
- 収束条件にあたるまで２と３を繰り返す
２と３をそれぞれ説明していきます。
breathing in
このステップはクラスターを追加するステップです。ランダムに追加するのではなく、精度（RMSE）の悪いクラスター付近に追加します。この時、既存のクラスターはステップ１で学習されたものを使います。
breathing out
このステップではクラスターを削除します。こちらは精度の良いクラスターから削除します。しかし、単純にRMSEで昇順ソートして先頭から削除することは問題です。下図の左端がbreathing inが終わった段階の各クラスターのRMSEを値を記載した図（論文より取得）です。値の小さい部分から取り除いてしまうと右端の図のように右上のクラスターがまとめて削除されてしまいます。そうなると結果として精度が悪化します。
そこで、以下のような手順で取り除きます。
- RMSEでクラスターをソート
- 最もRMSEの小さい（すぐ近くに異なるクラスターが存在する）クラスターを削除
- 削除されたクラスターに最も近いクラスターは取り除かないように除外フラグを立てる
- 除外フラグに立っていない最もRMSEの小さいクラスターを削除
- k個取り除けるまで繰り返す
典型的な適用例
人工的なデータを作成しいろいろ実験を行っているのでGitHubを見ていただいた方が良いと思います。以下一部抜粋
以上