Python: LightGBM の決定木を可視化して分岐を追ってみる

今回は、LightGBM が構築するブースターに含まれる決定木を可視化した上で、その分岐を追いかけてみよう。その過程を通して、LightGBM の最終的な出力がどのように得られているのかを確認してみよう。

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

$ sw_vers
ProductName:    macOS
ProductVersion: 11.2.3
BuildVersion:   20D91
$ python -V
Python 3.9.2

下準備

まずは動作に必要なパッケージをインストールする。

決定木の可視化のために graphviz を、並列計算のために OpenMP を入れておく。

$ brew install graphviz libomp

そして、Python のパッケージを入れる。

$ pip install lightgbm scikit-learn graphviz matplotlib

二値分類問題 (乳がんデータセット)

まずは乳がんデータセットを使って二値分類問題を扱ってみよう。

LightGBM には、lightgbm.Booster オブジェクトに含まれる決定木を可視化する API として lightgbm.plot_tree() という関数が用意されている。使うときは、tree_index オプションにイテレーション番号を指定することで、そのイテレーションで作成された決定木がグラフとして得られる。以下のサンプルコードでは、学習させた lightgbm.Booster の先頭にある決定木をグラフにプロットした。

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

from __future__ import annotations

from pprint import pprint

import lightgbm as lgb
from sklearn import datasets
from sklearn.model_selection import train_test_split
from matplotlib import pyplot as plt


def main():
    # 乳がんデータセットを使う
    dataset = datasets.load_breast_cancer()
    x, y = dataset.data, dataset.target
    # ホールドアウト
    train_x, eval_x, train_y, eval_y = train_test_split(x, y,
                                                        stratify=y,
                                                        shuffle=True,
                                                        random_state=42)
    # LightGBM のデータセット表現にする
    lgb_train = lgb.Dataset(train_x, train_y,
                            feature_name=list(dataset.feature_names))
    lgb_eval = lgb.Dataset(eval_x, eval_y, reference=lgb_train)
    # 学習パラメータ
    lgb_params = {
        'objective': 'binary',
        'metric': "binary_logloss",
        'verbose': -1,
        'seed': 42,
        'deterministic': True,
    }
    # 学習する
    booster = lgb.train(params=lgb_params,
                        train_set=lgb_train,
                        valid_sets=[lgb_train, lgb_eval],
                        num_boost_round=1_000,
                        early_stopping_rounds=50,
                        verbose_eval=10,
                        )

    # 検証用データの先頭の情報を出力する
    head_row = eval_x[0]
    pprint(dict(zip(dataset.feature_names, head_row)))
    # 1 本目の決定木だけを使って予測してみる
    single_tree_pred = booster.predict(data=[head_row],
                                       num_iteration=1)
    print(f'single tree pred: {single_tree_pred}')
    # 2 本目までの決定木を使って予測してみる
    double_tree_pred = booster.predict(data=[head_row],
                                       num_iteration=2)
    print(f'double tree pred: {double_tree_pred}')

    # 先頭の決定木を可視化してみる
    rows = 2
    cols = 1
    # 表示する領域を準備する
    fig = plt.figure(figsize=(12, 6))
    # 一本ずつプロットしていく
    for i in range(rows * cols):
        ax = fig.add_subplot(rows, cols, i + 1)
        ax.set_title(f'Booster index: {i}')
        lgb.plot_tree(booster=booster,
                      tree_index=i,
                      show_info='internal_value',
                      ax=ax,
                      )
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。実行すると、検証用データの先頭行の説明変数と、それを学習済みモデルで予測させたときのスコアが出力される。なお、以下は一例であり環境が変わると出力は異なる可能性がある。 (LightGBM の seed や deterministic オプションは、学習に CPU を使った上で完全に同一の環境でのみ結果が同じになることを仮定できる)

