k-meansクラスタリングをnumpyで実装してみる
k-means
k-meansとは
k-meansとは非階層型クラスタリングの手法です。
クラスタ重心との距離情報を使うことで、あらかじめ決めたクラスタ数に分けていきます。
よく使われる手法ですが、アルゴリズム自体は結構単純です。
各データのラベルを更新
すべてのデータに対して各クラスタ重心との距離を計算し、最も近いクラスタのラベルを割り当てます。 距離は何を使っても大丈夫ですが、一般的にはユークリッド距離が使われることが多いみたいです。
2, 3を収束するまで繰り返す
ラベル更新とクラスタ重心更新を何度も繰り返すことで近いデータに同じラベルが割り当てられます。 このとき、どれくらい繰り返し処理を行うかは実装次第ですが、一般的には「あらかじめ決められた繰り返し回数処理を行った」や「ラベルの更新が行われなくなった」あたりを条件にすることが多いようです。
numpyでの実装
scikit-learnとかを使えば簡単に使うこともできるのですが、勉強も兼ねてnumpyを使って自分で実装してみます。
実装は「機械学習のエッセンス」を参考にさせてもらいました。(特にデータとクラスタ重心の距離をまとめて計算する方法は参考になりました。)
import numpy as np import matplotlib.pyplot as plt class KMeans: def __init__(self, n_clusters, max_iter=100, random_state=0): self.n_clusters = n_clusters self.max_iter = max_iter self.cluster_centers = None np.random.seed(random_state) def fit_predict(self, X): # 1. クラスタ中心をランダムに決定する n = X.shape[0] pr = np.repeat(1/n,n) self.cluster_centers = X[np.random.choice(np.arange(n), self.n_clusters, p=pr, replace=False), :] labels = np.zeros(X.shape[0]) pre_labels = np.ones(X.shape[0]) iter_counter = 0 while iter_counter < self.max_iter: if (labels==pre_labels).all(): break # 2. 各データのラベルを最も近いクラスタのラベルに変更する dists = ((X[:, :, np.newaxis]-self.cluster_centers.T[np.newaxis, :, :])**2).sum(axis=1) pre_labels = labels labels = np.argmin(dists, axis=1) # 3. クラスタ重心を更新する for c in range(self.n_clusters): cluster_X = X[labels==c, :] self.cluster_centers[c, :] = cluster_X.mean(axis=0) iter_counter += 1 return labels #データの生成 np.random.seed(123) x1 = np.r_[np.random.normal(size=20,loc=1,scale=2),np.random.normal(size=20,loc=8,scale=2),np.random.normal(size=20,loc=15,scale=2),np.random.normal(size=20,loc=25,scale=2)] x2 = np.r_[np.random.normal(size=20,loc=15,scale=2),np.random.normal(size=20,loc=1,scale=2),np.random.normal(size=20,loc=20,scale=2),np.random.normal(size=20,loc=0,scale=2)] X = np.c_[x1,x2] # kmeansでのクラスタリング kmeans = KMeans(n_clusters=4) labels = kmeans.fit_predict(X) #可視化 plt.scatter(X[:,0],X[:,1],cmap="jet", c=labels, s=10,alpha=0.5) plt.show()
今回は4つの塊になるようなデータを人工的に生成して試してみました。
結果はこちらの通りです。
ちなみにクラスタ重心の初期値は、全データにランダムにラベルを割り振ってその重心を初期値にするのが本来のやり方のようなのですが、 その方法で初期化をすると更新時にどのデータにも割り当てられないラベルが発生して上手く動作しないことがありました。 なので今回はランダムにデータを選択して、それをクラスタ重心の初期値にする方法を取っています。(実装方法が悪いのかもしれませんが...)
実装のとき少し躓いた点としては、初期クラスタ重心を選択するときに同じデータが初期重心として選ばれて上手く動作しないことがありました。
np.random.choiceの引数にreplace=Falseを追加して重複を許さないようにすることで対策しています。
特異値分解をnumpyで実装してみる
特異値分解
特異値分解とは
固有値分解は正方行列にしか適用できませんでした。
それを任意のサイズの行列に適用できるように拡張したのが特異値分解と呼ばれるものです。
行列を下記のように3つの行列に分解することを特異値分解といいます。
ここで、とは直行行列、は行列の特異値を降順にi行i列に並べて残りを0にした行列です。
ここで特異値とは、の固有値の平方根で求められる値のことです。
つまり、正方行列にしか適用できなかった固有値分解を任意のサイズの行列で疑似的にできるように拡張したのが特異値分解と呼ばれる手法になります。
numpyでの実装
numpyのnp.linalg.svdメソッドを使うことで特異値分解をすることができます。
ちなみには対角成分(特異値)のみが返されるので、元の行列に復元したいときなどはnp.diagを使うなどして対角行列に戻す必要があります。
import numpy as np A = np.array([[1, 3, 4, 2], [1, -1, 2, -5], [-2, 2, -1, 3]]) U, S, VT = np.linalg.svd(A) print("U") print(U) print("Sigma") print(S) print("V^T") print(VT) print("") print("AA^Tの固有値の平方根") eigen_val, eigen_vec = np.linalg.eig(np.dot(A, A.T)) print(np.sqrt(eigen_val))
実行結果はこちら。
U [[-0.33922653 -0.93622648 -0.09168066] [ 0.7463714 -0.32718841 0.57954937] [-0.57258632 0.1281707 0.80976366]] Sigma [ 6.99488133 5.27073893 1.5135872 ] V^T [[ 0.22192192 -0.4159075 0.10127734 -0.87607905] [-0.28833837 -0.42217034 -0.85897888 0.02807978] [-0.7476666 0.5053795 -0.01148765 -0.43064399] [-0.55552344 -0.62720389 0.50176311 0.21504133]] AA^Tの固有値の平方根 [ 6.99488133 5.27073893 1.5135872 ]
固有値分解の定義とPython実装
固有値分解
固有値と固有ベクトル
ある正方行列と零ベクトルではないベクトル、スカラーが下記の条件を満たすとき、はの固有ベクトル、はの固有値といいます。
固有値の求め方
定義を少し変形させてみます。(ただし、は単位行列とする)
は零ベクトルではないので、が非正則(逆行列を持たない)である必要があります。 つまり、行列の行列式を0にするを計算すれば固有値を求めることができます。
例題
よって
のとき、
のとき、
固有値分解
が成り立つ。
とすると、
となる。
このように正方行列を固有値、固有ベクトルの行列に分解することを固有値分解といいます。
Pythonでの実装
Pythonで固有値、固有ベクトルを求めたいときはnp.linalg.eigを使えば簡単に計算できます。
import numpy as np A = np.array([[5, -2], [9, -6]]) eigen_val, eigen_vec = np.linalg.eig(A) print("行列Aの固有値") print(eigen_val) print("") print("行列Aの固有ベクトル") print(eigen_vec) print("") print("PΛP^-1の計算結果") P = eigen_vec Lambda = np.diag(eigen_val) print(np.dot(np.dot(P, Lambda), np.linalg.inv(P)))
実行結果はこちらのとおり
行列Aの固有値 [ 3. -4.] 行列Aの固有ベクトル [[0.70710678 0.21693046] [0.70710678 0.97618706]] PΛP^-1の計算結果 [[ 5. -2.] [ 9. -6.]]
ちゃんと手計算で求めた結果と同じになってます。
ブログをはじめました
自己紹介
地方のソフトウェア開発会社で働く5年目のシステムエンジニアです。
マイペースにAIや数学、プログラミングについて勉強しています。
今までは調べて終わりのことが多かったのですが、備忘録やアウトプットの練習も兼ねて勉強したことをブログにまとめていきたいと思います。
疲れない程度にのんびり続けていきたいと思いますので、よろしくお願いします。