CUBE SUGAR CONTAINER

技術系のこと書きます。

Python: アンサンブル学習の Voting を試す

今回は機械学習におけるアンサンブル学習の一種として Voting という手法を試してみる。 これは、複数の学習済みモデルを用意して多数決などで推論の結果を決めるという手法。 この手法を用いることで最終的なモデルの性能を上げられる可能性がある。 実装については自分で書いても良いけど scikit-learn に使いやすいものがあったので、それを選んだ。

sklearn.ensemble.VotingClassifier — scikit-learn 0.20.2 documentation

使った環境は次の通り。

$ sw_vers                                                  
ProductName:    Mac OS X
ProductVersion: 10.14.1
BuildVersion:   18B75
$ python -V                    
Python 3.7.1

下準備

まずは今回使うパッケージをインストールしておく。

$ pip install scikit-learn tqdm 

とにかく混ぜてみる

とりあえず、最初は特に何も考えず複数のモデルを使って Voting してみる。

以下のサンプルコードでは乳がんデータセットを使って Voting を試している。 使ったモデルはサポートベクターマシン、ランダムフォレスト、ロジスティック回帰、k-最近傍法、ナイーブベイズの五つ。 モデルの性能は 5-Fold CV を使って精度 (Accuracy) について評価している。

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

from collections import defaultdict

import numpy as np
from tqdm import tqdm
from sklearn import datasets
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import VotingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import accuracy_score
from sklearn.svm import SVC
from sklearn.model_selection import StratifiedKFold
from sklearn.naive_bayes import GaussianNB


def main():
    # 乳がんデータセットを読み込む
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    # voting に使う分類器を用意する
    estimators = [
        ('svm', SVC(gamma='scale', probability=True)),
        ('rf', RandomForestClassifier(n_estimators=100)),
        ('logit', LogisticRegression(solver='lbfgs', max_iter=10000)),
        ('knn', KNeighborsClassifier()),
        ('nb', GaussianNB()),
    ]

    accs = defaultdict(list)

    skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    for train_index, test_index in tqdm(list(skf.split(X, y))):
        X_train, X_test = X[train_index], X[test_index]
        y_train, y_test = y[train_index], y[test_index]

        # 分類器を学習する
        voting = VotingClassifier(estimators)
        voting.fit(X_train, y_train)

        # アンサンブルで推論する
        y_pred = voting.predict(X_test)
        acc = accuracy_score(y_test, y_pred)
        accs['voting'].append(acc)

        # 個別の分類器の性能も確認してみる
        for name, estimator in voting.named_estimators_.items():
            y_pred = estimator.predict(X_test)
            acc = accuracy_score(y_test, y_pred)
            accs[name].append(acc)

    for name, acc_list in accs.items():
        mean_acc = np.array(acc_list).mean()
        print(name, ':', mean_acc)


if __name__ == '__main__':
    main()

上記に適当な名前をつけて実行してみる。 それぞれのモデルごとに計測した性能が出力される。

$ python voting.py 
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:03<00:00,  1.64it/s]
voting : 0.957829934590227
svm : 0.9385148133897653
rf : 0.9648787995382839
logit : 0.949134282416314
knn : 0.9314659484417083
nb : 0.9401923816852635

なんと Voting するよりもランダムフォレスト単体の方が性能が良いという結果になってしまった。 このように Voting するからといって必ずしも性能が上がるとは限らない。 例えば今回のように性能が突出したモデルがあるなら、それ単体で使った方が良くなる可能性はある。 あるいは、極端に性能が劣るモデルがあるならそれは取り除いた方が良いかもしれない。 それ以外には、次の項目で説明するモデルの重み付けという手もありそう。

モデルに重みをつける

性能が突出したモデルを単体で使ったり、あるいは劣るモデルを取り除く以外の選択肢として、モデルの重み付けがある。 これは、多数決などで推論結果を出す際に、特定のモデルの意見を重要視・あるいは軽視するというもの。 scikit-learn の VotingClassifier であれば weights というオプションでモデルの重みを指定できる。

次のサンプルコードでは、ランダムフォレストとロジスティック回帰の意見を重要視するように重みをつけてみた。

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

from collections import defaultdict

import numpy as np
from tqdm import tqdm
from sklearn import datasets
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import VotingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import accuracy_score
from sklearn.svm import SVC
from sklearn.model_selection import StratifiedKFold
from sklearn.naive_bayes import GaussianNB


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    estimators = [
        ('svm', SVC(gamma='scale', probability=True)),
        ('rf', RandomForestClassifier(n_estimators=100)),
        ('logit', LogisticRegression(solver='lbfgs', max_iter=10000)),
        ('knn', KNeighborsClassifier()),
        ('nb', GaussianNB()),
    ]

    accs = defaultdict(list)

    skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    for train_index, test_index in tqdm(list(skf.split(X, y))):
        X_train, X_test = X[train_index], X[test_index]
        y_train, y_test = y[train_index], y[test_index]

        # 分類器に重みをつける
        voting = VotingClassifier(estimators,
                                  weights=[1, 2, 2, 1, 1])
        voting.fit(X_train, y_train)

        y_pred = voting.predict(X_test)
        acc = accuracy_score(y_test, y_pred)
        accs['voting'].append(acc)

        for name, estimator in voting.named_estimators_.items():
            y_pred = estimator.predict(X_test)
            acc = accuracy_score(y_test, y_pred)
            accs[name].append(acc)

    for name, acc_list in accs.items():
        mean_acc = np.array(acc_list).mean()
        print(name, ':', mean_acc)


if __name__ == '__main__':
    main()

上記を実行してみよう。 先ほどよりも Voting したときの性能は向上している。 ただ、やはりランダムフォレスト単体での性能には届いていない。

$ python weight.py
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:03<00:00,  1.59it/s]
voting : 0.9613697575990766
svm : 0.9385148133897653
rf : 0.9666487110427088
logit : 0.949134282416314
knn : 0.9314659484417083
nb : 0.9401923816852635

Seed Averaging

先ほどの例では、モデルに重み付けしてみたものの結局ランダムフォレストを単体で使った方が性能が良かった。 とはいえ Voting は一つのアルゴリズムだけを使う場合にも性能向上につなげる応用がある。 それが、続いて紹介する Seed Averaging という手法。 これは、同じアルゴリズムでも学習に用いるシード値を異なるものにしたモデルを複数用意して Voting するというやり方。

次のサンプルコードでは、Voting で使うアルゴリズムはランダムフォレストだけになっている。 ただし、初期化するときのシード値がそれぞれ異なっている。

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

from collections import defaultdict

import numpy as np
from tqdm import tqdm
from sklearn import datasets
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import VotingClassifier
from sklearn.metrics import accuracy_score
from sklearn.model_selection import StratifiedKFold


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    # Seed Averaging
    estimators = [
        ('rf1', RandomForestClassifier(n_estimators=100, random_state=0)),
        ('rf2', RandomForestClassifier(n_estimators=100, random_state=1)),
        ('rf3', RandomForestClassifier(n_estimators=100, random_state=2)),
        ('rf4', RandomForestClassifier(n_estimators=100, random_state=3)),
        ('rf5', RandomForestClassifier(n_estimators=100, random_state=4)),
    ]

    accs = defaultdict(list)

    skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    for train_index, test_index in tqdm(list(skf.split(X, y))):
        X_train, X_test = X[train_index], X[test_index]
        y_train, y_test = y[train_index], y[test_index]

        voting = VotingClassifier(estimators)
        voting.fit(X_train, y_train)

        y_pred = voting.predict(X_test)
        acc = accuracy_score(y_test, y_pred)
        accs['voting'].append(acc)

        for name, estimator in voting.named_estimators_.items():
            y_pred = estimator.predict(X_test)
            acc = accuracy_score(y_test, y_pred)
            accs[name].append(acc)

    for name, acc_list in accs.items():
        mean_acc = np.array(acc_list).mean()
        print(name, ':', mean_acc)


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python sa.py 
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:03<00:00,  1.37it/s]
voting : 0.9683878414774914
rf1 : 0.9648787995382839
rf2 : 0.9666179299730666
rf3 : 0.9683570604078492
rf4 : 0.9666179299730665
rf5 : 0.9666179299730666

今回は、最も性能の良い三番目のモデルよりも、わずかながら Voting した結果の方が性能が良くなっている。 これは、各モデルの推論結果を平均することで、最終的なモデルの識別境界がなめらかになる作用が期待できるためと考えられる。

Soft Voting と Hard Voting

Voting と一口に言っても、推論結果の出し方には Soft Voting と Hard Voting という二つのやり方がある。 分かりやすいのは Hard Voting で、これは単純に各モデルの意見を多数決で決めるというもの。 もうひとつの Soft Voting は、それぞれのモデルの出した推論結果の確率を平均するというもの。 そこで、続いては、それぞれの手法について詳しく見ていくことにする。

