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 付近にあるようだ。

以下のサンプルコードでは、上記の関数の最大値を 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プログラミング: Pythonのより良い書き方を学ぶ

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プログラミング: Pythonのより良い書き方を学ぶ

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 個の特徴量が新たに生成されることになる。

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

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

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

いじょう。

下準備

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

$ pip install scikit-learn matplotlib

サンプルデータについて

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

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

#!/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 なので二次元の特徴量が新たに得られる。 それらを散布図にしてみると、ちゃんと線形分離できそうな形で特徴量が生成されている。

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

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

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

課題について

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

いじょう。

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

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 の計測方法になる。 ただ、文章だけだとちょっと分かりにくいかもと思ったので図示してみる。 最初のベースラインを測る部分は一般的な機械学習モデルの汎化性能を計測するやり方と同じ。

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

上記において勘違いしやすいポイントは 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

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

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

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

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

いじょう。

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

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)

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>

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

いじょう。

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

Linux の UNIX パスワード認証について調べた

ふと Linux ディストリビューションのユーザ認証周りについて気になって、その中でも特に UNIX パスワード認証について調べてみた。 UNIX パスワード認証というのは、Linux に限らず Unix 系のディストリビューションで広く採用されているパスワードを使った認証の仕組み。 特に、ログイン用のパスワードを暗号化 (ハッシュ化) した上でパスワードファイル (/etc/passwd) やシャドウパスワードファイル (/etc/shadow) に保存するところが特徴となっている。 UNIX パスワード認証は Unix 系の色々なディストリビューションでそれぞれ実装されている。 その中でも、今回は Linux ディストリビューションの Ubuntu 18.04 LTS について見ていく。 先に断っておくと、かなり長い。

注: この記事では、やっていることがハッシュ化でも manpage の "encrypted" や "encryption" という表現を尊重して「暗号化」と記載している箇所がある

動作確認に使った環境は次の通り。

$ cat /etc/lsb-release
DISTRIB_ID=Ubuntu
DISTRIB_RELEASE=18.04
DISTRIB_CODENAME=bionic
DISTRIB_DESCRIPTION="Ubuntu 18.04.1 LTS"
$ uname -r
4.15.0-29-generic

パスワードファイルとシャドウパスワードファイルについて

太古の昔、ユーザのログインに使うパスワードはパスワードファイル (/etc/passwd) に記述されていたらしい。 しかし、今ではここにパスワードに関する情報が書かれることはほとんどない。 過去に暗号化されたパスワードが記載されていたフィールドには、代わりに x といった文字が載っている。

$ cat /etc/passwd | grep $(whoami)
vagrant:x:1000:1000:vagrant,,,:/home/vagrant:/bin/bash

なぜなら、このファイルの内容は誰でも読めるため。 仮にパスワードが暗号化 (ハッシュ化) されていたとしても、誰からでも読めるのはセキュリティ的には致命的といえる。 そこで、ログインに必要な情報はシャドウパスワードファイル (例えば /etc/shadow) に分離された。

$ ls -l /etc/passwd
-rw-r--r-- 1 root root 1662 Aug 24 09:07 /etc/passwd
$ ls -l /etc/shadow
-rw-r----- 1 root shadow 979 Aug 24 09:07 /etc/shadow

上記を見て分かる通り /etc/shadow は root ユーザか shadow グループのメンバー以外からは閲覧できない。

シャドウパスワードファイルにおいて、カンマで区切られた二つ目のカラムがユーザがログインに用いるパスワードを暗号化した文字列になる。 これは、ソルト付きのハッシュ化されたパスワードになっている。

$ sudo cat /etc/shadow | grep $(whoami)
vagrant:$6$xUMKs3vr$r/rGFiHf6.B0xsKDpustA2G9m1zGcWJnjZLMlyQeh.VDEpDEMEQ3AKSvrRv.NQEc2OBbban3BAHj4Nhwpa0t4.:17767:0:99999:7:::

ユーザがシステムにログインするときは、ユーザが入力した内容と上記の内容を突合して一致するかを確認することになる。

シャドウパスワードファイルの詳細については man コマンドで調べることができる。 man のセクション 5 には設定ファイルに関する情報が記載されている。

