勾配ブースティング

直感 #

絵を描くとき、まず全体を 1 色で塗り、次に「暗すぎる場所」と「明るすぎる場所」だけを薄い色で重ね塗りし、それを何十回も繰り返して仕上げる。勾配ブースティングはこの重ね塗りです。最初の 1 色が定数の予測 \(F_0\)、「暗すぎる・明るすぎる」が各点の残差、薄い色が学習率 \(\eta\) を掛けた小さな木の出力に当たります。1 回の重ね塗りは短い木 1 本なので粗い階段状です。それでも 50 回、100 回と重ねると、波形のような複雑な形にも追従します。直すのはいつも「いま一番外れている場所」なので、残差という坂を一歩ずつ降りている、つまり関数空間での勾配降下になっています。

何が起きているか（ステップで見る） #

勾配ブースティングは「いまの予測の間違い（残差）だけを、次の小さな木で直す」を繰り返します。縦の線が各点の残差 \(y_i - F_{m-1}(x_i)\)、青い階段状の線が現在の予測 \(F_m\) です。段階を 1 つ進めるごとに残差の縦線が縮んでいきます。

見るポイント：1 段階で足されるのは、1 回だけ分岐する短い木（スタンプ）の出力に学習率 0.5 を掛けたもので、灰色の線が前の段階の予測なので、1 段階でどこがどれだけ動いたかが分かります。残差の縦線が長い場所ほど優先的に直されるのが「勾配に沿って降りる」という意味です。段階を増やすほど訓練データへの当てはまりは良くなりますが、ノイズまで追い始めると過学習するので、n_estimators と learning_rate の組み合わせでこの速度を調整します。

式で見ると何を足しているか #

関数空間での勾配降下 #

勾配ブースティングは、任意の微分可能な損失関数\(L(y, F(\mathbf{x}))\)に対して「関数空間での勾配降下法」を行います。パラメータの値を少しずつ更新する通常の勾配降下と違い、更新されるのは予測関数\(F(\mathbf{x})\)そのもので、木を 1 本足すことが 1 ステップに当たります。予測関数は次のように加法的に構築します。

ここで\(h_m\)は\(m\)番目の弱学習器（通常は浅い決定木）、\(\eta\)は学習率（縮小係数）です。

疑似残差 #

各ステップで弱学習器がフィットする対象は「疑似残差」と呼ばれる損失関数の負の勾配です。

二乗誤差\(L = \frac{1}{2}(y - F)^2\)の場合、疑似残差は単純な残差\(y_i - F_{m-1}(\mathbf{x}_i)\)と一致し、ステップ再生の縦線がそのまま次の木の学習対象になります。絶対誤差や Huber では下の表のように形が変わり、外れ値の縦線が長くても次の木に伝わる値は \(\pm 1\) や \(\pm\delta\) に抑えられます。

アルゴリズムの手順 #

手順は 3 段階です。まず初期モデルを定数\(F_0(\mathbf{x}) = \arg\min_\gamma \sum_{i=1}^{n} L(y_i, \gamma)\)に置きます。二乗誤差なら\(y\)の平均、絶対誤差なら中央値です。次に\(m = 1, \ldots, M\)について、全サンプルの疑似残差\(r_{im}\)を計算し、決定木\(h_m\)を\({r_{im}}\)にフィットさせ、各リーフ領域\(R_{jm}\)で損失を最小にするステップサイズ\(\gamma_{jm} = \arg\min_\gamma \sum_{x_i \in R_{jm}} L(y_i, F_{m-1}(\mathbf{x}i) + \gamma)\)を求め、予測を\(F_m(\mathbf{x}) = F{m-1}(\mathbf{x}) + \eta \sum_j \gamma_{jm} \mathbf{1}(\mathbf{x} \in R_{jm})\)と更新します。最後に\(M\)本そろった\(F_M(\mathbf{x})\)を最終モデルとして出力します。GradientBoostingRegressor では\(M\)が n_estimators、\(\eta\)が learning_rate です。

損失関数と疑似残差の対応 #

縮小（Shrinkage）の理論的根拠 #

学習率\(\eta < 1\)を掛けて各木の寄与を抑えると、汎化性能が向上することが実験的に知られています。Friedman (2001) は\(\eta \approx 0.1\)を推奨しています。直感的には、各ステップで控えめに更新することで、後続の木が異なる側面から誤差を補正する余地が残り、より多様なアンサンブルが構築されます。