Hard Voting

まずは Hard Voting から見ていく。

次のサンプルコードでは、結果を分かりやすいようにするために scikit-learn のインターフェースを備えたダミーの分類器を書いた。 この分類器は、インスタンスを初期化したときに指定された値をそのまま返すだけの分類器になっている。 つまり、fit() メソッドでは何も学習しない。 この分類器を三つ使って、これまたダミーの学習データに対して Hard Voting してみよう。 scikit-learn の VotingClassifier はデフォルトで Soft Voting なので、Hard Voting するときは明示的に指定する必要がある。

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

from sklearn.base import BaseEstimator
from sklearn.ensemble import VotingClassifier


class EchoBinaryClassifier(BaseEstimator):
    """インスタンス化したときに指定した値をオウム返しする二値分類器"""

    def __init__(self, answer_proba):
        """
        :param answer_proba: オウム返しする値 (0 ~ 1)
        """
        self.answer_proba = answer_proba

    def fit(self, X, y=None):
        # 何も学習しない
        return self

    def predict(self, X, y=None):
        # 指定された値が 0.5 以上なら 1 を、0.5 未満なら 0 を推論結果として返す
        return [1 if self.answer_proba >= 0.5 else 0 for _ in X]

    def predict_proba(self, X, y=None):
        # 指定された値を推論結果の確率としてそのまま返す
        return [(float(1 - self.answer_proba), float(self.answer_proba)) for _ in X]


def main():
    # ダミーの入力 (何も学習・推論しないため)
    dummy = [0, 0, 1]

    estimators = [
        ('ebc1', EchoBinaryClassifier(1.0)),
        ('ebc2', EchoBinaryClassifier(1.0)),
        ('ebc3', EchoBinaryClassifier(0.0)),
    ]

    # Hard voting する
    voting = VotingClassifier(estimators, voting='hard')
    voting.fit(dummy, dummy)

    # voting の結果を表示する
    y_pred = voting.predict(dummy)
    print('predict:', y_pred)

    # Hard voting は単純な多数決なので確率 (probability) は出せない
    # y_pred_proba = voting.predict_proba(dummy)
    # print(y_pred_proba)


if __name__ == '__main__':
    main()

上記のサンプルコードにおいて三つの分類器は、正反対の推論結果を返すことになる。 具体的には、1 を返すものが二つ、0 を返すものが一つある。 Voting による最終的な推論結果はどうなるだろうか。

上記を実行してみよう。

$ python hard.py      
predict: [1 1 1]

全て 1 と判定された。 これは多数派の判定結果として 1 が二つあるためだ。

一応、もうちょっと際どい確率でも試してみよう。 今度は、それぞれのモデルが 0.51, 0.51, 0.0 を返すようになっている。 もし、確率で平均したなら (0.51 + 0.51 + 0.0) / 3 = 0.34 < 0.5 となって 0 と判定されるはずだ。

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

from sklearn.base import BaseEstimator
from sklearn.ensemble import VotingClassifier


class EchoBinaryClassifier(BaseEstimator):
    """インスタンス化したときに指定した値をオウム返しする二値分類器"""

    def __init__(self, answer_proba):
        """
        :param answer_proba: オウム返しする値 (0 ~ 1)
        """
        self.answer_proba = answer_proba

    def fit(self, X, y=None):
        # 何も学習しない
        return self

    def predict(self, X, y=None):
        # 指定された値が 0.5 以上なら 1 を、0.5 未満なら 0 を推論結果として返す
        return [1 if self.answer_proba >= 0.5 else 0 for _ in X]

    def predict_proba(self, X, y=None):
        # 指定された値を推論結果の確率としてそのまま返す
        return [(float(1 - self.answer_proba), float(self.answer_proba)) for _ in X]


def main():
    # ダミーの入力 (何も学習・推論しないため)
    dummy = [0, 0, 1]

    estimators = [
        ('ebc1', EchoBinaryClassifier(0.51)),
        ('ebc2', EchoBinaryClassifier(0.51)),
        ('ebc3', EchoBinaryClassifier(0.00)),
    ]

    # Hard voting する
    voting = VotingClassifier(estimators, voting='hard')
    voting.fit(dummy, dummy)

    # voting の結果を表示する
    y_pred = voting.predict(dummy)
    print('predict:', y_pred)

    # Hard voting は単純な多数決なので確率 (probability) は出せない
    # y_pred_proba = voting.predict_proba(dummy)
    # print(y_pred_proba)


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python hard2.py 
predict: [1 1 1]

これは、やはり多数派として 1 があるため。

Soft Voting

続いては、先ほどのサンプルコードをほとんどそのまま流用して手法だけ Soft Voting にしてみよう。 Soft Voting では確率の平均を取るため、今度は (0.51 + 0.51 + 0.0) / 3 = 0.34 < 0.5 となって 0 に判定されるはず。

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

from sklearn.base import BaseEstimator
from sklearn.ensemble import VotingClassifier


class EchoBinaryClassifier(BaseEstimator):
    """インスタンス化したときに指定した値をオウム返しする二値分類器"""

    def __init__(self, answer_proba):
        """
        :param answer_proba: オウム返しする値 (0 ~ 1)
        """
        self.answer_proba = answer_proba

    def fit(self, X, y=None):
        # 何も学習しない
        return self

    def predict(self, X, y=None):
        # 指定された値が 0.5 以上なら 1 を、0.5 未満なら 0 を推論結果として返す
        return [1 if self.answer_proba >= 0.5 else 0 for _ in X]

    def predict_proba(self, X, y=None):
        # 指定された値を推論結果の確率としてそのまま返す
        return [(float(1 - self.answer_proba), float(self.answer_proba)) for _ in X]


def main():
    # ダミーの入力 (何も学習・推論しないため)
    dummy = [0, 0, 1]

    estimators = [
        ('ebc1', EchoBinaryClassifier(0.51)),
        ('ebc2', EchoBinaryClassifier(0.51)),
        ('ebc3', EchoBinaryClassifier(0.00)),
    ]

    # Soft voting する
    voting = VotingClassifier(estimators, voting='soft')
    voting.fit(dummy, dummy)

    # voting の結果を表示する
    y_pred = voting.predict(dummy)
    print('predict:', y_pred)

    # Soft voting は確率の平均を出す
    y_pred_proba = voting.predict_proba(dummy)
    print('proba:', y_pred_proba)


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python soft.py         
predict: [0 0 0]
proba: [[0.66 0.34]
 [0.66 0.34]
 [0.66 0.34]]

無事、今度は判定結果が 0 になることが確認できた。

めでたしめでたし。

統計的学習の基礎 ―データマイニング・推論・予測―

統計的学習の基礎 ―データマイニング・推論・予測―

  • 作者: Trevor Hastie,Robert Tibshirani,Jerome Friedman,杉山将,井手剛,神嶌敏弘,栗田多喜夫,前田英作,井尻善久,岩田具治,金森敬文,兼村厚範,烏山昌幸,河原吉伸,木村昭悟,小西嘉典,酒井智弥,鈴木大慈,竹内一郎,玉木徹,出口大輔,冨岡亮太,波部斉,前田新一,持橋大地,山田誠
  • 出版社/メーカー: 共立出版
  • 発売日: 2014/06/25
  • メディア: 単行本
  • この商品を含むブログ (6件) を見る

Python: Optuna で機械学習モデルのハイパーパラメータを選ぶ

今回は、ハイパーパラメータを最適化するフレームワークの一つである Optuna を使ってみる。 このフレームワークは国内企業の Preferred Networks が開発の主体となっていて、ほんの数日前にオープンソースになったばかり。

ハイパーパラメータ自動最適化ツール「Optuna」公開 | Preferred Research

先に使ってみた印象について話してしまうと、基本は Hyperopt にかなり近いと感じた。 実際のところ、使っているアルゴリズムの基本は変わらないし、定義できるパラメータの種類もほとんど同じになっている。 おそらく Hyperopt を使ったことがある人なら、すぐにでも Optuna に切り替えることができると思う。

その上で Hyperopt との違いについて感じたのは二点。 まず、Define-by-run という特性によって複雑なパラメータを構成しやすくなっていること。 そして、もうひとつが最適化を分散処理 (Distributed Optimization) するときの敷居がとても低いことだった。 これらの点に関して詳しくは後述する。 あと、これは使い勝手とは違うけどフレームワークのソースコードが Python のお手本のようにキレイ。

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

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.14.1
BuildVersion:   18B75
$ python -V              
Python 3.6.7

下準備

Optuna は PyPI に登録されているので pip を使ってインストールできる。

$ pip install optuna

単純な例