$ man 5 shadow

上記の内容を確認していくと、次のような記述がある。

       encrypted password
           Refer to crypt(3) for details on how this string is interpreted.

どうやらファイルに記載されている暗号化されたパスワードは crypt(3) に関連しているようだ。 man のセクション 3 は開発者向けのライブラリ関数に関する情報なので、それ用の manpage をインストールした上で確認する。

$ sudo apt-get -y install manpages-dev
$ man 3 crypt

上記を見ると暗号化されたパスワードは次のようなフォーマットになっていることが分かる。 id がハッシュアルゴリズムで salt がハッシュするときに使うソルト、そして最後の encrypted がハッシュ化されたパスワードになる。

$id$salt$encrypted

id については次のような記述がある。 使えるハッシュアルゴリズムは環境や glibc のバージョンによって異なるようだ。 glibc というのは C 言語の標準ライブラリである libc の、GNU Project による実装を指している。

       id identifies the encryption method used instead of DES and  this  then
       determines  how  the  rest  of the password string is interpreted.  The
       following values of id are supported:

              ID  | Method
              ─────────────────────────────────────────────────────────
              1   | MD5
              2a  | Blowfish (not in mainline glibc; added in some
                  | Linux distributions)
              5   | SHA-256 (since glibc 2.7)
              6   | SHA-512 (since glibc 2.7)

C 言語の標準ライブラリというのは、例えば多くの入門書で「おまじない」として紹介されている #include <stdio.h> なんかがそれ。

先ほど確認した自分自身のパスワードハッシュでは id が 6 になっていた。 つまり、ハッシュアルゴリズムとして SHA-512 が使われていることが分かる。

$ sudo cat /etc/shadow | grep $(whoami)
vagrant:$6$xUMKs3vr$r/rGFiHf6.B0xsKDpustA2G9m1zGcWJnjZLMlyQeh.VDEpDEMEQ3AKSvrRv.NQEc2OBbban3BAHj4Nhwpa0t4.:17767:0:99999:7:::

crypt(3) を使ってみよう

それでは、実際に crypt(3) を使ってみることにしよう。 ただ、C 言語で書くのではなく Python のラッパーである crypt モジュールを使うことにした。

まずは Python のインタプリタを起動する。

$ python3

先ほど確認した自分自身の暗号化されたパスワードを再掲する。 ハッシュアルゴリズムは SHA-512 でソルトが xUMKs3vr だった。 そして、使われているパスワードは Vagrant の環境なので vagrant になっている。

$ sudo cat /etc/shadow | grep $(whoami)
vagrant:$6$xUMKs3vr$r/rGFiHf6.B0xsKDpustA2G9m1zGcWJnjZLMlyQeh.VDEpDEMEQ3AKSvrRv.NQEc2OBbban3BAHj4Nhwpa0t4.:17767:0:99999:7:::

上記の内容にもとづいて crypt モジュールの関数を呼び出してみよう。 ソルトについては暗号化されたパスワードの $ で区切られた部分をそのまま使って構わない。

>>> import crypt
>>> crypt.crypt('vagrant', '$6$xUMKs3vr$')
'$6$xUMKs3vr$r/rGFiHf6.B0xsKDpustA2G9m1zGcWJnjZLMlyQeh.VDEpDEMEQ3AKSvrRv.NQEc2OBbban3BAHj4Nhwpa0t4.'

見事に自分自身の暗号化されたパスワードと同じ文字列が生成できた。

ちなみに、システムで使えるハッシュアルゴリズムは crypt.methods で確認できる。

>>> crypt.methods
[<crypt.METHOD_SHA512>, <crypt.METHOD_SHA256>, <crypt.METHOD_MD5>, <crypt.METHOD_CRYPT>]

上記を見て分かる通り id が 2a に対応するハッシュアルゴリズムの Blowfish はこのシステムでは使えないらしい。

crypt(3) の実装を見てみよう

次は、念のため crypt(3) を実装している場所についても調べておくことにしよう。 特に興味がなければ、次のセクションまで飛ばしてもらって構わない。

