CUBE SUGAR CONTAINER

技術系のこと書きます。

Python: MLflow Models の Custom Python Models でデータを Pickle 以外に永続化する

以前、このブログでは MLflow Models の使い方について以下のようなエントリを書いた。 この中では、Custom Python Models を作るときに、データを Python の Pickle 形式のファイルとして永続化していた。 今回は、それ以外のファイルにデータを永続化する方法について書いてみる。

blog.amedama.jp

使った環境は次のとおり。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.15.7
BuildVersion:   19H2
$ python -V
Python 3.8.5
$ pip list | grep -i mlflow
mlflow                    1.11.0

もくじ

下準備

まずは、下準備として MLflow をインストールしておく。

$ pip install mlflow

Pickle 形式のファイルにモデルを永続化する

まずは、おさらいとして Custom Python Models を作るときに、モデルのデータを Pickle 形式のファイルに永続化する方法から。

以下にサンプルコードを示す。 この中では、定数を入力に加えるだけのモデルとして AddN というクラスを定義している。 このクラスは mlflow.pyfunc.PythonModel を継承しているため、mlflow.pyfunc.save_model()mlflow.pyfunc.load_model() を使ってファイルに読み書きできる。 サンプルコードでは、実際に定数として 5 を加える設定にしたインスタンスをファイルに永続化している。

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

from mlflow import pyfunc


class AddN(pyfunc.PythonModel):
    """指定した数を入力に加えるモデル"""

    def __init__(self, n):
        self.n = n

    def predict(self, context, model_input):
        return model_input.apply(lambda column: column + self.n)


def main():
    # Python の Pickle ファイルとしてモデルを永続化する
    model_path = 'add-n-pickle'
    add5_model = AddN(5)
    pyfunc.save_model(path=model_path,
                      python_model=add5_model)



if __name__ == '__main__':
    main()

上記のポイントは、mlflow.pyfunc.save_model()python_model にインスタンスを指定しているだけというところ。 この場合、永続化されるのはインスタンスを Pickle 形式で直列化したファイルになる。

上記を実行してみよう。

$ python saveaddnpkl.py

実行すると、以下のように MLflow Models のフォーマットに沿ってディレクトリができる。 この中で python_model.pkl が前述した AddN クラスのインスタンスを表す Pickle 形式のファイルになる。

$ ls add-n-pickle 
MLmodel         conda.yaml      python_model.pkl
$ python -m pickle add-n-pickle/python_model.pkl
<__main__.AddN object at 0x1048cdd90>

上記を Python から読み込んで使うサンプルコードも次のように用意した。

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

import pandas as pd
from mlflow import pyfunc


def main():
    model_path = 'add-n-pickle'
    loaded_model = pyfunc.load_model(model_path)

    x = pd.DataFrame(list(range(10, 21)), columns=['n'])
    y = loaded_model.predict(x)

    print(f'Input: {x.n.values}')
    print(f'Output: {y.n.values}')


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python loadaddnpkl.py
Input: [10 11 12 13 14 15 16 17 18 19 20]
Output: [15 16 17 18 19 20 21 22 23 24 25]

ちゃんと、読み込んだモデルが入力データに定数として 5 を加えていることが確認できる。

テキストファイルにインスタンスのアトリビュートを永続化してみる

それでは、続いて Pickle 以外のフォーマットのファイルも使って永続化するパターンを扱う。 これは、たとえば永続化したいモデルが Pickle 以外のフォーマットでファイルに読み書きする API がある場合などに使い勝手が良い。 実際に MLflow Models と XGBoost や LightGBM のインテグレーションは、フレームワークが提供する永続化用 API を流用して書かれている。

以下にサンプルコードを示す。 今回は、入力に加算する定数をテキストファイルとして、AddN クラスのインスタンスとは別に永続化するやり方を取った。 また、先ほどはファイルへの書き込みと読み込みを別の Python モジュールに分けたのに対して、これはひとつで完結させている。

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

import tempfile
import os

import pandas as pd
from mlflow import pyfunc


class AddN(pyfunc.PythonModel):
    """指定した数を入力に加えるモデル"""

    def __init__(self, n=None):
        self.n = n

    def load_context(self, context: pyfunc.PythonModelContext):
        """インスタンスの状態を復元するのに使われるコールバック"""
        # アーティファクトのパスを取り出す
        artifact_path = context.artifacts['n_file']
        with open(artifact_path, mode='r') as fp:
            # ファイルからモデルのパラメータを復元する
            self.n = int(fp.read())

    def predict(self,
                context: pyfunc.PythonModelContext,
                model_input: pd.DataFrame):
        return model_input.apply(lambda column: column + self.n)

    def __repr__(self):
        """インスタンスの文字列表現を得るときに呼ばれる特殊メソッド"""
        return f'<AddN n:{self.n}>'


def save_model(n: int, path: str):
    """モデルをアーティファクトに永続化するためのユーティリティ関数"""
    # モデルが動作するのに必要なパラメータをテキストファイルなどにアーティファクトとして書き出す
    filename = 'n.txt'
    with tempfile.TemporaryDirectory() as d:
        # 一時ディレクトリを用意して、そこにファイルを作る
        artifact_path = os.path.join(d, filename)
        with open(artifact_path, mode='w') as fp:
            fp.write(str(n))
        # 上記で作ったファイルをアーティファクトとして記録する
        # NOTE: path 以下にファイルがコピーされる
        artifacts = {
            'n_file': artifact_path,
        }
        pyfunc.save_model(path=path,
                          # このインスタンスに load_context() メソッド経由でパラメータが読み込まれる
                          python_model=AddN(),
                          artifacts=artifacts,
                          )


def main():
    # モデルをアーティファクトに永続化する
    model_path = 'add-n-artifacts'
    save_model(5, model_path)

    # モデルをアーティファクトから復元する
    loaded_model = pyfunc.load_model(model_path)

    # 動作を確認する
    x = pd.DataFrame(list(range(10, 21)), columns=['n'])
    y = loaded_model.predict(x)

    print(f'Input: {x.n.values}')
    print(f'Output: {y.n.values}')


if __name__ == '__main__':
    main()

上記にはポイントがいくつかある。 まず、書き込みに関しては mlflow.pyfunc.save_model() を呼ぶ際に artifacts というオプションを指定している。 これには、モデルのデータなどを記録したファイルへのパスを Python の辞書として渡す。 ちなみに、ここに指定するパスは実際には Artifact URI なので、リモートのストレージにあっても構わない。 この点は、おそらく MLflow Tracking と組み合わせて使うことを想定しているのだと思う。 もちろん、ファイルはあらかじめ、そのパス (繰り返しになるけど、実際には Artifact URI) に存在している必要がある。 ここに指定したファイルは、MLflow Models が作成するディレクトリへとコピーされる。 そして、読み込みに関しては mlflow.pyfunc.PythonModel を継承したクラスに load_context() というメソッドを実装する。 このメソッドは、インスタンスを Pickle 形式のファイルから読み込んだ後に呼ばれるようだ。 つまり、Pickle 以外のファイルにデータを保存しているときに、モデルの状態をそれで更新するときに使うことができる。

さて、前置きが長くなったけど実際に上記のサンプルコードを実行してみよう。

$ python addnartifacts.py
Input: [10 11 12 13 14 15 16 17 18 19 20]
Output: [15 16 17 18 19 20 21 22 23 24 25]

ちゃんと想定どおりの入出力になっている。

MLflow Models が作成したディレクトリを確認してみよう。 すると、artifacts というサブディレクトリがあることに気づく。

$ ls add-n-artifacts 
MLmodel         artifacts       conda.yaml      python_model.pkl

中を見ると、インスタンスのアトリビュートがテキストファイルとして書き込まれている。

$ cat add-n-artifacts/artifacts/n.txt   
5

それ以外には、アトリビュートが None の状態の AddN クラスのインスタンスが Pickle 形式で永続化されていることがわかる。

$ python -m pickle add-n-artifacts/python_model.pkl
<AddN n:None>

いじょう。

Python: MLflow Models を使ってみる

MLflow は MLOps に関連した OSS のひとつ。 いくつかのコンポーネントに分かれていて、それぞれを必要に応じて独立して使うことができる。

その中でも、今回扱う MLflow Models は主に学習済みモデルやパイプラインの取り回しに関するコンポーネント。 MLflow Models を使うことで、たとえば学習済みモデルの Serving やシステムへの組み込みが容易になる可能性がある。

使った環境は次のとおり。

$ sw_vers                     
ProductName:    Mac OS X
ProductVersion: 10.15.6
BuildVersion:   19G73
$ python -V                                        
Python 3.8.5
$ pip list | egrep "(mlflow|lightgbm|scikit-learn)"
lightgbm                  3.0.0
mlflow                    1.11.0
scikit-learn              0.23.2

もくじ

下準備

はじめに、必要なパッケージをインストールしておく。

$ pip install mlflow lightgbm scikit-learn seaborn category_encoders

モデルを MLflow Models で永続化する

論よりコードということで、いきなりだけど以下にサンプルコードを示す。 このサンプルコードでは Boston データセットを LightGBM で学習するコードになっている。 そして、学習させたモデルを MLflow Models を使って永続化している。

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

import lightgbm as lgb
from mlflow import lightgbm as mlflow_lgb
from sklearn import datasets
from sklearn.model_selection import train_test_split


def main():
    # Boston データセットを読み込む
    dataset = datasets.load_boston()
    train_x, train_y = dataset.data, dataset.target
    feature_names = list(dataset.feature_names)

    # 学習用データと検証用データに分割する
    x_tr, x_eval, y_tr, y_eval = train_test_split(train_x, train_y,
                                                  test_size=0.33,
                                                  shuffle=True,
                                                  random_state=42,
                                                  )
    # LightGBM のデータ形式に直す
    lgb_train = lgb.Dataset(x_tr, y_tr, feature_name=feature_names)
    lgb_eval = lgb.Dataset(x_eval, y_eval, reference=lgb_train)

    # モデルを学習する
    lgbm_params = {
        'objective': 'regression',
        'metric': 'rmse',
        'first_metric_only': True,
        'verbose': -1,
        'random_state': 42,
    }
    booster = lgb.train(params=lgbm_params,
                        train_set=lgb_train,
                        valid_sets=[lgb_train, lgb_eval],
                        num_boost_round=1_000,
                        early_stopping_rounds=100,
                        verbose_eval=50,
                        )

    # モデルを MLflow Models の形式で永続化する
    mlflow_lgb.save_model(booster, path='mlflow-lgbm')
    """
    # MLflow Tracking に残すならこうする
    with mlflow.start_run():
        mlflow_lgb.log_model(booster,
                             artifact_path='mlflow-lgbm')
    """


if __name__ == '__main__':
    main()

上記のモジュールを実行してみよう。

$ python lgbmlf.py
Training until validation scores don't improve for 100 rounds
[50]  training's rmse: 2.23963    valid_1's rmse: 3.59161
[100] training's rmse: 1.55562    valid_1's rmse: 3.4141
[150] training's rmse: 1.22661    valid_1's rmse: 3.36079
[200] training's rmse: 1.0165 valid_1's rmse: 3.35222
[250] training's rmse: 0.860022   valid_1's rmse: 3.34358
[300] training's rmse: 0.735403   valid_1's rmse: 3.34575
[350] training's rmse: 0.630665   valid_1's rmse: 3.35982
Early stopping, best iteration is:
[254] training's rmse: 0.848352   valid_1's rmse: 3.33877
Evaluated only: rmse

すると、次のようにディレクトリができる。 このディレクトリと中身のファイルが、MLflow Models を使って永続化したモデルを表している。 要するに、決められたフォーマットに沿って学習済みモデルをパッケージングしている。

$ ls mlflow-lgbm 
MLmodel     conda.yaml  model.lgb

この中で特に重要なのが MLmodel という YAML フォーマットで書かれたファイル。 このファイルには、そのモデルがどのように永続化されたかといった情報が記録されている。

$ cat mlflow-lgbm/MLmodel 
flavors:
  lightgbm:
    data: model.lgb
    lgb_version: 3.0.0
  python_function:
    data: model.lgb
    env: conda.yaml
    loader_module: mlflow.lightgbm
    python_version: 3.8.5
utc_time_created: '2020-09-30 09:44:55.890106'

なお、上記のフォーマットの詳細は次のドキュメントに記載されている。

www.mlflow.org

www.mlflow.org

