Python: KMeans 法を実装してみる

KMeans 法は、機械学習における教師なし学習のクラスタリングという問題を解くためのアルゴリズム。教師なし学習というのは、事前に教師データというヒントが与えられないことを指している。その上で、クラスタリングというのは未知のデータに対していくつかのまとまりを作る問題をいう。

今回取り扱う KMeans 法は、比較的単純なアルゴリズムにも関わらず広く使われているものらしい。実際に書いてみても、基本的な実装であればたしかにとてもシンプルだった。ただし、データの初期化をするところで一点考慮すべき内容があることにも気づいたので、それについても書く。 KMeans 法の具体的なアルゴリズムについてはサンプルコードと共に後述する。

今回使った環境は次の通り。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.12.3
BuildVersion:   16D32
$ python --version
Python 3.5.3

依存パッケージをインストールする

あらかじめ、今回ソースコードで使う依存パッケージをインストールしておく。グラフ描画ライブラリの matplotlib を Mac OS X で動かすには、pip でのインストール以外にもちょっとした設定が必要になる。

$ pip install scipy scikit-learn matplotlib
$ mkdir -p ~/.matplotlib
$ cat << 'EOF' > ~/.matplotlib/matplotlibrc
backend: TkAgg
EOF

まずは scikit-learn のお手本から

本来はお手本の例については後回しにしたいんだけど、今回は構成の都合で先に示しておく。 Python の機械学習系ライブラリである scikit-learn には KMeans 法の実装も用意されている。自分で書いた実装をお披露目する前に、どんなものなのか見てもらいたい。

次のサンプルコードでは、ダミーのデータを生成した上でそれをクラスタリングしている。 scikit-learn にはデータセットを生成する機能も備わっている。その中でも make_blobs() という関数は、いくつかのまとまりを持ったダミーデータを作ってくれる。あとは、そのダミーデータを KMeans 法の実装である sklearn.cluster.KMeans で処理させている。クラスタリングした処理結果は matplotlib を使って可視化した。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from matplotlib import pyplot as plt
from sklearn import datasets
from sklearn.cluster import KMeans


def main():
    # クラスタ数
    N_CLUSTERS = 5

    # Blob データを生成する
    dataset = datasets.make_blobs(centers=N_CLUSTERS)

    # 特徴データ
    features = dataset[0]
    # 正解ラベルは使わない
    # targets = dataset[1]

    # クラスタリングする
    cls = KMeans(n_clusters=N_CLUSTERS)
    pred = cls.fit_predict(features)

    # 各要素をラベルごとに色付けして表示する
    for i in range(N_CLUSTERS):
        labels = features[pred == i]
        plt.scatter(labels[:, 0], labels[:, 1])

    # クラスタのセントロイド (重心) を描く
    centers = cls.cluster_centers_
    plt.scatter(centers[:, 0], centers[:, 1], s=100,
                facecolors='none', edgecolors='black')

    plt.show()


if __name__ == '__main__':
    main()

上記のサンプルコードでは、クラスタ数として適当に 5 を選んだ。ここで「クラスタ数を選ぶ」というのはポイントとなる。 KMeans 法では、あらかじめクラスタ数をハイパーパラメータとして指定するためだ。つまり、データをいくつのまとまりとして捉えるかを人間が教えてやる必要がある。クラスタリングの手法によっては、これをアルゴリズムがやってくれる場合もある。

それでは、上記のサンプルコードを実行してみよう。

$ python kmeans_scikit.py

すると、こんな感じのグラフが出力される。それぞれの特徴ベクトルが所属するクラスタが、別々の色で表現されている。また、それぞれのクラスタの真ん中らへんにある黒い円は、後述するセントロイドという概念を表している。 f:id:momijiame:20170319081628p:plain ちなみに、生成されるダミーデータは毎回異なるので、出力されるグラフも異なる。

自分で実装してみる

scikit-learn がクラスタリングしたお手本を見たところで、次は KMeans 法を自分で書いてみよう。

まず、KMeans 法ではセントロイド (重心) という概念が重要になる。セントロイドというのは、文字通りクラスタの中心に置かれるものだ。アルゴリズムの第一歩としては、このセントロイドを作ることになる。セントロイドは、前述した通り各クラスタの中心に置かれる。そのため、セントロイドはハイパーパラメータとして指定したクラスタ数だけ必要となる。また、中心というのは、クラスタに属する特徴ベクトルの各次元でそれらの平均値にいる、ということ。