先ほど見た通り、実装しているのは glibc と思われる。 そこで、まずは今使っているシステムの glibc のバージョンを調べよう。 検索パスに入っている共有ライブラリを ldconfig コマンドで調べたら、そこから libc を探す。

$ ldconfig -p | grep libc.so
    libc.so.6 (libc6,x86-64, OS ABI: Linux 3.2.0) => /lib/x86_64-linux-gnu/libc.so.6

あとは共有ライブラリをそのまま実行するとバージョンが確認できる。 どうやら、この環境ではバージョン 2.27 が使われているようだ。 先ほど crypt(3) のハッシュアルゴリズムとして SHA-512 が使えるのはバージョン 2.7 以降とあったので、その内容とも整合している。

$ /lib/x86_64-linux-gnu/libc.so.6
GNU C Library (Ubuntu GLIBC 2.27-3ubuntu1) stable release version 2.27.
Copyright (C) 2018 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.
There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.
Compiled by GNU CC version 7.3.0.
libc ABIs: UNIQUE IFUNC
For bug reporting instructions, please see:
<https://bugs.launchpad.net/ubuntu/+source/glibc/+bugs>.

続いて公式サイトの記述に沿ってソースコードを取得する。

$ sudo apt-get -y install git
$ git clone git://sourceware.org/git/glibc.git
$ cd glibc/
$ git checkout release/2.27/master
Branch 'release/2.27/master' set up to track remote branch 'release/2.27/master' from 'origin'.
Switched to a new branch 'release/2.27/master'

とりあえず crypt(3) で検索してみると、それっぽいものが見つかる。

$ grep -r "crypt(3)" .
./crypt/ufc.c: * UFC-crypt: ultra fast crypt(3) implementation
./crypt/crypt-entry.c: * UFC-crypt: ultra fast crypt(3) implementation
./crypt/cert.c: * This crypt(3) validation program shipped with UFC-crypt
./crypt/crypt.c: * UFC-crypt: ultra fast crypt(3) implementation
./crypt/crypt_util.c: * UFC-crypt: ultra fast crypt(3) implementation
./crypt/README.ufc-crypt:- Crypt implementation plugin compatible with crypt(3)/fcrypt.
./crypt/README.ufc-crypt:- On most machines, UFC-crypt runs 30-60 times faster than crypt(3) when
./crypt/README.ufc-crypt:  that of crypt(3).
./crypt/README.ufc-crypt:such as the one said to be used in crypt(3).
./crypt/crypt-private.h: * UFC-crypt: ultra fast crypt(3) implementation
./crypt/crypt.h: * UFC-crypt: ultra fast crypt(3) implementation

ヘッダファイルに crypt() 関数らしき宣言があったり、各ハッシュアルゴリズムの実装も見受けられる。 関数のエントリポイントは crypt-entry.c にあって、そこから各実装に分岐しているみたい。

$ grep extern crypt/crypt.h | head -n 1
extern char *crypt (const char *__key, const char *__salt)
$ ls crypt | egrep "(md5|sha256|sha512)\.c"
md5.c
sha256.c
sha512.c

たしかに crypt(3) は glibc に実装があることが分かった。

crypt(3) は何処から使われているのか？

続いて気になるのは、暗号化されたパスワードを生成する crypt(3) が何処から使われているのかという点。 おそらくユーザがログインするタイミングで入力されたパスワードに対して crypt(3) を使っているはず。 そして、使った結果とシャドウパスワードファイルの内容を比較して、一致すればログイン成功という処理が何処かにあると考えられる。

一般的なシステムであれば、ユーザがログインする経路として代表的なものが二つ考えられる。 まず一つ目が、コンピュータに直接接続された物理的なディスプレイやキーボードを通してユーザ名とパスワードを入力するパターン。 この場合、端末デバイスには制御端末 (tty: teletypewriter) が割り当てられる。 そして二つ目が、SSH クライアントを通してネットワーク越しにユーザ名やパスワード、あるいは公開鍵の情報を渡すパターン。 この場合、端末デバイスには擬似端末 (pty: pseudo terminal) が割り当てられる。

制御端末 (tty: teletypewriter) からログインするパターン