また、conda.yaml というファイルには Conda の仮想環境に関する情報が記録されている。 これはつまり、永続化したモデルを利用するために必要な Conda の環境を構築するたのもの。 MLflow Models では、デフォルトで Conda の仮想環境上に学習済みモデルをデプロイすることを想定している。

たとえば、中身を見ると LightGBM が依存パッケージとして追加されていることがわかる。

$ cat mlflow-lgbm/conda.yaml 
channels:
- defaults
- conda-forge
dependencies:
- python=3.8.5
- pip
- pip:
  - mlflow
  - lightgbm==3.0.0
name: mlflow-env

永続化したモデルを使って推論用の REST API を立ち上げる

ここからは MLflow Models を使うことで得られる嬉しさについて書いていく。 MLflow には、MLflow Models で永続化したモデルを扱うための機能がいくつか用意されている。

たとえば、MLflow には mlflow というコマンドラインが用意されている。 このコマンドの models serve サブコマンドを使うと、学習済みモデルを使った推論用の REST API が気軽に立てられる。

実際に使ってみよう。 コマンドを実行する際に、--model-uri オプションには、先ほど永続化したディレクトリを指定する。 また、今回は Conda を使っていないので --no-conda オプションをつけた。 これで、デフォルトでは localhost の 5000 番ポートで推論用の API が立ち上がる

$ mlflow models serve \
    --no-conda \
    --model-uri mlflow-lgbm
2020/09/30 18:49:46 INFO mlflow.models.cli: Selected backend for flavor 'python_function'
2020/09/30 18:49:46 INFO mlflow.pyfunc.backend: === Running command 'gunicorn --timeout=60 -b 127.0.0.1:5000 -w 1 ${GUNICORN_CMD_ARGS} -- mlflow.pyfunc.scoring_server.wsgi:app'
[2020-09-30 18:49:47 +0900] [22853] [INFO] Starting gunicorn 20.0.4
[2020-09-30 18:49:47 +0900] [22853] [INFO] Listening at: http://127.0.0.1:5000 (22853)
[2020-09-30 18:49:47 +0900] [22853] [INFO] Using worker: sync
[2020-09-30 18:49:47 +0900] [22855] [INFO] Booting worker with pid: 22855

上記に推論させたいデータを HTTP で投げ込んでみよう。 たとえば curl コマンドを使って以下のようにする。

$ curl -X POST \
    -H "Content-Type:application/json" \
    --data '{"columns": ["CRIM", "ZN", "INDUS", "CHAS", "NOX", "RM", "AGE", "DIS", "RAD", "TAX", "PTRATIO", "B", "LSTAT"], "data": [[10.233, 0.0, 18.1, 0.0, 0.614, 6.185, 96.7, 2.1705, 24.0, 666.0, 20.2, 379.7, 18.03]]}' \
    http://localhost:5000/invocations
[15.444706764627714]

すると、推論の結果として 15.44... という結果が得られた。

永続化したモデルを使って CSV ファイルを処理する

また、同様に CSV のファイルを処理することもできる。 さっきと同じ内容を CSV ファイルに記録してみよう。

$ cat << 'EOF' > data.csv
CRIM,ZN,INDUS,CHAS,NOX,RM,AGE,DIS,RAD,TAX,PTRATIO,B,LSTAT
10.233,0.0,18.1,0.0,0.614,6.185,96.7,2.1705,24.0,666.0,20.2,379.7,18.03
EOF

今度は models predict というサブコマンドを使う。 --content-type オプションには csv を指定する。 そして、--input-path オプションに先ほど保存した CSV ファイルを指定する。

$ mlflow models predict \
    --no-conda \
    --model-uri mlflow-lgbm/ \
    --content-type csv \
    --input-path data.csv
2020/09/30 18:51:50 INFO mlflow.models.cli: Selected backend for flavor 'python_function'
[15.444706764627714]

先ほどと同じように、推論結果として 15.44... という値が得られた。

ただ、現状のままだと上手くいかない場面もある。 たとえば、CSV のカラムを一部入れ替えてみよう。 以下では CRIM カラムと ZN カラムの順番が入れ替わっている。

$ cat << 'EOF' > data.csv
ZN,CRIM,INDUS,CHAS,NOX,RM,AGE,DIS,RAD,TAX,PTRATIO,B,LSTAT
0.0,10.233,18.1,0.0,0.614,6.185,96.7,2.1705,24.0,666.0,20.2,379.7,18.03
EOF

このファイルを使ってもう一度同じことをしてみよう。 ちゃんとカラム名まで認識していれば結果は変わらないはず。

$ mlflow models predict \
    --no-conda \
    --model-uri mlflow-lgbm/ \
    --content-type csv \
    --input-path data.csv
2020/09/30 18:52:33 INFO mlflow.models.cli: Selected backend for flavor 'python_function'
[16.96720478471085]

しかし、残念ながら結果は変わってしまった。 つまり、先ほどのサンプルコードではカラム名の情報までは永続化できていない。

Signature を追加する

カラム名まで認識してほしいときは、モデルを永続化する際に Signature という情報を追加する必要がある。

以下にサンプルコードを示す。 先ほどのサンプルコードに、Pandas の DataFrame から自動的に Signature を認識させるコードを追加している。

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

import pandas as pd
import lightgbm as lgb
from mlflow import lightgbm as mlflow_lgb
from mlflow.models.signature import infer_signature
from sklearn import datasets
from sklearn.model_selection import train_test_split


def main():
    dataset = datasets.load_boston()
    train_x, train_y = dataset.data, dataset.target
    feature_names = list(dataset.feature_names)

    x_tr, x_eval, y_tr, y_eval = train_test_split(train_x, train_y,
                                                  test_size=0.33,
                                                  shuffle=True,
                                                  random_state=42,
                                                  )

    lgb_train = lgb.Dataset(x_tr, y_tr, feature_name=feature_names)
    lgb_eval = lgb.Dataset(x_eval, y_eval, reference=lgb_train)

    lgbm_params = {
        'objective': 'regression',
        'metric': 'rmse',
        'first_metric_only': True,
        'verbose': -1,
        'random_state': 42,
    }
    booster = lgb.train(params=lgbm_params,
                        train_set=lgb_train,
                        valid_sets=[lgb_train, lgb_eval],
                        num_boost_round=1_000,
                        early_stopping_rounds=100,
                        verbose_eval=50,
                        )

    y_tr_pred = booster.predict(x_tr, num_iteration=booster.best_iteration)
    # 入力が DataFrame であれば、場合によってはカラム名とデータ型を自動で認識してくれる
    x_tr_df = pd.DataFrame(x_tr, columns=dataset.feature_names)
    signature = infer_signature(x_tr_df, y_tr_pred)
    # 渡すデータと推論の結果を Signature として付与する
    mlflow_lgb.save_model(booster, path='mlflow-lgbm-with-sig', signature=signature)


if __name__ == '__main__':
    main()

上記を実行しよう。

$ python lgbmlfsig.py
Training until validation scores don't improve for 100 rounds
[50]  training's rmse: 2.23963    valid_1's rmse: 3.59161
[100] training's rmse: 1.55562    valid_1's rmse: 3.4141
[150] training's rmse: 1.22661    valid_1's rmse: 3.36079
[200] training's rmse: 1.0165 valid_1's rmse: 3.35222
[250] training's rmse: 0.860022   valid_1's rmse: 3.34358
[300] training's rmse: 0.735403   valid_1's rmse: 3.34575
[350] training's rmse: 0.630665   valid_1's rmse: 3.35982
Early stopping, best iteration is:
[254] training's rmse: 0.848352   valid_1's rmse: 3.33877
Evaluated only: rmse

今度は保存された MLmodel ファイルに signature という情報が付与されている。 中身を見るとカラム名とデータ型が入っている。

$ ls mlflow-lgbm-with-sig
MLmodel     conda.yaml  model.lgb
$ cat mlflow-lgbm-with-sig/MLmodel 
flavors:
  lightgbm:
    data: model.lgb
    lgb_version: 3.0.0
  python_function:
    data: model.lgb
    env: conda.yaml
    loader_module: mlflow.lightgbm
    python_version: 3.8.5
signature:
  inputs: '[{"name": "CRIM", "type": "double"}, {"name": "ZN", "type": "double"},
    {"name": "INDUS", "type": "double"}, {"name": "CHAS", "type": "double"}, {"name":
    "NOX", "type": "double"}, {"name": "RM", "type": "double"}, {"name": "AGE", "type":
    "double"}, {"name": "DIS", "type": "double"}, {"name": "RAD", "type": "double"},
    {"name": "TAX", "type": "double"}, {"name": "PTRATIO", "type": "double"}, {"name":
    "B", "type": "double"}, {"name": "LSTAT", "type": "double"}]'
  outputs: '[{"type": "double"}]'
utc_time_created: '2020-09-30 09:58:18.952375'

それでは、Signature を追加したモデルで推論させてみよう。 CSV ファイルは先ほどと同じものを使う。 つまり、モデルの学習時と推論時でカラムの順番が入れかわっている。

$ cat data.csv                 
ZN,CRIM,INDUS,CHAS,NOX,RM,AGE,DIS,RAD,TAX,PTRATIO,B,LSTAT
0.0,10.233,18.1,0.0,0.614,6.185,96.7,2.1705,24.0,666.0,20.2,379.7,18.03

永続化したモデルを使って推論させてみる。

$ mlflow models predict \
    --no-conda \
    --model-uri mlflow-lgbm-with-sig/ \
    --content-type csv \
    --input-path data.csv
2020/09/30 19:00:54 INFO mlflow.models.cli: Selected backend for flavor 'python_function'
[15.444706764627714]

今度は、カラムを入れかえているにも関わらず、結果が変わらなくなった。 上手くいっているようだ。

ちなみに手動で Signature の情報を指定するときは次のようにすれば良い。

    # 手動で Signature を構築する場合
    from mlflow.models.signature import ModelSignature
    from mlflow.types.schema import Schema
    from mlflow.types.schema import ColSpec
    input_schema = Schema([
        ColSpec('double', 'CRIM'),
        ColSpec('double', 'ZN'),
        ColSpec('double', 'INDUS'),
        ColSpec('double', 'CHAS'),
        ColSpec('double', 'NOX'),
        ColSpec('double', 'RM'),
        ColSpec('double', 'AGE'),
        ColSpec('double', 'DIS'),
        ColSpec('double', 'RAD'),
        ColSpec('double', 'TAX'),
        ColSpec('double', 'PTRATIO'),
        ColSpec('double', 'B'),
        ColSpec('double', 'LSTAT'),
    ])
    output_schema = Schema([ColSpec('double', 'MEDV')])
    signature = ModelSignature(inputs=input_schema, outputs=output_schema)
    mlflow_lgb.save_model(booster, path='mlflow-lgbm-with-sig', signature=signature)

Input Example を追加する

また、永続化するモデルにはサンプルとなる入力データも Input Example として同梱させることができる。 次は Input Example も追加してみよう。

以下にサンプルコードを示す。 やっていることは簡単で、学習させたデータの先頭の何件かを与えているだけ。

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

import pandas as pd
import lightgbm as lgb
from mlflow import lightgbm as mlflow_lgb
from mlflow.models.signature import infer_signature
from sklearn import datasets
from sklearn.model_selection import train_test_split


def main():
    dataset = datasets.load_boston()
    train_x, train_y = dataset.data, dataset.target
    feature_names = list(dataset.feature_names)

    x_tr, x_eval, y_tr, y_eval = train_test_split(train_x, train_y,
                                                  test_size=0.33,
                                                  shuffle=True,
                                                  random_state=42,
                                                  )

    lgb_train = lgb.Dataset(x_tr, y_tr, feature_name=feature_names)
    lgb_eval = lgb.Dataset(x_eval, y_eval, reference=lgb_train)

    lgbm_params = {
        'objective': 'regression',
        'metric': 'rmse',
        'first_metric_only': True,
        'verbose': -1,
        'random_state': 42,
    }
    booster = lgb.train(params=lgbm_params,
                        train_set=lgb_train,
                        valid_sets=[lgb_train, lgb_eval],
                        num_boost_round=1_000,
                        early_stopping_rounds=100,
                        verbose_eval=50,
                        )

    y_tr_pred = booster.predict(x_tr, num_iteration=booster.best_iteration)
    x_tr_df = pd.DataFrame(x_tr, columns=dataset.feature_names)
    signature = infer_signature(x_tr_df, y_tr_pred)
    # サンプルの入力データをつける
    input_example = x_tr_df.iloc[:5]
    mlflow_lgb.save_model(booster,
                          path='mlflow-lgbm-with-sig-and-example',
                          input_example=input_example,
                          signature=signature)