機械学習モデルのハイパーパラメータを扱う前に、まずはもっと単純な例から始めてみよう。 手始めに、次のような一つの変数 x を取る関数の最大値を探す問題について考えてみる。 グラフを見ると、最大値は 2 付近にあるようだ。

f:id:momijiame:20180818132951p:plain

以下のサンプルコードでは、上記の関数の最大値を Optuna で探している。 ただし、Optuna は今のところ (v0.4.0) は最小化にのみ対応しているため、目的関数は符号を反転している。

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

import math

import optuna


def objective(trial):
    """最小化する目的関数"""
    # パラメータが取りうる範囲
    x = trial.suggest_uniform('x', -5, +15)
    # デフォルトで最小化かつ現在は最小化のみのサポートなので符号を反転する
    return - math.exp(-(x - 2) ** 2) + math.exp(-(x - 6) ** 2 / 10) + 1 / (x ** 2 + 1)


def main():
    # 最適化のセッションを作る
    study = optuna.create_study()
    # 100 回試行する
    study.optimize(objective, n_trials=100)
    # 最適化したパラメータを出力する
    print('params:', study.best_params)


if __name__ == '__main__':
    main()

コメントもあるけど、念のためソースコードについて補足する。 まず、Optuna では一つの最適化したい問題を Study というオブジェクト (または概念) で取り扱う。 このオブジェクトの optimize() というメソッドに最適化したい目的関数を渡して試行を繰り返す。 サンプルコードにおいて、目的関数は objective() という名前で定義している。 また、各試行は Trial というオブジェクトで、目的関数の引数として渡される。 この Trial は過去に探索した内容の履歴を扱っていて、それを元に次に試行するパラメータの値が決まる。 なお、パラメータの探索にはデフォルトで TPE (Tree-structured Parzen Estimator) というベイズ最適化の一種が用いられる。

あんまり説明ばかりが長くなっても仕方ないので、先ほどのサンプルコードを実際に実行してみる。 すると、各試行ごとにログが出力されて探索の様子が分かるようになっている。

$ python helloptuna.py 
[I 2018-12-05 21:34:11,776] Finished a trial resulted in value: 0.5459903986050455. Current best value is 0.5459903986050455 with parameters: {'x': -0.9268011342183158}.
[I 2018-12-05 21:34:11,777] Finished a trial resulted in value: -0.11535207904386202. Current best value is -0.11535207904386202 with parameters: {'x': 1.2859333865518803}.
[I 2018-12-05 21:34:11,778] Finished a trial resulted in value: 0.014796275560628312. Current best value is -0.11535207904386202 with parameters: {'x': 1.2859333865518803}.
...(snip)...
[I 2018-12-05 21:34:12,484] Finished a trial resulted in value: -0.46165044813597855. Current best value is -0.597104581903869 with parameters: {'x': 2.029056526434356}.
[I 2018-12-05 21:34:12,495] Finished a trial resulted in value: 0.020372692700788855. Current best value is -0.597104581903869 with parameters: {'x': 2.029056526434356}.
[I 2018-12-05 21:34:12,508] Finished a trial resulted in value: -0.24643573203869212. Current best value is -0.597104581903869 with parameters: {'x': 2.029056526434356}.
params: {'x': 2.029056526434356}

結果として、たしかに 2 前後にある最大値が見つかった。

機械学習のモデルに適用する

続いては、実際に機械学習のモデルを使ってハイパーパラメータを探してみよう。 その準備として scikit-learn をインストールしておく。

$ pip install scikit-learn

まずは、乳がんデータセットと SVM (Support Vector Machine) の組み合わせを試してみる。 これは、二値分類という問題が分かりやすいのと SVM がハイパーパラメータに対して割と敏感な印象があるため。 最小化する目的関数の返り値は、モデルを 5-Fold CV で Accuracy について評価したもの。 また、カーネル関数のようなカテゴリ変数や、応答の特性が対数に比例する類のパラメータもちゃんと扱えることが分かる。

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

from functools import partial

import optuna
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate
from sklearn.svm import SVC
from sklearn import datasets


def objective(X, y, trial):
    """最小化する目的関数"""
    params = {
        'kernel': trial.suggest_categorical('kernel', ['rbf', 'sigmoid']),
        'C': trial.suggest_loguniform('C', 1e+0, 1e+2),
        'gamma': trial.suggest_loguniform('gamma', 1e-2, 1e+1),
    }

    # モデルを作る
    model = SVC(**params)

    # 5-Fold CV / Accuracy でモデルを評価する
    kf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    scores = cross_validate(model, X=X, y=y, cv=kf)
    # 最小化なので 1.0 からスコアを引く
    return 1.0 - scores['test_score'].mean()


def main():
    # 乳がんデータセットを読み込む
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target
    # 目的関数にデータを適用する
    f = partial(objective, X, y)
    # 最適化のセッションを作る
    study = optuna.create_study()
    # 100 回試行する
    study.optimize(f, n_trials=100)
    # 最適化したパラメータを出力する
    print('params:', study.best_params)


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python svmopt.py
[I 2018-12-05 21:36:11,195] Finished a trial resulted in value: 0.37257406694882644. Current best value is 0.37257406694882644 with parameters: {'kernel': 'rbf', 'C': 2.0190784147815957, 'gamma': 0.8877063047664268}.
[I 2018-12-05 21:36:11,332] Finished a trial resulted in value: 0.37257406694882644. Current best value is 0.37257406694882644 with parameters: {'kernel': 'rbf', 'C': 2.0190784147815957, 'gamma': 0.8877063047664268}.
[I 2018-12-05 21:36:11,413] Finished a trial resulted in value: 0.37257406694882644. Current best value is 0.37257406694882644 with parameters: {'kernel': 'rbf', 'C': 2.0190784147815957, 'gamma': 0.8877063047664268}.
...(snip)...
[I 2018-12-05 21:36:22,924] Finished a trial resulted in value: 0.37077337437475955. Current best value is 0.37077337437475955 with parameters: {'kernel': 'rbf', 'C': 86.65350935649266, 'gamma': 0.01187563934327901}.
[I 2018-12-05 21:36:23,093] Finished a trial resulted in value: 0.37257406694882644. Current best value is 0.37077337437475955 with parameters: {'kernel': 'rbf', 'C': 86.65350935649266, 'gamma': 0.01187563934327901}.
[I 2018-12-05 21:36:23,254] Finished a trial resulted in value: 0.37077337437475955. Current best value is 0.37077337437475955 with parameters: {'kernel': 'rbf', 'C': 86.65350935649266, 'gamma': 0.01187563934327901}.
params: {'kernel': 'rbf', 'C': 86.65350935649266, 'gamma': 0.01187563934327901}

スコアの改善自体は 0.2% 程度のものだけど、ちゃんと最初よりも性能の良いパラメータが見つかっていることが分かる。

もう少し複雑な例を扱ってみる

さきほどの例では、パラメータの組み合わせが単純だったので、正直 Optuna の真価を示すことができていなかった。 そこで、次はもう少しだけ複雑な例を扱ってみることにする。

次のサンプルコードでは、先ほどと比較して使う分類器が RandomForest と SVM の二種類に増えている。 当然ながら RandomForest と SVM ではインスタンス化するときの引数からして異なる。 こういったパターンでも、Optuna なら Trial オブジェクトが返した値を元に処理を分岐することで目的関数が簡単に表現できる。 Optuna では、このように実行しながら試行するパターンを決定していく特定を指して Define-by-run と呼んでいる。

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

from functools import partial

import optuna
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate
from sklearn.svm import SVC
from sklearn import datasets


def objective(X, y, trial):
    """最小化する目的関数"""

    # 使う分類器は SVM or RF
    classifier = trial.suggest_categorical('classifier', ['SVC', 'RandomForestClassifier'])

    # 選ばれた分類器で分岐する
    if classifier == 'SVC':
        # SVC のとき
        params = {
            'kernel': trial.suggest_categorical('kernel', ['rbf', 'sigmoid']),
            'C': trial.suggest_loguniform('C', 1e+0, 1e+2),
            'gamma': trial.suggest_loguniform('gamma', 1e-2, 1e+1),
        }
        model = SVC(**params)
    else:
        # RF のとき
        params = {
            'n_estimators': int(trial.suggest_loguniform('n_estimators', 1e+2, 1e+3)),
            'max_depth': int(trial.suggest_loguniform('max_depth', 2, 32)),
        }
        model = RandomForestClassifier(**params)

    kf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    scores = cross_validate(model, X=X, y=y, cv=kf, n_jobs=-1)
    return 1.0 - scores['test_score'].mean()