制御端末 (tty) を通してログインするパターンから見ていく。 このパターンでは、まず Linux ディストリビューションの初期化プロセスから考える必要がある。 Ubuntu 18.04 LTS では init プロセス (pid に 1 を持つプロセス) として systemd が用いられる。

そして systemd は agetty(8) というプログラムを起動する。 このプログラムはログインプロンプトの表示部分を担当している。 システムのプロセスを確認すると、制御端末が割り当てられた状態で agetty(8) が起動していることが分かる。

$ ps auxww | grep [a]getty
root       678  0.0  0.1  16180  1656 tty1     Ss+  16:00   0:00 /sbin/agetty -o -p -- \u --noclear tty1 linux
$ pstree | head -n 3
systemd-+-VBoxService---7*[{VBoxService}]
        |-accounts-daemon---2*[{accounts-daemon}]
        |-agetty

このプログラムに関する詳細は man ページのセクション 8 の agetty に記載されている。

$ man 8 agetty

上記の内容を見ると、このプログラムは /bin/login を呼び出すとある。 どうやら、このプログラムが実際のログイン部分を担当しているようだ。

DESCRIPTION
       agetty  opens  a  tty port, prompts for a login name and invokes the /bin/login command.  It is normally invoked by
       init(8).

ということで、次はこの login(1) のソースコードが読みたい。 なので、まずはこのコマンドがどのパッケージに含まれているのかを dpkg コマンドで確認する。

$ dpkg -S /bin/login
login: /bin/login

どうやら、そのものズバリ login というパッケージに含まれているらしい。

パッケージ名が分かったので apt-get source コマンドを使ってソースコードを取り寄せる。 すると、どうやら login(1) の実体は shadow というプログラム群に含まれることが分かった。

$ sudo apt-get -y install dpkg-dev
$ sudo sed -i.back -e "s/^# deb-src/deb-src/g" /etc/apt/sources.list
$ sudo apt-get update
$ apt-get source login
Reading package lists... Done
Picking 'shadow' as source package instead of 'login'
NOTICE: 'shadow' packaging is maintained in the 'Git' version control system at:
https://anonscm.debian.org/git/pkg-shadow/shadow.git
Please use:
git clone https://anonscm.debian.org/git/pkg-shadow/shadow.git
to retrieve the latest (possibly unreleased) updates to the package.
...

ドキュメントを見ると、このプログラム群が Linux のシャドウパスワードを扱うためのものだということが分かる。 どうやら核心に近づいてきたようだ。

$ cd shadow-4.5/
$ head doc/HOWTO -n 6 | tail -n 2
  Linux Shadow Password HOWTO
  Michael H. Jackson, mhjack@tscnet.com

中身を確認すると src というディレクトリにコマンドラインから叩くプログラムがまとまっていることが分かった。 例えば login(1) もここにあるし、他には usermod(8) や groupmod(8) なんかも見つかる。

$ ls src/ | grep login
login.c
login_nopam.c
nologin.c
sulogin.c
$ ls src/ | grep usermod
usermod.c
$ ls src/ | grep groupmod
groupmod.c

(2018/11/5 追記: 以降の説明はシステムで PAM を用いない場合の説明になっているため、一般的な環境の説明と異なっている。例えば、後述する関数が login コマンドにおいて呼び出されるフローは、ビルドする際に USE_PAM フラグを無効にした場合にしか通らない。)

改めて src/login.c を読み進めていくと、ログインパスワードを pw_auth() という関数に突っ込んでいる。 この関数は lib/pwauth.c で定義されていて、さらに lib/encrypt.c の pw_encrypt() という関数を呼び出している。 そして、この pw_encrypt() 関数の中に、探し求めていた crypt(3) の呼び出しが見つかる！ ついに見つけることができた。

$ grep "crypt (" lib/encrypt.c 
/*@exposed@*//*@null@*/char *pw_encrypt (const char *clear, const char *salt)
    cp = crypt (clear, salt);

さらに src/login.c 周辺を読み進めると crypt(3) で暗号化したパスワードは getpwnam(3) や getspnam(3) と比較していることが分かった。 この関数は /etc/passwd や /etc/shadow を読み取って、そこに記載されている暗号化されたパスワードを読み取るためのものだ。 関数を提供しているのは、まさにこの shadow で、ソースコードの実体は lib 以下に存在する。

