NGBoostを読んで、実装する。
不確実性を考慮した予測が可能と噂の NGBoost の論文を読んでみたので、全体のながれをまとめて見ました。加えて自分でも NGBoost を実装して、その結果を載せています。
元の論文 NGBoost: Natural Gradient Boosting for Probabilistic Prediction はこちら https://arxiv.org/abs/1910.03225。
Introduction
一般的な教師あり学習を考えます。このとき予測モデルは入力データ $X$ に対して予測値 $y$ を出力するように学習していきますが、たいていのモデルではひとつのデータに対しては予測値はひとつしか得ることができません。
例えばウェブ上の行動履歴から、ユーザーの年齢を予測してください、という問題があったとすると、ユーザーの期待される年齢そのものを返すようなモデルがそれに当たります。
例えばテーブルデータでデファクトスタンダードとなっている XGBoost や LightGBM のような勾配ブースティングもそうですし一般的な Deep Learning もそのような設計です。
しかし実際には予測には確信度のようなものが存在しているはずです。ユーザの年齢だけ答えてくださいと言われると数字を出力するしかありませんが、実際には「このユーザーは絶対に若い」とわかるようなデータだったり「年取っていると思うけれど自信ない」というようなデータも存在しているでしょう。 このようなデータによる不確実性も同時にモデリングしましょうというのがモチベーションです。
不確実性を考慮した予測
一般的なモデリングでは確率分布それ自体ではなくて、代表するようなパラメータを予測します。 例えば RMSE を目的関数としている場合では平均値 $\mu$ を予測する関数を作ることが目的です。
不確実性を考慮する予測とは入力 $X \in \mathbb{R}^k$ にたいして対応する予測値 $y$ がどのような分布になっているかを関数として予測する必要があります。すなわち
$$ P(y | X) $$
で表される確率分布 $P$ を推論する作業です。 このままだと分布 $P$ は無限の形状をとれて埒が明かないので特定のパラメータ $\theta \in \mathbb{R}^p$で表現できるような分布 $P_{\theta}$ に絞って考えていくことにします。 例えば $P$ に正規分布を仮定して正規分布は平均値 $\mu$ と分散 $\sigma$ の2つで説明できる分布なので $\theta = (\mu, \sigma)$ とおく、というような操作に相当します。
この仮定によって分布自体を考える必要はなくなり「パラメータ $\theta$ をいかに予測するか」という問題に変えることができます。 以下では $X$ から分布のパラメータ $\theta$ を出力するような関数 $F$ を推論することを考えていきます。
$\theta$ の推論のための評価基準
$\theta$ を推論するためには $\theta$ のクオリティを評価できる必要があります。すなわち $\theta$ とそのときの実現値 $y$ を入力として数値を返すような関数 $S$ を決めなくてはいけません(以下では $S$ が小さいほどよい関数だとします)。 最も良い $\theta$ とは何かを考えてみると、これは間違いなく $P_{\theta}$ が真の分布 $Q$ に一致するような場合ですね。理想的には。 これを言い換えると $\theta$ を評価する関数 $S$ は $P$ が真の分布 $Q$ になったとき、もっとも小さくなっている必要があります。これを数学的に記述すると
$$ \mathbb{E}_{y \sim Q} [S (Q, y)] \le \mathbb{E}_{y \sim Q} [S(P, y)] \ \ \forall P, Q $$
となります。ここで $y \sim Q$ は確率分布 $Q$ からのサンプル $y$ での期待値を表します。すなわち
$$ \mathbb{E}_{y \sim Q} [S (Q, y)] = \int S(Q, y) Q(y) dy $$
です。このような性質を満たしている指標はいくつかあり最も有名なものが負の対数尤度関数 $\mathcal{L}$ (MLE) です。 $$ L(P, y) = - \ln P(y) $$
今回は $P$ がパラメータ $\theta$ で表わされるものと仮定したので実際には $\theta$ の関数となりますので以下のように記述できます。
$$ \mathcal{L}(\theta, y) = - \ln P_{\theta}(y) $$
この他にも CRPS (以下では $\mathcal{C}$ と表します) という指標もありこちらは以下のように記述される値です。 (一般に MLE よりもロバストであるとのことです)
$$ \mathcal{C} = \int_{- \infty}^{y} F_{\theta} (z)^2 dz + \int_{y}^{\infty} (1 - F_{\theta}(z))^2 dz $$
ここで $F_{\theta}$ は分布 $P_{\theta}$ の累積分布関数です。 (予測関数 $F$ と同じ記号ですけれど違います。念の為。)
$\theta$ の最適化
これらのスコア $S$ があればパラメータが最適化できそうです。ここでは一般的な勾配法を使った場合を考えてみます。 $S$ に対しての $\theta$ に対する最急降下方向 $\nabla S$ は以下のように記述できます。
$$ \nabla S (\theta, y) \sim \lim_{\epsilon \to \infty} \arg \min_{d: | d | = \epsilon} S (\theta + d, y) $$
一見これで問題なさそうに見えますが $S$ はもともと分布 $P$ に対して計算されるものですからパラメータ空間 $|d|$ でのユークリッド距離を用いた範囲内で min をとると $S$ による歪みを考慮することができません。
ここで $S$ が満たす必要がある式を思い出して右辺から左辺を引くと
$$ \mathbb{E}_{y \sim Q} [S (P, y)] - \mathbb{E}_{y \sim Q} [S (P, y)] \ge 0 $$
となります。これは非負の値ですから $Q$ から $P$ がどのぐらい離れているかの度合い (距離) とみなすことができます。 例えば負の対数尤度であれば KL Divergence に相当し, CRPS であれば $L^2$ Divegence となります。
これらの Divergence は Q / P がどのようなパラメータを持っているかに不偏です (Divergenceの式自体には $\theta$ が出現しないことに注意してください)。したがってこの Divergence で表されるような空間で議論することができれば $\theta$ ではなく分布 $P$ の距離の意味で最も小さくなる方向を見つけることができるため、パラメータ $\theta$ のとり方に依存しない方向を見つけることができます。
すなわち最急降下方向を求める際の移動量を測る空間であったユークリッド距離 $| d |$ を、もともとのパラメータ $\theta$ と移動した点 $\theta +d$ での確率分布 $P_{\theta}, P_{\theta + d}$ との距離を測る Divergence $D_S$ で置き換えて
$$ \nabla \hat{S} (\theta, y) \sim \lim_{\epsilon \to \infty} \operatorname{argmin}_{d: D_S (P_{\theta} | P_{\theta + d}) = \epsilon} S (\theta + d, y) $$
とします。こうして求められる勾配はパラメータに不偏です。これを自然勾配 (Natural Gradient) と呼びます。 これを解くとスコア $S$ を用いたときの点 $\theta$ でのリーマン計量行列 $\mathcal{I}_{S}(\theta)$ としたとき
$$ \nabla \hat{S} (\theta, y) \propto \mathcal{I}_{S} (\theta)^{-1} \nabla S(\theta, y) $$
と記述できます。例えば $S = \mathcal{L}$ (負の対数尤度) とするとこれはフィッシャー情報量行列です。
$$ \mathcal{I}_{\mathcal{L}} (\theta) = \mathbb{E}_{y \sim P_{\theta}} [ \nabla_\theta \mathcal{L}(\theta, y) \nabla_\theta \mathcal{L}(\theta, y) ^T] $$
自然勾配を使うとパラメータ変換に堅牢に最適化を行えます。 Figure3 では正規分布に対して平均値 $\mu$ と分散を対数変換した $\log \sigma$ 場合の勾配の比較があります。これを見ると通常の勾配 (右側) よりも自然勾配 (左側) のほうが安定して最適解への方向をしめしていることがわかります。
$F$ の最適化
では次に $\theta$ を求めるような関数 $F$ を作っていくのですが NGBoost ではこれを Natural Gradient Boosting で学習します。 Natural Gradient Boosting は Gradient Boosting をちょっといじったものなので、まずは Gradient Boosting から。
Gradient Boosting
Gradient Boosting は求めるモデル $F$ が複数のモデル $f_i$ の足し算 $F_{m} = \sum_{i=1}^m f_i$ で表される加法的モデルと仮定します。 これを $m$ 個同時に最適化することは難しいので $m$ 番目のステップではすでに作った $i = { 1, \cdots, m-1 }$ のモデルは固定して、次に加えるモデルをスコア $S (\theta, y)$ の $\theta$ に対する勾配と一致するように学習します。すなわち $i$ 番目のデータに対する予測結果が
$$ f_{m} (x_i) = - \frac{S ( F_{m-1}(x_i), y_i)}{\partial F_{m-1}(x_i)} $$
となるよう $f_m$ を求めます。 ($m-1$番目までの状態で得られる $\theta_{i}^{m-1}$ は $F_{m-1}$による予測値 $F_{m-1}(x_i)$ に一致していることに注意してください.) 具体的には $g_m^i = - {S ( F_{m-1}(x_i), y_i)} / {\partial F_{m-1}(x_i)}$ とおいて二乗誤差
$$ \sum_{i=1}^N \left( g_m^i - f_{m}(x_i) \right)^2 $$
を最小化するように $f_m$ を決定します。このとき $f$ に使うモデルを弱学習器と呼び、一般には深さを制限した決定木などが使われます。
Natural Gradient Boosting
Gradient Boosting 目標としていた再急降下方向を自然勾配に置き換えてしまおう、というのが Natural Gradient Boosting です。すなわち $f_m$ を近づける勾配を自然勾配に置き換えます。要するにリーマン計量行列の逆行列で補正した
$$ \hat{g}_m^i = - \mathcal{I}_{S}(F_{m-1}(x_i))^{-1} \frac{S ( F_{m-1}(x_i), y_i)}{\partial F_{m-1}(x_i)} $$
を使って
$$ \sum_{i=1}^N \left( \hat{g}_m^i - f_{m}(x_i) \right)^2 $$
を最小化するように $f_m$ を決定します。
アルゴリズム
それでは実際に Natural Gradient Boosting のアルゴリズムがどのように動いていくかをまとめてみます。 上記でのべた仮定をふまえるとアルゴリズムを動かすためには初めに
- 予測対象とする分布 $P_{\theta}$
- 各 Boosting で学習するためのモデル $f$
- $\theta$ の良さを評価するための指標 $S$
の3つを決めておく必要があります。 これらを使って、以下の 1 ~ 4 を $M$ 回繰り返すことで学習器とステップサイズ $\{ \rho^m, f^m \}_{m=1}^M$ の集合を作っていきます。 (以下の $m$ は現在の boosting の iteration 数です)
1. 探索方向 (Natural Gradient) の計算
今のパラメータ $\theta^{(m-1)}$ を用いて $i$ 番目のデータに対する自然勾配 $g_i^m$ を計算します。
2. $f_m$ の学習
次にデータ $x_i$ から勾配 $g_i^m$ を予測するモデルを作成します
$$ f^m = \arg \min_{f} \sum_{i=1}^N \left( f(x_i) - g_i^m \right)^2 $$
このとき求めているのは予測値 $y$ に関係するものでなく、単にパラメータ $\theta$ の勾配です。したがって学習のためのモデルは Regression を行なうものなら基本的に何でも OK です。
3. ステップサイズ $\rho$ の決定
モデル $f^m$ が決定すれば各データ $x_i$ ごとにパラメータ $\theta_i$ がどちらに移動した良いかを $f(x_i)$ によって推論することができます。 (実際に次のステップ 4 ではこの $f(x_i)$ によって $\theta_i^m$ を更新します。) そのままこの更新方法を使っても良いですが、最適解からとても遠い (or 近い) ことを考慮して、勾配を定数倍して最もスコアが下がる点を探します。
$$ \rho^m = \arg \min_{\rho} \sum_{i=1}^N S\left(\theta^{(m-1)} - \rho f^{m} (x_i) , y_i \right) $$
実際には linear search によって探索し最適なステップサイズを $\rho^m$ とします。
4. パラメータ $\theta$ の更新
ステップサイズ $\rho$ とモデルによって予測された更新方向 $f(x_i)$ を使って実際にパラメータを更新します。このとき値をそのまま使わず定数 $\eta$ をかけて縮めます。
$$ \theta_{i}^{m} = \theta^{(m-1)} - \eta \cdot \rho^m f^m(x_i) $$
(Gradient Boosting などの文脈では learning rate と呼ばれることもあるパラメータです。一気に学習が進んでしまうのを防ぐのが狙いです。)
NGBoost の利点
論文中にあげられている NGBoost の利点は大きく以下の3つです。
1. いろんな分布で使える
予測対象としている分布の制限は単にパラメータで表現できるというものだけですから、いろんなケースで使用することができます。(たとえば負の二項分布だったり、ガンマ分布だったり)
2. 複数のパラメータの同時推定でも最適解を見つけやすい
通常の勾配を用いると初期値の時点でたまたまうまく予測できている lucky なデータがある場合それらのデータが支配的になってしまい、分散推定がそれらのたまたまうまく予測できたデータのみ進んでしまい、反対にたまたま初期値が大きく外れている unlucky なデータで学習が足りないという現象が発生します。 (Figure4) 自然勾配を使うとフィッシャー情報量行列などで勾配を補正するため、初期点にせずにまんべんなく学習をすすめることができます。
3. パラメータ変換に不偏
自然勾配の定義からもわかるようにパラメータ変換に依存せずに勾配を選ぶことができます。これによってユーザーは推定する分布のパラメータを好きに変数変換しても問題ありません。例えば正規分布を推定する分布とすると、予測するために平均値と分散が必要になりますが、分散は正である必要があるためログ変換したパラメータ $\log \theta$ を推論するほうが数値的にも安定して便利です。
このとき自然勾配を使っていれば問題なく最適化できますが普通のニュートン法だともはやこの問題は凸ではないため解くことが難しくなります。
実験と結果
table1 では不確実性を考慮できるいくつかのモデルと NGBoost の性能比較がされています。 評価指標は負の対数尤度 (NLL) です。 これを見ると負けているデータセットもあるものの概ね良いスコアを出せていることがわかります。
table2 では Natural Gradient を使わない通常の Gradient Boosting (Multiparameter Boosting) と、分散を固定して不確実性を考慮しない Boosting (Homoscedatic Boosting) との比較がされています。 これも評価指標は負の対数尤度です。 こちらを見ると NGBoost は Gradient Boosting の場合よりも性能が高いことがわかります(Natural Gradient にした成果)。 また分散を固定したモデルとも匹敵した性能がでていて分散を個別に推定できるという利点を考えると NGBoost を使って損はない、という結果になっています。
負の対数尤度では分布全体を考慮する NGBoost が有利でしょう、ということで点推定の結果を出したものが table3 です。こちらは評価指標が RMSE となっています。 NGBoost 以外のモデルは RMSE の点推定を目的関数に持っており NGBoost は違うにもかかわらず、他のモデルに匹敵した性能がでていることがわかります。
結論
- NGBoost は不確実性を考慮したモデリングを出来るように Gradient Boosting を拡張したモデルです。
- 安定した探索方向を作成するため Gradient Boosting の各ステップでの最急降下方向を Natural Gradient に変更しました。
- 基本的には Gradient Boosting と同じなのでモンテカルロ法やベイズの知識がなくとも簡単に不確実性を考慮した予測をすることが可能です。
- 点推定モデルとくらべて、点推定タスクに置いても遜色ない性能が見込まれるため「コスト無しで (for free)」不確実性を推論できます。
実装
ここまでは論文と調べた内容のまとめです。 実際に NGBoost を使う場合 python package が公式から提供されています https://github.com/stanfordmlgroup/ngboost 。 ライブラリ自体は pip から install できますのでさっくり試してみたい方はこちらを使うことをおすすめします。 User Guide もあり自分で Score $S$ や分布 $P$ を カスタムしたい場合のガイド もあります。
pip install --upgrade git+https://github.com/stanfordmlgroup/ngboost.git
以下は勉強兼ねて、僕が自分で実装してみたよという内容です。 コード自体はちょっと違いますが script にまとめたものは以下から参照ください。
自然勾配の計算
アルゴリズムの部分にもあったようにコアになるのは現状のパラメータから Natural Gradient を返す部分です。 対数尤度をスコアとしている場合勾配とフィッシャー情報量行列を計算する必要がありますが、上記の議論でもあったように、もとめるパラメータの変数をどのように定義するか(変数変換の定義方法)によってこれらの形式も変動します。
今回は正規分布を予測分布として平均値 m, 分散をログ変換した log_var を推論対象としました。
このとき勾配とフィッシャー情報量行列は以下のように計算できます。*1
var = np.exp(log_var) # log -> var に変換 gradients = [ (m - t) / var, 1 / 2 - ((t - m) ** 2) / (2 * var) ] fisher_matrix = [ [1 / var, 0], [0, 1 / 2] ]
弱学習器としては決定木 sklearn.tree.DecisionTreeRegressor を用います。
このときあまり精緻に木を作ってしまうとオーバーフィットのもとになるため max_depth / min_sample_leaf だけパラメータ指定をしています。
from sklearn.tree import DecisionTreeRegressor def build_weak_learner(): return DecisionTreeRegressor(max_depth=3, min_samples_leaf=4)
学習する
あとは Boosting 部分を実装すればOKです。 全体を記述すると以下のようになります。
import numpy as np from scipy.stats import norm from sklearn.tree import DecisionTreeRegressor def build_weak_learner(): return DecisionTreeRegressor(max_depth=3, min_samples_leaf=4) class LogVarianceNorm: def calculate_gradient(self, t, m, log_var): var = np.exp(log_var) gradients = [ (m - t) / var, 1 - ((t - m) ** 2) / var ] fisher_matrix = [ [1 / var, 0], [0, 1] ] return np.array(gradients), np.array(fisher_matrix) def score(self, t, m, log_var): var = np.exp(log_var) return - norm.logpdf(t, loc=m, scale=var ** .5) def true_function(X): return np.sin(3 * X) n_samples = 200 np.random.seed(71) X = np.random.uniform(-3, 1, n_samples) y = true_function(X) + np.random.normal(scale=abs(X), size=n_samples) log_var_norm = LogVarianceNorm() # hyper parameters learning_rate = 1e-2 n_iterations = 500 estimators = [] step_sizes = [] X = X.reshape(-1, 1) N = len(X) # initialize parameters by whole dataset mean and variance m, log_var = np.zeros(shape=(N,)), np.zeros(shape=(N,)) m += np.mean(y) log_var += np.log(np.var(y) ** .5) for iteration in range(n_iterations): # 1. compute natural gradient grads = [] for i in range(N): score, fisher = log_var_norm.calculate_gradient(y[i], m[i], log_var[i]) grad = np.linalg.solve(fisher, score) grads.append(grad) grads = np.array(grads) # 2. fit weak learner # mean estimator clf_mean = build_weak_learner() clf_mean.fit(X, y=grads[:, 0]) # log_var estimator clf_var = build_weak_learner() clf_var.fit(X, y=grads[:, 1]) directions = np.zeros(shape=(N, 2)) directions[:, 0] = clf_mean.predict(X) directions[:, 1] = clf_var.predict(X) estimators.append( [clf_mean, clf_var] ) # 3. linear search scores = [] stepsize_choices = np.linspace(.5, 2, 21) for stepsize in stepsize_choices: d = directions * stepsize score_i = log_sigma_norm.score(y, m - d[:, 0], log_var - d[:, 1]).mean() scores.append(score_i) best_idx = np.argmin(scores) rho_i = stepsize_choices[best_idx] stepsize_i = learning_rate * rho_i step_sizes.append(stepsize_i) # 4. update parameters grad_parameters = directions * stepsize_i m -= grad_parameters[:, 0] log_var -= grad_parameters[:, 1] if iteration % 50 == 0: print(f'[iter: {iteration}]\tscore: {scores[best_idx]:.4f}\tstepsize: {rho_i:.3f}')
推論する
推論の際にはすべての弱学習器の結果を足しあわせてパラメータを生成します。
xx = np.linspace(-4, 2, 501) # 予測対象のデータ xx_mean = np.zeros_like(xx) + np.mean(y) xx_log_var = np.zeros_like(xx) + np.log(np.var(y) ** .5) # すべての弱学習器の予測結果を重み付けして足しあわせ for i in range(len(estimators)): xx_mean -= step_sizes[i] * estimators[i][0].predict(xx.reshape(-1, 1)) xx_log_var -= step_sizes[i] * estimators[i][1].predict(xx.reshape(-1, 1)) # 分散に関してはログ変換しているので元に戻す xx_var = np.exp(xx_log_var)
推論結果をグラフにすると以下のようになります。 ノイズ項を入力 $X$ の絶対値に比例するように設定していますが、それに追随するような予測値が得られていることがわかります。
Boosting Iteration の進行と予測の変化をアニメーションにすると以下のような感じになります。初め平均値を合わせたあと、段々と分散が正しい値に落ち着いていっていることがわかります。
弱学習器に木構造を使っているので当たり前ですが、外挿される値は平均も分散も定数になるのが面白いなと思いました。Dropout MCMC の結果 と比べると大きく違う傾向ですので、タスクによってどのように差がでるのかなども見てみたいです。
上記コードですが、Kaggle Notebook: https://www.kaggle.com/nyk510/simple-ngboost-implementation にも載せてみました。 ご参考までに。
感想
- XGBoost のような GBDT を拡張するという話かと思っていたらちょっとちがった
- 点推定の性能評価で Gradient Boosting 比較しかなかったので GBDT と比較したときどの程度性能が落ちるのか見てみたい。
- ベースは Gradient Boosting なので GBDT のような葉の重みに対する正則化 (reg_alpha / reg_beta) が入っておらず weak learner が overfit しないようにうまく調整する必要がありそうだと感じた。実装していても
max_depthを指定しなかったりmin_samples_leafを大きくすると予測結果がかなりギザギザして木構造っぽい予測になりがちなので、データに合わせたチューニングが必要そう。 - table1~3にあるNavalってデータセットは簡単すぎんか? (全部予測誤差 0)。逆にどんなのか見てみたい気もする。
- linear search でステップ幅を探索する部分があるのが気になった。
- これらは各 boosting で選ばれたデータに対して学習したモデルの方向で探索するためたまたまいい方向が選ばれたり・逆があったりで大きく $\rho$ の値が変わる場合がありそう。
- これによって boosting 回数の合計を制御しても, colsample / subsample の度合いでステップ合計が変わるので例えば CV で n_estimator を決定して再学習させたとき期待した量だけパラメータが更新されない可能性があるような気がしている。 (GBDTの場合は
learning_rateのみなのでn_estimatorsさえ決めれば更新ステップ量は固定できる) そんな変なこと考えなくてもよい?
参考文献
- NGBoost: https://github.com/stanfordmlgroup/ngboost
- 公式の実装です。普通に使う分にはこちらを使いましょう。
- 統計的学習の基礎 (第十章)
- 通称カステラ本。久しぶりに開きましたが Boosting の章は得にわかりやすくてとても好きです。
- Gradient Boosting と XGBoost · ZABURO
- Pure な Boosting から Gradient Boosting ~ XGBoost までの流れを丁寧に説明されているブログです。
*1:ヘシアンが対角行列ですので, gradient と事前に割り算をして計算量を減らせますが定式化通りやっています
