GBDT（勾配ブースティング決定木）の理論とフルスクラッチ実装とその解説

はじめに

引き続き、今度はGBDTに関しての記事を書いていこうと思います。GBDTといえば、最近も機械学習の第一手段として人気の、XGboostやLightGBMの元となっているモデルです。（LightGBMはXGboostを基に改善したモデルですが）概念を理解している人は多くいるものの、図解的な概念での理解にとどまっており、理解しながら使用してる人は実際は少ないのではないかと思っているモデルでもあります。バギングに比べ、ブースティングは少々理解するのが難しい部分でもあるので、そんなGBDTに関して理論を中心に、フルスクラッチのコードを添えて記事を書いていこうと思います。
元々Xgboostまでは触れようと思っていたのですが、自分のなまけ癖が限界を迎えようとしており、扱えないかもしれませんが許してください。前回のランダムフォレストみたいなグダグダ解説にならないように気を付けながら進めたいと思います。

GBDT（勾配ブースティング決定木）とは

1. 概要

勾配ブースティング決定木（GBDT）は、弱学習器（通常は決定木）を逐次的に組み合わせて強力な予測モデルを構築する機械学習アルゴリズムです。各ステップで、前のモデルの残差（実際の値と予測値の差）を予測する新しい決定木を追加することで、モデルを改善していきます。

2.数学的な説明

目的関数：

$L(\phi) = \sum_{i=1}^n l(y_i, \phi(x_i))$

この目的関数は、モデル全体の性能を評価するために使用されます。GBDTの目標は、この関数を最小化することです。
構成要素：
1.

$L(\phi)$ : 全体の損失関数：これは、すべてのデータポイントにわたる個々の損失の合計を表します。 $\phi$ はモデル全体を表す関数です。

$l(y_i, \phi(x_i))$ : 各データポイントの損失：これは、個々のデータポイント $i$ に対する予測値と実際の値の間の誤差を測定します。

$y_i$ : 実際の値（目標変数）：これは、モデルが予測しようとしている実際の値です。

$\phi(x_i)$ : モデルの予測値：これは、入力 $x_i$ に対するモデルの予測出力です。

$n$ : データポイントの数：訓練データセット内の全サンプル数です。

　　
一般的に使用される損失関数：
1. 回帰問題の場合

二乗誤差（Mean Squared Error, MSE）: $l(y_i, \phi(x_i)) = (y_i - \phi(x_i))^2$

絶対誤差（Mean Absolute Error, MAE）: $l(y_i, \phi(x_i)) = |y_i - \phi(x_i)|$

Huber損失：MSEとMAEのハイブリッド。外れ値に対してより頑健です。

　　
2. 二値分類の問題の場合

対数損失（Log Loss）または二値交差エントロピー： $l(y_i, \phi(x_i)) = -[y_i \log(\phi(x_i)) + (1-y_i) \log(1-\phi(x_i))$ ]

ヒンジ損失（Hinge Loss）： $l(y_i, \phi(x_i)) = \max(0, 1 - y_i\phi(x_i))$

　　
3. 多クラス分類問題の場合

多クラス対数損失（Multiclass Log Loss）： $l(y_i, \phi(x_i)) = -\sum_{k=1}^K y_{ik} \log(\phi_k(x_i))$ ここで、 $K$ はクラスの数、 $y_{ik}$ は $i$ 番目のサンプルが $k$ クラスに属する場合1、そうでなければ0です。

　　
GBDTの良いところは、問題の性質や要求される特性に対応して、様々な損失関数を使い、それを最小化するように、モデルが作成される点にあります。
目的関数を最小化することで、GBDTは予測と実際の値の差を最小限に抑えるようにモデルを調整します。各イテレーションで、モデルは前のステップの残差（実際の値と予測値の差）を予測する新しい決定木を追加することで、徐々にこの目的を達成していきます。

各ステップの解説

Step1

初期モデル $F_0(x)$ の定義

$F_0(x) = \arg\min_{\gamma} \sum_{i=1}^n L(y_i, \gamma)$

この段階では、最も単純な予測モデルを作成します。通常、全データポイントの平均値や中央値を使用します。これは、損失関数 $L$ を最小化する定数 $\gamma$ を見つけることに相当します。例えば、二乗誤差を損失関数として使用する場合、この初期モデルは全てのデータポイントの平均値になります。

Step2

モデルの逐次的改善（For $m = 1$ to $M$ ）

　　
a)負の勾配の計算

$r_{im} = -\left[\frac{\partial L(y_i, F(x_i))}{\partial F(x_i)}\right]{F(x)=F{m-1}(x)}$

