のんびり亀エンジニアの勉強日記

勉強したこと調べたことをゆるくまとめていきます。

k-meansクラスタリングをnumpyで実装してみる

k-means

k-meansとは

k-meansとは非階層型クラスタリングの手法です。
クラスタ重心との距離情報を使うことで、あらかじめ決めたクラスタ数に分けていきます。
よく使われる手法ですが、アルゴリズム自体は結構単純です。

  1. ランダムにクラスタ重心を決める
    あらかじめ決めたクラスタ数分、ランダムに重心となる座標を決定します。

    f:id:chamekame:20210111232326p:plain

  2. 各データのラベルを更新
    すべてのデータに対して各クラスタ重心との距離を計算し、最も近いクラスタのラベルを割り当てます。 距離は何を使っても大丈夫ですが、一般的にはユークリッド距離が使われることが多いみたいです。
    f:id:chamekame:20210111232720p:plain

  3. クラスタ重心を更新
    クラスタの重心を、それぞれのラベルが割り当てられたデータの平均に更新します。
    f:id:chamekame:20210111233022p:plain

  4. 2, 3を収束するまで繰り返す
    ラベル更新とクラスタ重心更新を何度も繰り返すことで近いデータに同じラベルが割り当てられます。 このとき、どれくらい繰り返し処理を行うかは実装次第ですが、一般的には「あらかじめ決められた繰り返し回数処理を行った」や「ラベルの更新が行われなくなった」あたりを条件にすることが多いようです。
    f:id:chamekame:20210111235544p:plain

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つの塊になるようなデータを人工的に生成して試してみました。
結果はこちらの通りです。
f:id:chamekame:20210112001754p:plain

ちなみにクラスタ重心の初期値は、全データにランダムにラベルを割り振ってその重心を初期値にするのが本来のやり方のようなのですが、 その方法で初期化をすると更新時にどのデータにも割り当てられないラベルが発生して上手く動作しないことがありました。 なので今回はランダムにデータを選択して、それをクラスタ重心の初期値にする方法を取っています。(実装方法が悪いのかもしれませんが...)

実装のとき少し躓いた点としては、初期クラスタ重心を選択するときに同じデータが初期重心として選ばれて上手く動作しないことがありました。
np.random.choiceの引数にreplace=Falseを追加して重複を許さないようにすることで対策しています。

特異値分解をnumpyで実装してみる

特異値分解

特異値分解とは

固有値分解は正方行列にしか適用できませんでした。
それを任意のサイズの行列に適用できるように拡張したのが特異値分解と呼ばれるものです。

 m\times n行列 {\bf A}を下記のように3つの行列に分解することを特異値分解といいます。

 \displaystyle
{\bf A}{\bf x} = {\bf U}{\bf \Sigma}{\bf T}

ここで、 {\bf U} {\bf V}は直行行列、 {\bf \Sigma}は行列 {\bf A}の特異値 \sigma_iを降順にi行i列に並べて残りを0にした行列です。

 \displaystyle
{\bf \Sigma} = 
\begin{bmatrix}
\sigma_0 & 0 & \cdots & 0 \\
0 & \sigma_1 & \cdots & 0 \\
\vdots & \vdots &  & \vdots \\
0 & 0 & \cdots & \sigma_n \\
\vdots & \vdots &  & \vdots \\
0 & 0 & \cdots & 0 \\
\end{bmatrix}

ここで特異値とは、 {\bf A}{\bf A}^T固有値平方根で求められる値のことです。
つまり、正方行列にしか適用できなかった固有値分解を任意のサイズの行列で疑似的にできるように拡張したのが特異値分解と呼ばれる手法になります。

numpyでの実装

numpyのnp.linalg.svdメソッドを使うことで特異値分解をすることができます。
ちなみに \Sigmaは対角成分(特異値)のみが返されるので、元の行列に復元したいときなどは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 ]

特異値がちゃんと {\bf A}{\bf A}^T固有値平方根になっていることが分かります。

固有値分解の定義とPython実装

固有値分解

固有値固有ベクトル

ある正方行列 {\bf A} と零ベクトルではないベクトル {\bf x}スカラー {\lambda}が下記の条件を満たすとき、 {\bf x} {\bf A} 固有ベクトル {\lambda} {\bf A} 固有値といいます。

 \displaystyle
{\bf A}{\bf x} = \lambda{\bf x}
固有値の求め方

定義を少し変形させてみます。(ただし、 {\bf I}単位行列とする)

 \displaystyle
\begin{align}
{\bf A}{\bf x} &= \lambda{\bf x} \\
{\bf A}{\bf x} &= \lambda{\bf I}{\bf x} \\
({\bf A}-\lambda{\bf I}){\bf x} &= {\bf 0}
\end{align}

 {\bf x}は零ベクトルではないので、 {\bf A}-\lambda{\bf I}が非正則(逆行列を持たない)である必要があります。 つまり、行列 {\bf A}-\lambda{\bf I}行列式を0にする \lambdaを計算すれば固有値を求めることができます。

例題

   {\bf A} = \begin{bmatrix}5 & -2 \\9 & -6 \end{bmatrix}固有値 \lambda固有ベクトル {\bf x}を求める。

 \displaystyle
\begin{align}
{\bf A}-\lambda{\bf I} &= \begin{bmatrix}5-\lambda & -2 \\9 & -6-\lambda \end{bmatrix}
\end{align}

  よって

 \displaystyle