$ python binary.py 
Training until validation scores don't improve for 50 rounds
[10]  training's binary_logloss: 0.247243 valid_1's binary_logloss: 0.266197
[20]  training's binary_logloss: 0.116632 valid_1's binary_logloss: 0.158874
[30]  training's binary_logloss: 0.0581821    valid_1's binary_logloss: 0.113181
[40]  training's binary_logloss: 0.0286961    valid_1's binary_logloss: 0.0965949
[50]  training's binary_logloss: 0.0140411    valid_1's binary_logloss: 0.0985209
[60]  training's binary_logloss: 0.00667688   valid_1's binary_logloss: 0.10083
[70]  training's binary_logloss: 0.00317889   valid_1's binary_logloss: 0.104945
[80]  training's binary_logloss: 0.00160051   valid_1's binary_logloss: 0.115742
[90]  training's binary_logloss: 0.00082228   valid_1's binary_logloss: 0.129502
Early stopping, best iteration is:
[45]  training's binary_logloss: 0.0201096    valid_1's binary_logloss: 0.0943608
{'area error': 28.62,
 'compactness error': 0.01561,
 'concave points error': 0.009199,
 'concavity error': 0.01977,
 'fractal dimension error': 0.003629,
 'mean area': 493.8,
 'mean compactness': 0.1117,
 'mean concave points': 0.02995,
 'mean concavity': 0.0388,
 'mean fractal dimension': 0.06623,
 'mean perimeter': 82.51,
 'mean radius': 12.75,
 'mean smoothness': 0.1125,
 'mean symmetry': 0.212,
 'mean texture': 16.7,
 'perimeter error': 2.495,
 'radius error': 0.3834,
 'smoothness error': 0.007509,
 'symmetry error': 0.01805,
 'texture error': 1.003,
 'worst area': 624.1,
 'worst compactness': 0.1979,
 'worst concave points': 0.08045,
 'worst concavity': 0.1423,
 'worst fractal dimension': 0.08557,
 'worst perimeter': 93.63,
 'worst radius': 14.45,
 'worst smoothness': 0.1475,
 'worst symmetry': 0.3071,
 'worst texture': 21.74}
single tree pred: [0.66326872]
double tree pred: [0.69607225]

今回は、以下のようなグラフが得られた。

f:id:momijiame:20210319004235p:plain — 乳がんデータセットを学習したモデルの先頭にある決定木

検証用データの先頭行が、どのリーフに落ちるのかを決定木から確認してみよう。まずは、Booster index: 0 から。

最初の条件は worst perimeter <= 112.800 になっていて、データは 93.63 なので yes に分岐する。次は worst concave points <= 0.146 で、0.08045なので yes に分岐する。以下、同様に area error <= 34.750 で 28.62 なので yes。 worst texture <= 30.045 で 21.74 なので yes。 mean radius <= 13.875 で 12.75 なので yes。 mean radius <= 12.310 で 12.75 なので no。最終的に leaf 7 に落ちて、内部的なスコアは 0.678 になった。

さて、この 0.678 という値は、先頭の決定木だけを使って予測した 0.66326872 というスコアとは少し乖離がある。これは当然のことで、実際には内部的なスコアにシグモイド関数がかかるため。

Python のインタプリタを別で起動して確認してみよう。

$ python

次のようにシグモイド関数を定義する。

>>> import numpy as np
>>> def sigmoid(x):
...     return 1. / (1. + np.exp(-x))
...

leaf 7 の内部的なスコアをシグモイド関数にかけると、最終的な予測とほぼ同じ値が得られる。微妙にズレているのはグラフに出力するときの値がデフォルトだと小数点 3 桁で丸められているから。

>>> sigmoid(0.678)
0.6632921720482895

同じように 2 本 (イテレーション) 目の決定木も確認してみよう。 2 本目は分岐の詳細は省略するけど、最終的に leaf 0 に落ちて 0.151 というスコアになる。 2 本目までの決定木を使った予測は、各決定木から得られる内部的なスコアを足してシグモイド関数にかければ良い。