この段階では、現在のモデル $F_{m-1}(x)$ に対する損失関数の勾配を計算します。負の勾配を使用することで、損失を減少させる方向を特定します。これは、現在のモデルの「誤差」または「残差」と解釈できます。

b)決定木 $h_m(x)$ のフィッティング

決定木 $h_m(x)$ を負の勾配 $r_{im}$ にフィットさせます。つまり、前のステップで計算した「誤差」を予測する決定木を作成します。この決定木は、現在のモデルの弱点を補うように設計されています。

c) 最適な係数 $\gamma_m$ の決定

$\gamma_m = \arg\min_{\gamma} \sum_{i=1}^n L(y_i, F_{m-1}(x_i) + \gamma h_m(x_i))$

新しく作成した決定木 $h_m(x)$ をモデルに追加する際の最適な重み（係数） $\gamma_m$ を決定します。これは、新しい決定木を追加した後の全体の損失を最小化するように選択されます。

　　
d)モデルの更新

$F_m(x) = F_{m-1}(x) + \gamma_m h_m(x)$

最後に、最適な係数を用いて新しい決定木をモデルに追加します。これにより、モデルは徐々に改善されていきます。

このプロセス（a〜d）を $M$ 回繰り返すことで、モデルは段階的に複雑になり、予測精度が向上していきます。

Step3

最終モデル $F_M(x)$ の出力

$M$ 回の反復が完了した後、最終的なモデル $F_M(x)$ が得られます。このモデルは、初期モデルと $M$ 個の重み付けされた決定木の和として表現されます

$F_M(x) = F_0(x) + \sum_{m=1}^M \gamma_m h_m(x)$

この最終モデルは、元の問題に対する予測を行うために使用されます。
GBDTの特徴は、各ステップで前のモデルの誤差を修正する新しい決定木を追加することで、徐々にモデルを改善していく点にあります。これにより、複雑な非線形関係を捉えつつ、過学習を抑制することができます。

重要な特徴

a) 正則化：過学習を防ぐために、学習率 $\eta$ を導入します

$F_m(x) = F_{m-1}(x) + \eta \cdot \gamma_m h_m(x)$

正則化は過学習（オーバーフィッティング）を防ぐための重要な技術です。GBDTにおいて、学習率 $\eta$ （0 < $\eta$ ≤ 1）を導入することで正則化を実現します。

これは各イテレーションでのモデルの更新量を制御する概念です。

小さな $\eta$ を使用することで、モデルの更新を慎重に行い、急激な変化を抑制します。これには、過学習のリスク低減やモデルの汎化性能の向上、アンサンブル全体の安定性の向上などに寄与します。

　　
b)特徴量の重要度：モデル全体での各特徴量の重要度を計算できます。

$\text{Importance}k = \sum{j=1}^J I_j^2 \cdot \mathbf{1}(v_j = k)$

$\text{Importance}_k$ : 特徴量 $k$ の重要度

$I_j$ : $j$ 番目の分割（ノード）の重要度

$v_j$ : $j$ 番目の分割に使用された変数（特徴量）

$\mathbf{1}(v_j = k)$ : 指示関数。 $v_j = k$ の場合1、それ以外は0

計算プロセス

モデル内の全ての決定木を走査します。
各分割（ノード）で、その分割の重要度 $I_j$ を計算します。
特徴量 $k$ が使用された分割のみを考慮し、それらの重要度の二乗和を取ります。

分割の重要度 $I_j$ の計算

通常、分割によって減少した不純度（例：ジニ不純度や情報利得）を使用します。
分割によって大きく不純度が減少した場合、その分割（および使用された特徴量）は重要と見なされます。
重要なハイパーパラメータ
木の深さ
木の数（反復数）
学習率
サンプリング率
正則化パラメータ以上が重要なハイパーパラメータになります。
GBDTの問題点
ハイパーパラメータのチューニングが必要
解釈性が低い（個々の木は理解しやすいが、全体としては複雑な構造になっているため）
計算コストが高い
フルスクラッチ実装

GBDTのフルスクラッチ実装は以下です。

import numpy as np

class DecisionNode:
    def __init__(self, feature_index=None, threshold=None, left=None, right=None, value=None):
        # 決定木の各ノードを表現するクラス
        self.feature_index = feature_index  # 分割に使用する特徴量のインデックス
        self.threshold = threshold  # 分割の閾値
        self.left = left  # 左の子ノード（閾値未満のデータ）
        self.right = right  # 右の子ノード（閾値以上のデータ）
        self.value = value  # 葉ノードの場合の予測値（回帰問題では平均値）