\begin{align}
\det({\bf A}-\lambda{\bf I}) &= 0 \\
(5-\lambda)(-6-\lambda)- (-2) \times9 &= 0 \\
\lambda^2 + \lambda - 12 &= 0 \\
(\lambda+4)(\lambda-3) = 0 \\
\lambda = -4, 3
\end{align}

   \lambda=-4のとき、

 \displaystyle
\begin{align}
\begin{bmatrix}9 & -2 \\9 & -2 \end{bmatrix} {\bf x} = {\bf 0} \\
{\bf x} = \begin{bmatrix}2 \\9 \end{bmatrix}
\end{align}

   \lambda=3のとき、

 \displaystyle
\begin{align}
\begin{bmatrix}2 & -2 \\9 & -9 \end{bmatrix} {\bf x} = {\bf 0} \\
{\bf x} = \begin{bmatrix}1 \\1 \end{bmatrix}
\end{align}

固有値分解

正方行列 {\bf A} 固有ベクトル {\bf x}_1, {\bf x}_2, \cdots , {\bf x}_n固有値 \lambda_1, \lambda_2, \cdots , \lambda_nのとき

 \displaystyle
\begin{align}
{\bf A}\begin{bmatrix} {\bf x}_1 & {\bf x}_2 & \cdots & {\bf x}_n \end{bmatrix} &= \begin{bmatrix} \lambda_1{\bf x}_1 & \lambda_2{\bf x}_2 & \cdots & \lambda_n{\bf x}_n \end{bmatrix}  \\
&= \begin{bmatrix} {\bf x}_1 & {\bf x}_2 & \cdots & {\bf x}_n \end{bmatrix} 
       \begin{bmatrix} 
       \lambda_1 & 0 & \cdots & 0\\
       0 & \lambda_2 & \cdots & 0\\
       \vdots & \vdots &  & \vdots\\
       0 & 0 & \cdots & \lambda_n\\
       \end{bmatrix} 
\end{align}

が成り立つ。

 \displaystyle
\begin{align}
{\bf P} =  \begin{bmatrix} {\bf x}_1 & {\bf x}_2 & \cdots & {\bf x}_n \end{bmatrix} , 
\Lambda =  \begin{bmatrix} 
       \lambda_1 & 0 & \cdots & 0\\
       0 & \lambda_2 & \cdots & 0\\
       \vdots & \vdots &  & \vdots\\
       0 & 0 & \cdots & \lambda_n\\
       \end{bmatrix} 
\end{align}

とすると、

 \displaystyle
\begin{align}
{\bf A} = {\bf P}\Lambda{\bf P}^{-1}
\end{align}

となる。
このように正方行列 {\bf A}固有値固有ベクトルの行列に分解することを固有値分解といいます。

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.]]

ちゃんと手計算で求めた結果と同じになってます。

行列、ベクトルと多項式の関係

機械学習分野では、多項式を行列・ベクトルで表現していることがよくあります。
落ち着いて考えれば何てことないのですが、"なんでこんな式変換になるの?"と混乱することがあったので改めてまとめてみます。

一次式のベクトル表現


\begin{align}
w_{1}x_{1} +w_{2}x_{2} + \cdots + w_{n}x_{n}&=\sum_{i=1}^{n} w_{i}x_{i} \\
&=
\begin{bmatrix}
w_{1} & w_{2} & \cdots & w_{n} \\
\end{bmatrix}
\begin{bmatrix}
x_{1} \\ x_{2} \\ \vdots \\ x_{n} \\
\end{bmatrix}
 \\
&={\bf w}^T {\bf x}
\end{align}

重み付き総和などは、よくこのような式で表現されます。
また、このように表現される式を線形結合と呼びます。

二次式の行列・ベクトル表現


\begin{align}
{\bf x} = 
\begin{bmatrix}
x_1 \\x_2
\end{bmatrix}, 

{\bf A} = 
\begin{bmatrix}
a & b \\
c & d
\end{bmatrix}
\end{align}

としたとき、


\begin{align}
ax_1^2+(b+c)x_{1}x_{2}+dx_2^2 &= 

\begin{bmatrix}
x_1 & x_2
\end{bmatrix}
\begin{bmatrix}
a & b \\
c & d
\end{bmatrix}
\begin{bmatrix}
x_1 \\ x_2
\end{bmatrix}\\

 &= {\bf x}^T{\bf A}{\bf x}

\end{align}

このように表現される式を二次形式と呼びます。

二乗和のベクトル表現


\begin{align}
x_{1}^2 + x_{2}^2 + \cdots + x_{n}^2&=\sum_{i=1}^{n} x_{i}^2 \\
&=
\begin{bmatrix}
x_{1} & x_{2} & \cdots & x_{n} \\
\end{bmatrix}
\begin{bmatrix}
x_{1} \\ x_{2} \\ \vdots \\ x_{n} \\
\end{bmatrix}
 \\
&={\bf x}^T {\bf x}\\
&=\begin{Vmatrix}
{\bf x}
\end{Vmatrix}^2\\
\end{align}

ブログをはじめました

自己紹介

地方のソフトウェア開発会社で働く5年目のシステムエンジニアです。
マイペースにAIや数学、プログラミングについて勉強しています。
今までは調べて終わりのことが多かったのですが、備忘録やアウトプットの練習も兼ねて勉強したことをブログにまとめていきたいと思います。

疲れない程度にのんびり続けていきたいと思いますので、よろしくお願いします。