if __name__ == '__main__':
    main()

上記を実行しよう。

$ python lgbmlfeg.py
Training until validation scores don't improve for 100 rounds
[50]  training's rmse: 2.23963    valid_1's rmse: 3.59161
[100] training's rmse: 1.55562    valid_1's rmse: 3.4141
[150] training's rmse: 1.22661    valid_1's rmse: 3.36079
[200] training's rmse: 1.0165 valid_1's rmse: 3.35222
[250] training's rmse: 0.860022   valid_1's rmse: 3.34358
[300] training's rmse: 0.735403   valid_1's rmse: 3.34575
[350] training's rmse: 0.630665   valid_1's rmse: 3.35982
Early stopping, best iteration is:
[254] training's rmse: 0.848352   valid_1's rmse: 3.33877
Evaluated only: rmse

見ると、今度は input_example.json というファイルがディレクトリに追加されている。

$ ls mlflow-lgbm-with-sig-and-example 
MLmodel         conda.yaml      input_example.json  model.lgb
$ cat mlflow-lgbm-with-sig-and-example/input_example.json 
{"columns": ["CRIM", "ZN", "INDUS", "CHAS", "NOX", "RM", "AGE", "DIS", "RAD", "TAX", "PTRATIO", "B", "LSTAT"], "data": [[10.233, 0.0, 18.1, 0.0, 0.614, 6.185, 96.7, 2.1705, 24.0, 666.0, 20.2, 379.7, 18.03], [0.67191, 0.0, 8.14, 0.0, 0.538, 5.813, 90.3, 4.682, 4.0, 307.0, 21.0, 376.88, 14.81], [0.14455, 12.5, 7.87, 0.0, 0.524, 6.172, 96.1, 5.9505, 5.0, 311.0, 15.2, 396.9, 19.15], [0.11132, 0.0, 27.74, 0.0, 0.609, 5.983, 83.5, 2.1099, 4.0, 711.0, 20.1, 396.9, 13.35], [0.12802, 0.0, 8.56, 0.0, 0.52, 6.474, 97.1, 2.4329, 5.0, 384.0, 20.9, 395.24, 12.27]]}

試しに、このサンプルを推論させてみよう まずは REST API を立ち上げる。

$ mlflow models serve \
    --no-conda \
    --model-uri mlflow-lgbm-with-sig-and-example
2020/09/30 19:05:39 INFO mlflow.models.cli: Selected backend for flavor 'python_function'
2020/09/30 19:05:39 INFO mlflow.pyfunc.backend: === Running command 'gunicorn --timeout=60 -b 127.0.0.1:5000 -w 1 ${GUNICORN_CMD_ARGS} -- mlflow.pyfunc.scoring_server.wsgi:app'
[2020-09-30 19:05:39 +0900] [23035] [INFO] Starting gunicorn 20.0.4
[2020-09-30 19:05:39 +0900] [23035] [INFO] Listening at: http://127.0.0.1:5000 (23035)
[2020-09-30 19:05:39 +0900] [23035] [INFO] Using worker: sync
[2020-09-30 19:05:39 +0900] [23037] [INFO] Booting worker with pid: 23037

サンプルの JSON ファイルを使って REST API を叩く。

$ curl -X POST \
    -H "Content-Type:application/json" \
    --data "$(cat mlflow-lgbm-with-sig-and-example/input_example.json)" \
    http://localhost:5000/invocations
[15.444706764627714, 16.79758862860849, 25.64257218297901, 19.626464010328057, 20.184689951658456]

ちゃんと推論できているようだ。 今のところクライアント側からサンプルの情報は得られないのかな。 とはいえ、モデルがどんな入力を受け取るかソースコードを見て調べることってよくある。 なので、管理する上で助かるといえば助かるのかな。

前処理が必要なデータセットで試す

ところで、ここまでのサンプルコードには前処理が入っていなかった。 しかし、実際には前処理が存在しない機械学習のコードなんて考えられないだろう。 続いては前処理を含んだコードを MLflow Models で扱う方法について考えている。

たとえば、以下のサンプルコードでは前処理と推論の処理を scikit-learn の Pipeline としてまとめている。 Pipeline にまとめるには、関連するオブジェクトが scikit-learn のインターフェースに準拠している必要がある。 そこで LightGBM の分類器としては LGBMClassifier を使った。 また、ラベルエンコードには category_encoders の実装を使っている。 分類するデータには Titanic データセットを使った。

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

import seaborn as sns
import category_encoders as ce
from mlflow import sklearn as mlflow_sklearn
from mlflow.models.signature import ModelSignature
from mlflow.types.schema import Schema
from mlflow.types.schema import ColSpec
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
import lightgbm as lgb


def main():
    # Titanic データを読み込む
    df = sns.load_dataset('titanic')

    # 使う特徴量
    feature_names = [
        'class',
        'sex',
        'age',
        'sibsp',
        'parch',
        'fare',
        'embark_town',
        'deck',
    ]
    train_x = df[feature_names]
    train_y = df['survived']

    # 前処理 (ラベルエンコード)
    cols = ['class', 'sex', 'embark_town', 'deck']
    encoder = ce.OrdinalEncoder(cols)
    encoded_train_x = encoder.fit_transform(train_x)

    # 学習用データと検証用データに分割する
    x_tr, x_eval, y_tr, y_eval = train_test_split(encoded_train_x,
                                                  train_y,
                                                  test_size=0.33,
                                                  shuffle=True,
                                                  random_state=42,
                                                  )
    # 分類器
    clf = lgb.LGBMClassifier(n_estimators=1_000,
                             first_metric_only=True,
                             random_state=42,
                             )
    # 学習させる
    clf.fit(x_tr, y_tr,
            early_stopping_rounds=100,
            eval_set=[(x_eval, y_eval)],
            verbose=50,
            eval_metric='binary_logloss',
            )

    # 学習させたエンコーダーとモデルをパイプラインにまとめる
    steps = [
        ('preprocessing', encoder),
        ('classification', clf)
    ]
    pipeline = Pipeline(steps)

    # パイプラインを MLflow Models で保存する
    # NOTE: Categorical な型があると MLflow がスキーマをうまく推測できない
    input_schema = Schema([
        ColSpec('string', 'class'),
        ColSpec('string', 'sex'),
        ColSpec('double', 'age'),
        ColSpec('long', 'sibsp'),
        ColSpec('long', 'parch'),
        ColSpec('double', 'fare'),
        ColSpec('string', 'embark_town'),
        ColSpec('string', 'deck'),
    ])
    output_schema = Schema([ColSpec('double', 'survived')])
    signature = ModelSignature(inputs=input_schema,
                               outputs=output_schema)
    input_example = train_x.iloc[:5]
    mlflow_sklearn.save_model(pipeline,
                              path='mlflow-sklearn-pipeline',
                              signature=signature,
                              input_example=input_example,
                              )


if __name__ == '__main__':
    main()

上記を実行しよう。

$ python skpipemlf.py
Training until validation scores don't improve for 100 rounds
[50]  valid_0's binary_logloss: 0.457466
[100]  valid_0's binary_logloss: 0.510931
Early stopping, best iteration is:
[25]  valid_0's binary_logloss: 0.427704
Evaluated only: binary_logloss

永続化したモデルで推論してみよう。 データは次のように CSV のファイルとして記録しておく。 見ると分かるとおり、扱う上ではエンコードが必要となるカラムが複数含まれている。

$ cat << 'EOF' > data.csv 
class,sex,age,sibsp,parch,fare,embark_town,deck
Third,male,22.0,1,0,7.25,Southampton,
First,female,38.0,1,0,71.2833,Cherbourg,C
Third,female,26.0,0,0,7.925,Southampton,
First,female,35.0,1,0,53.1,Southampton,C
Third,male,35.0,0,0,8.05,Southampton,
EOF

しかし、先ほどのサンプルコードでは前処理を含めたパイプラインを MLflow Models で永続化している。 そのため、前処理が必要なデータをそのまま放り込んでも推論できる。 DeprecationWarning は出ているところは愛嬌ということで。

$ mlflow models predict \
    --no-conda \
    --model-uri mlflow-sklearn-pipeline/ \
    --content-type csv \
    --input-path data.csv
2020/09/30 19:19:29 INFO mlflow.models.cli: Selected backend for flavor 'python_function'
/Users/amedama/.virtualenvs/py38/lib/python3.8/site-packages/patsy/constraint.py:13: DeprecationWarning: Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated since Python 3.3, and in 3.9 it will stop working
  from collections import Mapping
[0, 1, 0, 1, 0]

ところで上記を見てわかるとおり、結果はバイナリの整数に丸められてしまっている。 これは MLflow Models の sklearn モジュールでは、モデルの predict() メソッドを呼ぶように作られているため。 scikit-learn のインターフェースでは分類器の predict() が整数に丸めた結果を返してしまう。

モデルが確率 (predict_proba) を返すようにする

ただ、丸めた結果だけでは困るケースが多いはず。 なので、試しに predict_proba() の結果を返すようにしてみよう。 やり方は簡単で LGBMClassifier を継承して predict()predict_proba() にすりかえるクラスを用意する。

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

import seaborn as sns
import category_encoders as ce
from mlflow import sklearn as mlflow_sklearn
from mlflow.models.signature import ModelSignature
from mlflow.types.schema import Schema
from mlflow.types.schema import ColSpec
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
import lightgbm as lgb


class LGBMClassifierWrapper(lgb.LGBMClassifier):
    """predict() の処理を predict_proba() にリダイレクトするラッパー"""

    def predict(self, *args, **kwargs):
        # 処理をリダイレクトする
        proba = super().predict_proba(*args, **kwargs)
        # Positive の確率を返す
        return proba[:, 1]

def main():
    df = sns.load_dataset('titanic')

    feature_names = [
        'class',
        'sex',
        'age',
        'sibsp',
        'parch',
        'fare',
        'embark_town',
        'deck',
    ]
    train_x = df[feature_names]
    train_y = df['survived']

    cols = ['class', 'sex', 'embark_town', 'deck']
    encoder = ce.OrdinalEncoder(cols)
    encoded_train_x = encoder.fit_transform(train_x)

    x_tr, x_eval, y_tr, y_eval = train_test_split(encoded_train_x,
                                                  train_y,
                                                  test_size=0.33,
                                                  shuffle=True,
                                                  random_state=42,
                                                  )
    # 処理をラップした分類器を使う
    clf = LGBMClassifierWrapper(n_estimators=1_000,
                                first_metric_only=True,
                                random_state=42,
                                )
    clf.fit(x_tr, y_tr,
            early_stopping_rounds=100,
            eval_set=[(x_eval, y_eval)],
            verbose=50,
            eval_metric='binary_logloss',
            )

    steps = [
        ('preprocessing', encoder),
        ('classification', clf)
    ]
    pipeline = Pipeline(steps)

    input_schema = Schema([
        ColSpec('string', 'class'),
        ColSpec('string', 'sex'),
        ColSpec('double', 'age'),
        ColSpec('long', 'sibsp'),
        ColSpec('long', 'parch'),
        ColSpec('double', 'fare'),
        ColSpec('string', 'embark_town'),
        ColSpec('string', 'deck'),
    ])
    output_schema = Schema([ColSpec('double', 'survived')])
    signature = ModelSignature(inputs=input_schema,
                               outputs=output_schema)
    input_example = train_x.iloc[:5]
    mlflow_sklearn.save_model(pipeline,
                              path='mlflow-sklearn-pipeline-with-proba',
                              signature=signature,
                              input_example=input_example,
                              )


if __name__ == '__main__':
    main()

上記を実行しよう。

$ python skpipemlfp.py
Training until validation scores don't improve for 100 rounds
[50]  valid_0's binary_logloss: 0.457466
[100]  valid_0's binary_logloss: 0.510931
Early stopping, best iteration is:
[25]  valid_0's binary_logloss: 0.427704
Evaluated only: binary_logloss

推論させてみると、今度はちゃんと浮動小数点の結果になっている。

$ mlflow models predict \                   
    --no-conda \
    --model-uri mlflow-sklearn-pipeline-with-proba/ \
    --content-type csv \
    --input-path data.csv