そして、データセットに含まれる特徴ベクトルは、必ず最寄りのセントロイドに属する。ここでいう最寄りとは、特徴ベクトル同士のユークリッド距離が近いという意味だ。属したセントロイドが、それぞれのクラスタを表している。

上記の基本的な概念を元に KMeans 法のアルゴリズムは次のようなステップに分けることができる

最初のセントロイドを何らかの方法で決める
特徴ベクトルを最寄りのセントロイドに所属させる
各クラスタごとにセントロイドを再計算する
上記 2, 3 番目の操作を繰り返す
もしイテレーション前後で所属するセントロイドが変化しなければ、そこで処理を終了する

ここで、一番最初のステップである「最初のセントロイドを何らかの方法で決める」のがポイントだった。例えば、次のページではあらかじめ各特徴ベクトルをランダムにクラスタリングしてから、セントロイドを決めるやり方を取っている。

tech.nitoyon.com

上記のやり方であれば、次のようなサンプルコードになる。ここでは KMeans というクラスで KMeans 法を実装している。インターフェースは、先ほどお手本として提示した scikit-learn と揃えた。そのため main() 関数は基本的に変わっておらず、見るべきは KMeans クラスだけとなっている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from matplotlib import pyplot as plt
from sklearn import datasets
import numpy as np


class KMeans(object):
    """KMeans 法でクラスタリングするクラス"""

    def __init__(self, n_clusters=2, max_iter=300):
        """コンストラクタ

        Args:
            n_clusters (int): クラスタ数
            max_iter (int): 最大イテレーション数
        """
        self.n_clusters = n_clusters
        self.max_iter = max_iter

        self.cluster_centers_ = None

    def fit_predict(self, features):
        """クラスタリングを実施する

        Args:
            features (numpy.ndarray): ラベル付けするデータ

        Returns:
            numpy.ndarray: ラベルデータ
        """
        # 初期データは各要素に対してランダムにラベルをつける
        pred = np.random.randint(0, self.n_clusters, len(features))

        # クラスタリングをアップデートする
        for _ in range(self.max_iter):

            # 各クラスタごとにセントロイド (重心) を計算する
            self.cluster_centers_ = np.array([features[pred == i].mean(axis=0)
                                              for i in range(self.n_clusters)])

            # 各特徴ベクトルから最短距離となるセントロイドを基準に新しいラベルをつける
            new_pred = np.array([
                np.array([
                    self._euclidean_distance(p, centroid)
                    for centroid in self.cluster_centers_
                ]).argmin()
                for p in features
            ])

            if np.all(new_pred == pred):
                # 更新前と内容を比較して、もし同じなら終了
                break

            pred = new_pred

        return pred

    def _euclidean_distance(self, p0, p1):
        return np.sum((p0 - p1) ** 2)


def main():
    # クラスタ数
    N_CLUSTERS = 5

    # Blob データを生成する
    dataset = datasets.make_blobs(centers=N_CLUSTERS)

    # 特徴データ
    features = dataset[0]
    # 正解ラベルは使わない
    # targets = dataset[1]

    # クラスタリングする
    cls = KMeans(n_clusters=N_CLUSTERS)
    pred = cls.fit_predict(features)

    # 各要素をラベルごとに色付けして表示する
    for i in range(N_CLUSTERS):
        labels = features[pred == i]
        plt.scatter(labels[:, 0], labels[:, 1])

    centers = cls.cluster_centers_
    plt.scatter(centers[:, 0], centers[:, 1], s=100,
                facecolors='none', edgecolors='black')

    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python kmeans_random.py

実行すると、次のようなグラフが表示される。 f:id:momijiame:20170319153805p:plain 一見すると、上手くいっているようだ。

しかし、何度か実行すると、次のようなエラーが出ることがある。

$ python kmeans_random.py 
kmeans_random.py:41: RuntimeWarning: Mean of empty slice.
  for i in range(self.n_clusters)])
/Users/amedama/.virtualenvs/kmeans/lib/python3.5/site-packages/numpy/core/_methods.py:73: RuntimeWarning: invalid value encountered in true_divide
  ret, rcount, out=ret, casting='unsafe', subok=False)

表示されるグラフは、次のようなものになる。 f:id:momijiame:20170319153914p:plain 何やら、ぜんぜん上手いことクラスタリングできていない。

上記が起こる原因について調べたところ、セントロイドの初期化に問題があることが分かった。初期化時に特徴ベクトルに対してランダムにラベルをつけるやり方では、一つの特徴ベクトルも属さないセントロイドが生じる恐れがあるためだ。