class DecisionTreeRegressor:
    def __init__(self, max_depth=None, min_samples_split=2):
        # 決定木回帰器の初期化
        self.max_depth = max_depth  # 木の最大深さ（過学習防止のため）
        self.min_samples_split = min_samples_split  # ノードを分割するための最小サンプル数
        self.root = None  # ルートノード（木の開始点）
        self.feature_importances_ = None  # 各特徴量の重要度を格納する配列

    def fit(self, X, y):
        # 決定木のトレーニングを行うメソッド
        self.n_features = X.shape[1]  # 特徴量の数を取得
        self.feature_importances_ = np.zeros(self.n_features)  # 特徴量重要度の初期化
        self.root = self._grow_tree(X, y)  # 再帰的に木を成長させる

    def _grow_tree(self, X, y, depth=0):
        # 再帰的に決定木を成長させるメソッド
        n_samples, n_features = X.shape
        
        # 停止条件のチェック
        if (self.max_depth is not None and depth >= self.max_depth) or \
           n_samples < self.min_samples_split:
            return DecisionNode(value=np.mean(y))  # 葉ノードを返す（予測値は平均）

        best_feature, best_threshold = self._best_split(X, y)  # 最適な分割点を見つける
        
        # 分割できない場合（例：すべての特徴量が同じ値の場合）
        if best_feature is None:
            return DecisionNode(value=np.mean(y))  # 葉ノードを返す

        # データを左右に分割
        left_indices = X[:, best_feature] < best_threshold
        right_indices = ~left_indices

        # 左右の子ノードを再帰的に成長させる
        left = self._grow_tree(X[left_indices], y[left_indices], depth + 1)
        right = self._grow_tree(X[right_indices], y[right_indices], depth + 1)

        # 特徴量の重要度を更新（二乗誤差の減少量を加算）
        self.feature_importances_[best_feature] += (
            self._calculate_mse_reduction(y, y[left_indices], y[right_indices])
        )

        # 新しい内部ノードを返す
        return DecisionNode(feature_index=best_feature, threshold=best_threshold,
                            left=left, right=right)

    def _best_split(self, X, y):
        # 最適な分割点を見つけるメソッド
        m = y.size
        if m <= 1:
            return None, None  # 分割できない場合

        # 親ノードの二乗誤差を計算
        parent_var = np.var(y) * m  # 分散 * サンプル数 = 二乗誤差の合計

        best_var_reduction = 0  # 最大の分散減少量（初期値）
        best_feature = None  # 最適な特徴量（初期値）
        best_threshold = None  # 最適な閾値（初期値）

        for feature in range(self.n_features):
            thresholds = np.unique(X[:, feature])  # ユニークな閾値候補を取得
            for threshold in thresholds:
                # データを左右に分割
                left = y[X[:, feature] < threshold]
                right = y[X[:, feature] >= threshold]
                
                if len(left) == 0 or len(right) == 0:
                    continue  # 片方のノードにデータがない場合はスキップ

                # 子ノードの二乗誤差を計算し、分散の減少量を求める
                var_reduction = parent_var - (np.var(left) * len(left) + 
                                              np.var(right) * len(right))

                # より良い分割が見つかった場合、更新
                if var_reduction > best_var_reduction:
                    best_var_reduction = var_reduction
                    best_feature = feature
                    best_threshold = threshold

        return best_feature, best_threshold

    def _calculate_mse_reduction(self, parent, left, right):
        # 二乗誤差の減少量を計算するメソッド
        m = len(parent)
        m_l, m_r = len(left), len(right)

        # 各ノードの二乗誤差を計算
        mse_parent = np.sum((parent - parent.mean()) ** 2)
        mse_left = np.sum((left - left.mean()) ** 2)
        mse_right = np.sum((right - right.mean()) ** 2)

        # 二乗誤差の減少量を計算
        mse_reduction = mse_parent - (mse_left + mse_right)
        return mse_reduction / m  # サンプル数で正規化

    def predict(self, X):
        # 新しいデータに対する予測を行うメソッド
        return np.array([self._traverse_tree(x, self.root) for x in X])

    def _traverse_tree(self, x, node):
        # 決定木を走査して予測値を得るメソッド
        if node.value is not None:
            return node.value  # 葉ノードの場合、その値を返す

        if x[node.feature_index] < node.threshold:
            return self._traverse_tree(x, node.left)  # 左の子ノードへ
        return self._traverse_tree(x, node.right)  # 右の子ノードへ

