MLエンジニアへの道 #10 - 勾配ブースト

ML

Last Edited: 8/13/2024

このブログ記事では、勾配ブーストの仕組みについて紹介します。

ML

決定木が陥りがちな過学習を防ぐのに多くのメカニズム（剪定など）がありますが、それでも決定木は高次元で複雑なデータに苦労する傾向があります。 これは、適切にデータをモデル化するためには木がより複雑になる必要があり、その結果として過学習が発生する可能性があるためです。

バギングとブースティング

モデルのシンプルさを維持しながら、高次元で複雑なデータにも適合させるにはどうすればよいでしょうか？この問題の1つの解決策はバギングです。 バギングでは、多くの単純な木を作成し、その結果を組み合わせて予測を行います。最も直感的な方法は、ランダムサンプルとランダムな特徴を使用して多く の木を訓練し、その後、過半数の投票に基づいて予測を行うことです。これをランダムフォレストと呼び、非常に効果的な手法として知られています （ここでは詳しく説明しません）。

Bag vs Boost

もう一つの小さな木を組み合わせる方法として、小さな木の連鎖を作成し、それぞれの木が前の木が犯した誤りを修正しようとする方法があげられます。 この手法はブースティングと呼ばれ、高次元で複雑なデータの学習に非常に効果的であることが分かっています。今回の記事では、このブースティング を理解するために、勾配ブーストと呼ばれるブースティングの例を見ていきたいと思います。

勾配ブースト (Gradient Boost)

勾配ブーストのアルゴリズムは次のように定義されます。

入力: データ (xi,yi)I=1n{(x_i, y_i)}_{I=1}^{n} と微分可能なコスト関数 L(y,F(x))L(y, F(x))。

ステップ1: 定数値でモデルを初期化する: F0(x)=arg min⁡γ∑I=1nL(yi,γ)F_0(x) = \argmin_{\gamma} \sum_{I=1}^{n} L(y_i, \gamma)

ステップ2: m=1m=1 から MM まで繰り返す:

  • rim=−[∂L(yi,F(xi))∂F(xi)]F(x)=Fm−1(x)r_{im} = -[\frac{\partial L(y_i, F(x_i))}{\partial F(x_i)}]_{F(x) = F{m-1}(x)} を i=1,…,ni = 1,…,n に対して計算する。
  • rimr_{im} の値に適合する木を構築し、端末領域(リーフノード) RjmR_{jm} を j=1,…,Jmj = 1, … , J_m に対して作成する。
  • j=1,…,Jmj = 1, … , J_m に対して γjm=arg min⁡γ∑xi∈RijL(yi,Fm−1(xi)+γ)\gamma_{jm} = \argmin_{\gamma} \sum_{x_i \in R_{ij}} L(y_i, F_{m-1} (x_i) + \gamma) を計算する。
  • Fm(x)=Fm−1(x)+ν∑j=1JmγmI(x∈Rjm)F_m (x) = F_{m-1} (x) + \nu \sum_{j=1}^{J_m} \gamma_m I(x \in R_{jm}) を更新する。

ここで、xix_i と yiy_i はそれぞれデータポイント ii の説明変数と予測変数を表します。最初に定数初期予測 F0(x)F_0(x) から始め、前回の予測の残差 rimr_{im} に適合する MM 個の木を生成します。（木の適合のメカニズムは、決定木に関する以前の記事でカバーされています。）

木は、木の出力値と学習率 ν\nu の積を加えることで連鎖します。通常、MM は約100に設定され、ν\nu は0から1の値に設定され、木は8から32のリーフノードで 構成されます。多くの小さな木がこの方法で組み合わされるとき、勾配ブーストが最も効果的であることが実験的に示されています。

このアルゴリズムが、真の値と予測値の間の勾配/残差を取り、それを学習率で掛けることでモデルのパラメータを徐々に更新するという点で、勾配降下法に似ているこ とに気付くかもしれません。このような類似したメカニズムが決定木にも効果的であることがわかります。