def main():
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target
    f = partial(objective, X, y)
    study = optuna.create_study()
    study.optimize(f, n_trials=100)
    print('params:', study.best_params)


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python clfopt.py
[I 2018-12-05 21:38:48,588] Finished a trial resulted in value: 0.04390919584455555. Current best value is 0.04390919584455555 with parameters: {'classifier': 'RandomForestClassifier', 'n_estimators': 733.1380680046733, 'max_depth': 2.3201133434932117}.
[I 2018-12-05 21:38:49,972] Finished a trial resulted in value: 0.045648326279338236. Current best value is 0.04390919584455555 with parameters: {'classifier': 'RandomForestClassifier', 'n_estimators': 733.1380680046733, 'max_depth': 2.3201133434932117}.
[I 2018-12-05 21:38:50,084] Finished a trial resulted in value: 0.37257406694882644. Current best value is 0.04390919584455555 with parameters: {'classifier': 'RandomForestClassifier', 'n_estimators': 733.1380680046733, 'max_depth': 2.3201133434932117}.
...(snip)...
[I 2018-12-05 21:42:30,883] Finished a trial resulted in value: 0.033382070026933386. Current best value is 0.028133897652943496 with parameters: {'classifier': 'RandomForestClassifier', 'n_estimators': 453.56949779973104, 'max_depth': 30.86562241046072}.
[I 2018-12-05 21:42:31,710] Finished a trial resulted in value: 0.03866102347056566. Current best value is 0.028133897652943496 with parameters: {'classifier': 'RandomForestClassifier', 'n_estimators': 453.56949779973104, 'max_depth': 30.86562241046072}.
[I 2018-12-05 21:42:31,826] Finished a trial resulted in value: 0.37077337437475955. Current best value is 0.028133897652943496 with parameters: {'classifier': 'RandomForestClassifier', 'n_estimators': 453.56949779973104, 'max_depth': 30.86562241046072}.
params: {'classifier': 'RandomForestClassifier', 'n_estimators': 453.56949779973104, 'max_depth': 30.86562241046072}

先ほどよりも実行に時間はかかかるものの、着実にスコアが上がる様子が見て取れる。

Storage と最適化の分散処理 (Distributed Optimization)

続いては Optuna の重要なオブジェクトとして Storage について扱う。 これは Study や各 Trial といった探索に関する情報を記録しておくための場所を示している。

Storage は普段オンメモリにあるため、実行ごとに揮発してしまう。 しかし、指定すれば実装を RDB (Relational Database) などに切り替えて使うことができる。 こうすれば探索した履歴をディスクに永続化できるので、途中までの計算結果が何かのはずみにパーという事態も減らせる。 Optuna では SQLAlchemy を使って RDB を抽象化して扱っているため SQLite3 や MySQL などを切り替えながら利用できる。 この点はバックエンドが基本的に MongoDB しかない Hyperopt に対して優位な点だと思う。

説明が長くなってきたので、そろそろ実際に使ってみる。 最初は Storage の実装として SQLite3 を使うパターンで試してみよう。 まずは optuna コマンドを使ってデータベースを初期化する。 このとき Study に対して分かりやすい名前をつけておく。

$ optuna create-study --study 'distributed-helloworld' --storage 'sqlite:///example.db'
[I 2018-12-05 22:04:07,038] A new study created with name: distributed-helloworld
distributed-helloworld

これで必要な SQLite3 のデータベースが作られる。

$ file example.db 
example.db: SQLite 3.x database, last written using SQLite version 3025003
$ sqlite3 example.db '.tables'
studies                  trial_params             trial_values           
study_system_attributes  trial_system_attributes  trials                 
study_user_attributes    trial_user_attributes    version_info

指定した名前で Study に関する情報も記録されている。

$ sqlite3 example.db 'SELECT * FROM studies'
1|distributed-helloworld|MINIMIZE

この時点では、まだ全くパラメータを探索していないので Trial に関する情報はない。

$ sqlite3 example.db 'SELECT COUNT(1) AS count FROM trials'
0

上記の Storage を使って探索するサンプルコードが次の通り。 Study のオブジェクトを作るときに、名前とバックエンドに関する情報を指定してやる。 試行回数は意図的に抑えてある。

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

import math

import optuna


def objective(trial):
    """最小化する目的関数"""
    x = trial.suggest_uniform('x', -5, +15)
    return - math.exp(-(x - 2) ** 2) + math.exp(-(x - 6) ** 2 / 10) + 1 / (x ** 2 + 1)


def main():
    # NOTE: 実行前のデータベースを作ること
    # $ optuna create-study --study 'distributed-helloworld' --storage 'sqlite:///example.db'

    # study_name, storage は作成したデータベースの内容と合わせる
    study = optuna.Study(study_name='distributed-helloworld', storage='sqlite:///example.db')

    study.optimize(objective, n_trials=5)
    print('params:', study.best_params)


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python dbstudy.py 
[I 2018-12-05 22:05:20,037] Finished a trial resulted in value: 0.6938982089558796. Current best value is 0.6938982089558796 with parameters: {'x': 7.969822697335578}.
[I 2018-12-05 22:05:20,075] Finished a trial resulted in value: 0.9714561387259608. Current best value is 0.6938982089558796 with parameters: {'x': 7.969822697335578}.
[I 2018-12-05 22:05:20,113] Finished a trial resulted in value: -0.515546413663528. Current best value is -0.515546413663528 with parameters: {'x': 1.7262219554205434}.
[I 2018-12-05 22:05:20,151] Finished a trial resulted in value: 0.5172823829449016. Current best value is -0.515546413663528 with parameters: {'x': 1.7262219554205434}.
[I 2018-12-05 22:05:20,191] Finished a trial resulted in value: 0.9779376537980013. Current best value is -0.515546413663528 with parameters: {'x': 1.7262219554205434}.
params: {'x': 1.7262219554205434}

5 回だけ試行が実行されて、ベストスコアは -0.515546413663528 となった。

データベースを確認すると Trial の回数が更新されていることが分かる。 つまり、探索した内容を永続化しつつ処理が進められている。

$ sqlite3 example.db 'SELECT COUNT(1) AS count FROM trials'
5

そして、サンプルコードをもう一度実行してみよう。 すると、初回の探索時に示されているベストスコアが前回から引き継がれていることがわかる。 つまり、過去にデータベースに永続化した内容を使いながら探索が進められている。 これはソースコードレベルでも確認していて、過去に試行した結果を全て引っ張ってきて、その内容を元に次の試行内容を決定していた。

$ python dbstudy.py 
[I 2018-12-05 22:05:29,048] Finished a trial resulted in value: 0.9358553655678095. Current best value is -0.515546413663528 with parameters: {'x': 1.7262219554205434}.
[I 2018-12-05 22:05:29,079] Finished a trial resulted in value: 0.04114448607195525. Current best value is -0.515546413663528 with parameters: {'x': 1.7262219554205434}.
[I 2018-12-05 22:05:29,118] Finished a trial resulted in value: 0.7127626383671348. Current best value is -0.515546413663528 with parameters: {'x': 1.7262219554205434}.
[I 2018-12-05 22:05:29,157] Finished a trial resulted in value: -0.5533154442856236. Current best value is -0.5533154442856236 with parameters: {'x': 2.2005954763853053}.
[I 2018-12-05 22:05:29,197] Finished a trial resulted in value: 0.4829461021707021. Current best value is -0.5533154442856236 with parameters: {'x': 2.2005954763853053}.
params: {'x': 2.2005954763853053}

確認すると、ちゃんと試行回数が増えている。

$ sqlite3 example.db 'SELECT COUNT(1) AS count FROM trials'
10

ちなみに、この仕組みは最適化の分散処理とも深く関わりがある。 というのも、このデータベースを複数のマシンから参照できる場所に置いた場合、どうなるだろうか。 各マシンは、データベースの内容を参照しつつ分散処理でパラメータの探索が可能になる。

また、複数のマシンを使わなくても同じマシンから複数のプロセスで探索するコードを実行すればマルチプロセスで処理が可能になる。 まあ、ただこの用途についてはあまり現実的ではない。 というのも Study 自体が n_jobs という並列度を制御するオプションを持っているのと、一般的にモデルの学習・推論部分が並列化されているので。

バックエンドを MySQL にしてみる

最適化の分散処理という観点でいえば本番環境で SQLite3 を使うことは、おそらくないはず。 そこで、次は Storage のバックエンドに MySQL を試してみる。

まずは MySQL をインストールする。

$ brew install mysql
$ brew services start mysql

Optuna 用のデータベースを作る。

$ mysql -u root -e "CREATE DATABASE IF NOT EXISTS optuna"

SQLAlchemy が接続に使うデータベースドライバのパッケージをインストールする。

$ pip install mysqlclient

あとは Storage の指定を MySQL が使われるように変えるだけ。

$ optuna create-study --study 'distributed-mysql' --storage 'mysql://root@localhost/optuna'
[I 2018-12-05 23:44:09,129] A new study created with name: distributed-mysql
distributed-mysql