ランダムにラベルをつけるやり方では、セントロイドはおのずとデータセット全体の平均あたりに生じやすくなる。場合によっては、セントロイドが別のセントロイドに囲まれるなどして、一つの特徴ベクトルも属さないセントロイドがでてくる。サンプルコードでは、そのような事態を想定していなかったのでセントロイドが計算できなくなって上手く動作しなかった。

先ほどの KMeans 法を図示しているページでも、特徴ベクトルの数 N に対してクラスタの数 k を増やして実行してみよう。一つの特徴ベクトルも属さないクラスタが生じることが分かる。

初期化を改良してみる

一つの特徴ベクトルも属さないクラスタが生じることについて、問題がないと捉えることもできるかもしれない。その場合には、一つも属さないクラスタのセントロイドを、前回の位置から動かなさないようにすることでケアできそうだ。

とはいえ、ここでは別のアプローチで問題を解決してみることにする。それというのは、初期のセントロイドをデータセットに含まれるランダムな特徴ベクトルにする、というものだ。つまり、最初に作るセントロイドの持つ特徴ベクトルがランダムに選んだ要素の特徴ベクトルと一致する。こうすれば各セントロイドは、少なくとも一つ以上の特徴ベクトルが属することは確定できる。ユークリッド距離がゼロになる特徴ベクトルが、必ず一つは存在するためだ。

次のサンプルコードでは、初期のセントロイドの作り方だけを先ほどと変更している。ここでポイントとなるのは、最初にセントロイドとして採用する特徴ベクトルは重複しないように決める必要があるということだ。そこで、初期のセントロイドは特徴ベクトルをシャッフルした上で先頭から k 個を取り出して決めている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from matplotlib import pyplot as plt
from sklearn import datasets
import numpy as np


class KMeans(object):
    """KMeans 法でクラスタリングするクラス"""

    def __init__(self, n_clusters=2, max_iter=300):
        """コンストラクタ

        Args:
            n_clusters (int): クラスタ数
            max_iter (int): 最大イテレーション数
        """
        self.n_clusters = n_clusters
        self.max_iter = max_iter

        self.cluster_centers_ = None

    def fit_predict(self, features):
        """クラスタリングを実施する

        Args:
            features (numpy.ndarray): ラベル付けするデータ

        Returns:
            numpy.ndarray: ラベルデータ
        """
        # 要素の中からセントロイド (重心) の初期値となる候補をクラスタ数だけ選び出す
        feature_indexes = np.arange(len(features))
        np.random.shuffle(feature_indexes)
        initial_centroid_indexes = feature_indexes[:self.n_clusters]
        self.cluster_centers_ = features[initial_centroid_indexes]

        # ラベル付けした結果となる配列はゼロで初期化しておく
        pred = np.zeros(features.shape)

        # クラスタリングをアップデートする
        for _ in range(self.max_iter):
            # 各要素から最短距離のセントロイドを基準にラベルを更新する
            new_pred = np.array([
                np.array([
                    self._euclidean_distance(p, centroid)
                    for centroid in self.cluster_centers_
                ]).argmin()
                for p in features
            ])

            if np.all(new_pred == pred):
                # 更新前と内容が同じなら終了
                break

            pred = new_pred

            # 各クラスタごとにセントロイド (重心) を再計算する
            self.cluster_centers_ = np.array([features[pred == i].mean(axis=0)
                                              for i in range(self.n_clusters)])

        return pred

    def _euclidean_distance(self, p0, p1):
        return np.sum((p0 - p1) ** 2)


def main():
    # クラスタ数
    N_CLUSTERS = 5

    # Blob データを生成する
    dataset = datasets.make_blobs(centers=N_CLUSTERS)

    # 特徴データ
    features = dataset[0]
    # 正解ラベルは使わない
    # targets = dataset[1]

    # クラスタリングする
    cls = KMeans(n_clusters=N_CLUSTERS)
    pred = cls.fit_predict(features)

    # 各要素をラベルごとに色付けして表示する
    for i in range(N_CLUSTERS):
        labels = features[pred == i]
        plt.scatter(labels[:, 0], labels[:, 1])

    centers = cls.cluster_centers_
    plt.scatter(centers[:, 0], centers[:, 1], s=100,
                facecolors='none', edgecolors='black')

    plt.show()


if __name__ == '__main__':
    main()