回帰

次に、勾配ブーストが回帰でどのように機能するかを見てみましょう。アルゴリズムは複雑に見えるかもしれませんが、回帰問題に適用するのは驚くほど簡単です。 まず、コスト関数を二乗残差の半分に設定します。

L(y,F(x))=12(y−F(x))2 L(y, F(x)) = \frac{1}{2} (y-F(x))^2

残差の二乗の半分を使用するのは、γ\gamma を最小化するための微分を取るときに定数2がキャンセルできるからです。次に、ステップ1に進み、次の条件を満たすように初期予測 F0(x)F_0(x) を設定します。

F0(x)=arg min⁡γ∑I=1nL(yi,γ) F_0(x) = \argmin_{\gamma} \sum_{I=1}^{n} L(y_i, \gamma)

上記の式を γ\gamma に関して微分し、それを0に設定して、コストの合計を最小化する γ\gamma を求めます。

∂∂γL(y,γ)=−(y−γ)∂∂γ∑I=1nL(yi,γ)=∑I=1nγ−yi=nγ−∑I=1nyinγ−∑I=1nyi=0γ=∑I=1nyin \frac{\partial}{\partial \gamma} L(y, \gamma) = -(y-\gamma) \\ \frac{\partial}{\partial \gamma} \sum_{I=1}^{n} L(y_i, \gamma) = \sum_{I=1}^{n} \gamma - y_i \\ = n \gamma - \sum_{I=1}^{n} y_i \\ n \gamma - \sum_{I=1}^{n} y_i = 0 \\ \gamma = \frac{\sum_{I=1}^{n} y_i}{n}

したがって、初期予測 F0(x)F_0(x) は予測変数の平均値であることが分かります。

F0(x)=∑I=1nyin(=yˉ) F_0(x) = \frac{\sum_{I=1}^{n} y_i}{n} (= \bar{y})

これで、ステップ2に進む準備ができました。まず、残差 rimr_{im} を計算する必要があります。

rim=−[∂L(yi,F(xi))∂F(xi)]F(x)=Fm−1(x)=−(Fm−1(xi)−yi)=yi−Fm−1(xi) r_{im} = -[\frac{\partial L(y_i, F(x_i))}{\partial F(x_i)}]_{F(x) = F_{m-1}(x)} \\ = -(F_{m-1}(x_i) - y_i) = y_i - F_{m-1}(x_i)

上記で計算された残差に適合する新しい木 mm を作成します。次に、各リーフノードについて、次の式を使用して γ\gamma を計算します。

γjm=arg min⁡γ∑xi∈RijL(yi,Fm−1(xi)+γ) \gamma_{jm} = \argmin_{\gamma} \sum_{x_i \in R_{ij}} L(y_i, F_{m-1} (x_i) + \gamma)

コストの合計を最小にする γ\gamma を求めるために、コストの合計を γ\gamma に関して微分します。

L(y,Fm−1(x)+γ)=12(y−(Fm−1(x)+γ))2∂∂γL=(Fm−1(x)+γ)−y∂∂γ∑xi∈RijL(yi,Fm−1(xi)+γ)=∑xi∈Rij(Fm−1(xi)+γ)−yi=I(xi∈Rij)γ+∑xi∈RijFm−1(xi)−yiI(xi∈Rij)γjm+∑xi∈RijFm−1(xi)−yi=0γjm=∑xi∈Rijyi−Fm−1(xi)I(xi∈Rij) L(y, F_{m-1} (x) + \gamma) = \frac{1}{2} (y-(F_{m-1}(x)+\gamma))^2 \\ \frac{\partial}{\partial \gamma} L = (F_{m-1}(x)+\gamma) - y \\ \frac{\partial}{\partial \gamma} \sum_{x_i \in R_{ij}} L(y_i, F_{m-1} (x_i) + \gamma) = \sum_{x_i \in R_{ij}} (F_{m-1}(x_i)+\gamma) - y_i \\ = I(x_i \in R_{ij})\gamma + \sum_{x_i \in R_{ij}} F_{m-1}(x_i) - y_i \\ I(x_i \in R_{ij}) \gamma_{jm} + \sum_{x_i \in R_{ij}} F_{m-1}(x_i) - y_i = 0 \\ \gamma_{jm} = \frac{\sum_{x_i \in R_{ij}} y_i - F_{m-1}(x_i)}{I(x_i \in R_{ij}) }