ちゃんと MySQL にテーブルができていることを確認する。

$ mysql -u root -D optuna -e "SHOW TABLES"
+-------------------------+
| Tables_in_optuna        |
+-------------------------+
| studies                 |
| study_system_attributes |
| study_user_attributes   |
| trial_params            |
| trial_system_attributes |
| trial_user_attributes   |
| trial_values            |
| trials                  |
| version_info            |
+-------------------------+

次のサンプルコードもバックエンドに MySQL が使われるように指定している。

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

import math

import optuna


def objective(trial):
    """最小化する目的関数"""
    x = trial.suggest_uniform('x', -5, +15)
    return - math.exp(-(x - 2) ** 2) + math.exp(-(x - 6) ** 2 / 10) + 1 / (x ** 2 + 1)


def main():
    # NOTE: 実行前のデータベースを作ること
    # $ optuna create-study --study 'distributed-mysql' --storage 'mysql://root@localhost/optuna'

    # study_name, storage は作成したデータベースの内容と合わせる
    study = optuna.Study(study_name='distributed-mysql', storage='mysql://root@localhost/optuna')

    study.optimize(objective, n_trials=100)
    print('params:', study.best_params)


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python mysqlstudy.py 
[I 2018-12-05 23:45:07,051] Finished a trial resulted in value: 0.021628217984214517. Current best value is 0.0216282 with parameters: {'x': 12.4684}.
[I 2018-12-05 23:45:07,093] Finished a trial resulted in value: 0.9874823396086119. Current best value is 0.0216282 with parameters: {'x': 12.4684}.
[I 2018-12-05 23:45:07,140] Finished a trial resulted in value: 0.8407030713413823. Current best value is 0.0216282 with parameters: {'x': 12.4684}.
...(snip)...
[I 2018-12-05 23:45:12,811] Finished a trial resulted in value: 0.7561641667979796. Current best value is -0.507712 with parameters: {'x': 2.28833}.
[I 2018-12-05 23:45:12,925] Finished a trial resulted in value: 1.0257911319634438. Current best value is -0.507712 with parameters: {'x': 2.28833}.
[I 2018-12-05 23:45:13,003] Finished a trial resulted in value: 0.989037566945567. Current best value is -0.507712 with parameters: {'x': 2.28833}.
params: {'x': 2.28833}

データベースに記録されている試行回数が増えている。

$ mysql -u root -D optuna -e "SELECT COUNT(1) AS count FROM trials"
+-------+
| count |
+-------+
|   100 |
+-------+

ばっちり。

Python: Annoy の近似最近傍探索 (ANN) を試す

今回は Spotify の作った近似最近傍探索 (ANN: Approximate Nearest Neighbor algorithms search) ライブラリの Annoy を試してみる。 ANN は k-NN (k-Nearest Neighbor algorithms search) の一種で、厳密な解を追い求めない代わりに高いスループットが得られる。

ちなみに ANN のライブラリごとのベンチマークを公開している Web サイトがある。 その中でいうと Annoy は大体のベンチマークで真ん中くらいの位置にある。 その代わり Annoy はインストールが簡単という利点がある。

ANN-Benchmarks

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14
BuildVersion:   18A391
$ python -V                     
Python 3.7.1

下準備

まずは Annoy 本体と、データセットの読み込みなどに使うために scikit-learn をインストールしておく。

$ pip install annoy scikit-learn

上記のように pip から一発でインストールできるのは Annoy の良いところだと思う。

Annoy を試してみる

早速だけど Annoy のサンプルコードを以下に示す。 データセットには Iris データセットを使って、分類精度を 3-fold CV を使って計測した。

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

from sklearn import datasets
from annoy import AnnoyIndex
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import accuracy_score


def main():
    # Iris データセットを読み込む
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target

    # 3-fold CV
    skf = StratifiedKFold(n_splits=3, shuffle=True, random_state=42)
    for train_index, test_index in skf.split(X, y):
        # データを学習用と検証用に分ける
        X_train, X_test = X[train_index], X[test_index]
        y_train, y_test = y[train_index], y[test_index]

        # ユークリッド距離で探索するモデルを用意する
        annoy_index = AnnoyIndex(X.shape[1], metric='euclidean')

        # モデルを学習させる
        for i, x in enumerate(X_train):
            annoy_index.add_item(i, x)

        # k-d tree をビルドする
        annoy_index.build(n_trees=10)

        # 近傍 1 点を探索して、そのクラスを返す
        y_pred = [y_train[annoy_index.get_nns_by_vector(v, 1)[0]] for v in X_test]

        # 精度 (Accuracy) を計算する
        score = accuracy_score(y_test, y_pred)
        print('acc:', score)


if __name__ == '__main__':
    main()

上記の実行結果は次の通り。 だいたい 96% 前後の精度が出ている。

$ python helloannoy.py 
acc: 0.9607843137254902
acc: 0.9607843137254902
acc: 0.9583333333333334

もう少し詳しく見ていく

ここからはもう少し詳しく Annoy の API を眺めてみる。 まずは Python のインタプリタを起動する。

$ python

とりあえず動作確認用に Iris データセットを読み込んで、それをホールドアウト検証用に分割しておく。

>>> from sklearn import datasets
>>> dataset = datasets.load_iris()
>>> X, y = dataset.data, dataset.target
>>> from sklearn.model_selection import train_test_split
>>> X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)

続いて肝心の Annoy のモデルを用意する。 これには AnnoyIndex というクラスを用いる。 インスタンス化するときは、データの次元数と、距離の計算に用いるメトリックを指定する。 先ほどはユークリッド距離 (euclidean) を使っていたけど、ここではコサイン距離 (angular) を使っている。 これは正規化したユークリッド距離を表しているらしい。

>>> from annoy import AnnoyIndex
>>> annoy_index = AnnoyIndex(X.shape[1], metric='angular')

続いてモデルに学習データを登録していく。 ここで指定したインデックスが後ほど探索結果として得られる。

>>> for i, x in enumerate(X_train):
...     annoy_index.add_item(i, x)
... 

Annoy は k-d tree を元に探索を高速化している。 そこで続いては学習データを元に木をビルドする。

>>> annoy_index.build(n_trees=10)
True

これで分類の準備が整った。 あとは AnnoyIndex#get_nns_by_vector() メソッドで探索結果が得られる。 これには、近傍を見つけたいデータと、見つけたい近傍数を与える。 例えば近傍数が 1 ならこんな感じ。

>>> annoy_index.get_nns_by_vector(X_test[0], 1)
[85]

近傍数が 3 なら次のようになる。。

>>> annoy_index.get_nns_by_vector(X_test[0], 3)
[85, 71, 26]

近傍数と共に距離も手に入れたいときは include_distances オプションに True を渡す。 すると結果がタプルになって近傍のインデックスと共に距離が得られる。

>>> annoy_index.get_nns_by_vector(X_test[0], 3, include_distances=True)
([85, 71, 26], [0.021935366094112396, 0.025716988369822502, 0.034560393542051315])

Annoy が返すのは近傍のインデックスなので、ラベルに変換するには次のようにする。 分類させたデータのラベルを見ると 1 で、近傍 3 点のラベルも 1 になっているので、これはちゃんと分類できている。

>>> y_train[85], y_train[71], y_train[26]
(1, 1, 1)
>>> y_test[0]
1

scikit-learn の API に対応させてみる

続いて、試しに Annoy の実装を使って scikit-learn の分類器を書いてみた。 次のサンプルコードの AnnoyClassifier がそれに当たる。 scikit-learn の API に対応したことで sklearn.model_selection.cross_validate() がそのまま使えるようになっている。

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

from collections import Counter

from sklearn import datasets
from annoy import AnnoyIndex
from sklearn.base import BaseEstimator
from sklearn.base import ClassifierMixin
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate
from sklearn.utils import check_X_y
from sklearn.utils import check_array