Python で確かめる #

以下では sin と cos を重ねた波形に GradientBoostingRegressor を当て、まず基本の当てはまりを見てから、loss（外れ値への反応）、n_estimators（本数）、learning_rate（歩幅）を 1 つずつ動かして訓練データへの当てはまりがどう変わるかを図にします。

1
2
3

import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor

訓練データに回帰モデルを当てはめる #

実験用に、$\sin x$ と $\cos 4x$ を足して 10 倍し、傾きと一様乱数のノイズを加えた 500 点の波形データを用意します。木を 50 本、学習率 0.5 で当てて、赤い予測線がどこまで波形を追えるかを見ます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

# 訓練データ
X = np.linspace(-10, 10, 500)[:, np.newaxis]
noise = np.random.rand(X.shape[0]) * 10

# 目的変数
y = (
    (np.sin(X).ravel() + np.cos(4 * X).ravel()) * 10
    + 10
    + np.linspace(-10, 10, 500)
    + noise
)

# 回帰モデルを作成
reg = GradientBoostingRegressor(
    n_estimators=50,
    learning_rate=0.5,
)
reg.fit(X, y)
y_pred = reg.predict(X)

# 訓練データへのフィッティング具合を確認する
plt.figure(figsize=(10, 5))
plt.scatter(X, y, c="k", marker="x", label="訓練データ")
plt.plot(X, y_pred, c="r", label="最終的に作成されたモデルの予測", linewidth=1)
plt.xlabel("x")
plt.ylabel("y")
plt.title("訓練データへのフィッティング具合")
plt.legend()
plt.show()

赤い線は 12 個以上ある山と谷のほぼすべてを追っています。既定の max_depth=3 の木は葉が最大 8 枚なので 1 本では階段 8 段分しか表せませんが、50 本を学習率 0.5 で重ねることで細かい波形になっています。

損失関数の結果への影響 #

loss を "squared_error"、"absolute_error"、"huber"、"quantile" と変えて、訓練データへの当てはまりがどう変わるかを比べます。今回は外れ値入りです。80 点ごとに 1 点、値を 60〜80 ほど上に飛ばした点を 7 個混ぜてあります。二乗誤差は外れ値の残差を 2 乗して評価するので、予測線が外れ値に向かって跳ねます。一方 absolute_error と huber は外れ値の疑似残差が符号や \(\delta\) に抑えられるため、7 点の外れ値を予測しに行かず、周囲の正常な点に沿った線のままです。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

# 訓練データ
X = np.linspace(-10, 10, 500)[:, np.newaxis]

# 外れ値を用意
noise = np.random.rand(X.shape[0]) * 10
for i, ni in enumerate(noise):
    if i % 80 == 0:
        noise[i] = 70 + np.random.randint(-10, 10)

# 目的変数
y = (
    (np.sin(X).ravel() + np.cos(4 * X).ravel()) * 10
    + 10
    + np.linspace(-10, 10, 500)
    + noise
)

for loss in ["squared_error", "absolute_error", "huber", "quantile"]:
    # 回帰モデルを作成
    reg = GradientBoostingRegressor(
        n_estimators=50,
        learning_rate=0.5,
        loss=loss,
    )
    reg.fit(X, y)
    y_pred = reg.predict(X)

    # 訓練データへのフィッティング具合を確認する
    plt.figure(figsize=(10, 5))
    plt.scatter(X, y, c="k", marker="x", label="訓練データ")
    plt.plot(X, y_pred, c="r", label="最終的に作成されたモデルの予測", linewidth=1)
    plt.xlabel("x")
    plt.ylabel("y")
    plt.title(f"訓練データへのフィッティング具合, loss={loss}", fontsize=18)
    plt.legend()
    plt.show()

squared_error の図では外れ値 7 点の真下で赤線が上に跳ね、周囲の正常な点から離れています。absolute_error と huber では同じ場所でも線が波形に留まり、quantile（alpha=0.9）では線全体が点の帯の上端寄りに走ります。外れ値を含むデータでは huber が二乗誤差の滑らかさと絶対誤差の頑健さの中間になるので、まず試す候補です。

n_estimators の結果への影響 #