上記からわかるように、γjm\gamma_{jm} は常に RjmR_{jm} における平均残差となります。これらの γ\gamma を各木のリーフノードの出力値として使用します。次に、 新しい木の結果に学習率 ν\nu を掛けて新しい Fm(x)F_m(x) を作成します。

したがって、手順を次のように書き直すことができます。

入力: データ (xi,yi)I=1n{(x_i, y_i)}_{I=1}^{n} と微分可能なコスト関数（二乗残差の半分）。

ステップ1: 予測変数 yy の平均値として、定数値 F0(x)F_0(x) でモデルを初期化する。

ステップ2: m=1m=1 から MM まで繰り返す:

  • 疑似残差 rim=yi−Fm−1(x)r_{im} = y_i - F_{m-1}(x) を i=1,…,ni = 1,…,n に対して計算する。
  • rimr_{im} の値に適合する木を構築し、端末領域 RjmR_{jm} を j=1,…,Jmj = 1, … , J_m に対して作成する。
  • j=1,…,Jmj = 1, … , J_m に対して、RjmR_{jm} の平均残差として γjm\gamma_{jm} を計算する。
  • Fm(x)=Fm−1(x)+ν∑j=1JmγmI(x∈Rjm)F_m (x) = F_{m-1} (x) + \nu \sum_{j=1}^{J_m} \gamma_m I(x \in R_{jm}) を更新する。

バイナリー分類

バイナリー分類における勾配ブーストは回帰の場合と似ていますが、いくつか重要な違いがあります。ロジスティック回帰と同様に、 木は負の対数尤度を最小化するようにロジットに適合します。

p=eF(x)1+eF(x)F(x)=log(p1−p)L(yi,F(x))=−(yilog(p)+(1−yi)log(1−p)) p = \frac{e^{F(x)}}{1+e^{F(x)}} \\ F(x) = log(\frac{p}{1-p}) \\ L(y_i, F(x)) = - (y_i log(p) + (1-y_i) log(1-p))

負の対数尤度を pp の代わりに予測されたロジット F(x)F(x) に関して表しましょう。

L(yi,F(x))=−(yilog(p)+(1−yi)log(1−p))=−(yi(log(p)−log(1−p))+log(1−p))=−(yiF(x)+log(1−p))=−(yiF(x)+log(1−eF(x)1+eF(x)))=−(yiF(x)+log(1−11+eF(x)))=−(yiF(x)−log(1+eF(x)))=−yiF(x)+log(1+eF(x)) L(y_i, F(x)) = - (y_i log(p) + (1-y_i) log(1-p)) \\ = - (y_i (log(p) - log(1-p)) + log(1-p)) \\ = - (y_i F(x) + log(1-p)) \\ = - (y_i F(x) + log(1-\frac{e^{F(x)}}{1+e^{F(x)}})) \\ = - (y_i F(x) + log(1-\frac{1}{1+e^{F(x)}})) \\ = - (y_i F(x) - log(1+e^{F(x)})) \\ = - y_i F(x) + log(1+e^{F(x)})

ステップ1として、コスト関数 L(yi,γ)L(y_i, \gamma) の γ\gamma に関する微分を求める必要があります。

