はじめに

突然ですが、パターン認識についていろいろ書いていこうと思います。

用語について

パターン認識という言葉を使います。いろいろな流派？があり、機械学習、人工知能などいろいろありますが、ここではパターン認識として考えます。それぞれ、突っ込みどころのある用語ですが、まぁそれはおいておいて。

流れ

パターン認識の全てについてですが、モデルを用いて行いことをします。

モデルを用意します。モデルのパラメータを求める必要があります。これを「学習」と呼びます。

モデルを用いて、観測が何であるかを推定します。

K-means algorithm

たくさんのデータがあり、それらが分類されていたとします。
新しく入ってきたデータが、それまでに分類されていたどのクラスに属するかを考える、というタスクを考えます。

そうすると、学習として求めたいモデルは各クラスの代表値、認識は観測値が属するクラスを求める、ことに当たります。

使い方として、

たくさんのサンプルを分類する（クラスタリング）
分類したクラスタをもとに、新しいサンプルがどのクラスタに属するか求める

というのを考えてみたいと思います。

サンプルデータ

ここでは、次のように分布している2次元データを考えます。

この例だと、何となく４つのクラスタに分類されそうだなー、と思うと思います。
そういう分類を行いたいと思います。

実際には多次元のデータですし、クラスタの分類が明確ではない場合もあります。さらに、クラスタ数がいくつかも普通は分からない状態で使います。（そういう場合の「学習」を「教師無し学習」と呼ぶことがあります。）

これらの学習用のＮ個の観測を $x_n$ と書くことにします。

今、クラスタ数をＫ＝４とします。各クラスタは代表値をがあるとします。各学習サンプルは、４個あるクラスタのどれかに属します。最も自然な考え方として、代表値に対して最も近くにあるクラスタに属すると考えます。n番目のサンプルが属するクラスタのインデックスk[n] は

k[n] = \arg \min_k ||x_n -  \mu_k||^2

とする。自乗は無くてもあってここでは同じであるが、後の目的関数の最小化で説明しやすいのでつけることにする。

もし、クラスタの代表値が分かっていれば、各サンプルがどこに属するかを求めることができます。一方で、各サンプルがどのクラスタに属するかが分かれば、各クラスタに属するサンプルを用いて代表値を更新することができます。これを繰り返すことが、k-means clustering です。

以下は実行結果です。
最初に乱数で初期の代表値を与えました。初期値の与え方はいろいろありますし、実はかなりデリケートな問題です。それは実際に動かすと出てきますが、ここではとりあえず置いておきます。

☆が各クラスタの代表値の推定値です。繰り返しのステップをくりかえすごとに各クラスタの中心に収束していくことが分かります。

実装は、いろいろあるでしょうが、ここでは変数を全て保持して計算を行う、以下のようにしました。

def kmeans_update_alignment(x:ndarray, mu:ndarray) -> ndarray:
    """
    サンプルがどのクラスタに属するかを求める。
    parameters
    ----------
    x[N, D]: samples
    mu[K, D]: centroids
    """
    N, D = x.shape
    r = np.zeros((N, K))
    for n in range(N):
        r[n, argmin([ sum(v*v for v in x[n, :] - mu[k, :]) for k in range(K)])] = 1
    return r

def kmeans_update_centroid(x:ndarray, r:ndarray) -> ndarray:
    """
    各クラスタの平均値を求める
    Parameters
    ----------
    x, training data (N, D)
    r, alignment information(K, D)
    """
    mu = dot(r.transpose(), x)
    for k in range(K):
        if sum(r[:,k]) < 0.01:
            continue
        mu[k,:] = mu[k,:] / sum(r[:,k])
    return mu

各クラスタの平均を求めるとき、もしそのクラスタに属するサンプルが一つもない場合があり得ます。そのときは更新しないような実装になっています。

###　目的関数の最小化

このようにad hoc に更新を行い、無事に収束するところを紹介しましたが、このクラスタリングはいつでも収束します。そのことを確認するために、定式化して考えます。

J = \sum_{n=1}^N \sum_{k=1}^K r_{n,k} || x_n - \mu_k ||^2

を最小化します。

(1) 各サンプルが属するクラスタを求める

$\mu$ を固定したときの$r_{n,k}$の最適化です。かならず目的関数が減少します。

(2) クラスタの中心の更新

$r_{n,k}$固定したときに最適な目的関数を最小にする$\mu_k$を求めます。これは、各クラスタごとに独立に最適化することができます。

J = \sum_{x_n\in C_1} || x_n - \mu_k||^2 + \cdots + \sum_{x_n\in C_1} || x_n - \mu_k||^2

自乗誤差を最小にするベクトルは平均ベクトルです。形をみると、各クラスタを正規分布と思い、その正規分布の平均の最尤推定値を求めることと同じ計算をしていますね。

\hat{\mu_k} = \frac{1}{N_k}\sum_{x_n \in C_k} x_n

以上から、単調に減少し、かつ、これは下限があるので、結果として必ず収束することが分かります。（「収束」の定義が何かはおいておいて＾＾；）先の計算での目的関数は、下記のように収束しています。最初は代表値が遠くにあるので目的関数は大きいですが、あとは単調に収束していきます。k-means ではサンプルがどのクラスタに属するかという０１なので、全てのサンプルのクラスタへの帰属が更新により変化しなくなった時が、クラスタリングの処理が終了したときになります。

行列による表現

実は、これらの計算を行列で書くことができます。さきに書いたpython コードは行列表現で書いています。これらは何となく昔からある慣習？なのかな、と思っていますが、ここでは一応紹介だけしました。（時間があれば後で式も追加）

まとめ

K-means 法の概要と実装を復習しました。いろいろ思い出してきました。

緒言

アルゴリズムについては：

初期値の設定について：初期値はデリケートな問題です。クラスタを徐々に増やしていく方法に、ＬＢＧアルゴリズムというのがあります。余力があれば紹介したい
制約条件がある場合：分布について、ノルムが１など、制約がある場合は応用例として考えらえれる。spherical k-means と呼ばれるものがこれに当たるのか？

実際にいろいろ動かすと分かることですが：

サンプル数が多いと収束までのステップが長くなる。（境界のサンプルの取り合いでセントロイドが微妙に更新される）

サンプルが割り当てられないクラスタもある。
シミュレーションの場合、ガウス分布の平均値に収束するわけではない。（サンプルの割り当ての最適性とは異なる。生成モデルが異なるので。）
実装を公開
グラフの余白の調整

等々。

今日はここまで。次回はGMMの予定。
(2021/01/23)