class AnnoyClassifier(BaseEstimator, ClassifierMixin):
    """scikit-learn API に準拠した Annoy の分類器"""

    def __init__(self, n_trees, metric='angular', n_neighbors=1, search_k=-1):
        # k-d tree の数
        self.n_trees_ = n_trees
        # 計算に用いる距離
        self.metric_ = metric
        # 近傍数
        self.n_neighbors_ = n_neighbors
        # 精度とスループットのトレードオフに用いるパラメータ
        self.search_k_ = search_k
        # ビルドしたモデル
        self.clf_ = None
        # 学習データのクラスラベルを入れる
        self.train_y_ = None

    def fit(self, X, y):
        # 入力をバリデーションする
        check_X_y(X, y)

        # 学習データのクラスラベルを保存しておく
        self.train_y_ = y

        # Annoy のモデルを用意する
        self.clf_ = AnnoyIndex(X.shape[1], metric=self.metric_)

        # モデルを学習させる
        for i, x in enumerate(X):
            self.clf_.add_item(i, x)

        # k-d tree を構築する
        self.clf_.build(n_trees=self.n_trees_)

        return self

    def predict(self, X):
        # 入力をバリデーションする
        check_array(X)

        # 分類結果を返す
        y_pred = [self._predict(x) for x in X]
        return y_pred

    def _predict(self, x):
        # 近傍を見つける
        neighbors = self.clf_.get_nns_by_vector(x, self.n_neighbors_, search_k=self.search_k_)
        # インデックスをクラスラベルに変換する
        neighbor_classes = self.train_y_[neighbors]
        # 最頻値を取り出す
        counter = Counter(neighbor_classes)
        most_common = counter.most_common(1)
        # 最頻値のクラスラベルを返す
        return most_common[0][0]

    def get_params(self, deep=True):
        # 分類器のパラメータ
        return {
            'n_trees': self.n_trees_,
            'metric': self.metric_,
            'n_neighbors': self.n_neighbors_,
            'search_k': self.search_k_,
        }


def main():
    # Iris データセットを読み込む
    dataset = datasets.load_iris()
    X, y = dataset.data, dataset.target

    # Annoy を scikit-learn API でラップした分類器
    clf = AnnoyClassifier(n_trees=10)
    # 3-fold CV
    skf = StratifiedKFold(n_splits=3, shuffle=True, random_state=42)
    # 精度を評価指標にして汎化性能を計測する
    score = cross_validate(clf, X, y, cv=skf, scoring='accuracy')
    mean_test_score = score.get('test_score').mean()
    print('acc:', mean_test_score)


if __name__ == '__main__':
    main()

上記の実行結果が次の通り。

$ python scikitannoy.py  
acc: 0.9669117647058824

めでたしめでたし。

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

Python: k-NN Feature Extraction について

k-NN Feature Extraction (k-近傍法を用いた特徴量抽出) という手法があるらしい。 これは、文字通り k-NN (k-Nearest Neighbor algorithm: k-近傍法) を特徴量の抽出に応用したもの。 興味深かったので、今回は自分でも Python を使って実装してみた。 手法について知ったのは、以下のブログを目にしたのがきっかけ。

upura.hatenablog.com

また、上記は以下のブログに記載のある R の実装を参考にしているとのことだった。

Feature Extraction with KNN • fastknn

ただ、先ほどの Python API では、特徴量を付与する対象のデータをパラメータとして指定できない点が気になった。 具体的には、以下のような交差検証を使った性能の計測が難しいのではと感じた。

  • データセットを学習用と検証用に分割する
  • 学習用データに対して k-NN Feature Extraction で特徴量を付与する
    • 付与するデータと付与に使うデータに分けながら順番に付与していく
  • 特徴量が付与された学習用データを使ってモデルを学習させる
  • 学習用データを使って検証用データに k-NN Feature Extraction で特徴量を付与する
  • 特徴量が付与された検証用データを使ってモデルの性能を計測する

R の実装では、特徴量を付与するデータと付与に使うデータを別々に指定できる。 そこで、今回書いてみたものは R の実装にならって別々に指定できるようにしてみた。

使った環境は次の通り。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.13.6
BuildVersion:   17G3025
$ python -V
Python 3.7.1

k-NN Feature Extraction について

ここからは k-NN Feature Extraction の具体的な手法について解説していく。 まず、k-NN Feature Extraction では、名前の通り k-NN (k-近傍法) を特徴量の抽出に応用している。 元々の k-NN は分類問題を解くためのアルゴリズムだった。

最初に k-NN では、分類したい対象のデータから最も近い k 点の学習データを探し出す。 k 点の k には 1 や 3 など、任意の数を用いる。 最も近い k 点の学習データが見つかったら、その学習データについているクラスを使って対象のデータを分類する。 k が 1 より大きいときは、多数決などを使って分類先のクラスを決める。

k-NN Feature Extraction では、上記のアルゴリズムのうち前半部分だけを利用する。 具体的には、最も近い k 点の学習データを探し出す部分。 このとき、最も近い k 点の学習データを探すために、データ間の距離を計算することになる。 k-NN Feature Extraction では、このデータ間の距離を特徴量として用いる。

少し k-NN と異なるのは、k 点の近傍データを探すのがクラスごとという点。 つまり、二値分類問題であれば、学習データの中のクラス 0 とクラス 1 ごとに近傍データをそれぞれ探す。 もちろん多値分類問題であれば、クラス数 c の数だけ近傍データを探すことになる。 さらに、k が 1 を越える場合には、1 ~ k までの全ての数で近傍データを探すようだ。 これによって k * c 個の特徴量が新たに生成されることになる。

文章だけだとなかなか分かりにくいと思うので図示してみる。 まず、以下のような感じで特徴量を付与したいデータとクラスが三つに分かれた学習データがあるとする。

f:id:momijiame:20181111134125p:plain

特徴量を付与したいデータから、それぞれのクラスの近傍データまでの距離を計算して特徴量にする。

f:id:momijiame:20181111134354p:plain

もし k > 1 であれば、複数の近傍点までの距離を足し合わせたものを特徴量にする。

f:id:momijiame:20181111134457p:plain

いじょう。

下準備

実際に実装してみる前に、必要なパッケージを環境にインストールしておく。

$ pip install scikit-learn matplotlib

サンプルデータについて

特徴量の付与に使うサンプルデータについては、以下のような線形分離できないものを用いた。

f:id:momijiame:20181111135431p:plain

上記を生成したのは次のようなサンプルコードになる。

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

import numpy as np
from matplotlib import pyplot as plt


def main():
    # サンプルデータを生成する
    x0 = np.random.rand(500) - 0.5
    x1 = np.random.rand(500) - 0.5
    X = np.array(list(zip(x0, x1)))
    y = np.array([1 if i0 * i1 > 0 else 0 for i0, i1 in X])

    # データを可視化する
    plt.scatter(X[y == 0, 0], X[y == 0, 1])
    plt.scatter(X[y == 1, 0], X[y == 1, 1])
    plt.show()


if __name__ == '__main__':
    main()

適当な名前をつけて実行すると、先ほどのグラフが出力される。

$ python visualize.py

k-NN Feature Extraction で特徴量を抽出してみる

それでは、実際に k-NN Feature Extraction を試してみよう。 以下のサンプルコードの knn_extract() 関数が特徴量を生成する関数になる。 基本的な使い方としては、学習データとしてxtr と、そのクラスデータとして ytr を渡す。 そして、特徴量を生成する対象のデータとして xte を渡す。

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

import functools
from concurrent.futures import ProcessPoolExecutor
import multiprocessing

import numpy as np
from sklearn.model_selection import KFold
from matplotlib import pyplot as plt


def _distance(a, b):
    """2 点間の距離を計算する"""
    distance = np.linalg.norm(a - b)
    return distance


def _knn_distance(xtr_c, target, k):
    """特定クラスの最近傍な教師データへの距離を計算する"""
    # 各要素への距離を計算する
    distances = np.array([_distance(target, x) for x in xtr_c])
    # 距離で昇順ソートする
    sorted_distances = np.sort(distances)
    # 先頭 k 件を取り出す (最近傍への距離 k 点)
    nearest_distances = sorted_distances[:k]
    # 距離を足し算して返す
    sum_distances = np.sum(nearest_distances)
    return sum_distances


def _knn_distance_fold(xtr_c, k, folds, target):
    """過学習を防ぐためにサンプリングしながら最近傍への距離の平均を計算する"""
    # 途中までの計算結果を代入するための入れ物
    distances = np.empty([folds])

    # K-分割することでサンプリングする
    kf = KFold(shuffle=True, n_splits=folds)
    for i, (train_index, _) in enumerate(kf.split(xtr_c)):
        # サンプリングした教師信号
        xtr_c_sampled = xtr_c[train_index]
        # ターゲットと教師信号との距離を計算する
        distance = _knn_distance(xtr_c_sampled, target, k)
        distances[i] = distance

    # 各サンプルでの距離の平均を返す
    average_distance = distances.mean()
    return average_distance


def knn_extract(xtr, ytr, xte, k=1, folds=5, nprocesses=-1):
    if nprocesses == -1:
        # デフォルトでは CPU のコア数で並列計算する
        nprocesses = multiprocessing.cpu_count()

    # 教師データに含まれているクラス (ラベル) の種類
    # NOTE: CV するときは含まれるクラスが偏らないようにすること
    classes = np.unique(ytr)

    # 生成する特徴量の入れ物を用意する
    features = np.zeros([len(xte), len(classes) * k])

    for i, class_ in enumerate(classes):
        # 処理対象のクラス
        xtr_c = xtr[ytr == class_]

        for j, k_n in enumerate(range(1, k + 1), 1):
            # クラス class_ で近傍数 k の距離を計算する関数に実引数を部分適用する
            f = functools.partial(_knn_distance_fold, xtr_c, k_n, folds)

            # クラス class_ で近傍数 k で計算した特徴ベクトル
            with ProcessPoolExecutor(max_workers=nprocesses) as executor:
                feature = executor.map(f, xte)
                features[:, i * j] = list(feature)

    return features


