A Blog Entry on Bayesian Computation by an Applied Mathematician
$$
$$
より図が見やすい PDF 版は こちら .
導入
歴史
ハード -平均法はモデルフリーのクラスタリングアルゴリズムである.Voronoi 分割による競争学習の一形態とも見れる.
一方で原論文 (Lloyd, 1982 ) では,パルス符号変調 の文脈で,アナログ信号の量子化の方法として提案している.
実際, -平均法は(非可逆)データ圧縮にも用いられる.クラスター中心での画像の値と,それ以外では帰属先のクラスター番号のみを保存すれば良いというのである.このようなアプローチを ベクトル量子化 (vector quantization) という.
ソフト -平均法とは,このようなデータ点のクラスターへの一意な割り当てを,ソフトマックス関数を用いて軟化したアルゴリズムであり,多少アルゴリズムとしての振る舞いは改善するとされている.
最適化アルゴリズムとしての見方
個のデータ の クラスへの -平均クラスタリングアルゴリズムは,ハードとソフトの二種類存在するが,いずれも という損失関数の逐次最小化アルゴリズムとみなせる.
この見方は,EM アルゴリズム への一般化の軸となる.
しかしこの目的関数 が非凸関数であることが,第 5 節で示す実験結果で見る通り,アルゴリズムに強い初期値依存性をもたらす.
これをクラスタがより互いに離れるように,クラス割り当て法を修正することで対策した 最遠点クラスタリング (Gonzalez, 1985 ) が提案された.これは k-means++ (Arthur and Vassilvitskii, 2007 ) とも呼ばれており,データ圧縮法として用いた場合は復元誤差が でバウンド出来ることも示されている.
用いるデータ
実際のコードとデータを用いて -平均法を解説する.
まずは,解説にために作られた,次のような3つのクラスタからなる2次元のデータを考え,これの正しいクラスタリングを目指す.
ハード -平均法
アルゴリズムの説明
head -means algorithm はデータ とクラスタ数 ,そして初期クラスター中心 の3組をパラメータに持つ.
soft -means algorithm 3.1 はさらに硬度パラメータ を持つ.
numpy
の提供する行列積を利用して,これを Python により実装した例を以下に示す.ソフト -平均法の実装と対比できるように,負担率を通じた実装を意識した例である.
アノテーションを付してあるので,該当箇所(右端の丸囲み数字)をクリックすることで適宜解説が読めるようになっている.
def hkmeans_2d(data, K, init, max_iter= 100 ):
"""
2次元データに対するハード K-平均法の実装例.
Parameters:
- data: (N,2)-numpy.ndarray
- K: int クラスター数
- init: (2,K)-numpy.ndarray 初期値
Returns:
- clusters: (N,)-numpy.ndarray クラスター番号
"""
1 N = data.shape[0 ]
2 I = data.shape[1 ]
3 m = init
4 r = np.zeros((K, N), dtype= float )
for _ in range (max_iter):
# Assignment Step
for i in range (N):
5 distances = np.array([d(data[i], m[:,k]) for k in range (K)])
6 k_hat = np.argmin(distances)
7 r[:,i] = 0
r[k_hat,i] = 1
# Update Step
8 new_m = np.zeros_like(m, dtype= float )
9 numerator = np.dot(r, data)
10 denominator = np.sum (r, axis= 1 )
11 for k in range (K):
if denominator[k] > 0 :
new_m[:,k] = numerator[k] / denominator[k]
else :
new_m[:,k] = m[:,k]
12 if np.allclose(m, new_m):
break
m = new_m
13 return np.argmax(r, axis= 0 )
1
データ数を取得している.
2
データの次元を取得している.今回はすべて2次元データを用いる.
3
クラスター中心に引数として受け取った初期値を代入. -行列であることに注意.
4
負担率を -行列として格納している.その理由は後ほど行列積を通じた計算を行うためである.dtype=float
の理由は後述.
5
この distances
変数は (K,)-numpy.ndarray
になる.すなわち,第 成分が,第 クラスター中心との距離となっているようなベクトルである.ただし,d
は Euclid 距離を計算する関数として定義済みとした.
6
距離が最小となるクラスター番号 を, 番目のデータについて求める.
7
に基づいて負担率を更新するが,ループ内で前回の結果をリセットする必要があることに注意.
8
ここで dtype=float
と指定しないと,初め引数 init
が整数のみで構成されていた場合に,Python の自動型付機能が int
型だと判定し,クラスター中心 m
の値が整数に限られてしまう.すると,アルゴリズムがすぐに手頃な格子点に収束してしまう.
9
numpy
の行列積を計算する関数 np.dot
を使用している.更新式 の分子を行列積と見たのである.
10
分母 (denominator) は -行列 r
の行和として得られる.
11
ゼロによる除算が起こらないように場合わけをしている.
12
クラスター中心がもはや変わらない場合はアルゴリズムを終了する.
13
負担率の最も大きいクラスター番号を返す.今回は hat_k
の列をそのまま返せば良いが,soft -means アルゴリズムにも通じる形で実装した.
ただし,本記事の背後では次の実装を用いる.
クラスター中心の推移のヒストリーを保存して図示に利用したり,負担率 r
の中身を見たりすることが出来るようにするため,assignment step と update step とに分けてクラスメソッドとして実装し,run
メソッドでそれらを呼び出すようにしている.これに fetch_cluster
と fetch_history
メソッドを加えることで,クラスター番号とクラスター中心の推移を取得することが出来る.フィールド .r
から(最終的な)負担率を見ることもできる.
Code
class kmeans_2d:
"""
2次元データに対するソフト K-平均法の実装.
Usage:
kmeans = kmeans_2d(data, K, init, beta)
kmeans.run()
Parameters:
- data: (N,2)-numpy.ndarray
- K: int クラスター数
- init: (2,K)-numpy.ndarray 初期値
- beta: float 硬度パラメータ
"""
def __init__ (self , data, K, init, beta, max_iter= 100 ):
self .data = np.array(data, dtype= float )
self .K = K
self .init = np.array(init, dtype= float )
self .beta = float (beta)
self .max_iter = max_iter
self .N = data.shape[0 ] # データ数
self .I = data.shape[1 ] # 次元数 今回は2
self .m = init # クラスター中心の初期化.2×K行列.
self .r = np.zeros((K, self .N), dtype= float ) # 負担率.K×N行列.
self .history = [init.copy()] # クラスター中心の履歴.2×K行列.
def soft_assigment(self ):
"""soft K-means の場合の負担率の更新"""
for i in range (self .N):
distances = np.array([d(self .data[i], self .m[:,j]) ** 2 for j in range (self .K)]) # (N,)-numpy.ndarray
denominator_ = np.sum (np.exp(- self .beta * distances)) # 分母
self .r[:,i] = np.exp(- self .beta * distances) / denominator_
def hard_assigment(self ):
"""hard K-means の場合の負担率の更新"""
for i in range (self .N):
distances = np.array([d(self .data[i], self .m[:,j]) for j in range (self .K)]) # (N,)-numpy.ndarray
k_hat = np.argmin(distances) # 最小距離のクラスター番号
self .r[:,i] = 0 # 前のループの結果をリセット
self .r[k_hat,i] = 1
def update(self ):
"""クラスター中心の更新"""
new_m = np.zeros_like(self .m, dtype= float ) # ここで float にしないと,クラスター中心が整数に限られてしまう.
numerator = np.dot(self .r, self .data) # (K,2)-numpy.ndarray
denominator = np.sum (self .r, axis= 1 ) # 各クラスターの負担率の和
for k in range (self .K):
if denominator[k] > 0 :
new_m[:,k] = numerator[k] / denominator[k]
else :
new_m[:,k] = self .m[:,k]
self .m = new_m
def fetch_cluster(self ):
"""最終的なクラスター番号を格納した (N,)-array を返す"""
return np.argmax(self .r, axis= 0 )
def fetch_history(self ):
"""クラスター中心の履歴を格納したリストを,3次元の np.array に変換して返す"""
return np.stack(self .history, axis= 0 )
def run_soft(self ):
"""soft K-means アルゴリズムの実行"""
for _ in range (self .max_iter):
self .soft_assigment()
self .update()
self .history.append(self .m.copy())
if np.allclose(self .history[- 1 ], self .history[- 2 ]):
break
def run_hard(self ):
"""hard K-means アルゴリズムの実行"""
for _ in range (self .max_iter):
self .hard_assigment()
self .update()
self .history.append(self .m.copy())
if np.allclose(self .history[- 1 ], self .history[- 2 ]):
break
なお,この実装は などの場合にオーバーフローが起こることに注意.これへの対処は logsumexp
の使用などが考えられる.
初期値依存性
次の2つの初期値を与えてみる. と, は変えずに の -座標を だけ下げたもの とを初期値として与えてみる.
正解数: 51 正解率: 56.7 % 反復数: 9 回
別の初期値を与えてみる(右下の点 を だけ下に下げただけ):
正解数: 85 正解率: 94.4 % 反復数: 7 回
結果が全く変わり, を与えた方が,大きく正解に近づいている.具体的には,右下の初期値 は右上の島に行くが, は左下の島に行ってくれる.
ハード -平均アルゴリズムは初期値に敏感である ことがよく分かる.
局所解への収束
直前の結果2 ではクラスター2と3の境界線で4つのミスを犯しており,これを修正できないか試したい.
そこで,答えに近いように, を初期値として与えてみて,正答率の変化を観察する.
正解数: 85 正解率: 94.4 % 反復数: 5 回
もはや初期値から殆ど動いていないが,目標のクラスター3に分類された3つの点が,相変わらず3のままであり,加えてクラスター2の中心がこれらから逃げているようにも見えるので,クラスター2の初期値をよりクラスター3に近いように誘導し,クラスター3の中心をより右側から開始する:
正解数: 85 正解率: 94.4 % 反復数: 6 回
こんなに誘導をしても,正しく分類してくれない.
実は,以上2つの初期値では,最終的に3つのクラスター中心は同じ値に収束している.よって,これ以上どのように初期値を変更しても,正答率は上がらないシナリオが考えられる.
以上の観察から,ハード -平均法はある種の 局所解に収束する ようなアルゴリズムであると考えられる.
ソフト -平均法
アルゴリズムの説明
ハード -平均法2 では,負担率 は のいずれかの値しか取らなかった.この振る舞いを, で定まる ソフトマックス関数 を用いて,「軟化」する.
ここでは, として, の形で用い, の代わりに とする.ただし, は 上の Euclid 距離とした.
硬度パラメータ
は 硬度 (stiffness) または逆温度と呼ぶ. は距離の次元を持つ.
のときは温度が無限大の場合にあたり,常に負担率は一様になる.絶対零度に当たる の極限が hard -means アルゴリズムに相当する.
逆温度 を連続的に変化させることで,クラスタ数に分岐が起こる,ある種の相転移現象を見ることができる.
実装
実装は例えば hard -means アルゴリズム2 から,負担率計算の部分のみを変更すれば良い:
for i in range (N):
1 distances = np.array([d(data[i], m[:,k]) for k in range (K)])
2 denominator_ = np.sum (np.exp(- beta * distances))
3 r[:,i] = np.exp(- beta * distances) / denominator_
1
データ とクラスター中心 との距離を計算し,ベクトル を distances
に格納している.
2
負担率の計算 を2段階に分けて行なっており,分母を先に計算して変数 denominator_
に格納している.
3
すでに計算してある分母 denominator_
を用いてデータ の負担率 を計算し, -行列 r
の各列に格納している.
挙動の変化の観察
逆温度をはじめに としてみる.図 1 と全く同様な初期値 を与えてみると,次の通りの結果を得る:
正解数: 44 vs. 51 正解率: 48.9 % vs. 56.7 % 反復数: 28 回 vs. 9 回
クラスターの境界が変化しており,正解率は悪化している.さらに,反復数が9回であったところから,3倍に増えている(28回).
また,右上の2つのクラスター中心の収束先は,微妙にずれているが ほとんど一致している 点も注目に値する.
centers = history[- 1 , :, :]
df = pd.DataFrame(centers, columns= ['Cluster1' , 'Cluster2' , 'Cluster3' ])
print (df)
Cluster1 Cluster2 Cluster3
x 2.397456 2.397535 -0.036071
y 2.047565 2.047580 -1.448288
図 2 で与えた初期値 も与えてみる.
正解数: 85 vs. 85 正解率: 94.4 % vs. 94.4 % 反復数: 59 回 vs. 7 回
クラスター境界と正答率は変わらないが,反復数がやはり7回から大きく増えている.
結果はやはり 図 3 とは大きく異なっており,ハード -平均法で観察された初期値鋭敏性が,変わらず残っている.
加えてこの場合も 図 3 のクラスター1と2と同様に,クラスター2と3の中心がほぼ一致している.
Cluster1 Cluster2 Cluster3
x 2.466833 -0.369537 0.447958
y 2.124961 -1.076874 -1.543758
の場合のソフト -平均法は,この例では クラスター中心が融合する傾向にある ようである.
一般に, が小さく,温度が大きいほど,エネルギーランドスケープに極小点が少なくなり,クラスターは同じ場所へ収束しやすくなると予想される.
高温になるほどクラスター数は減少する
初期値を直前で用いた で固定とし,さらに温度を上げて,逆温度を としてみる.
正解数: 68 vs. 85 正解率: 75.6 % vs. 94.4 % 反復数: 101 回 vs. 59 回
反復数はさらに増加し,全てがほとんど同じクラスターに属する結果となってしまった.
Cluster1 Cluster2 Cluster3
x 1.715903 0.862511 1.329066
y 1.012398 -0.099845 0.511186
温度が大変に高い状態では,全てが乱雑で,3つのクラスターが一様・公平に負担率を持つようになった.そのため,第一歩からほとんど全体の中心へと移動し,反復数が減る.
次に,温度を少し下げて,逆温度を としてみる.
正解数: 85 vs. 85 正解率: 94.4 % vs. 94.4 % 反復数: 17 回 vs. 59 回
初めて soft -means アルゴリズムを用いた場合で,3つのクラスター中心がはっきりと別れた.反復回数は, の場合と比べればやはり落ち着いている.
しかし,正解率は head -means の場合( 図 2 など)と全く同じである.実は,最終的なクラスター中心も 図 2 の最終的なクラスター中心とほとんど同じになっている.
今回のソフト -平均法の最終的なクラスター中心
Cluster1 Cluster2 Cluster3
x 2.416113 0.881629 -1.338782
y 2.086327 -1.934090 -0.816316
図 2 のハード -平均法の最終的なクラスター中心
Cluster1 Cluster2 Cluster3
x 2.426102 0.868333 -1.323353
y 2.091429 -1.948458 -0.765176
以上より,ソフト -平均法は温度を上げるほどクラスター数が少なくなり,温度を下げるほどクラスター数は上がり,十分に温度を下げるとハード -平均法に挙動が似通う .
最適な硬度の選択
ではクラスターが2つに縮退し, では hard -means アルゴリズムの結果とほとんど変わらなくなる.その中間では次のように挙動が変わる:
正解数: 50 vs. 86
正解率: 55.6 % vs. 95.6 %
正解数: 85 vs. 85
正解率: 94.4 % vs. 94.4 %
やはり,温度が高い場合はクラスター中心が合流・融合してしまいやすいが,冷却することでクラスター数は大きい状態で安定する,と言えるだろう.
本番データセットでの実験
今まで使っていたデータ1.3 はクラスターのオーバーラップはなかったため,いわば優しいデータであった.ここからはよりデータ生成過程が複雑なデータを用いて,ソフト -平均法の挙動を観察する.
データの概観
今度は,次の4クラスのデータを用いる.
実は,これは4つの Gauss 分布から生成されたデータである.
最適な温度の選択
正解数: 377 vs. 366 正解率: 83.8 % vs. 81.3 % 反復数: 49 回 vs. 62 回
正解数: 386 vs. 375 正解率: 85.8 % vs. 83.3 % 反復数: 101 回 vs. 70 回
正解数: 378 vs. 378 正解率: 84.0 % vs. 84.0 % 反復数: 39 回 vs. 14 回
実験結果まとめ
データ 1.3 に対して,(初期値 で)ソフト -平均法を適用すると,
の場合で結果はハード -平均法と変わらなくなる.
の場合で結果はクラスターがほとんど2つになり, では計算機上では実際に2つになってしまう.
正答率は で最大であった.
を大きくするほど,反復回数は減少していった.
データ 4.1 に対しても,以上の4点について同様の傾向が確認できた.
こうしてソフト -平均法とハード -平均法の性質は分かった.主に
初期値依存性
クラスタ数 の選択法
の問題が未解決であり,恣意性が残る.
これを,ベイズ混合モデリングにより解決する方法を次稿で紹介する:
References
Arthur, D., and Vassilvitskii, S. (2007). K-means++: The advantages of careful seeding. In Proceedings of the eighteenth annual ACM-SIAM symposium on discrete algorithms , pages 1027–1035. USA: Society for Industrial; Applied Mathematics.
Bishop, C. M. (2006).
Pattern recognition and machine learning . Springer New York.
Gonzalez, T. F. (1985).
Clustering to minimize the maximum intercluster distance .
Theoretical Computer Science ,
38 , 293–306.
Lloyd, S. (1982).
Least squares quantization in PCM .
IEEE Transactions on Information Theory ,
28 (2), 129–137.
MacKay, D. J. C. (2003).
Information theory, inference and learning algorithms . Cambridge University Press.
MacQueen, J. (1967).
Some methods for classification and analysis of multivariate observations . In
Proceedings of the fifth berkeley symposium on mathematical statistics and probability ,Vol. 1, pages 281–297.