>>> sigmoid(0.678 + 0.151)
0.6961434435735563

サンプルコードから得られた 0.69607225 というスコアと、ほぼ同じ結果が得られることがわかる。以下、同様にすべてのイテレーションで作られた決定木のスコアを足し合わせていくことで最終的な結果が得られる。

回帰問題 (ボストンデータセット)

続いてはボストンデータセットを使って回帰問題を扱ってみる。以下のサンプルコードは先ほどとやっていることはほとんど同じ。問題が二値分類から回帰になっているだけ。

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

from __future__ import annotations

from pprint import pprint

import lightgbm as lgb
from sklearn import datasets
from sklearn.model_selection import train_test_split
from matplotlib import pyplot as plt


def main():
    # ボストンデータセットを使う
    dataset = datasets.load_boston()
    x, y = dataset.data, dataset.target
    # ホールドアウト
    train_x, eval_x, train_y, eval_y = train_test_split(x, y,
                                                        shuffle=True,
                                                        random_state=42)
    # LightGBM のデータセット表現にする
    lgb_train = lgb.Dataset(train_x, train_y,
                            feature_name=list(dataset.feature_names))
    lgb_eval = lgb.Dataset(eval_x, eval_y, reference=lgb_train)
    # 学習パラメータ
    lgb_params = {
        'objective': 'regression',
        'metric': "rmse",
        'verbose': -1,
        'seed': 42,
        'deterministic': True,
    }
    # 学習する
    booster = lgb.train(params=lgb_params,
                        train_set=lgb_train,
                        valid_sets=[lgb_train, lgb_eval],
                        num_boost_round=1_000,
                        early_stopping_rounds=50,
                        verbose_eval=10,
                        )

    # 検証用データの先頭の情報を出力する
    head_row = eval_x[0]
    pprint(dict(zip(dataset.feature_names, head_row)))
    # 1 本目の決定木だけを使って予測してみる
    single_tree_pred = booster.predict(data=[head_row],
                                       num_iteration=1)
    print(f'single tree pred: {single_tree_pred}')
    # 2 本目までの決定木を使って予測してみる
    double_tree_pred = booster.predict(data=[head_row],
                                       num_iteration=2)
    print(f'double tree pred: {double_tree_pred}')

    # 先頭の決定木を可視化してみる
    rows = 2
    cols = 1
    # 表示する領域を準備する
    fig = plt.figure(figsize=(12, 6))
    # 一本ずつプロットしていく
    for i in range(rows * cols):
        ax = fig.add_subplot(rows, cols, i + 1)
        ax.set_title(f'Booster index: {i}')
        lgb.plot_tree(booster=booster,
                      tree_index=i,
                      show_info='internal_value',
                      ax=ax,
                      )
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python regression.py 
Training until validation scores don't improve for 50 rounds
[10]  training's rmse: 4.71003    valid_1's rmse: 4.87178
[20]  training's rmse: 3.18599    valid_1's rmse: 3.9085
[30]  training's rmse: 2.68799    valid_1's rmse: 3.66692
[40]  training's rmse: 2.3669 valid_1's rmse: 3.55354
[50]  training's rmse: 2.13835    valid_1's rmse: 3.41701
[60]  training's rmse: 1.96456    valid_1's rmse: 3.38303
[70]  training's rmse: 1.81511    valid_1's rmse: 3.35055
[80]  training's rmse: 1.68986    valid_1's rmse: 3.34631
[90]  training's rmse: 1.59394    valid_1's rmse: 3.34022
[100] training's rmse: 1.49454    valid_1's rmse: 3.30722
[110] training's rmse: 1.41423    valid_1's rmse: 3.30035
[120] training's rmse: 1.33056    valid_1's rmse: 3.28779
[130] training's rmse: 1.25246    valid_1's rmse: 3.26555
[140] training's rmse: 1.19406    valid_1's rmse: 3.25197
[150] training's rmse: 1.13264    valid_1's rmse: 3.24115
[160] training's rmse: 1.07332    valid_1's rmse: 3.23656
[170] training's rmse: 1.02584    valid_1's rmse: 3.22715
[180] training's rmse: 0.983137   valid_1's rmse: 3.2189
[190] training's rmse: 0.940608   valid_1's rmse: 3.22034
[200] training's rmse: 0.898673   valid_1's rmse: 3.22073
[210] training's rmse: 0.862312   valid_1's rmse: 3.22325
[220] training's rmse: 0.827644   valid_1's rmse: 3.22175
[230] training's rmse: 0.795422   valid_1's rmse: 3.22017
Early stopping, best iteration is:
[184] training's rmse: 0.966589   valid_1's rmse: 3.21349
{'AGE': 84.1,
 'B': 395.5,
 'CHAS': 0.0,
 'CRIM': 0.09178,
 'DIS': 2.6463,
 'INDUS': 4.05,
 'LSTAT': 9.04,
 'NOX': 0.51,
 'PTRATIO': 16.6,
 'RAD': 5.0,
 'RM': 6.416,
 'TAX': 296.0,
 'ZN': 0.0}