def main():
    # サンプルデータを生成する
    x0 = np.random.rand(500) - 0.5
    x1 = np.random.rand(500) - 0.5
    X = np.array(list(zip(x0, x1)))
    y = np.array([1 if i0 * i1 > 0 else 0 for i0, i1 in X])

    classes = np.unique(y)
    k = 1
    extracted_features = np.empty([0, len(classes) * k])

    # データを 5 分割して、それぞれで特徴量を作っていく
    skf = KFold(n_splits=5)
    for train_index, test_index in skf.split(X):
        # 学習用データと検証用データに分ける
        X_train, X_test = X[train_index], X[test_index]
        y_train, y_test = y[train_index], y[test_index]

        # 検証用データに対して、学習用データを使って特徴量を抽出する
        extracted_feature = knn_extract(X_train, y_train, X_test, k)
        extracted_features = np.append(extracted_features, extracted_feature, axis=0)

    # k-NN Feature Extraction の結果を可視化する
    plt.scatter(extracted_features[:, 0], extracted_features[:, 1], c=y)
    plt.show()


if __name__ == '__main__':
    main()

上記に適当な名前をつけて実行しよう。

$ python knnfe.py

実行結果として得られるグラフは次の通り。 今回の問題は k = 1, c = 2 なので二次元の特徴量が新たに得られる。 それらを散布図にしてみると、ちゃんと線形分離できそうな形で特徴量が生成されている。

f:id:momijiame:20181111135927p:plain

続いて Iris データセットに使ってみた場合についても紹介する。 まず、データセットの読み込み部分を以下のように変更する。

# Iris データセットを読み込む
from sklearn import datasets
iris = datasets.load_iris()
X, y = iris.data, iris.target

こちらも、四次元の特徴量から以下のような二次元の特徴量が得られる。 それなりに上手く分離できそうな雰囲気がある。

f:id:momijiame:20181111140135p:plain

課題について

実際に試した上で気になった点としては、やはり計算量があげられる。 今回使ったようなサイズの小さいデータセットでも、実行してみると結構な時間がかかる。 これは最近傍の探索が素朴な実装になっているから。 そのため、実用上は ANN (Approximate Nearest Neighbor: 近似最近傍) を使う必要があると考えられる。

いじょう。

Python: 特徴量の重要度を Permutation Importance で計測する

学習させた機械学習モデルにおいて、どの特徴量がどれくらい性能に寄与しているのかを知りたい場合がある。 すごく効く特徴があれば、それについてもっと深掘りしたいし、あるいは全く効かないものがあるなら取り除くことも考えられる。 使うフレームワークやモデルによっては特徴量の重要度を確認するための API が用意されていることもあるけど、そんなに多くはない。 そこで、今回はモデルやフレームワークに依存しない特徴量の重要度を計測する手法として Permutation Importance という手法を試してみる。 略称として PIMP と呼ばれたりすることもあるようだ。

この手法を知ったのは、以下の Kaggle のノートブックを目にしたのがきっかけだった。

Permutation Importance | Kaggle

あんまりちゃんと読めてないけど、論文としては Altmann et al. (2010) になるのかな?

Permutation importance: a corrected feature importance measure | Bioinformatics | Oxford Academic

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.13.6
BuildVersion:   17G3025
$ python -V     
Python 3.7.1

Permutation Importance の考え方について

Permutation Importance の考え方はとてもシンプル。 ある特徴量が使い物にならないとき、モデルの性能がどれだけ落ちるかで重要度を計測する。 もし性能が大きく落ちるなら、その特徴量はモデルにおいて重要だと考えられるし、その逆もまたしかり。 下手をすると性能がむしろ上がって、ない方がマシということが分かったりするかもしれない。

この説明だけだと頭の中にクエスチョンマークがたくさん浮かぶと思うので、順を追って説明していく。 Permutation Importance では、まずデータセットを学習用のデータと検証用のデータに分割する。 その上で、学習用のデータを使ってモデルを学習させる。 そして、何も手を加えていない検証用のデータを使ってまずは性能を計測する。 性能の計測には、解決したい問題に対して任意の適切な評価指標を用いる。 ここまでの一連の流れは、一般的な機械学習モデルの汎化性能の計測方法と何も変わらない。 このときに測った性能を、後ほどベースラインとして用いる。

ベースラインが計測できたところで、続いては検証用データについて特定の特徴量だけを使い物にならない状態にする。 これには、手法の名前にもある通り "Permutation (交換・置換)" を用いる。 具体的には、特徴量をランダムにシャッフルすることで、目的変数に対して特徴量の相関が取れない状態にする。 特定の特徴量がシャッフルされた検証用データが完成したら、それを学習済みのモデルに与えて性能を計測する。 このとき、どれだけベースラインから性能が変化するかを確認する。

上記を全ての特徴量に対して実施して、ベースラインに対する性能の変化を特徴量ごとに比較する。 以上が Permutation Importance の計測方法になる。 ただ、文章だけだとちょっと分かりにくいかもと思ったので図示してみる。 最初のベースラインを測る部分は一般的な機械学習モデルの汎化性能を計測するやり方と同じ。

f:id:momijiame:20181118112928p:plain

ベースラインが計測できたら、続いて検証用データの特定の特徴量を Permutation で使い物にならない状態にする。 そして、その検証用データを使って性能を改めて検証する。 このとき、ベースラインからの性能の変化を観測する。

f:id:momijiame:20181118112944p:plain

上記において勘違いしやすいポイントは Permutation Importance でモデルの再学習は必要ないという点。 なんとなく大元のデータセットをシャッフルした上でモデルを再学習させて...と考えてしまうんだけど、実は必要ない。 モデルを再学習すると、それに計算コストも必要になるし、ベースラインで使ったモデルとの差異も大きくなってしまう。 そこで、学習済みのモデルと検証用データだけを使って性能を計測する。

実際に計測してみる

続いては実際に Permutation Importance を確認してみよう。

下準備として、使うパッケージをあらかじめインストールしておく。

$ pip install pandas scikit-learn matplotlib

今回はデータセットにみんな大好き Iris データセットを、モデルには Random Forest を選んだ。 結果がなるべく安定するように 5-fold CV でそれぞれ計算した結果を出している。 以下は手作業で計算処理を書いてるけど eli5 というパッケージに PIMP を計算するためのモジュールがあったりもする。

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

from collections import defaultdict

import numpy as np
import pandas as pd
from sklearn import datasets
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score
from sklearn.model_selection import StratifiedKFold
from matplotlib import pyplot as plt


def permuted(df):
    """特定のカラムをシャッフルしたデータフレームを返す"""
    for column_name in df.columns:
        permuted_df = df.copy()
        permuted_df[column_name] = np.random.permutation(permuted_df[column_name])
        yield column_name, permuted_df


def pimp(clf, X, y, cv=None, eval_func=accuracy_score):
    """PIMP (Permutation IMPortance) を計算する"""
    base_scores = []
    permuted_scores = defaultdict(list)

    if cv is None:
        cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

    for train_index, test_index in cv.split(X, y):
        # 学習用データと検証用データに分割する
        X_train, y_train = X.iloc[train_index], y.iloc[train_index]
        X_test, y_test = X.iloc[test_index], y.iloc[test_index]

        # 学習用データでモデルを学習する
        clf.fit(X_train, y_train)

        # まずは何もシャッフルしていないときのスコアを計算する
        y_pred_base = clf.predict(X_test)
        base_score = eval_func(y_test, y_pred_base)
        base_scores.append(base_score)

        # 特定のカラムをシャッフルした状態で推論したときのスコアを計算する
        permuted_X_test_gen = permuted(X_test)
        for column_name, permuted_X_test in permuted_X_test_gen:
            y_pred_permuted = clf.predict(permuted_X_test)
            permuted_score = eval_func(y_test, y_pred_permuted)
            permuted_scores[column_name].append(permuted_score)

    # 基本のスコアとシャッフルしたときのスコアを返す
    np_base_score = np.array(base_scores)
    dict_permuted_score = {name: np.array(scores) for name, scores in permuted_scores.items()}
    return np_base_score, dict_permuted_score


def score_difference_statistics(base, permuted):
    """シャッフルしたときのスコアに関する統計量 (平均・標準偏差) を返す"""
    mean_base_score = base.mean()
    for column_name, scores in permuted.items():
        score_differences = scores - mean_base_score
        yield column_name, score_differences.mean(), score_differences.std()