L(yi,γ)=−yiγ+log(1+eγ)∂∂γL(yi,γ)=−yi+eγ1+eγ=−yi+pγ L(y_i, \gamma) = - y_i \gamma + log(1+e^{\gamma}) \\ \frac{\partial}{\partial \gamma} L(y_i, \gamma) = -y_i + \frac{e^{\gamma}}{1+e^{\gamma}} \\ = -y_i + p_{\gamma}

pγp_\gamma はロジット γ\gamma に対応する確率です。次に、コストの合計を最小化するために以下を求めます。

∂∂γ∑i=1nL(yi,γ)=∑i=1n−yi+pγ=npγ−∑i=1nyinpγ−∑i=1nyi=0pγ=∑i=1nyin \frac{\partial}{\partial \gamma} \sum_{i=1}^{n} L(y_i, \gamma) = \sum_{i=1}^{n} -y_i + p_{\gamma} \\ = n p_{\gamma} - \sum_{i=1}^{n} y_i \\ n p_{\gamma} - \sum_{i=1}^{n} y_i = 0 \\ p_{\gamma} = \frac{\sum_{i=1}^{n} y_i}{n}

この結果から、コストの合計を最小化する確率はクラスが1である確率によって与えられることが分かります。 次に、ロジット関数を使用して pγp_{\gamma} から γ\gamma を計算しましょう。

γ=log(pγ1−pγ)=log(∑i=1nyin1−∑i=1nyin)=log(∑i=1nyi∑i=1n1−yi) \gamma = log(\frac{p_{\gamma}}{1-p_{\gamma}}) \\ = log(\frac{\frac{\sum_{i=1}^{n} y_i}{n}}{1-\frac{\sum_{i=1}^{n} y_i}{n}}) \\ = log(\frac{\sum_{i=1}^{n} y_i}{\sum_{i=1}^{n} 1-y_i})

上記の結果から、クラスが1であるロジットは、負の対数尤度コストを最小化する初期のロジット予測 F0(x)F_0(x) または γ\gamma と同じであることが分かります。 次に、ステップ2では rjmr_{jm} を計算する必要があります。

rim=−[∂L(yi,F(xi))∂F(xi)]F(x)=Fm−1(x)=−(−yi+eFm−1(xi)1+eFm−1(xi))=yi−pFm−1(xi) r_{im} = -[\frac{\partial L(y_i, F(x_i))}{\partial F(x_i)}]_{F(x) = F_{m-1}(x)} \\ = -(-y_i + \frac{e^{F_{m-1}(x_i)}}{1+e^{F_{m-1}(x_i)}}) \\ = y_i - p_{F_{m-1}(x_i)}

ここで、クラスと前回の予測からの予測確率との間の疑似残差に対して木を適合させる必要があります。次に、予測を更新するために出力値 γ\gammas を計算します。出力値は、コストの合計を最小化する値である必要があります。

γjm=arg min⁡γ∑xi∈RijL(yi,Fm−1(xi)+γ)=arg min⁡γ∑xi∈Rij−yi(Fm−1(xi)+γ)+log(1+eFm−1(xi)+γ) \gamma_{jm} = \argmin_{\gamma} \sum_{x_i \in R_{ij}} L(y_i, F_{m-1} (x_i) + \gamma) \\ = \argmin_{\gamma} \sum_{x_i \in R_{ij}} - y_i (F_{m-1} (x_i) + \gamma) + log(1+e^{F_{m-1} (x_i) + \gamma})

しかし、上記の式の微分を取り γ\gamma を解くことは非常に難しいため、テイラー級数の二次多項式でその微分を近似します。（テイラー級数による近似が うまくいく理由はこの記事の範囲を超えています。ここでは、近似がインターセプト、一次導関数、二次導関数（値、傾き、凹凸）を同じにすることで機能するこ とだけ知っていれば十分です。より詳しく学びたい場合は、3Blue1Brownの Taylor series | Chapter 11, Essence of calculusを参照してください。）