2020/09/30 19:24:09 INFO mlflow.models.cli: Selected backend for flavor 'python_function'
/Users/amedama/.virtualenvs/py38/lib/python3.8/site-packages/patsy/constraint.py:13: DeprecationWarning: Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated since Python 3.3, and in 3.9 it will stop working
  from collections import Mapping
[0.20596751692506657, 0.9204031036579741, 0.4137860515635112, 0.9286529085847584, 0.0976711088927368]

永続化した内容を Python から読み込んで使う

ここまでの内容は、永続化した内容を常に mlflow コマンドから読み込んで使ってきた。 しかし、Python のコードから読み込んで使いたいケースも当然あるはず。

以下のサンプルコードでは先ほど永続化したモデルを読み込んで使っている。 具体的には mlflow.pyfunc.load_model() を使えばモデルが読み込める。

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

import seaborn as sns
from mlflow import pyfunc


def main():
    # Titanic データを読み込む
    df = sns.load_dataset('titanic')

    # Categorical 型は文字列に直す
    df = df.astype({
        'class': str,
        'deck': str,
    })

    # 使う特徴量の名前
    feature_names = [
        'class',
        'sex',
        'age',
        'sibsp',
        'parch',
        'fare',
        'embark_town',
        'deck',
    ]
    train_x = df[feature_names]
    train_y = df['survived']

    # 保存したモデルを読み込む
    model_path = 'mlflow-sklearn-pipeline-with-proba'
    # 汎用な pyfunc モジュールから読み出せる
    loaded_model = pyfunc.load_model(model_path)
    """
    # あるいは sklearn 向けモジュールから読んでも良い
    from mlflow import sklearn as mlflow_sklearn
    loaded_model = mlflow_sklearn.load_model(model_path)
    """

    # 保存したモデルで予測する
    # NOTE: ここで予測しているのはモデルが見たことのあるデータなので、あくまでデモとして
    train_y_pred = loaded_model.predict(train_x)

    # 先頭を表示してみる
    print(f'Inference: {train_y_pred[:5]}')
    # 正解
    print(f'GroundTruth: {train_y.values[:5]}')


if __name__ == '__main__':
    main()

ポイントとしては、永続化に使ったモジュールが何であれ、この統一されたインターフェースから読み出せるということ。 ようするに mlflow.sklearnmlflow.lightgbm などのモジュールを使って永続化したモデルであっても、ひとつの API で読める。 MLmodel ファイルには loader_module という項目に、モデルの復元に使うモジュールが指定されているため、このようなことが実現できる。 復元したモデルには predict() メソッドがあるので、あとはこれを使って推論すれば良い。

上記を実行してみよう。

$ python load.py                      
/Users/amedama/.virtualenvs/py38/lib/python3.8/site-packages/patsy/constraint.py:13: DeprecationWarning: Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated since Python 3.3, and in 3.9 it will stop working
  from collections import Mapping
Inference: [0.20596752 0.9204031  0.41378605 0.92865291 0.09767111]
GroundTruth: [0 1 1 1 0]

ちゃんと推論できている。

Custom Python Models を作る

先ほど扱ったサンプルコードでは、前処理とモデルが scikit-learn のインターフェースを備えていることを前提としていた。 しかし、扱うコードによっては scikit-learn のインターフェースがない場合もあるはず。 続いては、そんな場合にどうすれば良いかを扱う。

最も簡単なやり方は mlflow.pyfunc.PythonModel を継承したクラスを作るというもの。 継承したクラスの predict() メソッドに、生データから推論するまでに必要な処理のパイプラインを詰め込む。 そして、このクラスのインスタンスを mlflow.pyfunc.save_model() で永続化してやれば良い。

以下にサンプルコードを示す。 今度は LightGBM の標準 API を使っているため scikit-learn のインターフェースに準拠していない。 つまり、scikit-learn の Pipeline にまとめる作戦が使えない状況を意図的に作り出している。

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

import seaborn as sns
import category_encoders as ce
from mlflow import pyfunc
from mlflow.models.signature import ModelSignature
from mlflow.pyfunc.model import get_default_conda_env
from mlflow.types.schema import Schema
from mlflow.types.schema import ColSpec
from sklearn.model_selection import train_test_split
import lightgbm as lgb


class InferencePipeline(pyfunc.PythonModel):
    """推論に使うパイプライン"""

    def __init__(self, preprocessor, estimator):
        self.preprocessor = preprocessor
        self.estimator = estimator

    def predict(self, context, model_input):
        """入力を推論結果に変換する過程"""
        transformed_input = self.preprocessor.transform(model_input)
        prediction = self.estimator.predict(transformed_input)
        return prediction


def main():
    df = sns.load_dataset('titanic')

    feature_names = [
        'class',
        'sex',
        'age',
        'sibsp',
        'parch',
        'fare',
        'embark_town',
        'deck',
    ]
    train_x = df[feature_names]
    train_y = df['survived']

    cols = ['class', 'sex', 'embark_town', 'deck']
    encoder = ce.OrdinalEncoder(cols)
    encoded_train_x = encoder.fit_transform(train_x)

    x_tr, x_eval, y_tr, y_eval = train_test_split(encoded_train_x,
                                                  train_y,
                                                  test_size=0.33,
                                                  shuffle=True,
                                                  random_state=42,
                                                  )

    # lightgbm.train() を使う
    # 返ってくる Booster オブジェクトには scikit-learn インターフェースがない
    lgb_train = lgb.Dataset(x_tr, y_tr)
    lgb_eval = lgb.Dataset(x_eval, y_eval, reference=lgb_train)

    lgb_params = {
        'objective': 'binary',
        'metrics': 'binary_logloss',
        'first_metric_only': True,
        'random_state': 42,
        'verbose': -1,
    }
    booster = lgb.train(lgb_params, lgb_train,
                        valid_sets=lgb_eval,
                        num_boost_round=1_000,
                        early_stopping_rounds=100,
                        verbose_eval=50,
                        )

    # 前処理と推論の処理を mlflow.pyfunc.PythonModel を継承したクラスのインスタンスにまとめる
    pipeline = InferencePipeline(encoder, booster)

    input_schema = Schema([
        ColSpec('string', 'class'),
        ColSpec('string', 'sex'),
        ColSpec('double', 'age'),
        ColSpec('long', 'sibsp'),
        ColSpec('long', 'parch'),
        ColSpec('double', 'fare'),
        ColSpec('string', 'embark_town'),
        ColSpec('string', 'deck'),
    ])
    output_schema = Schema([ColSpec('double', 'survived')])
    signature = ModelSignature(inputs=input_schema,
                               outputs=output_schema)

    input_example = train_x.iloc[:5]

    # 動作に必要な依存ライブラリを追加する
    conda_env = get_default_conda_env()
    deps = conda_env['dependencies']
    other_deps = deps[-1]  # XXX: ちょっと決め打ちすぎ
    other_deps['pip'].append('category_encoders')
    other_deps['pip'].append('scikit-learn')
    other_deps['pip'].append('lightgbm')

    # 永続化する
    pyfunc.save_model(path='mlflow-custom-pyfunc-model',
                      python_model=pipeline,
                      signature=signature,
                      input_example=input_example,
                      conda_env=conda_env,
                      )


if __name__ == '__main__':
    main()

MLmodel は次のように記録されている。 モデルの本体は Pickle オブジェクトとして python_model.pkl にある

$ cat mlflow-custom-pyfunc-model/MLmodel   
flavors:
  python_function:
    cloudpickle_version: 1.6.0
    env: conda.yaml
    loader_module: mlflow.pyfunc.model
    python_model: python_model.pkl
    python_version: 3.8.5
saved_input_example_info:
  artifact_path: input_example.json
  pandas_orient: split
  type: dataframe
signature:
  inputs: '[{"name": "class", "type": "string"}, {"name": "sex", "type": "string"},
    {"name": "age", "type": "double"}, {"name": "sibsp", "type": "long"}, {"name":
    "parch", "type": "long"}, {"name": "fare", "type": "double"}, {"name": "embark_town",
    "type": "string"}, {"name": "deck", "type": "string"}]'
  outputs: '[{"name": "survived", "type": "double"}]'
utc_time_created: '2020-09-30 09:11:38.717424'

永続化した内容を使って推論させてみよう。

$ mlflow models predict \
    --no-conda \
    --model-uri mlflow-custom-pyfunc-model/ \
    --content-type csv \
    --input-path data.csv
2020/09/30 20:00:10 INFO mlflow.models.cli: Selected backend for flavor 'python_function'
/Users/amedama/.virtualenvs/py38/lib/python3.8/site-packages/patsy/constraint.py:13: DeprecationWarning: Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated since Python 3.3, and in 3.9 it will stop working
  from collections import Mapping
[0.20596751692506657, 0.9204031036579741, 0.4137860515635112, 0.9286529085847584, 0.0976711088927368]

ちゃんと動いていることがわかる。

ちなみに、今回は扱わなかったけどモデルの情報を Pickle 以外のファイルに artifacts として保存することもできるようだ。 また、さらに複雑なモデルや Python 以外の言語を使う場合には、自分で Custom Flavor を書くこともできる。

とりあえず、そんな感じで。

Python: LightGBM の cv() 関数から得られるモデルの特徴量の重要度を可視化してみる

今回は LightGBM の cv() 関数から得られる複数の学習済み Booster から特徴量の重要度を取り出して可視化してみる。 それぞれの Booster 毎のバラつきなどから各特徴量の傾向などが確認できるかもしれない。

使った環境は次のとおり。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.15.6
BuildVersion:   19G2021
$ python -V         
Python 3.8.5

下準備

あらかじめ必要なパッケージをインストールしておく。 なお、LightGBM のバージョン 3.0 以降 (2020-08-22 現在 RC1) をインストールしておくと学習済みモデルを取り出すのが楽になる。

$ pip install "lightgbm>=3.0.0rc1" scikit-learn seaborn

バージョン 3.0 以前の場合には次の記事を参照のこと。 ちなみに、以下のコールバックを使うやり方はバージョン 3.0 以降でも利用できる。

blog.amedama.jp

複数の学習済み Booster の特徴量の重要度を可視化する

早速だけど以下にサンプルコードを示す。 このサンプルコードでは、擬似的に生成した二値分類のデータセットを使っている。 状況としては、特徴量は 100 次元あるものの、その中で本当に有益なものは先頭の 5 次元しかない。 複数の学習済み Booster から得られる特徴量の重要度を可視化するには箱ひげ図を使った。 箱ひげ図の項目は、重要度の平均値を使ってソートしている。

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

import lightgbm as lgb
from matplotlib import pyplot as plt
import seaborn as sns
import pandas as pd
from sklearn.datasets import make_classification
from sklearn.model_selection import StratifiedKFold


def main():
    # 疑似的な教師データを作るためのパラメータ
    dist_args = {
        # データ点数
        'n_samples': 10_000,
        # 次元数
        'n_features': 100,
        # その中で意味のあるもの
        'n_informative': 5,
        # 重複や繰り返しはなし
        'n_redundant': 0,
        'n_repeated': 0,
        # タスクの難易度
        'class_sep': 0.65,
        # 二値分類問題
        'n_classes': 2,
        # 生成に用いる乱数
        'random_state': 42,
        # 特徴の順序をシャッフルしない (先頭の次元が informative になる)
        'shuffle': False,
    }
    # 教師データを作る
    train_x, train_y = make_classification(**dist_args)
    # 擬似的な特徴量の名前
    feature_names = [f'col{i}' for i in range(dist_args['n_features'])]

    # LightGBM が扱うデータセットの形式に直す
    lgb_train = lgb.Dataset(train_x, train_y,
                            feature_name=feature_names)
    # 学習用のパラメータ
    lgbm_params = {
        'objective': 'binary',
        'metric': 'binary_logloss',
        'first_metric_only': True,
        'verbose': -1,
    }
    # データの分割方法
    folds = StratifiedKFold(n_splits=5,
                            shuffle=True,
                            random_state=42,
                            )
    # 交差検証
    cv_result = lgb.cv(lgbm_params,
                       lgb_train,
                       folds=folds,
                       num_boost_round=1_000,
                       early_stopping_rounds=100,
                       verbose_eval=100,
                       # 学習済みモデルを取り出す (v3.0 以降)
                       return_cvbooster=True,
                       )

    # 学習済みモデルから特徴量の重要度を取り出す
    cvbooster = cv_result['cvbooster']
    raw_importances = cvbooster.feature_importance(importance_type='gain')
    feature_name = cvbooster.boosters[0].feature_name()
    importance_df = pd.DataFrame(data=raw_importances,
                                 columns=feature_name)
    # 平均値でソートする
    sorted_indices = importance_df.mean(axis=0).sort_values(ascending=False).index
    sorted_importance_df = importance_df.loc[:, sorted_indices]
    # 上位をプロットする
    PLOT_TOP_N = 20
    plot_cols = sorted_importance_df.columns[:PLOT_TOP_N]
    _, ax = plt.subplots(figsize=(8, 8))
    ax.grid()
    ax.set_xscale('log')
    ax.set_ylabel('Feature')
    ax.set_xlabel('Importance')
    sns.boxplot(data=sorted_importance_df[plot_cols],
                orient='h',
                ax=ax)
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python lgbcvimp.py
[100]  cv_agg's binary_logloss: 0.280906 + 0.0141782
[200] cv_agg's binary_logloss: 0.280829 + 0.018406

