雑な要約

　本記事では，K-meansクラスタリングのアルゴリズムの紹介と，アルゴリズムの各ステップで目的関数が単調減少することの粗い証明を書いています．また，pythonでK-meansクラスタリングのアルゴリズムをスクラッチ(?)で実装しています．

K-meansクラスタリングとは

　K-meansクラスタリングとは，教師なし学習(正解ラベルなしで学習させる)の一種で，非階層的クラスタリングの手法の最も有名なものの１つです．主な目的は，データの構造からグループを作ることです．例えば，以下のようなデータの散布図(注:極端な例です)が得られたとします．

パッと見た感じ，このデータは２つの集団から構成されるデータだと推測できます．このデータに対してpythonのライブラリsklearnのKMeans()を適用して，クラスのラベルを推定し，ラベルごとに色分けした結果が以下のようになります．

以上のように，K-meansクラスタリングはデータのグループを推定する手法，もしくは，クラスラベルを推定する手法とみなせます．

各種記号と目的関数

　以下の節で利用する記号と目的関数を先にまとめておきます．

X: データ行列($n\times p$)
g: クラスター(グループ)の数
E: メンバーシップ行列($n\times g$)．i番目の個体がどのグループに属するかを表す行列

\begin{eqnarray}
e_{ij} = \left\{
\begin{array}
1 1&\ if x_i \in G_j\\
0 &otherwize
\end{array}
\right.
\end{eqnarray}

C: セントロイド行列($g\times p$)．各グループの中心を格納する行列
$||A||_F^2=trA^TA=trAA^T$: Frobeniusノルム

　K-meansクラスタリングのモデル式は，データの各点が各クラスター中心とメンバーシップ行列で表されるという考え方から，$$X = EC + error$$と表されるため，目的関数は以下のようになります．

f(E,C)=||error||_F^2=||X-EC||_F^2

アルゴリズム

　以下K-meansクラスタリングのアルゴリズムです．

step1. クラスター中心を初期化する．
step2. i番目のデータ点と各クラスター中心の距離を計算し，最も近いクラスター中心に対応するメンバーシップ行列の列に1を代入する．それ以外の列には0を代入する．これを全てのiについて行う($i=1,2,\cdots,n$)．
step3. 新しく計算されたメンバーシップ行列を用いてクラスター中心を計算する．
step4. 収束するまでstep2,3を反復する．

　K-meansクラスタリングのアルゴリズムは以上のように，メンバーシップ行列Eとセントロイド行列Cを交互に更新していくアルゴリズムです．step2,3の反復で目的関数が単調減少することが知られています．

単調減少することの証明

　各ステップで目的関数が減少することを示せばokです．まずは，step2(メンバーシップ行列の更新)に関してです．目的関数を要素ごとに書き換えてみると，

\begin{eqnarray}
f(E,C) &=& ||X-EC||_F^2\\
&=& tr(X-EC)^T(X-EC)\\
&=& \sum_{i=1}^n\sum_{j=1}^p(x_{ij}-\sum_{k=1}^ge_{ik}c_{kj})^2\\
&=& \sum_{i=1}^n\sum_{k=1}^ge_{ij}\sum_{j=1}^p(x_{ij}-c_{kj})^2
\end{eqnarray}

となります．以上の式展開から目的関数を最小化したい場合，データ点$x_i$を最も近いクラスター中心のクラスターに割り当てれば良いということがわかります．よってstep2の操作で目的関数は減少します．ちなみにstep2での更新は，以下のようにも書き換えられます．

\begin{eqnarray}
e_{ik} = \left\{
\begin{array}
11\ &if\ \sum_{j=1}^p(x_{ij}-c_{kj})^2\le\sum_{j=1}^p(x_{ij}-c_{lj})^2\ (l=1,\cdots,g,l\ne k)\\
0& otherwise
\end{array}
\right.
\end{eqnarray}

　次にstep3についてです．step3では，step2で得られたメンバーシップ行列(データ点がどのクラスターに所属するかの情報を持つ行列)を用いて各クラスターの重心を計算しています．目的関数をセントロイド行列Cに関して微分すると，

\begin{eqnarray}
\frac{\partial}{\partial C}f(C) &=&\frac{\partial}{\partial C}tr(X-EC)^T(X-EC)\\
&=&\frac{\partial}{\partial C}tr(X^TX-2X^TEC + C^TE^TEC)\\
&=&-2E^TX+2(E^TE)C
\end{eqnarray}

これを$\bf{O}$解くことで，メンバーシップ行列Eを固定した状態の目的関数を最小にするセントロイド行列Cは$$C = (E^TE)^{-1}E^TX$$と表すことができます．ここで，$E^TE$は各クラスターに属するデータ点の数を対角成分に持つ対角行列で，$(E^TE)^{-1}$はそれらの逆数を対角成分に持つ対角行列です．$E^TX$はその$(k,j)$成分($:=ex_{kj}$とする)が，

ex_{kj} = \sum_{i=1}^ne_{ik}x_{ij}

つまり，k番目のクラスターに属する個体の総和を表す行列を表しています．以上の話から$C=(E^TE)^{-1}E^TX$の$(k,j)$成分は$$c_{kj} = \frac{1}{n_k}\sum_{i=1}^ne_{ik}x_{ij}$$(=重心の計算式!)となり，各クラスターの重心を計算することで目的関数が減少することがわかりました．
　以上から，step2,3を反復することで目的関数は単調減少します．

申し訳コード(スクラッチでK-means実装しました)

相変わらず汚い．

# -*- coding: utf8 -*-

import numpy as np 
import random
import matplotlib.pyplot as plt 

from sklearn.cluster import KMeans 

# 例で用いた散布図
sigma1 = np.array([[1,2],[2,1.7]]); mu1 = np.array([-1,-1]); mu2 = np.array([10,10])

np.random.seed(222)#今日2/22だから
ran1 = np.random.multivariate_normal(mu1,sigma1,100)
np.random.seed(222)
ran2 = np.random.multivariate_normal(mu2,sigma1,100)
ranran = np.vstack([ran1,ran2])[random.sample(list(np.arange(200)),200)]#インデックスごちゃ混ぜにする
## 色分け前
plt.scatter(ran1[:,0],ran1[:,1],color="gray")
plt.scatter(ran2[:,0],ran2[:,1],color="gray")
# plt.savefig("適当なディレクトリ")
plt.show()

## K-meansやってみる(sklearn.clusterのKMeansを使う)
model_km = KMeans(n_clusters=2)
pred_clus = model_km.fit_predict(ranran) #分類結果のクラス

# 色分け後
for i in range(2):
	label = ranran[pred_clus == i]
	plt.scatter(label[:,0],label[:,1])
# plt.savefig("適当なディレクトリ")
plt.show()






# K-means(スクラッチ)
def kmeans(data, n_clus, max_iter=10,tol=10**(-6),nstart=10):
	X = data.astype("float64"); n,p = X.shape;loss_best = np.Inf;loss_values = []
	count = 0;E_hat = np.zeros((n,n_clus));C_hat = np.zeros((n_clus,p))
	while count < nstart:
		print("------nstart: %d -------" % count)
		count += 1
		#クラスター中心の初期値の発生
		idx = random.sample(list(np.arange(n)),n_clus)
		C = X[idx,:]

		inner_cnt = 0
		loss = np.Inf;losses = []
		while True:
			print("inner count %d: %f" % (inner_cnt,loss))
			inner_cnt += 1

			#メンバーシップ行列Eの更新
			E = np.zeros((n,n_clus))
			for i in range(n):
				dist = []
				for c in range(n_clus):
					dist.append(np.sum((X[i,:]-C[c,:])**2))#各セントロイドとの距離を計算
				mindex = np.where(dist==np.min(dist))
				E[i,mindex] = 1#一番近いセントロイドのクラスターに1を代入

			#セントロイド行列Cの更新(各クラスの重心計算)
			D = np.dot(E.T,E);ns = np.diag(D)#各クラスターに属する点の個数の配列
			C = (np.dot(X.T,E)/ns).T

			new_loss = np.trace(np.dot((X-np.dot(E,C)).T,(X-np.dot(E,C))))#目的関数値の計算
			losses.append(new_loss)
			diff = loss - new_loss
			if diff < tol:
				break
			if inner_cnt > max_iter:
				print("Not converge!!!")
				break
			loss = new_loss
		if new_loss < loss_best:
			C_hat = C; E_hat = E;
			loss_best = new_loss; loss_values = losses

	return C_hat, E_hat, loss_values

if __name__ == "__main__":
	sigma2 = np.array([[2,3.7],[3.7,-7]]);mu3 = np.array([13,-20])
	#各分布からデータ点生成
	np.random.seed(1);g1 = np.random.multivariate_normal(mu1,sigma1,100)
	np.random.seed(2);g2 = np.random.multivariate_normal(mu2,sigma1,100)
	np.random.seed(3);g3 = np.random.multivariate_normal(mu3,sigma2,100)
	#色分け前
	plt.scatter(g1[:,0],g1[:,1],color="black")
	plt.scatter(g2[:,0],g2[:,1],color="black")
	plt.scatter(g3[:,0],g3[:,1],color="black")
	plt.title("raw data")
	plt.show()

	data = np.vstack([g1,g2,g3])[random.sample(list(np.arange(300)),300)]#インデックスごちゃ混ぜ

	C, E, loss = kmeans(data,n_clus=3)

	#結果の可視化
	cols = np.array(["red","blue","green"])
	for i in range(300):
		cl_idx = np.where(E[i,]==1)
		plt.scatter(data[i,0],data[i,1],color=cols[cl_idx])
	plt.scatter(C[:,0],C[:,1],color="black",marker="*")
	plt.show()

このスクリプトをターミナル等で実行すると，
これが，

こうなります．ちょっとプログラムいじってもっとクラスター構造が不明瞭なデータを発生させてみても面白いかもしれません．また，データの発生を乱数のシード固定していますが，お好みで好きなの使ってください．

次なにやる？

未定．K-meansやったから混合ガウス分布のパラメータ推定とかもありかも．