single tree pred: [23.08757859]
double tree pred: [23.24927528]

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

f:id:momijiame:20210319010222p:plain — ボストンデータセットを学習したモデルの先頭にある決定木

先ほどと同じように決定木の分岐を追いかけてみる。分岐を辿ると、1 本目の決定木は leaf 10 に落ちて 23.088 というスコアが得られる。これは先頭 1 本だけを使った予測と同じ値になっており、回帰では内部的なスコアがそのまま最終的な出力となることがわかる。

同様に 2 本目の分岐を辿ると leaf 11 に落ちて 0.162 というスコアになる。 2 本目までを使った予測は、両方の決定木のスコアを足し合わせることで得られる。

>>> 23.088 + 0.162
23.25

多値分類問題 (あやめデータセット)

続いてはあやめデータセットを使って多値分類問題を扱う。基本的にはこれまでと変わらないけど、多値分類問題は内部的に作られる決定木の数が多い。

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

from __future__ import annotations

from pprint import pprint

import lightgbm as lgb
from sklearn import datasets
from sklearn.model_selection import train_test_split
from matplotlib import pyplot as plt


def main():
    # あやめデータセットを使う
    dataset = datasets.load_iris()
    x, y = dataset.data, dataset.target
    # ホールドアウト
    train_x, eval_x, train_y, eval_y = train_test_split(x, y,
                                                        stratify=y,
                                                        shuffle=True,
                                                        random_state=42)
    # LightGBM のデータセット表現にする
    lgb_train = lgb.Dataset(train_x, train_y,
                            feature_name=list(dataset.feature_names))
    lgb_eval = lgb.Dataset(eval_x, eval_y, reference=lgb_train)
    # 学習パラメータ
    lgb_params = {
        'objective': 'multiclass',
        'metric': 'softmax',
        'num_class': 3,
        'verbose': -1,
        'seed': 42,
        'deterministic': True,
    }
    # 学習する
    booster = lgb.train(params=lgb_params,
                        train_set=lgb_train,
                        valid_sets=[lgb_train, lgb_eval],
                        num_boost_round=1_000,
                        early_stopping_rounds=50,
                        verbose_eval=10,
                        )

    # 検証用データの先頭の情報を出力する
    head_row = eval_x[0]
    pprint(dict(zip(dataset.feature_names, head_row)))
    # 1 本目の決定木だけを使って予測してみる
    single_tree_pred = booster.predict(data=[head_row],
                                       num_iteration=1)
    print(f'single tree pred: {single_tree_pred}')
    # 2 本目までの決定木を使って予測してみる
    double_tree_pred = booster.predict(data=[head_row],
                                       num_iteration=2)
    print(f'double tree pred: {double_tree_pred}')

    # 先頭の決定木を可視化してみる
    rows = 2
    cols = 3
    # 表示する領域を準備する
    fig = plt.figure(figsize=(14, 6))
    # 一本ずつプロットしていく
    for i in range(rows * cols):
        ax = fig.add_subplot(rows, cols, i + 1)
        ax.set_title(f'Booster index: {i}')
        lgb.plot_tree(booster=booster,
                      tree_index=i,
                      show_info='internal_value',
                      ax=ax,
                      )
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python multiclass.py
Training until validation scores don't improve for 50 rounds
[10]  training's multi_logloss: 0.287552  valid_1's multi_logloss: 0.366785
[20]  training's multi_logloss: 0.119475  valid_1's multi_logloss: 0.232667
[30]  training's multi_logloss: 0.0678466 valid_1's multi_logloss: 0.234263
[40]  training's multi_logloss: 0.0346539 valid_1's multi_logloss: 0.270086
[50]  training's multi_logloss: 0.016588  valid_1's multi_logloss: 0.350929
[60]  training's multi_logloss: 0.00939384    valid_1's multi_logloss: 0.384578
[70]  training's multi_logloss: 0.00567975    valid_1's multi_logloss: 0.414652
Early stopping, best iteration is:
[26]  training's multi_logloss: 0.0854995 valid_1's multi_logloss: 0.216513
{'petal length (cm)': 1.3,
 'petal width (cm)': 0.2,
 'sepal length (cm)': 4.4,
 'sepal width (cm)': 3.2}