それぞれの関数の情報は man コマンドから参照できる。

$ man 3 getspnam

その後、認証に成功するとパスワードファイルから読み取ったログインシェルを起動することもソースコードから確認できる。

$ grep -r shell login.c
    if (pwd->pw_shell[0] == '*') { /* subsystem root */
        pwd->pw_shell++;   /* skip the '*' */
        err = shell (tmp, pwd->pw_shell, newenvp); /* fake shell */
        /* exec the shell finally */
        err = shell (pwd->pw_shell, (char *) 0, newenvp);

これで Ubuntu 18.04 LTS におけるログインの流れをソースコードから確認できた。

• ユーザからの入力パスワードを暗号化する
• 上記をシャドウパスワードファイルに記載されている暗号化されたパスワードと比較する
• 一致するときはログインシェルを起動する

実際に制御端末 (tty) からログインした状態で pstree を確認すると agetty(8) がいなくなって、代わりに login(1) 経由で bash(1) が起動している。

$ pstree
systemd─┬─VBoxService───7*[{VBoxService}]
...
        ├─login───bash

ちなみにシステムの /sbin/agetty や /bin/login をリネームすると、なかなか面白い挙動になるので試してみるのもおすすめ。 例えば agetty が見つからない状態でシステムが起動するとログインプロンプトが表示されないことを確認できる。

擬似端末 (pty: pseudo terminal) からログインするパターン

続いては SSH クライアントを通してシステムにログインするパターン。 これは単純で、ようするに SSH デーモンに認証部分があるだろうことは容易に想像がつく。

なので、まずは sshd(8) の含まれるパッケージを確認する。

$ dpkg -S $(which sshd)
openssh-server: /usr/sbin/sshd

先ほどと同じ要領で OpenSSH のソースコードを取り寄せる。

$ sudo apt-get source openssh-server
Reading package lists... Done
Picking 'openssh' as source package instead of 'openssh-server'
NOTICE: 'openssh' packaging is maintained in the 'Git' version control system at:
https://salsa.debian.org/ssh-team/openssh
Please use:
git clone https://salsa.debian.org/ssh-team/openssh
to retrieve the latest (possibly unreleased) updates to the package.
...
$ cd openssh-7.6p1/

この中では、例えば crypt(3) を xcrypt() 関数というプラットフォーム非依存な形でラップしていたりする。

$ grep " crypt(" */*.c
openbsd-compat/xcrypt.c:        crypted = crypt(password, salt);
openbsd-compat/xcrypt.c:        crypted = crypt(password, salt);
openbsd-compat/xcrypt.c:    crypted = crypt(password, salt);

あるいは各所で getspnam(3) や getpwnam(3) が使われている。

$ grep "getspnam" *.c
auth.c:     spw = getspnam(pw->pw_name);
auth-shadow.c:  if ((spw = getspnam((char *)user)) == NULL) {
$ grep "getspnam" */*.c
openbsd-compat/xcrypt.c:    struct spwd *spw = getspnam(pw->pw_name);

ということで、こちらの経路に関しても crypt(3) や getspnam(3) を駆使してユーザを認証しているらしいことが分かった。 プロセスを確認しても、ユーザがログインした状態では sshd(8) 経由で bash(1) が起動していることが分かる。

$ pstree | grep sshd
        |-sshd---sshd---sshd---bash-+-grep

まとめ

分かったことは次の通り。

• シャドウパスワードファイル (例えば /etc/shadow) には暗号化されたユーザのログインパスワードが書き込まれている
• 暗号化されたパスワードは crypt(3) で生成する
• シャドウパスワードファイルに書き込まれている情報は getspnam(3) で読み取る
• ユーザの入力内容から生成したパスワードハッシュと、シャドウパスワードファイルに書き込まれた内容を比較することで認証する

いじょう。

ふつうのLinuxプログラミング 第2版 Linuxの仕組みから学べるgccプログラミングの王道

詳解UNIXプログラミング 第3版