本記事のサマリ

　勾配ブースティング木をPythonでスクラッチ実装しました。本記事では、スクラッチ実装したコードとそれを実行した結果について書いていこうと思います。なお、勾配ブースティング木の理論については、下記の記事に書いたので興味がある方はご覧ください。

aisinkakura-datascientist.hatenablog.com

勾配ブースティング木のアルゴリズム

　勾配ブースティング木のアルゴリズムを図１に載せます。

　今回は図１に載せたアルゴリズムを実装していきます。なお、損失関数は最小二乗誤差 $L(y, x) = \dfrac{1}{2}(y-x) ^ 2$ とします。

実装するクラスの内容

　勾配ブースティング木のモデルを作成するメソッド（.fit）と予測をするメソッド（.predict）が含まれたクラスGradientBoostingHandmadeを実装していきます。GradientBoostingHandmadeの引数を表１にまとめました。

　また、GradientBoostingHandmadeの外から使われる関数の一覧も表２に記載しておきます。

　これらを実装した結果は下記です。

class GradientBoostingHandmade:
    
    def __init__(self, m, WeakLearner, **weaklearner_params):
        self._m = m
        self._WeakLearner = WeakLearner
        self._weaklearner_params = weaklearner_params
        
    @property
    def m(self):
        return self._m
    
    @property
    def WeakLearner(self):
        return self._WeakLearner
    
    @property
    def weaklearner_params(self):
        return self._weaklearner_params
                
    def _calculate_pseudo_residuals(self, y, fmx):
        return y-fmx
    
    def _calculate_model_predict_class(self, X, pseudo_residuals):
        model = self.WeakLearner(**self.weaklearner_params)
        model.fit(X, pseudo_residuals)
        predict_class = model.apply(X_train)
        
        return (model, predict_class)
    
    def _calculate_gamma(self, y, fmx, predict_class):
        data = {
            "y":y,
            "fmx":fmx,
            "predict_class":predict_class
        }
        
        tmp_df = pd.DataFrame(data)
        
        # ここからの処理が損失関数によって変わる
        # 今回はMSEの場合
        tmp_df["for_calculate_gamma"] = tmp_df["y"]-tmp_df["fmx"]
        
        # 所属するクラスごとに残差の平均を計算する
        df_gamma = (
            tmp_df
            .groupby("predict_class", as_index=False)
            .agg({"for_calculate_gamma":"mean"})
            .rename(columns={"for_calculate_gamma":"gamma"})
        )
        
        return df_gamma
    
    def _initialize_prediction_method_on_MSE(self, X, y):
        return np.full(len(X), np.mean(y))
    
    def _minus_gradient_on_MSE(self, target, feature, prediction_method):
        return target-prediction_method(feature)
    
    
    def fit(self, X, y):
        # インデックスの振り直し(これをしないとindexがぐちゃぐちゃになって動かない)
        X = X.reset_index(drop=True)
        y = y.reset_index(drop=True)
        
        self.f0x = fmx = self._initialize_prediction_method_on_MSE(X, y)
        
        self.model_gamma_list = []
        
        for _ in range(self.m):
            pseudo_residuals = self._calculate_pseudo_residuals(y, fmx)
            model, predict_class = self._calculate_model_predict_class(X, pseudo_residuals)
            df_gamma = self._calculate_gamma(y, fmx, predict_class)
            self.model_gamma_list.append([model, df_gamma])
            fmx = self.predict(X)
    
    def predict(self, X):
        gamma_result_df_list = []
        for model, df_predict_class_gamma in self.model_gamma_list:
            result_class = model.apply(X)
            df_result = pd.DataFrame(result_class, columns=["predict_class"])
            df_merged = df_result.merge(df_predict_class_gamma, on="predict_class", how="left")
            df_merged_drop_predict_class = df_merged.drop(columns=["predict_class"])
            gamma_result_df_list.append(df_merged_drop_predict_class)
        
        gamma_result_df = pd.concat(gamma_result_df_list, axis=1, join="outer")
        # np.full(len(X), self.f0x[0])はXに大きさを揃えるための処理
        pred = gamma_result_df.sum(axis=1)+np.full(len(X), self.f0x[0])
        
        return pred

　各関数は、下記のように使います。

from sklearn.tree import DecisionTreeRegressor
from sklearn.model_selection import train_test_split

# 訓練データとテストデータに分割
X_train, X_test, y_train, y_test = train_test_split(df[features_name], df["Y"], test_size=0.2, random_state=1)

gbh = GradientBoostingHandmade(m=100, WeakLearner=DecisionTreeRegressor, max_depth=1)
gbh.fit(X_train, y_train)
pred = gbh.predict(X_test)

# 予測結果の表示（テストデータのyと予測結果の散布図を表示）
## MSEを使った評価は本記事の最後に述べる
import matplotlib.pyplot as plt
plt.scatter(y_test, pred)

実装の肝

　図１中の２－bから２－dを実装することが非常に難しく感じました。以下のことに気が付いてなんとか実装できたという感じです。２－dでは漸化式

$\begin{eqnarray} f_m(x) = f_{m-1}(x)+T(x; \hat{\Theta}_m) \end{eqnarray}$

を計算しています。しかし、この式はもともとは

$\begin{eqnarray} f_m(x) &=& \displaystyle{\sum_{i=1}^N}T(x; \hat{\Theta}_m) \\ &=& \displaystyle{\sum_{k=1}^m} \displaystyle{\sum_{j=1}^J} \gamma_{jk}I(x \in R_{jk}) \tag{1} \end{eqnarray}$