def main():
    # Iris データセットを読み込む
    dataset = datasets.load_iris()
    X = pd.DataFrame(dataset.data, columns=dataset.feature_names)
    y = pd.Series(dataset.target)

    # 計測に使うモデルを用意する
    clf = RandomForestClassifier(n_estimators=100)

    # Permutation Importance を計測する
    base_score, permuted_scores = pimp(clf, X, y)

    # 計測結果から統計量を計算する
    diff_stats = list(score_difference_statistics(base_score, permuted_scores))

    # カラム名、ベーススコアとの差、95% 信頼区間を取り出す
    sorted_diff_stats = sorted(diff_stats, key=lambda x: x[1])
    column_names = [name for name, _, _ in sorted_diff_stats]
    diff_means = [diff_mean for _, diff_mean, _ in sorted_diff_stats]
    diff_stds_95 = [diff_std * 1.96 for _, _, diff_std in sorted_diff_stats]

    # グラフにプロットする
    plt.plot(column_names, diff_means, marker='o', color='r')
    plt.errorbar(column_names, diff_means, yerr=diff_stds_95, ecolor='g', capsize=4)

    plt.title('Permutation Importance')
    plt.grid()
    plt.xlabel('column')
    plt.ylabel('difference')
    plt.show()


if __name__ == '__main__':
    main()

上記に適当な名前をつけて実行してみよう。

$ python pimp.py

すると、次のようなグラフが得られる。

f:id:momijiame:20181110230620p:plain

上記はベースラインに対する各特徴量の性能の平均的な落ち具合を可視化している。 Sepal length や Sepal width に比べると Petal length や Petal width の方が性能が大きく落ちていることが分かる。 ここから、このモデルでは Petal length や Petal width の方が識別に有効に使われていることが推測できる。

緑色のヒゲは、性能の落ち具合のバラつきが正規分布になるという仮定の元での 95% 信頼区間を示している。 実際に正規分布になるかは確認していないし、モデルにもよるだろうけど分散を確認するために描いてみた。

ちなみに、実際に色んなデータやモデルで計測してみると、意外と毎回結果が変わったりすることがある。 これは、そもそもモデルが決定論的に学習しないことだったり、試行回数が少ないことが原因として考えられる。 例えば計測する n-fold CV を何回も繰り返すなど、試行回数を増やすとより結果が安定しやすいかもしれない。

いじょう。

OpenSSH のコネクションが切れにくいように Keepalive を送る

たまに SSH のコネクションが頻繁に切れる環境があるので、定期的にデータを送受信することで切断されるのを防ぎたい。 これは OpenSSH のクライアントであれば ServerAliveInterval を設定することで実現できる。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.13.6
BuildVersion:   17G3025
$ ssh -V
OpenSSH_7.6p1, LibreSSL 2.6.2

コマンドラインからであれば、次のように -o オプションで ServerAliveInterval に秒数を指定しながらリモートに接続する。

$ ssh -o ServerAliveInterval=60 <username>@<remotehost>

毎回コマンドラインで指定するのは面倒なので OpenSSH のコンフィグファイル (典型的には ~/.ssh/config) に記述しても良い。

Host *
  ServerAliveInterval 60

このオプションの詳細については ssh_config のセクション 5 に記載がある。

$ man 5 ssh_config

内容を以下に引用する。 この値が設定されていると、一定の期間にサーバからデータを受信しなかったときは SSH のコネクションを通してメッセージを送る、とある。

     ServerAliveInterval
             Sets a timeout interval in seconds after which if no data has been
             received from the server, ssh(1) will send a message through the
             encrypted channel to request a response from the server.  The default is
             0, indicating that these messages will not be sent to the server.

ちゃんとサーバまでメッセージが届けば、それに対する返答がくる。 それによってコネクションが維持されるという寸法。 ちなみに、メッセージを送信する回数に ServerAliveCountMax で上限を設けることもできるようだ。

いじょう。

OpenSSH[実践]入門 (Software Design plus)

OpenSSH[実践]入門 (Software Design plus)

Python: 実行環境が Jupyter Notebook か判定する

今回は Python の実行環境が Jupyter Notebook か否かを判定する方法について。 ベースのアイデアは以下の Stackoverflow からお借りした。

stackoverflow.com

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.13.6
BuildVersion:   17G3025
$ python -V
Python 3.6.7
$ pip list --format=columns | grep notebook
notebook          5.7.0

下準備

まずは今回使用するパッケージをインストールしておく。

$ pip install notebook ipywidgets tqdm

実行環境が Jupyter Notebook か否かをチェックする関数

早速だけど実行環境をチェックする関数は次の通り。

def is_env_notebook():
    """Determine wheather is the environment Jupyter Notebook"""
    if 'get_ipython' not in globals():
        # Python shell
        return False
    env_name = get_ipython().__class__.__name__
    if env_name == 'TerminalInteractiveShell':
        # IPython shell
        return False
    # Jupyter Notebook
    return True

この関数を使うと、例えば次のように環境ごとにインポートするパッケージや関数を切り替えることができる。

if is_env_notebook():
    # Jupyter Notebook environment
    from tqdm import tqdm_notebook as tqdm
else:
    # Other shells
    from tqdm import tdqm

関数の動作原理について

ここからは関数がどのように動作するかを確認していく。

まずは Jupyter Notebook を起動しておこう。

$ jupyter-notebook

新しい Notebook を作ったら、以下のコードを使ってグローバルスコープのオブジェクトを全て表示してみよう。

import pprint
pprint.pprint(globals())

すると、例えば次のような出力になる。 この中に get_ipython() という関数があることが確認できる。

{'In': ['', 'import pprint\npprint.pprint(globals())'],
 'Out': {},
 '_': '',
 '__': '',
 '___': '',
 '__builtin__': <module 'builtins' (built-in)>,
 '__builtins__': <module 'builtins' (built-in)>,
 '__doc__': 'Automatically created module for IPython interactive environment',
 '__loader__': None,
 '__name__': '__main__',
 '__package__': None,
 '__spec__': None,
 '_dh': ['/Users/amedama/Documents/temporary'],
 '_i': '',
 '_i1': 'import pprint\npprint.pprint(globals())',
 '_ih': ['', 'import pprint\npprint.pprint(globals())'],
 '_ii': '',
 '_iii': '',
 '_oh': {},
 'exit': <IPython.core.autocall.ZMQExitAutocall object at 0x104733c18>,
 'get_ipython': <bound method InteractiveShell.get_ipython of <ipykernel.zmqshell.ZMQInteractiveShell object at 0x104733dd8>>,
 'pprint': <module 'pprint' from '/Users/amedama/.pyenv/versions/3.6.7/lib/python3.6/pprint.py'>,
 'quit': <IPython.core.autocall.ZMQExitAutocall object at 0x104733c18>}

では、次に別のターミナルから通常の Python インタプリタを起動しよう。

$ python

先ほどと同じコードを実行してみる。 すると、こちらには get_ipython() 関数が見当たらない。 その他にも、いくつか Jupyter Notebook の環境ではあったオブジェクトが無いことも分かる。

>>> from pprint import pprint
>>> pprint(globals())
{'__annotations__': {},
 '__builtins__': <module 'builtins' (built-in)>,
 '__doc__': None,
 '__loader__': <class '_frozen_importlib.BuiltinImporter'>,
 '__name__': '__main__',
 '__package__': None,
 '__spec__': None,
 'pprint': <function pprint at 0x10b08de18>}

この特徴を元に、環境を一次切り分けできる。 少なくとも get_ipython() 関数がグローバルスコープになければ Jupyter Notebook の環境ではなさそう。

$ python

...

>>> 'get_ipython' in globals()
False

ただし、これだけだと上手くいかない場合がある。 それは IPython の環境で、IPython のインタプリタを起動するとグローバルスコープに get_ipython() 関数が見つかる。

$ ipython

...

In [1]: 'get_ipython' in globals()                                                                                                          
Out[1]: True

ただし、一つ救いがあって get_ipython() 関数の返り値が Jupyter Notebook とは異なっている。 IPython のインタプリタで返されるのは TerminalInteractiveShell というクラスのインスタンス。

In [2]: get_ipython()                                                                                                                       
Out[2]: <IPython.terminal.interactiveshell.TerminalInteractiveShell at 0x1080a7710>

それに対し Jupyter Notebook では ZMQInteractiveShell というクラスのインスタンスになっている。

$ jupyter-notebook

...

In [1]: 'get_ipython' in globals()
Out[1]: True

In [2]: get_ipython()
Out[2]: <ipykernel.zmqshell.ZMQInteractiveShell at 0x1081f0198>

そのため、クラスの違いを元に環境を識別できる。

いじょう。