n_estimators を 5 本から 100 本まで 5 本刻みで増やし、訓練データの平均二乗誤差（MSE）がどう下がるかを見ます。木を足すたびに残差は縮みます。したがって MSE は単調に下がりますが、ある本数を超えると改善幅が小さくなり頭打ちになります。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

from sklearn.metrics import mean_squared_error as MSE

# 訓練データ
X = np.linspace(-10, 10, 500)[:, np.newaxis]
noise = np.random.rand(X.shape[0]) * 10

# 目的変数
y = (
    (np.sin(X).ravel() + np.cos(4 * X).ravel()) * 10
    + 10
    + np.linspace(-10, 10, 500)
    + noise
)

# n_estimatorsを変えてモデルを作成してみる
n_estimators_list = [(i + 1) * 5 for i in range(20)]
mses = []
for n_estimators in n_estimators_list:
    # 回帰モデルを作成
    reg = GradientBoostingRegressor(
        n_estimators=n_estimators,
        learning_rate=0.3,
    )
    reg.fit(X, y)
    y_pred = reg.predict(X)
    mses.append(MSE(y, y_pred))

# n_estimatorsを変えた時のmean_squared_errorをプロット
plt.figure(figsize=(10, 5))
plt.plot(n_estimators_list, mses, "x")
plt.xlabel("n_estimators")
plt.ylabel("Mean Squared Error(training data)")
plt.title(f"訓練データへのフィッティング具合", fontsize=18)
plt.grid()
plt.show()

learning_rate の結果への影響 #

n_estimators=30 に固定し、learning_rate を 0.1 から 2.0 まで 0.1 刻みで変えて、訓練データの MSE の対数を比べます。0.1〜0.3 のように小さすぎると 30 本では残差を直しきれず MSE が高止まりし、1.0 を超えると 1 歩が大きすぎて残差を行き過ぎ、本数を増やしても MSE が下がらず発散することが確認できます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

# n_estimatorsを変えてモデルを作成してみる
learning_rate_list = [np.round(0.1 * (i + 1), 1) for i in range(20)]
mses = []
for learning_rate in learning_rate_list:
    # 回帰モデルを作成
    reg = GradientBoostingRegressor(
        n_estimators=30,
        learning_rate=learning_rate,
    )
    reg.fit(X, y)
    y_pred = reg.predict(X)
    mses.append(np.log(MSE(y, y_pred)))

# n_estimatorsを変えた時のmean_squared_errorをプロット
plt.figure(figsize=(10, 5))
plt_index = [i for i in range(len(learning_rate_list))]
plt.plot(plt_index, mses, "x")
plt.xticks(plt_index, learning_rate_list, rotation=90)
plt.xlabel("learning_rate", fontsize=15)
plt.ylabel("log(Mean Squared Error) (training data)", fontsize=15)
plt.title(f"訓練データへのフィッティング具合", fontsize=18)
plt.grid()
plt.show()

対数 MSE は 0.1 から 0.5 あたりまで下がり、そこから 1.0 付近までは横ばい、1.2〜1.5 を超えると急に上がります。小さすぎる学習率は 30 本では足りず、大きすぎる学習率は残差を行き過ぎて発散する、という両側の失敗が 1 枚の図に出ています。

まとめ #

勾配ブースティングは、定数の初期予測から始めて「残差（損失の負の勾配）に浅い木を当て、学習率 \(\eta\) を掛けて足す」を n_estimators 回繰り返す手法で、50 本程度でも sin と cos を重ねた波形に追従できました。learning_rate と n_estimators はトレードオフで、学習率 0.1〜0.3 なら本数を 100 本以上、0.5 なら 50 本程度が目安になり、1.0 を超えると発散します。loss を huber や absolute_error にすると外れ値の疑似残差が頭打ちになり、squared_error のように予測線が外れ値へ跳ねることを防げます。過学習対策は max_depth を 3〜5 に留めることです。本数は n_iter_no_change による早期終了で決めるのが実践的で、木が 1 本ずつ足されていく様子を図で追いたい場合は 勾配ブースティングの可視化 に続きます。

推定器数と予測曲線 #

スライダーで n_estimators を 1 本から 200 本まで動かすと、階段状の予測線が滑らかになっていく一方で、train_mse は下がり続けるのに test_mse が 50 本あたりから上がり始める様子を確認できます。