k-means法を解説して実装してみる

こんにちは、ぐぐりら(@guglilac)です。
春休みに入って卒論の研究室から追い出されたため、 GPUが使えなくなりました( ;  ; )

ので、この機会に他の機会学習の手法をさらっと復習しておこう!と思い、解説とその実装をする記事をこれからいくつか書いていこうかなと思っています。

なるべくフルスクラッチしたいけどどうせ途中で挫折してscikit-learnに逃げる未来が見えるwww

今回は最初なのできつくなさそうなとこから。
教師なし学習のクラスタリング手法としてポピュラーなk-means法を解説して実装してみます。

こんな感じのきらきらーな結果ができた!
かきごおりみたい。



k-means法とは

k平均法とも言います。
教師なし学習の一種であるクラスタリングを行う手法であり、アルゴリズムも単純なことから実装もしやすく、最初にはちょうど良さそうです。がんばろう。

k-means法はこんな感じのアルゴリズムです。
これによってN個のデータをk個のグループにクラスタリングします。
  1. 適当にk個のグループに対応する代表点を決める(initialize)
  2. N個のデータを代表点が一番近いグループにそれぞれ割り当てる
  3. 割り当てたデータをグループごとに平均することで代表点を計算し直す
以降、更新量が小さくなるまで2,3のステップを繰り返す
k平均法の平均はステップ3の処理に対応している感じですね。
直感的にもそれっぽいアルゴリズムで理解しやすいと思います。

最初にこのアルゴリズムを知った時、「めちゃ適当につくった感がすごいなww」をと思ったんですが、導出する過程もちゃんとあるので、次は導出過程を解説します。

アルゴリズムの導出

損失関数を定義して、エラー度合いを計算できるようにします。
この損失関数が小さくなるように代表点を修正していくアルゴリズムになっています。

まず損失関数$J$は次のように定義されます。

\[J=\sum_{n=1}^{N}\sum_{k=1}^{K}r_{nk}(x_n-\mu_k)^2\]

xがデータで、N個あります。
kはグループに対応するインデックスで、全部でK個あります。
$\mu_k$はk番目のグループの代表点を表すベクトルです。
$r_{nk}$は、n個目のデータがk番目のグループに属していたら1、そうでなければ0をとる変数とします。

日本語で書くと、データが属しているグループの代表点との距離の二乗の和が損失関数ということですね。

この損失関数を小さくするように$\mu_k$を修正していけばいいわけですが、$\mu_k$を修正すると所属する点も変わる(すなわち$r_{nk}$も変わる)ので、$\mu_k$と$r_{nk}$を交互に固定し修正する処理を行います。

数学的には、

  1. $\mu_k$を固定してJを$r_{nk}$で偏微分して最小化する
  2. $r_{nk}$を固定してJを$\mu_k$で偏微分して最小化する
という二つのステップを交互に行います。

1ステップ目は大丈夫だと思います。
$r_{nk}$は各$n$に対して一つの$k$でのみ1で他の$k$では0になるので、各$n$について一番距離の小さい代表点の所を1にしてあげれば$J$は小さくなります。

2ステップ目はちゃんと偏微分しなくてはなりません。
偏微分して0になるところが最小になるので、そのような$\mu_k$を求めます。

\begin{eqnarray*}\frac{\partial J}{\partial \mu_k}&=&\frac{\partial}{\partial \mu_k}\sum_{n=1}^{N}\sum_{k=1}^{K}r_{nk}(x_n-\mu_k)^2\\&=&\frac{\partial}{\partial \mu_k}\sum_{n=1}^{N}\sum_{k=1}^{K}r_{nk}(-2x_n\mu_k+\mu_k^2)\\&=&-2\sum_{n=1}^{N}\sum_{k=1}^{K}r_{nk}(x_n-\mu_{k})\\&=&0\end{eqnarray*}

あとはこれを$\mu_k$について解いてあげることで、

\[\mu_k=\frac{\sum_{n=1}^{N}r_{nk}x_n}{\sum_{n=1}^{N}r_{nk}}\]

となり、最適な$\mu_k$は所属するデータの平均をとったものであるという結果が導かれます。ばんざい。

実験

平均ベクトルの異なる4つの2次元正規分布から1分布につき250個のデータをサンプリングしてデータを生成します。なので全部で1000点のデータが得られます。

この1000個のデータをk-means法で4つのグループにクラスタリングします。
果たして発生元の分布ごとに正しくクラスタリングしてくれるのでしょうか。

ソースコードにも書いてありますが、平均ベクトルは(2,2),(2,-2),(-2,2),(-2,-2)にしてあります。

プロットしてみるとこんな感じになりました。

なんとなく4つのグループから作られてるのが人の目からもわかります。
ごましおみたい

ソースコード

k-means法の部分だけ載せます。
データ生成は省略。まあ好きなデータでやってみるといいでしょう。
例によってPython3による実装です。

Colors = ["r", "g", "m", "b"]
max_iter = 100
Means = [[2, 2], [-2, -2], [2, -2], [-2, 2]]
n = 250
classes = len(Means)
X, Y = genarateData(Means, n)
plot(X, Y, n)
N = n * classes


mu = [np.random.uniform(min(X), max(X), 2) for _ in range(classes)]

print("initialized mu:{}".format(mu))

for _ in range(max_iter):
    r = np.zeros(N, dtype=int)
    # 一番近い代表点に割り当てる
    for i in range(N):
        r[i] = np.argmin([np.linalg.norm([X[i], Y[i]] - mu[k])
                          for k in range(classes)])
    # 代表点を修正する
    mu_prev = np.array(copy.copy(mu))
    L = [[] for _ in range(classes)]

    for i in range(N):
        L[r[i]].append(np.array([X[i], Y[i]]))
    mu = np.array([np.mean(L[k], axis=0) for k in range(classes)])
    assert mu.shape == mu_prev.shape, "{},{}".format(mu.shape, mu_prev.shape)
    diff = mu - mu_prev

    # グラフ作成
    plt.figure()
    for i in range(N):
        plt.scatter(X[i], Y[i], s=30,
                    c=Colors[r[i]], alpha=0.5, marker="x")

    for i in range(classes):
        plt.scatter([mu[i, 0]], [mu[i, 1]], c=Colors[i],
                    marker='o', edgecolors='k', linewidths=1)
    plt.title("iter:{}".format(_))
    plt.show()

    # 収束判定
    if np.abs(np.sum(diff)) < 0.0001:
        print("converged")
        break

結果

イテレーションは7回で収束判定に引っかかりました。
1回目と4回目と7回目のクラスタリングの様子を載せます。

円のマークが各グループの代表点です。

ピンク優勢。
もうすでにキレイに分かれてます。



収束後の各グループの代表点の座標はこんな感じ。


xy
オレンジ1.91.88
2.01-2.03
-2.012.08
ピンク-2.09-1.95

設定した平均ベクトルと近くなっています。よかった。


まとめ

教師なし学習のクラスタリングの手法の一つであるk-means法について解説し、実装するとこまでやりました。

アルゴリズム自体は大したことないけれど、データ作ったり図示したりするとこまでいれると結構時間かかるなあ。もっと自由にpythonを使えるようになりたいですね。

読んでくれてありがとうございました。

コメント