から導出されたものです。実装では（１）式を使った方が私としてはわかりやすかったので（１）式を使って実装しています。

　（１）式を使った実装の手順を説明していきます。シンプルに考えるために、弱分類器が３個、つまり $m=3$ の場合を考えてみます。ある $x$ について、１個目、２個目、３個目の弱分類器それぞれで、 $x \in R _ {j _ 1 1}$ 、 $x \in R _ {j _ 2 2}$ 、 $x \in R _ {j _ 3 3}$ に割り振られたとします。このとき、（１）式を $\sum$ を使わずに書いてみると、

$\begin{eqnarray} f_m(x) &=& \gamma_{j_1 1}+\gamma_{j_2 2}+\gamma_{j_3 3} \tag{2} \end{eqnarray}$

となります。（２）式の右辺を見てもらえるとわかるように $f _ m(x)$ は、各弱分類器が割り振った解領域のパラメータ $\gamma$ の和になっています。図２はこれをイメージ化した図です。解領域と $\gamma$ が対応しています。

　以上のことを利用して、メソッドpredictを実装しました。

実験

　ここでは、私が実装したGradientBoostingHandmadeとライブラリXGBRegressorの結果を比較してみます。当然、ライブラリXGBRegressorの方が実行時間が短く、精度も高いと予想されます。しかし、比較してみることで、自分が実装したGradientBoostingHandmadeが適切に実装されているかどうかを確かめることができると考えます。

サンプルデータ

　今回は、サンプルデータを自分で用意することにしました。サンプルデータは以下の条件で作成します。

特徴量
- 標準正規分布にしたがう10個の独立な変数。つまり、一つのサンプルに対して、 $X _ i \sim N(0, 1)\ \ (i=1, \cdots, 10)$ とする。
目的変数
- 各特徴量を二乗して和を取った。つまり、 $Y = \displaystyle{\sum _ {i=1} ^ {10}} X _ i ^ 2$ とする。
サンプルサイズ
- 2000個とした。

サンプルデータを作成するスクリプトは、下記の通りです。

# データの作成に関するパラメータ
mu, sigma = 0, 1 # 平均と標準偏差を設定
sample = 2000 # サンプルサイズ
J = 10 # 特徴量の個数
features_name = ["X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9", "X10"]

# データの作成
## 特徴量の作成
import numpy as np

X = np.random.normal(mu, sigma, size=(sample, J)) # 正規分布から乱数を生成

import pandas as pd
df = pd.DataFrame(X, columns=features_name)

# データの作成
### 目的変数の作成
X2 = np.apply_along_axis(lambda x: np.sum(x ** 2), axis=1, arr=X) # 自由度カイ二乗分布にしたがう乱数の作成

### カイ二乗分布の中央値を計算
from scipy.stats import chi2
median = chi2.ppf(0.5, J)  # 中央値

### カイ二乗分布の中央値より大きいと1, 小さいと0とする
Y = X2
df["Y"] = Y.astype("float")

モデル作成

　弱分類器を1個から100個まで増やしながらMSEの変化を見ていきます。自作のGradientBoostingHandmadeは下記のように書きました。

from sklearn.tree import DecisionTreeRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

# データの準備
X_train, X_test, y_train, y_test = train_test_split(df[features_name], df["Y"], test_size=0.2, random_state=1)

# MSEの結果を格納するリスト
MSE_by_GradientBoostingHandmade_list = []

for m in range(1, 100):
    gbh = GradientBoostingHandmade(m, WeakLearner=DecisionTreeRegressor, max_depth=1)
    gbh.fit(X_train, y_train)
    pred = gbh.predict(X_test)
    MSE_by_GradientBoostingHandmade_list.append(mean_squared_error(y_test, pred))

同様にxgboostを使ってみると下記のように実装できました。パラメータは基本的にはdefaultを使い、max_depthとn_estimatorsは上の条件と合わせた設定にしました。

import xgboost as xgb

MSE_by_xgboost_list = []
for m in range(1, 100):
    # xgboostモデルの作成
    reg = xgb.XGBRegressor(max_depth=1, n_estimators=m, objective="reg:squarederror")
    reg.fit(X_train, y_train)
    pred = reg.predict(X_test)
    MSE_by_xgboost_list.append(mean_squared_error(y_test, pred))

結果

　まず計算時間についてですが、xgboostの方が圧倒的に早かったです。私の持っているPC（CPU：intel core i3 10TH GEN）でxgboostは3秒でした。一方で自作のGradientBoostingHandmadeは11分21秒かかりました。圧倒的な差ですね、、、。　また、MSEの変化は図３のようになります。

GradientBoostinghandmadeの方が全体的にMSEが小さいです。ただし、これはxgboostでハイパーパラメータのチューニングをしていないからだと思います。とりあえず、弱分類器の個数を増やせばMSEが小さくなっていることを確認できて良かったです*1。

まとめ

　今回の記事では、勾配ブースティング木を実装してみました。スクラッチ実装することで勾配ブースティング木の仕組みの基本を理解できました。しかし、計算時間については自力で実装した場合、ライブラリよりも遥かに多くの時間を必要としてしまうこともわかりました。今後の課題としては、なぜライブラリxgboostは高速に計算できるのか調べたいと思っています。

*1:今回のデータセットであれば、弱分類器100個程度じゃ過学習しないと理解

初心者データサイエンティストの備忘録

調べたことは全部ここに書いて自分の辞書を作る

勾配ブースティング木を深く理解する～実装編～