single tree pred: [[0.4084366 0.2957817 0.2957817]]
double tree pred: [[0.47189453 0.26405274 0.26405274]]

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

f:id:momijiame:20210319010334p:plain — あやめデータセットを学習したモデルの先頭にある決定木

LightGBM では、多値分類問題を扱う際に「クラス数 x イテレーション数」本の決定木が作られる。先頭にある「クラス数」本の決定木が、各クラスの出力を得るのに使われる。今回の例でいえばあやめの品種は 3 種類なので、先頭の 3 本が 1 イテレーション目のそれぞれの品種に対応することになる。同様に、4 ~ 6 本目が 2 イテレーション目のそれぞれの品種に対応する。ようするに、上記のグラフでいうと縦に並んでいる決定木がそれぞれの品種 (クラス) に対応しているということ。

今回も決定木の分岐を追いかけてみよう。 1 本目の決定木は leaf 0 に落ちて -0.884 というスコアになる。同様に、2 本目は leaf 0 に落ちて -1.207 になる。 3 本目は leaf 2 に落ちて -1.207 になる。

さて、1 イテレーション目の内部的なスコアは [-0.884, -1.207, -1.207] になった。これを 1 イテレーション目の最終的な出力である [0.4084366 0.2957817 0.2957817] にするにはソフトマックス関数にかける。

以下のように、ソフトマックス関数を定義する。

>>> def softmax(x):
...     return np.exp(x) / np.sum(np.exp(x))
...

内部的なスコアをソフトマックス関数にかけてみよう。

>>> softmax([-0.884, -1.207, -1.207])
array([0.40850546, 0.29574727, 0.29574727])

すると、最終的な出力とほぼ同じ値になった。 2 イテレーション目以降の処理は、これまでと同じなので省略する。足すだけ。

まとめ

LightGBM の決定木を可視化して分岐を追いかけることで最終的な予測がどのように得られるのかを確認できた。

Linuxで動かしながら学ぶTCP/IPネットワーク入門

作者:もみじあめ
発売日: 2020/02/29
メディア: Kindle版

CUBE SUGAR CONTAINER

技術系のこと書きます。

Python: LightGBM の決定木を可視化して分岐を追ってみる

もくじ

下準備

二値分類問題 (乳がんデータセット)

回帰問題 (ボストンデータセット)

多値分類問題 (あやめデータセット)

まとめ