L(yi,Fm−1(xi)+γ)≈L(yi,Fm−1(xi))+ddFL(yi,Fm−1(xi))γ+12ddF2L(yi,Fm−1(xi))γ2ddγL(yi,Fm−1(xi)+γ)≈ddFL(yi,Fm−1(xi))+ddF2L(yi,Fm−1(xi))γ L(y_i, F_{m-1} (x_i) + \gamma) \approx L(y_i, F_{m-1} (x_i)) + \frac{d}{dF} L(y_i, F_{m-1} (x_i))\gamma + \frac{1}{2} \frac{d}{dF^2} L(y_i, F_{m-1} (x_i))\gamma^2 \\ \frac{d}{d\gamma} L(y_i, F_{m-1} (x_i) + \gamma) \approx \frac{d}{dF} L(y_i, F_{m-1} (x_i)) + \frac{d}{dF^2} L(y_i, F_{m-1} (x_i))\gamma \\

上記の式をゼロに設定し、γ\gamma を解きます。

ddFL(yi,Fm−1(xi))+ddF2L(yi,Fm−1(xi))γ=0γ=−ddFL(yi,Fm−1(xi))ddF2L(yi,Fm−1(xi))=yi−pFm−1(xi)ddF(−yi+eFm−1(xi)1+eFm−1(xi)) \frac{d}{dF} L(y_i, F_{m-1} (x_i)) + \frac{d}{dF^2} L(y_i, F_{m-1} (x_i))\gamma = 0 \\ \gamma = \frac{-\frac{d}{dF} L(y_i, F_{m-1} (x_i))}{\frac{d}{dF^2} L(y_i, F_{m-1} (x_i))} \\ = \frac{y_i - p_{F_{m-1}(x_i)}}{\frac{d}{dF} (-y_i + \frac{e^{F_{m-1}(x_i)}}{1+e^{F_{m-1}(x_i)}})}

ここで、上の分子が計算した前回の決定木の擬似残差と同じであることが分かります。あとは分母を計算すれば良いです。

ddF(−yi+eFm−1(xi)1+eFm−1(xi))=ddF(−yi+eFm−1(xi)(1+eFm−1(xi))−1)=−(1+eFm−1(xi))−2e2Fm−1(xi)+(1+eFm−1(xi))−1eFm−1(xi)=−e2Fm−1(xi)(1+eFm−1(xi))2+eFm−1(xi)(1+eFm−1(xi))=−e2Fm−1(xi)(1+eFm−1(xi))2+eFm−1(xi)+e2Fm−1(xi)(1+eFm−1(xi))2=eFm−1(xi)(1+eFm−1(xi))2=eFm−1(xi)(1+eFm−1(xi))1(1+eFm−1(xi))=eFm−1(xi)(1+eFm−1(xi))(1−eFm−1(xi)(1+eFm−1(xi)))=pFm−1(xi)(1−pFm−1(xi)) \frac{d}{dF} (-y_i + \frac{e^{F_{m-1}(x_i)}}{1+e^{F_{m-1}(x_i)}}) = \frac{d}{dF} (-y_i + e^{F_{m-1}(x_i)}(1+e^{F_{m-1}(x_i)})^{-1}) \\ = -(1+e^{F_{m-1}(x_i)})^{-2}e^{2F_{m-1}(x_i)}+(1+e^{F_{m-1}(x_i)})^{-1}e^{F_{m-1}(x_i)} \\ = -\frac{e^{2F_{m-1}(x_i)}}{(1+e^{F_{m-1}(x_i)})^{2}}+\frac{e^{F_{m-1}(x_i)}}{(1+e^{F_{m-1}(x_i)})} \\ = -\frac{e^{2F_{m-1}(x_i)}}{(1+e^{F_{m-1}(x_i)})^{2}}+\frac{e^{F_{m-1}(x_i)}+e^{2F_{m-1}(x_i)}}{(1+e^{F_{m-1}(x_i)})^{2}} \\ = \frac{e^{F_{m-1}(x_i)}}{(1+e^{F_{m-1}(x_i)})^{2}} \\ = \frac{e^{F_{m-1}(x_i)}}{(1+e^{F_{m-1}(x_i)})} \frac{1}{(1+e^{F_{m-1}(x_i)})} \\ = \frac{e^{F_{m-1}(x_i)}}{(1+e^{F_{m-1}(x_i)})} (1 - \frac{e^{F_{m-1}(x_i)}}{(1+e^{F_{m-1}(x_i)})}) \\ = p_{F_{m-1}(x_i)} (1-p_{F_{m-1}(x_i)})