class GBDT:
    def __init__(self, n_estimators=100, learning_rate=0.1, max_depth=3):
        # GBDTモデルの初期化
        self.n_estimators = n_estimators  # ブースティングの反復回数
        self.learning_rate = learning_rate  # 学習率（ステップサイズ）
        self.max_depth = max_depth  # 各決定木の最大深さ
        self.trees = []  # 決定木のリスト

    def fit(self, X, y):
        # GBDTモデルのトレーニングを行うメソッド
        self.initial_prediction = np.mean(y)  # 初期予測値（全体の平均）
        F = np.full(len(y), self.initial_prediction)  # 現在の予測値

        for m in range(self.n_estimators):
            residuals = y - F  # 残差（実際の値 - 予測値）を計算
            
            # 残差に対して決定木を学習
            tree = DecisionTreeRegressor(max_depth=self.max_depth)
            tree.fit(X, residuals)
            
            # 予測値を更新（学習率を考慮）
            update = self.learning_rate * tree.predict(X)
            F += update
            
            self.trees.append(tree)  # 学習した木を保存

    def predict(self, X):
        # 新しいデータに対する予測を行うメソッド
        predictions = np.full(len(X), self.initial_prediction)  # 初期予測値で初期化
        for tree in self.trees:
            # 各木の予測を加算（学習率を考慮）
            predictions += self.learning_rate * tree.predict(X)
        return predictions

    def feature_importance(self):
        # 特徴量の重要度を計算するメソッド
        importances = np.zeros(self.trees[0].n_features)  # 重要度を格納する配列
        for tree in self.trees:
            importances += tree.feature_importances_  # 各木の重要度を加算
        return importances / len(self.trees)  # 木の数で正規化

決定木の成長プロセス（_grow_tree メソッド）
最適な分割点の選択方法（_best_split メソッド）
二乗誤差の減少量の計算（_calculate_mse_reduction メソッド）
GBDTの学習プロセス（fit メソッド）
予測の方法（predict メソッド）
特徴量重要度の計算（feature_importance メソッド）
実際にサンプルデータを使って予測するコードは以下です。

import numpy as np
from sklearn.datasets import make_regression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score

# 前述のGBDTとDecisionTreeRegressorクラスの定義がここにあると仮定します

# サンプルデータの生成
X, y = make_regression(n_samples=1000, n_features=10, noise=0.1, random_state=42)

# データを訓練セットとテストセットに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# GBDTモデルの初期化
gbdt = GBDT(n_estimators=100, learning_rate=0.1, max_depth=3)

# モデルの学習
gbdt.fit(X_train, y_train)

# テストデータで予測
y_pred = gbdt.predict(X_test)

# 評価指標の計算
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)

print(f"Mean Squared Error: {mse:.4f}")
print(f"R-squared Score: {r2:.4f}")

# 特徴量の重要度を表示
feature_importance = gbdt.feature_importance()
for i, importance in enumerate(feature_importance):
    print(f"Feature {i}: {importance:.4f}")

# 新しいデータポイントでの予測例
new_data_point = np.random.rand(1, 10)  # 10個の特徴を持つ1つのデータポイント
prediction = gbdt.predict(new_data_point)
print(f"\n新しいデータポイントの予測値: {prediction[0]:.4f}")

# 学習曲線の可視化（オプション）
import matplotlib.pyplot as plt

# 訓練データのサイズを変えて学習と評価を行う
train_sizes = np.linspace(0.1, 0.9, 5)  # 0.1から0.9までの5点を生成
train_errors = []
test_errors = []

for size in train_sizes:
    X_train_subset, _, y_train_subset, _ = train_test_split(X_train, y_train, train_size=size, random_state=42)
    gbdt = GBDT(n_estimators=100, learning_rate=0.1, max_depth=3)  # 新しいモデルを作成
    gbdt.fit(X_train_subset, y_train_subset)
    
    train_pred = gbdt.predict(X_train_subset)
    test_pred = gbdt.predict(X_test)
    
    train_errors.append(mean_squared_error(y_train_subset, train_pred))
    test_errors.append(mean_squared_error(y_test, test_pred))

# 学習曲線のプロット
plt.figure(figsize=(10, 6))
plt.plot(train_sizes, train_errors, 'r-', label='Training error')
plt.plot(train_sizes, test_errors, 'b-', label='Test error')
plt.xlabel('Training set size ratio')
plt.ylabel('Mean Squared Error')
plt.title('Learning Curve for GBDT')
plt.legend()
plt.show()

結果は以下のようになります。

最後に

今回の記事では、あまりとっ散らかった説明にならず、理解しやすい説明ができたのではなきかと思っております。意外とGBDTに関して、概念的な理解にとどまり、理論を理解していた人は少なかったのではないでしょうか？理論を理解できている状態でツールを使用するのと、理解していない状態でツールを使用するのは大きく違ってくると個人的には思っています。また、特徴量重要度含め、SHAPで解釈する場合にも多重共線性には注意してください（予測や分類タスクの精度には関係しませんが）フルスクラッチのコード含め、初学者の方の理解の一助になれば幸いです。次回はXGboostをやるのか？本当にやりたかった因果推論方面にもう行くのかで悩んでいますが、ご期待ください。