すると、次のようなグラフが得られる。 重要度の平均値を上位 20 件について表示している。

f:id:momijiame:20200822195419p:plain
LightGBM の cv() 関数から得られた複数の Booster から可視化した特徴量の重要度

やはり、先頭の 5 次元は特徴量の重要度が高い。 一方で、それ以外の特徴量も多少なりとも重要であるとモデルが考えていることが見て取れる。 ここらへんは Null Importance などを使うことで判断できるだろう。

blog.amedama.jp

補足

元々は棒グラフとエラーバーを使って可視化することを考えていた。 しかし、エラーバーに表示する内容に何を使うのが適切か悩んで Twitter につぶやいたところ、以下のような助言をいただくことができた。

ありがとうございます。 全人類は Kaggle 本を買おう。

Python: CatBoost を GPU で学習させる

勾配ブースティング決定木を扱うフレームワークの CatBoost は、GPU を使った学習ができる。 GPU を使うと、CatBoost の特徴的な決定木の作り方 (Symmetric Tree) も相まって、学習速度の向上が見込める場合があるようだ。 今回は、それを試してみる。

使った環境は次のとおり。

$ cat /etc/*-release
DISTRIB_ID=Ubuntu
DISTRIB_RELEASE=20.04
DISTRIB_CODENAME=focal
DISTRIB_DESCRIPTION="Ubuntu 20.04.1 LTS"
NAME="Ubuntu"
VERSION="20.04.1 LTS (Focal Fossa)"
ID=ubuntu
ID_LIKE=debian
PRETTY_NAME="Ubuntu 20.04.1 LTS"
VERSION_ID="20.04"
HOME_URL="https://www.ubuntu.com/"
SUPPORT_URL="https://help.ubuntu.com/"
BUG_REPORT_URL="https://bugs.launchpad.net/ubuntu/"
PRIVACY_POLICY_URL="https://www.ubuntu.com/legal/terms-and-policies/privacy-policy"
VERSION_CODENAME=focal
UBUNTU_CODENAME=focal
$ python3 -V
Python 3.8.2
$ python3 -m pip list | grep -i catboost
catboost               0.24

使用するハードウェアは Google Compute Engine の N1 Standard 2 インスタンスに NVIDIA Tesla T4 を 1 台アタッチしている。

$ grep "model name" /proc/cpuinfo | head -n 1
model name   : Intel(R) Xeon(R) CPU @ 2.30GHz
$ grep processor /proc/cpuinfo | wc -l
2
$ nvidia-smi
Sat Aug 22 07:05:51 2020       
+-----------------------------------------------------------------------------+
| NVIDIA-SMI 450.51.06    Driver Version: 450.51.06    CUDA Version: 11.0     |
|-------------------------------+----------------------+----------------------+
| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |
| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
|                               |                      |               MIG M. |
|===============================+======================+======================|
|   0  Tesla T4            On   | 00000000:00:04.0 Off |                    0 |
| N/A   47C    P8    10W /  70W |     70MiB / 15109MiB |      0%      Default |
|                               |                      |                  N/A |
+-------------------------------+----------------------+----------------------+
                                                                               
+-----------------------------------------------------------------------------+
| Processes:                                                                  |
|  GPU   GI   CI        PID   Type   Process name                  GPU Memory |
|        ID   ID                                                   Usage      |
|=============================================================================|
|    0   N/A  N/A       889      G   /usr/lib/xorg/Xorg                 59MiB |
|    0   N/A  N/A       983      G   /usr/bin/gnome-shell               10MiB |
+-----------------------------------------------------------------------------+

下準備

はじめに、CatBoost とその他使うパッケージをインストールしておく。

$ python3 -m pip install catboost scikit-learn

なお、CatBoost では GPU のリソースを CUDA で扱うので、使用するマシンにはあらかじめ CUDA をインストールしておく。 今回は CUDA をインストールする部分については手順から省略する。 次のスニペットを実行して、結果が 0 でなければ GPU のリソースが CatBoost から見えていることがわかる。

$ python3 -c "from catboost.utils import get_gpu_device_count; print(get_gpu_device_count())"
1

CatBoost を GPU を使って学習する

以下のサンプルコードでは、擬似的に作った二値分類のデータセットを CatBoost で学習させている。 ポイントは、学習するときに渡す辞書のパラメータに task_type というキーで GPU を指定するところ。 CatBoost から GPU のリソースが認識できていれば、これだけで GPU を使った学習ができる。

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

import sys
import time
import logging
from contextlib import contextmanager

from catboost import CatBoost
from catboost import Pool
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.metrics import log_loss


LOGGER = logging.getLogger(__name__)


@contextmanager
def timeit():
    """処理にかかった時間を計測してログに出力するコンテキストマネージャ"""
    start = time.time()
    yield
    end = time.time()
    elapsed = end - start
    LOGGER.info(f'Elapsed Time: {elapsed:.2f} sec')


def main():
    logging.basicConfig(level=logging.INFO,
                        stream=sys.stderr,
                        )

    # 疑似的な教師信号を作るためのパラメータ
    dist_args = {
        # データ点数
        'n_samples': 100_000,
        # 次元数
        'n_features': 1_000,
        # その中で意味のあるもの
        'n_informative': 100,
        # 重複や繰り返しはなし
        'n_redundant': 0,
        'n_repeated': 0,
        # タスクの難易度
        'class_sep': 0.65,
        # 二値分類問題
        'n_classes': 2,
        # 生成に用いる乱数
        'random_state': 42,
        # 特徴の順序をシャッフルしない (先頭の次元が informative になる)
        'shuffle': False,
    }
    # 教師データを作る
    train_x, train_y = make_classification(**dist_args)
    # データセットを学習用と検証用に分割する
    x_tr, x_val, y_tr, y_val = train_test_split(train_x, train_y,
                                                test_size=0.3,
                                                shuffle=True,
                                                random_state=42,
                                                stratify=train_y)
    # CatBoost が扱うデータセットの形式に直す
    train_pool = Pool(x_tr, label=y_tr)
    valid_pool = Pool(x_val, label=y_val)
    # 学習用のパラメータ
    params = {
        # タスク設定と損失関数
        'loss_function': 'Logloss',
        # 学習率
        'learning_rate': 0.02,
        # 学習ラウンド数
        'num_boost_round': 5_000,
        # 検証用データの損失が既定ラウンド数減らなかったら学習を打ち切る
        # NOTE: ラウンド数を揃えたいので今回は使わない
        # 'early_stopping_rounds': 100,
        # 乱数シード
        'random_state': 42,
        # 学習に GPU を使う場合
        'task_type': 'GPU',
    }
    # モデルを学習する
    model = CatBoost(params)
    with timeit():
        model.fit(train_pool,
                  eval_set=valid_pool,
                  verbose_eval=100,
                  use_best_model=True,
                  )
    # 検証用データを分類する
    y_pred = model.predict(valid_pool,
                           prediction_type='Probability')
    # ロジスティック損失を確認する
    metric = log_loss(y_val, y_pred)
    LOGGER.info(f'Validation Metric: {metric}')


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python3 catgpubench.py
0: learn: 0.6917674   test: 0.6918382  best: 0.6918382 (0)   total: 41.1ms  remaining: 3m 25s
100:   learn: 0.5880012   test: 0.5928815  best: 0.5928815 (100) total: 2.34s   remaining: 1m 53s
200:   learn: 0.5159286   test: 0.5238258  best: 0.5238258 (200) total: 4.52s   remaining: 1m 47s

...

4800:  learn: 0.0799767   test: 0.1156692  best: 0.1156685 (4799)    total: 1m 36s   remaining: 4.01s
4900:  learn: 0.0790359   test: 0.1150905  best: 0.1150903 (4899)    total: 1m 38s   remaining: 1.99s
4999:  learn: 0.0780408   test: 0.1144641  best: 0.1144641 (4999)    total: 1m 40s   remaining: 0us
bestTest = 0.1144640706
bestIteration = 4999
INFO:__main__:Elapsed Time: 109.65 sec
INFO:__main__:Validation Metric: 0.11446408301027777

ちゃんと GPU を使って学習できた。 GPU が使われているかは nvidia-smi などでリソースの状態を確認すると良いと思う。

CPU を使った学習と比べてみる

一応、CPU を使って学習したときと比べてみよう。 ただ、先ほど使ったインスタンスは CPU のコア数があまりにも少ない。 そのため、似たような値段で借りられる N1 Standard 16 インスタンスを使って比較する。

ハードウェアの環境的には次のとおり。

$ grep "model name" /proc/cpuinfo | head -n 1
model name   : Intel(R) Xeon(R) CPU @ 2.30GHz
$ grep processor /proc/cpuinfo | wc -l
16

あらかじめ前置きしておくと、GPU を使って高速化が見込めるかは、データセットの特性や学習時のオプション、ハードウェアなど様々なパラメータに依存する。 なので、今回の内容はあくまで「特定の環境で試したときにこうなった」という結果に過ぎない。

ソースコードを編集して、学習時のパラメータで GPU を使わないようにする。

    # 学習用のパラメータ
    params = {
        # タスク設定と損失関数
        'loss_function': 'Logloss',
        # 学習率
        'learning_rate': 0.02,
        # 学習ラウンド数
        'num_boost_round': 5_000,
        # 検証用データの損失が既定ラウンド数減らなかったら学習を打ち切る
        # NOTE: ラウンド数を揃えたいので今回は使わない
        # 'early_stopping_rounds': 100,
        # 乱数シード
        'random_state': 42,
        # 学習に GPU を使う場合
        # 'task_type': 'GPU',
    }

そして、実行しよう。

$ python3 catcpubench.py
0: learn: 0.6916098   test: 0.6916659  best: 0.6916659 (0)   total: 239ms    remaining: 19m 56s
100:   learn: 0.5917145   test: 0.5961182  best: 0.5961182 (100) total: 7.93s   remaining: 6m 24s
200:   learn: 0.5218843   test: 0.5286355  best: 0.5286355 (200) total: 15.5s   remaining: 6m 10s

...

4800:  learn: 0.0643858   test: 0.1035075  best: 0.1035075 (4800)    total: 5m 40s   remaining: 14.1s
4900:  learn: 0.0629871   test: 0.1023799  best: 0.1023799 (4900)    total: 5m 47s   remaining: 7.01s
4999:  learn: 0.0618029   test: 0.1015037  best: 0.1015037 (4999)    total: 5m 53s   remaining: 0us

bestTest = 0.1015037231
bestIteration = 4999

INFO:__main__:Elapsed Time: 356.12 sec
INFO:__main__:Validation Metric: 0.10150372305811575

すると、今回の環境では GPU を学習に使った場合と比較して約 3 倍の時間がかかった。

補足

Google Compute Engine のインスタンスタイプごとの料金設定は次のとおり。

cloud.google.com

cloud.google.com

今回は us-central1-a ゾーンのインスタンスを使用した。 利用したインスタンスの料金は、現時点 (2020-08-22) で次のとおり。

  • GPU

    • n1-standard-2 + NVIDIA Tesla T4
      • $0.445/h ($0.0950 + $0.35)
  • CPU

    • n1-standard-16
      • $0.7600/h

Python: SHAP (SHapley Additive exPlanations) を LightGBM と使ってみる

SHAP は協力ゲーム理論にもとづいて機械学習モデルを解釈する手法と、その実装を指している。 今回は、あまり理論の部分には踏み込むことなく、使い方を中心として書いていく。

使った環境は次のとおり。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.15.6
BuildVersion:   19G73
$ python -V            
Python 3.8.5

下準備

はじめに、利用するパッケージをインストールしておく。

$ pip install shap lightgbm scikit-learn matplotlib jupyterlab

また、SHAP は Jupyter 経由で使った方がインタラクティブな表示ができる。 そのため、今回は Jupyter Lab 上で操作することを想定する。

$ jupyter-lab

Jupyter Lab を立ち上げたら実験に使うノートブックを作成しておこう。

Boston データセットと LightGBM で SHAP を使う

はじめに、SHAP のパッケージをインポートする。 そして、Jupyter Lab 上でグラフを表示できるように %matplotlib inlineshap.initjs() を実行しておく。

%matplotlib inline
import shap
shap.initjs()

今回は題材として Boston データセットを使う。

train_x, train_y = shap.datasets.boston()

今回はモデルとして LightGBM を用いる。 そこで、学習用のデータと検証用のデータに分割しておこう。

from sklearn.model_selection import train_test_split
tr_x, val_x, tr_y, val_y = train_test_split(train_x, train_y, shuffle=True, random_state=42)

分割したデータを LightGBM のデータセット表現にする。

import lightgbm as lgb
lgb_train = lgb.Dataset(tr_x, tr_y)
lgb_val = lgb.Dataset(val_x, val_y, reference=lgb_train)

あとは回帰タスクとして学習しよう。

lgbm_params = {
    'objective': 'regression',
    'metric': 'rmse',
    'verbose': -1,
}
booster = lgb.train(lgbm_params,
                    lgb_train,
                    valid_sets=lgb_val,
                    num_boost_round=1000,
                    early_stopping_rounds=100,
                    verbose_eval=50,
                    )

ここからは学習したモデルが、どのように推論をするのか SHAP を使って解釈していく。 決定木系のモデルを解釈するためには shap.TreeExplainer というクラスを使おう。 モデルと学習に使ったデータを渡してオブジェクトを作る。

explainer = shap.TreeExplainer(booster, data=tr_x)

そして、TreeExplainer を使って、モデルがどのように推論するか解釈したいデータについて SHAP Value を計算しよう。 この SHAP Value は、入力したのと同じ次元と要素数で得られる。 そして、値が大きいほど推論において影響が大きいと見なすことができる。

tr_x_shap_values = explainer.shap_values(tr_x)

つまり、行方向に見れば「特定の予測に、それぞれの特徴量がどれくらい寄与したか」と解釈できる。 同様に、列方向に見れば「予測全体で、その特徴量がどれくらい寄与したか」と解釈できる。

SHAP Value は自分で可視化しても良いけど、組み込みでいくつかグラフを描画する仕組みが用意されている。 ここからは、それらを使い分けなどと共に見ていこう。

Summary Plot

はじめに、Summary Plot から。 このグラフは、デフォルトでは特徴量ごとに SHAP Value を一軸の散布図として描画する。

shap.summary_plot(shap_values=tr_x_shap_values,
                  features=tr_x,
                  feature_names=tr_x.columns)

得られるグラフは次のとおり。

f:id:momijiame:20200813184109p:plain
Summary Plot

横軸が SHAP Value で、0 から離れているほど推論において影響を与えていることになる。 縦軸の特徴量は SHAP Value の絶対値の平均値でソートされている。 それぞれの要素の色は、その特徴量の値の大小を表している。

ちなみに、Summary Plot は SHAP Value の絶対値の平均値を使った棒グラフにもできる。

shap.summary_plot(shap_values=tr_x_shap_values,
                  features=tr_x,
                  feature_names=tr_x.columns,
                  plot_type='bar')

そして、このグラフは、特徴量の重要度と解釈することもできる。

f:id:momijiame:20200813184241p:plain
Summary (Bar) Plot

試しに LightGBM に組み込まれている特徴量の重要度と比較してみよう。

lgb.plot_importance(booster, importance_type='gain')

上記で得られるグラフは次のとおり。

f:id:momijiame:20200813184413p:plain
LightGBM Feature Importance (Gain)

微妙に特徴量の順番や値の割合は異なるものの、概ね似たような順序になっていることがわかる。

Dependence Plot

続いては Dependence Plot に移る。 これは、特定の特徴量について、取りうる値と SHAP Value の関係を散布図にしたもの。 Summary Plot では色で表現されていた特徴量の値の大小に軸を割り当てている。 たとえば LSTAT (給与の低い職業に従事する人口の割合) についてプロットしてみよう。

shap.dependence_plot(ind='LSTAT',
                     interaction_index=None,
                     shap_values=tr_x_shap_values,
                     features=tr_x,
                     feature_names=tr_x.columns)

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

f:id:momijiame:20200813184813p:plain
Dependence Plot

このグラフからは LSTAT が低いほど SHAP Value も高く、推論の結果に大きく影響を与えることが見て取れる。

また、Dependence Plot では interaction_index に別の特徴量を指定できる。 たとえば RM (その地域における住居の平均的な部屋数) を使ってみよう。

shap.dependence_plot(ind='LSTAT',
                     interaction_index='RM',
                     shap_values=tr_x_shap_values,
                     features=tr_x,
                     feature_names=tr_x.columns)

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

f:id:momijiame:20200813185126p:plain
Dependence (Interaction) Plot

このグラフからは、LSTAT が高くなると住居の平均的な部屋数も減る傾向が見て取れる。

Waterfall Plot

続いては Waterfall Plot について。 Waterfall Plot は、これまでと違って特定の予測に着目して可視化する。 たとえば教師データの先頭行を指定してみよう。

shap.waterfall_plot(expected_value=explainer.expected_value,
                    shap_values=tr_x_shap_values[0],
                    features=tr_x.iloc[0],
                    feature_names=tr_x.columns)

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

f:id:momijiame:20200813185550p:plain
Waterfall Plot

このグラフでは、目的変数の平均値からスタートして最終的に推論した値に各特徴量がどのように寄与しているかを可視化している。 (LightGBM の学習は平均値を基準として各条件にもとづいて分岐していく) たとえば LSTAT は最終的に推論した 37.9 という値の内訳において +5.91 する効果があった、ということになる。

Force Plot

Waterfall Plot は各推論を個別に見ていくものなので、全体を見たいときは次の Force Plot を使うことになる。

shap.force_plot(base_value=explainer.expected_value,
                shap_values=tr_x_shap_values,
                features=tr_x,
                feature_names=tr_x.columns)

このグラフでは、複数の推論について内訳を一度に確認できる。

f:id:momijiame:20200813190620p:plain
Force Plot

ただし、特定の予測だけに絞ってデータを与えれば、個別に見ることもできる。

shap.force_plot(base_value=explainer.expected_value,
                shap_values=tr_x_shap_values[0],
                features=tr_x.iloc[0],
                feature_names=tr_x.columns)

その場合、次のようなグラフになる。

f:id:momijiame:20200813190832p:plain
Force Plot (個別)

Decision Plot

Force Plot と同じように、複数の予測について内訳を見たいときは Devision Plot も使える。

shap.decision_plot(base_value=explainer.expected_value,
                shap_values=tr_x_shap_values,
                features=tr_x,
                feature_names=list(tr_x.columns))

Decision Plot は Waterfall Plot を折れ線グラフとして同じ描画領域に重ね合わせたような感じ。 あんまりデータ数が多いとわけがわからない。

f:id:momijiame:20200813191139p:plain
Devision Plot

試しに先頭 5 件だけに絞って可視化してみよう。

shap.decision_plot(base_value=explainer.expected_value,
                shap_values=tr_x_shap_values[:5],
                features=tr_x.iloc[:5],
                feature_names=list(tr_x.columns))

特定の特徴量によって最終的な結果は大きく異なる様子が見て取れる。

f:id:momijiame:20200813191303p:plain
Devision Plot (2)

グループ単位などで可視化するときに、傾向を確認しやすいかもしれない。

そんなかんじで。

参考

shap.readthedocs.io

dropout009.hatenablog.com

Python: LightGBM の cv() 関数の実装について

今回は LightGBM の cv() 関数について書いてみる。 LightGBM の cv() 関数は、一般的にはモデルの性能を評価する交差検証に使われる。 一方で、この関数から取り出した学習済みモデルを推論にまで使うユーザもいる。 今回は、その理由やメリットとデメリットについて書いてみる。

cv() 関数から取り出した学習済みモデルを使う理由とメリット・デメリットについて

一部のユーザの間では有名だけど、LightGBM の cv() 関数は各 Fold の決定木の増やし方に特色がある。 まず、LightGBM では決定木の集まりを Booster というオブジェクトで管理している。 Booster が内包する決定木の本数は、ラウンド (イテレーション) 数として認識できる。

https://github.com/microsoft/LightGBM/blob/v3.0.0rc1/python-package/lightgbm/basic.py#L1930

ちなみに、train() 関数を使って得られるのは、この Booster というオブジェクト。 一般的に、train() 関数を使って自前で交差検証をするときは、この Booster を Fold 毎にひとつずつ学習させることになる。

一方で、cv() 関数では全 Fold を並列で、複数の Booster を一度に学習させる。 具体的には、すべての Fold で歩調を合わせながら、それぞれの Booster のラウンド数をひとつずつ増やしている。 このとき、検証用のデータに対するメトリックも、ラウンド (イテレーション) 毎に「全 Fold の平均」で計算される。 つまり、全 Fold の平均的なメトリックが悪化するタイミングで Early Stopping がかかる。

言いかえると、cv() 関数から得られる学習済みモデルは Booster が内包する決定木の本数がすべて同じに揃う。 それに対して、train() 関数を使ってひとつずつ Booster を学習する方法では、ラウンド数が Fold によってバラつくことになる。 バラつきが小さいときは良いけど、ときには大きく偏ることもあって、その際は性能の見積もりや推論に悪影響があると考えられる。 この点から、cv() 関数では Fold 毎の偏りを考慮した、ようするに無難なモデルを得られることが期待できる。 なお、各 Fold から複数の Booster が得られるので、推論するときは Averaging などで対応する。

また、ターゲットの情報を使った特徴量抽出やスタッキングをするときも、この点は都合が良い。 これらのユースケースでは、一般的にはリークを防ぐために Out-of-Fold で処理することになる。 となると、データの全体を使って学習することが難しいので、各 Fold ごとに学習したモデルを使えると使い勝手が良い。

と、ここまでメリットばかり説明してきたけど、もちろんデメリットもある。 前述したとおり、cv() 関数では各 Fold の Booster を同時に並列で学習させていく。 そのため、学習に使うデータやモデルを一度にメモリに載せることになる。 つまり、train() 関数を使って Booster をひとつずつ学習するときよりも、相対的にメモリの制約は厳しくなると考えられる。 また、他の Fold を使って補える部分もあるとはいえ Out-of-Fold したデータは学習に使えない点もデメリットとして挙げられる。

cv() 関数の実装について

ここからは LightGBM のコードを軽く追いかけてみよう。

はじめに、LightGBM のコアといえる部分は C++ で書かれている。 Python では、それを ctypes モジュールを使った Binding として呼び出している。

自身の Python 実行環境で LightGBM のインストール先パスがわかっているときは LightGBM の共有ライブラリを探してみると良い。 上記でいうコアは Python 実行環境の中に「lib_lightgbm.so」として存在している。

$ python -c "import site; print (site.getsitepackages())"
['/Users/amedama/.virtualenvs/py38/lib/python3.8/site-packages']
$ file  ~/.virtualenvs/py38/lib/python3.8/site-packages/lightgbm/lib_lightgbm.so
/Users/amedama/.virtualenvs/py38/lib/python3.8/site-packages/lightgbm/lib_lightgbm.so: Mach-O 64-bit dynamically linked shared library x86_64

上記の共有ライブラリは Python Binding の Booster クラスから呼ばれている。

https://github.com/microsoft/LightGBM/blob/v3.0.0rc1/python-package/lightgbm/basic.py#L1930

ctypes モジュールで読み込んだライブラリを _LIB として呼び出している部分がそれ。

https://github.com/microsoft/LightGBM/blob/v3.0.0rc1/python-package/lightgbm/basic.py#L1988,L1991

そして、train() 関数や cv() 関数は、上記の Booster を学習させるためのラッパーになっている。

https://github.com/microsoft/LightGBM/blob/v3.0.0rc1/python-package/lightgbm/engine.py#L18

https://github.com/microsoft/LightGBM/blob/v3.0.0rc1/python-package/lightgbm/engine.py#L394

cv() 関数に着目して読んでいくと、以下で全 Fold の Booster を同時に更新していることがわかる。

https://github.com/microsoft/LightGBM/blob/v3.0.0rc1/python-package/lightgbm/engine.py#L592

また、Early Stopping は全 Fold の平均的なメトリックを元に発火することが確認できる。

https://github.com/microsoft/LightGBM/blob/v3.0.0rc1/python-package/lightgbm/engine.py#L593,L609

いじょう。

Python: Null Importance を使った特徴量選択について

今回は特徴量選択 (Feature Selection) の手法のひとつとして使われることのある Null Importance を試してみる。 Null Importance というのは、目的変数をシャッフルして意味がなくなった状態で学習させたモデルから得られる特徴量の重要度を指す。 では、それを使ってどのように特徴量選択をするかというと、シャッフルしなかったときの重要度との比率をスコアとして計算する。 もし、シャッフルしたときの重要度が元となった重要度よりも小さくなっていれば、スコアは大きくなって特徴量に意味があるとみなせる。 一方で、シャッフルしたときの重要度が元とさほど変わらなければ、スコアは小さくなってその特徴量は単なるノイズに近い存在と判断できる。 あとはスコアに一定の閾値を設けたり、上位 N 件を取り出すことで特徴量選択ができるようだ。

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

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.15.6
BuildVersion:   19G73
$ python -V                             
Python 3.8.3

下準備

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

$ pip install scikit-learn lightgbm matplotlib pandas tqdm

題材とするデータについて

題材としては scikit-learn が生成するダミーデータを用いる。 sklearn.datasets.make_classification() 関数を使うと、データの次元数や推論の難易度などを調整してダミーデータが作れる。

以下のサンプルコードでは二値分類のダミーデータを生成している。 このダミーデータは全体で 100 次元あるものの、先頭の 5 次元だけが推論する上で意味のあるデータとなっている。 試しにダミーデータを RandomForest で分類して、モデルに組み込みの特徴量の可視化してみよう。

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

import numpy as np
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate
from matplotlib import pyplot as plt


def main():
    # 疑似的な教師信号を作るためのパラメータ
    args = {
        # データ点数
        'n_samples': 1_000,
        # 次元数
        'n_features': 100,
        # その中で意味のあるもの
        'n_informative': 5,
        # 重複や繰り返しはなし
        'n_redundant': 0,
        'n_repeated': 0,
        # タスクの難易度
        'class_sep': 0.65,
        # 二値分類問題
        'n_classes': 2,
        # 生成に用いる乱数
        'random_state': 42,
        # 特徴の順序をシャッフルしない (先頭の次元が informative になる)
        'shuffle': False,
    }
    # ノイズを含んだ教師データを作る
    X, y = make_classification(**args)

    # 分類器にランダムフォレストを使う
    clf = RandomForestClassifier(n_estimators=100,
                                 random_state=42)

    # Stratified 5-Fold CV + ROC AUC でモデルを検証する
    folds = StratifiedKFold(n_splits=5,
                            shuffle=True,
                            random_state=42,
                            )
    # すべての次元を使った場合
    cv_all_result = cross_validate(clf, X, y,
                                  cv=folds,
                                  return_estimator=True,
                                  scoring='roc_auc',
                                  )
    # 意味のある特徴量だけを使った場合
    cv_ideal_result = cross_validate(clf,
                                     # 使う次元を先頭だけに絞り込む
                                     X[:, :args['n_informative']],
                                     y,
                                     return_estimator=False,
                                     scoring='roc_auc',
                                     )

    # それぞれの状況でのメトリックのスコア
    print('All used AUC:', np.mean(cv_all_result['test_score']))
    print('Ideal AUC:', np.mean(cv_ideal_result['test_score']))

    # 学習済みモデルを取り出す
    clfs = cv_all_result['estimator']

    # 特徴量の重要度を取り出す
    importances = [clf.feature_importances_ for clf in clfs]
    mean_importances = np.mean(importances, axis=0)
    std_importances = np.std(importances, axis=0)
    sorted_indices = np.argsort(mean_importances)[::-1]

    # 重要度の高い特徴量を表示する
    MAX_TOP_N = 10
    rank_n = min(X.shape[1], MAX_TOP_N)
    print('Feature importance ranking (TOP {rank_n})'.format(rank_n=rank_n))
    for rank, idx in enumerate(sorted_indices[:rank_n], start=1):
        params = {
            'rank': rank,
            'idx': idx,
            'importance': mean_importances[idx]
        }
        print('{rank}. feature {idx:02d}: {importance}'.format(**params))

    # 特徴量の重要度を可視化する
    plt.figure(figsize=(6, 8))
    plt.barh(range(rank_n),
             mean_importances[sorted_indices][:rank_n][::-1],
             color='g',
             ecolor='r',
             yerr=std_importances[sorted_indices][:rank_n][::-1],
             align='center')
    plt.yticks(range(rank_n), sorted_indices[:rank_n][::-1])
    plt.ylabel('Features')
    plt.xlabel('Importance')
    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。 はじめに、すべての次元を使って学習した場合と、実際に意味のある次元 (先頭 5 次元) だけを使った場合の ROC AUC を出力している。 そして、それに続いて値の大きさで降順ソートした特徴量の重要度が表示される。

$ python rfimp.py
All used AUC: 0.8748692786445511
Ideal AUC: 0.947211653938503
Feature importance ranking (TOP 10)
1. feature 00: 0.0843898282531266
2. feature 03: 0.05920082172020726
3. feature 01: 0.05713024480668112
4. feature 04: 0.05027373664878948
5. feature 02: 0.017385737849956347
6. feature 55: 0.010379957110544775
7. feature 52: 0.009663802073627043
8. feature 87: 0.00907575944742235
9. feature 34: 0.008999985108968961
10. feature 82: 0.008999816679827738

先頭の 5 次元だけを使った場合の方が ROC AUC のスコアが高くなっていることがわかる。 また、RandomForest のモデルに組み込まれている特徴量の重要度を見ても先頭の 5 次元が上位にきていることが確認できる。 上記の重要度を棒グラフにプロットしたものは以下のようになる。

f:id:momijiame:20200805181014p:plain
RandomForest から得られる特徴量の重要度

この特徴量の重要度をそのまま使って特徴量選択をすることもできる。 たとえば上位 N 件を取り出した場合で交差検証のスコアを比較することが考えられる。 しかし、スコアの増加が小さな特徴量は単なるノイズかどうかの判断が難しい。 たとえば上記であれば「2」の特徴量は事前知識がなければノイズかどうか怪しいところ。

特徴量の重要度を Null Importance と比較する

そこで、特徴量の重要度をより細かく検証するために Null Importance と比較する。 前述したとおり Null Importance は目的変数をシャッフルして学習させたモデルから計算した特徴量の重要度となる。 Null Importance を基準として、元の特徴量の重要度がどれくらい大きいか比べることでノイズかどうかの判断がしやすくなる。

以下のサンプルコードでは Null Importance と元の特徴量の重要度を上位 10 件についてヒストグラムとしてプロットしている。 Null Importance は何回も計算するのが比較的容易なので、このようにヒストグラムで比較できる。

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

import logging
import sys
from itertools import chain

from tqdm import tqdm
import numpy as np
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from matplotlib import pyplot as plt
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate


LOGGER = logging.getLogger(__name__)


def feature_importance(X, y):
    """特徴量の重要度を交差検証で計算する"""
    clf = RandomForestClassifier(n_estimators=100,
                                 random_state=42)

    folds = StratifiedKFold(n_splits=5,
                            shuffle=True,
                            random_state=42,
                            )

    cv_result = cross_validate(clf, X, y,
                               cv=folds,
                               return_estimator=True,
                               scoring='roc_auc',
                               n_jobs=-1,
                               )
    clfs = cv_result['estimator']

    importances = [clf.feature_importances_ for clf in clfs]
    mean_importances = np.mean(importances, axis=0)

    return mean_importances


def main():
    # 結構時間がかかるのでログを出す
    logging.basicConfig(level=logging.INFO, stream=sys.stderr)

    args = {
        'n_samples': 1_000,
        'n_features': 100,
        'n_informative': 5,
        'n_redundant': 0,
        'n_repeated': 0,
        'class_sep': 0.65,
        'n_classes': 2,
        'random_state': 42,
        'shuffle': False,
    }
    X, y = make_classification(**args)

    # ベースとなる特徴量の重要度を計算する
    LOGGER.info('Starting base importance calculation')
    base_importance = feature_importance(X, y)

    # データのシャッフルに再現性をもたせる
    np.random.seed(42)

    # Null Importance を何度か計算する
    LOGGER.info('Starting null importance calculation')
    TRIALS_N = 20
    null_importances = []
    for _ in tqdm(range(TRIALS_N)):
        # 目的変数をシャッフルする
        y_permuted = np.random.permutation(y)
        # シャッフルした状態で特徴量の重要度を計算する
        null_importance = feature_importance(X, y_permuted)
        null_importances.append(null_importance)
    null_importances = np.array(null_importances)
    # 列と行を入れ替える
    transposed_null_importances = null_importances.transpose(1, 0)

    sorted_indices = np.argsort(base_importance)[::-1]
    sorted_base_importance = base_importance[sorted_indices]
    sorted_null_importance = transposed_null_importances[sorted_indices]

    # ベースとなる特徴量の重要度の上位と Null Importance を可視化する
    LOGGER.info('Starting visualization')
    HIST_ROWS = 3
    HIST_COLS = 3
    fig, axes = plt.subplots(HIST_ROWS, HIST_COLS,
                             figsize=(8, 18))
    for index, ax in enumerate(chain.from_iterable(axes)):
        # Null Importance をヒストグラムにする
        col_null_importance = sorted_null_importance[index]
        ax.hist(col_null_importance, label='Null Importance', color='b')
        # ベースとなる特徴量の重要度の場所に縦線を引く
        col_base_importance = sorted_base_importance[index]
        ax.axvline(col_base_importance, label='Base Importance', color='r')
        # グラフの体裁
        ax.set_xlabel('Importance')
        ax.set_ylabel('Frequency')
        ax.set_title(f'Feature: {sorted_indices[index]}')
        ax.legend()
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python nullimphist.py

すると、次のようなヒストグラムが得られる。

f:id:momijiame:20200805185435p:plain
Null Importance のヒストグラムとの比較

先頭の 5 次元については Null Importance の分布から大きく離れている。 一方で、それ以降の特徴量の重要度は Null Importance の分布の中に入ってしまっていることがわかる。 ここから、先頭の 5 次元以降はノイズに近い特徴量であると判断する材料になる。

Null Importance を使って特徴量を選択してみる

それでは、実際に特徴量の選択までやってみよう。 以下のサンプルコードでは、Null Importance と元の特徴量の重要度の比率からスコアを計算している。 そして、スコアの大きさを元に上位 N% の特徴量を取り出して交差検証のスコアを確認している。

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

import logging
import sys

from tqdm import tqdm
import numpy as np
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from matplotlib import pyplot as plt
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate


LOGGER = logging.getLogger(__name__)


def _cross_validate(X, y):
    clf = RandomForestClassifier(n_estimators=100,
                                 random_state=42)

    folds = StratifiedKFold(n_splits=5,
                            shuffle=True,
                            random_state=42,
                            )

    cv_result = cross_validate(clf, X, y,
                               cv=folds,
                               return_estimator=True,
                               scoring='roc_auc',
                               n_jobs=-1,
                               )
    return cv_result


def cv_mean_feature_importance(X, y):
    """特徴量の重要度を交差検証で計算する"""
    cv_result = _cross_validate(X, y)
    clfs = cv_result['estimator']
    importances = [clf.feature_importances_ for clf in clfs]
    mean_importances = np.mean(importances, axis=0)

    return mean_importances


def cv_mean_test_score(X, y):
    """OOF Prediction のメトリックを交差検証で計算する"""
    cv_result = _cross_validate(X, y)
    mean_test_score = np.mean(cv_result['test_score'])
    return mean_test_score


def main():
    logging.basicConfig(level=logging.INFO, stream=sys.stderr)

    n_cols = 100
    args = {
        'n_samples': 1_000,
        'n_features': n_cols,
        'n_informative': 5,
        'n_redundant': 0,
        'n_repeated': 0,
        'class_sep': 0.65,
        'n_classes': 2,
        'random_state': 42,
        'shuffle': False,
    }
    X, y = make_classification(**args)
    columns = np.arange(n_cols)

    LOGGER.info('Starting base importance calculation')
    base_importance = cv_mean_feature_importance(X, y)

    np.random.seed(42)

    LOGGER.info('Starting null importance calculation')
    TRIALS_N = 20
    null_importances = []
    for _ in tqdm(range(TRIALS_N)):
        y_permuted = np.random.permutation(y)
        null_importance = cv_mean_feature_importance(X, y_permuted)
        null_importances.append(null_importance)
    null_importances = np.array(null_importances)

    # Null Importance を用いた特徴量選択の一例
    # Base Importance と Null Importance の値の比をスコアとして計算する
    criterion_percentile = 50  # 基準となるパーセンタイル (ここでは中央値)
    percentile_null_imp = np.percentile(null_importances,
                                        criterion_percentile,
                                        axis=0)
    null_imp_score = base_importance / (percentile_null_imp + 1e-6)
    # スコアの大きさで降順ソートする
    sorted_indices = np.argsort(null_imp_score)[::-1]

    # 上位 N% の特徴量を使って性能を比較してみる
    use_feature_importance_top_percentages = [100, 75, 50, 25, 10, 5, 1]

    mean_test_scores = []
    selected_cols_len = []
    for percentage in use_feature_importance_top_percentages:
        # スコアが上位のカラムを取り出す
        sorted_columns = columns[sorted_indices]
        num_of_features = int(n_cols * percentage / 100)
        selected_cols = sorted_columns[:num_of_features]
        X_selected = X[:, selected_cols]
        LOGGER.info(f'Null Importance score TOP {percentage}%')
        LOGGER.info(f'selected features: {selected_cols}')
        LOGGER.info(f'selected feature len: {len(selected_cols)}')
        selected_cols_len.append(len(selected_cols))

        mean_test_score = cv_mean_test_score(X_selected, y)
        LOGGER.info(f'mean test score: {mean_test_score}')
        mean_test_scores.append(mean_test_score)

    # 結果を可視化する
    _, ax1 = plt.subplots(figsize=(8, 4))
    ax1.plot(mean_test_scores, color='b', label='mean test score')
    ax1.set_xlabel('percentile')
    ax1.set_ylabel('mean test score')
    ax1.legend()
    ax2 = ax1.twinx()
    ax2.plot(selected_cols_len, color='r', label='selected features len')
    ax2.set_ylabel('selected features len')
    ax2.legend()
    plt.xticks(range(len(use_feature_importance_top_percentages)),
               use_feature_importance_top_percentages)
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python nullimpselect.py
INFO:__main__:Starting base importance calculation
INFO:__main__:Starting null importance calculation
100%|███████████████████████████████████████████| 20/20 [00:56<00:00,  2.82s/it]
INFO:__main__:Null Importance score TOP 100%
INFO:__main__:selected features: [ 0  3  1  4  2 55 13 52 34 87 67 18 29 56 82 94 23 75 58 89 54 47 97 61
 85 93 64 44 43 45 63  7 37 38 27 41 26 12 72 10 51 49 96 50 22 53 81 66
 78 69 60 79 16 90 19 80 24 17 36 71 15  8 39 20 98 40 11 86 84  6 74 73
  9 65 31 35 92 91 68 48 25 99 88 76 42 28  5 70 30 95 33 83 14 59 57 32
 21 46 77 62]
INFO:__main__:selected feature len: 100
INFO:__main__:mean test score: 0.8658391949639143
INFO:__main__:Null Importance score TOP 75%
INFO:__main__:selected features: [ 0  3  1  4  2 55 13 52 34 87 67 18 29 56 82 94 23 75 58 89 54 47 97 61
 85 93 64 44 43 45 63  7 37 38 27 41 26 12 72 10 51 49 96 50 22 53 81 66
 78 69 60 79 16 90 19 80 24 17 36 71 15  8 39 20 98 40 11 86 84  6 74 73
  9 65 31]
INFO:__main__:selected feature len: 75
INFO:__main__:mean test score: 0.8774203740902301
INFO:__main__:Null Importance score TOP 50%
INFO:__main__:selected features: [ 0  3  1  4  2 55 13 52 34 87 67 18 29 56 82 94 23 75 58 89 54 47 97 61
 85 93 64 44 43 45 63  7 37 38 27 41 26 12 72 10 51 49 96 50 22 53 81 66
 78 69]
INFO:__main__:selected feature len: 50
INFO:__main__:mean test score: 0.909876116663287
INFO:__main__:Null Importance score TOP 25%
INFO:__main__:selected features: [ 0  3  1  4  2 55 13 52 34 87 67 18 29 56 82 94 23 75 58 89 54 47 97 61
 85]
INFO:__main__:selected feature len: 25
INFO:__main__:mean test score: 0.932729764573096
INFO:__main__:Null Importance score TOP 10%
INFO:__main__:selected features: [ 0  3  1  4  2 55 13 52 34 87]
INFO:__main__:selected feature len: 10
INFO:__main__:mean test score: 0.9513031915436441
INFO:__main__:Null Importance score TOP 5%
INFO:__main__:selected features: [0 3 1 4 2]
INFO:__main__:selected feature len: 5
INFO:__main__:mean test score: 0.9629145987828075
INFO:__main__:Null Importance score TOP 1%
INFO:__main__:selected features: [0]
INFO:__main__:selected feature len: 1
INFO:__main__:mean test score: 0.6537074505769904

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

f:id:momijiame:20200805185850p:plain
Null Importance を用いて特徴量選択の可視化

今回は先頭 5 次元が重要とわかっている。 上記のグラフをみても上位 5% (= 5 次元) が選択されるまではスコアが順調に伸びている様子が確認できる。 しかし上位 1% (= 1 次元) までいくと推論に使える特徴量まで削ってしまいスコアが落ちていることがわかる。 ちなみに、実際には上記のように比較するのではなくスコアに一定の閾値を設けたり上位 N 件を決め打ちで選択してしまう場合も多いようだ。

LightGBM + Pandas を使った場合のサンプル

続いては、よくある例としてモデルに LightGBM を使って、データが Pandas のデータフレームの場合も確認しておこう。 サンプルコードは次のとおり。

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

import logging
import sys

import numpy as np
import pandas as pd
import lightgbm as lgb
from sklearn.metrics import roc_auc_score
from tqdm import tqdm
from sklearn.datasets import make_classification
from sklearn.model_selection import StratifiedKFold
from matplotlib import pyplot as plt


LOGGER = logging.getLogger(__name__)


class ModelExtractionCallback(object):
    """see: https://blog.amedama.jp/entry/lightgbm-cv-model"""

    def __init__(self):
        self._model = None

    def __call__(self, env):
        self._model = env.model

    def _assert_called_cb(self):
        if self._model is None:
            raise RuntimeError('callback has not called yet')

    @property
    def cvbooster(self):
        self._assert_called_cb()
        return self._model


def _cross_validate(train_x, train_y, folds):
    """LightGBM で交差検証する関数"""
    lgb_train = lgb.Dataset(train_x, train_y)

    model_extraction_cb = ModelExtractionCallback()
    callbacks = [
        model_extraction_cb,
    ]

    lgbm_params = {
        'objective': 'binary',
        'metric': 'auc',
        'first_metric_only': True,
        'verbose': -1,
    }
    lgb.cv(lgbm_params,
           lgb_train,
           folds=folds,
           num_boost_round=1_000,
           early_stopping_rounds=100,
           callbacks=callbacks,
           verbose_eval=20,
           )
    return model_extraction_cb.cvbooster


def _predict_oof(cv_booster, train_x, train_y, folds):
    """学習済みモデルから Out-of-Fold Prediction を求める"""
    split = folds.split(train_x, train_y)
    fold_mappings = zip(split, cv_booster.boosters)
    oof_y_pred = np.zeros_like(train_y, dtype=float)
    for (_, val_index), booster in fold_mappings:
        val_train_x = train_x.iloc[val_index]
        y_pred = booster.predict(val_train_x,
                                 num_iteration=cv_booster.best_iteration)
        oof_y_pred[val_index] = y_pred
    return oof_y_pred


def cv_mean_feature_importance(train_x, train_y, folds):
    """交差検証したモデルを使って特徴量の重要度を計算する"""
    cv_booster = _cross_validate(train_x, train_y, folds)
    importances = cv_booster.feature_importance(importance_type='gain')
    mean_importance = np.mean(importances, axis=0)
    return mean_importance


def cv_mean_test_score(train_x, train_y, folds):
    """交差検証で OOF Prediction の平均スコアを求める"""
    cv_booster = _cross_validate(train_x, train_y, folds)
    # OOF Pred を取得する
    oof_y_pred = _predict_oof(cv_booster, train_x, train_y, folds)
    test_score = roc_auc_score(train_y, oof_y_pred)
    return test_score


def main():
    logging.basicConfig(level=logging.INFO, stream=sys.stderr)

    n_cols = 100
    args = {
        'n_samples': 1_000,
        'n_features': n_cols,
        'n_informative': 5,
        'n_redundant': 0,
        'n_repeated': 0,
        'class_sep': 0.65,
        'n_classes': 2,
        'random_state': 42,
        'shuffle': False,
    }
    X, y = make_classification(**args)

    # Pandas のデータフレームにする
    col_names = [f'col{i}' for i in range(n_cols)]
    train_x = pd.DataFrame(X, columns=col_names)
    # インデックスが 0 から始まる連番とは限らないのでこういうチェックを入れた方が良い
    train_x.index = train_x.index * 10 + 10
    train_y = pd.Series(y, name='target')

    folds = StratifiedKFold(n_splits=5,
                            shuffle=True,
                            random_state=42,
                            )

    LOGGER.info('Starting base importance calculation')
    base_importance = cv_mean_feature_importance(train_x, train_y, folds)

    LOGGER.info('Starting null importance calculation')
    TRIALS_N = 20
    null_importances = []
    for _ in tqdm(range(TRIALS_N)):
        train_y_permuted = np.random.permutation(train_y)
        null_importance = cv_mean_feature_importance(train_x,
                                                     train_y_permuted,
                                                     folds)
        null_importances.append(null_importance)
    null_importances = np.array(null_importances)

    criterion_percentile = 50
    percentile_null_imp = np.percentile(null_importances,
                                        criterion_percentile,
                                        axis=0)
    null_imp_score = base_importance / (percentile_null_imp + 1e-6)
    sorted_indices = np.argsort(null_imp_score)[::-1]

    # 上位 N% の特徴量を使って性能を比較してみる
    use_feature_importance_top_percentages = [100, 75, 50, 25, 10, 5, 1]

    mean_test_scores = []
    percentile_selected_cols = []
    for percentage in use_feature_importance_top_percentages:
        sorted_columns = train_x.columns[sorted_indices]
        num_of_features = int(n_cols * percentage / 100)
        selected_cols = sorted_columns[:num_of_features]
        selected_train_x = train_x[selected_cols]
        LOGGER.info(f'Null Importance score TOP {percentage}%')
        LOGGER.info(f'selected features: {list(selected_cols)}')
        LOGGER.info(f'selected feature len: {len(selected_cols)}')
        percentile_selected_cols.append(selected_cols)

        mean_test_score = cv_mean_test_score(selected_train_x, train_y, folds)
        LOGGER.info(f'mean test_score: {mean_test_score}')
        mean_test_scores.append(mean_test_score)

    _, ax1 = plt.subplots(figsize=(8, 4))
    ax1.plot(mean_test_scores, color='b', label='mean test score')
    ax1.set_xlabel('Importance TOP n%')
    ax1.set_ylabel('mean test score')
    ax1.legend()
    ax2 = ax1.twinx()
    ax2.plot([len(cols) for cols in percentile_selected_cols],
             color='r', label='selected features len')
    ax2.set_ylabel('selected features len')
    ax2.legend()
    plt.xticks(range(len(use_feature_importance_top_percentages)),
               use_feature_importance_top_percentages)
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python lgbmnullimp.py

すると、次のようなヒストグラムが得られる。 こちらも、上位 5% を選択するまで徐々にスコアが上がって、上位 1% までいくとスコアが落ちていることがわかる。

f:id:momijiame:20200805190732p:plain
データに Pandas モデルに LightGBM を使った場合

めでたしめでたし。