分母は、pFm−1(xi)p_{F_{m-1}(x_i)} と 1−pFm−1(xi)1 - p_{F_{m-1}(x_i)} の積となります。したがって、領域RjmR_{jm} においてコストの合計を最小化する γjm\gamma_{jm} の値は次のようになります:

γjm=∑yi−pFm−1(xi)∑pFm−1(xi)(1−pFm−1(xi)) \gamma_{jm} = \frac{\sum y_i - p_{F_{m-1}(x_i)}}{\sum p_{F_{m-1}(x_i)} (1-p_{F_{m-1}(x_i)})}

この値を各領域で計算し、学習率 ν\nu で Fm(x)F_m(x) を更新します。これらに基づいて、ステップは次のように書き換えることができます。

入力: データ (xi,yi)i=1n{(x_i, y_i)}_{i=1}^{n} と微分可能なコスト関数、負の対数尤度。

ステップ 1: モデルを定数値 F0(x)F_0(x) をyy のロジットとして初期化する。

ステップ 2: m=1m=1 から MM まで繰り返す:

  • 疑似残差 rim=yi−pFm−1(x)r_{im} = y_i - p_{F_{m-1}(x)} を i=1,…,ni = 1,…,n に対して計算する。
  • rimr_{im} の値に木を適合させ、j=1,…,Jmj = 1, … , J_m に対して端末領域 RjmR_{jm} を作成する。
  • j=1,…,Jmj = 1, … , J_m に対して、RjmR_{jm} 内の残差の合計を pFm−1(x)(1−pFm−1(x))p_{F_{m-1}(x)}(1-p_{F_{m-1}(x)}) の合計で割った値 γjm\gamma_{jm} を計算する。
  • Fm(x)=Fm−1(x)+ν∑j=1JmγjmI(x∈Rjm)F_m (x) = F_{m-1} (x) + \nu \sum_{j=1}^{J_m} \gamma_{jm} I(x \in R_{jm}) を更新する。

多クラス分類における勾配ブーストも同じように行われます。ロジスティック関数の代わりにソフトマックス関数を使用して確率をロジットにマッピングし、 同じ方法でロジットに木を適合させます。しかし記事がかなり長くなっているので、多クラス分類の詳細とコード実装については省略します。上記を実装するには、 sklearn ライブラリの GradientBoostRegressor と GradientBoostClassifier を使用すると便利です。

バイナリー分類における疑似残差

バイナリー分類において疑似残差を計算したとき、残差を以前の予測確率とその補数の積で割った値が得られました。数学的には理にかなっていますが、これが疑似残差の値にどのように影響するかを理解するために、 関数 p(1−p)p(1-p) をプロットしてみましょう。

上記のグラフからわかるように、関数はp=0p=0 または 11 のときに0になり、p=0.5p=0.5 のときに最大点を持つ美しい曲線であるということが分かります。 p(1−p)p(1-p)は疑似残差の分母にあたるため、疑似残差は以前の予測が0または1に近いときに高くなり、以前の予測が0.5に近いときには比較的小さくなることが 伺えます。これは、予測が極端な値にあり、間違えているときは迅速に調整し、現在の予測が0.5に近いときには過学習を避けるために慎重になることを促進しており、 モデルがデータから学習する際に望ましい動作と一致しています。このようにp(1−p)p(1-p)を用いることは学習に役に立つ大事な工程である訳です。

リソース

TOPに戻る↑