python: loggingの出力値を文字列として取得したい
python の logging で出力した info とかを文字列として取得したい! という場合の方法についてのメモです。
下準備
- 単純な logger と stream handler (コンソールへの出力のハンドラ) を用意します。
- 詳しくは https://docs.python.org/ja/3/howto/logging.html#logging-advanced-tutorial など参考にしてみてください。
今回は logger / hander 両方に INFO level をつけましたので info よりも重要度が高いものだけ console に output されるようになっています。
from logging import getLogger, StreamHandler, Formatter handler = StreamHandler() handler.setLevel('INFO') logger = getLogger('nyk.510') logger.setLevel('INFO') logger.addHandler(handler)
https://docs.python.org/ja/3/howto/logging.html#loggers
組み込みの深刻度の中では DEBUG が一番低く、 CRITICAL が一番高くなります。たとえば、深刻度が INFO と設定されたロガーは INFO, WARNING, ERROR, CRITICAL のメッセージしか扱わず、 DEBUG メッセージは無視します。
なるほど。というわけで、一旦試してみましょう。
logger.warning('warn') logger.info('foo') logger.debug('debug') # debug はでないよ warn foo
確かに debug は出ないようになっていますね。
loggingの出力値を文字列として取得
さて本題の logging の出力をテキストとして取得する、です。これは要するに上記の例で言うと warn / foo / debug みたいな文字列を取得したい、ということです。 結論をいうと StringIO を stream にもつような handler を作成して logger に付与すればOKです。 テキストとして取得するっていうのはだいたい log をどこかに保存したいとかいう気持ちがあると思いますので、ちょっとおしゃれな formatter にして時間等も取得できるようにしています。
log_capture_io = io.StringIO() stream_handler = StreamHandler(stream=log_capture_io) # オシャに formatting formatter = Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') stream_handler.setFormatter(formatter) stream_handler.setLevel('INFO') logger.addHandler(stream_handler)
この状態で logger に先と同じように log を記録します。
logger.warning('warn') logger.info('foo') logger.debug('debug') # debug は console にでないよ
log の取得
作成した StreamIO から getvalue
すればOKです。
- 通常のコンソールアウトプットは単に文字列だったが,
formatter
をリッチにしているので取得される文字列には何時 log が記録されたかなどの情報も入っている - コンソールの方も普通の
StreamHandler
に formatter を設定すれば時間も表示できる.
s = log_capture_io.getvalue() s.splitlines() ['2020-08-01 07:59:41,132 - nyk.510 - WARNING - warn', '2020-08-01 07:59:41,135 - nyk.510 - INFO - foo']
後片付け
ずっと handler が付いていると記録され続けるので、いらなくなったら消しましょう
- io の close
- handler のひも付けを logger から削除
removeHandler
# 終わったら消しましょう log_capture_io.close() logger.removeHandler(stream_handler) log_capture_io.closed # True
close してしまうと value はもう取れませんので注意
log_capture_io.getvalue() --------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-10-b5f4bf9c6e8d> in <module> ----> 1 log_capture_io.getvalue() ValueError: I/O operation on closed file
上記コードは gist にもありますので参考にしてください ;). Logging Recording · GitHub
機械学習なdockerfileを書くときに気をつけとくと良いこと
みなさん機械学習系の環境構築はどうやってますか? 僕は最近は Docker を使った管理を行っています。
特に師匠も居なかったので、ぐぐったり人のイメージを見たり手探りで docker をつかいつかいしている中で、最初からやっとけばよかったなーということがいくつかあるのでメモとして残しておきます。
大きく2つです。
- キャッシュは消す
- テストを書く
キャッシュは消す
ライブラリをいろいろと install すると大抵の場合ダウンロードしたファイルを保存されている場合が多いです。何かのタイミングで再びそのライブラリをインストールする際にはダウンロードしたファイルを使って、素早くインストールすることができます (この仕組みがキャッシュです)。
キャッシュがあると容量が重くなるという欠点があります。重たいイメージは pull に単に時間がかかりますから、システムとしてデプロイする時にトラフィックが多く発生してしまう + 更新に時間がかかるといったデメリットがあります。
ですのでライブラリを入れたあとにはキャッシュを削除すると吉です。
機械学習系だと基本的に ubuntu を base image として使うことが多いと思いますので、以下ではすべて ubuntu を使用していることを想定しています。 また、機械学習に必要なライブラリを入れるために使うコマンドは基本的に
- apt-get (or apt 絶対使う)
- pip (基本使う. conda only なら別.)
- conda (場合によって使う)
3つだと思いますので、それぞれ述べていきます。
apt-get の場合
ubuntu では apt-get
(もしくは単に apt
) でライブラリを入れますが、コマンドとしてキャッシュを消す機能がついていますので、それを使うのが便利です。
RUN apt-get autoremove -y &&\ apt-get clean &&\ rm -rf /usr/local/src/*
順番に
autoremove
で必要がないライブラリを削除してclean
でキャッシュされているすべてのライブラリを消して- インストールに使ったソースコードがある場合消す
ということをやっています。
pip の場合
pip には残念ながらキャッシュディレクトリを消すコマンドが用意されていません。一番簡単なのは pip install する時には cache を作らないようにする方法です。pip>=6 であれば no-cache-dir option をつけるとキャッシュせずに install が行われます。
pip install --no-cache-dir scipy
pip を最新版にしてからインストールすることが多いと思いますのでワンライナーで以下のように記述するといい感じです。
RUN pip install -U pip &&\ pip install --no-cache-dir scipy
インストールするライブラリはたくさんでてきくるので pip freeze で作れる requirements.txt
で管理する場合も多いと思います。その場合にも単に --no-cache-dir
をつけるだけです。
COPY /path/to/requirements.txt requirements.txt # あらかじめ requirements.txt を ADD or COPY などで # docker 内部に送っておくことが必要 RUN pip install -U pip &&\ pip install --no-cache-dir -r requrements.txt
もしくはインストールしたあとに cache directory を消すという方法もあります。Ubuntu なら場所は ~/.cache/pip
ですので全部ふっ飛ばせばよいです。
RUN pip install -r requirements.txt &&\ rm -rf ~/.cache/pip
Conda の場合
conda にはキャッシュ削除の clean
が用意されています。https://docs.conda.io/projects/conda/en/latest/commands/clean.html
いろいろインストールした最後に呼び出しましょう。
RUN conda install -y \ numpy \ scipy &&\ conda clean -i -t -y
その他
バイナリからインストールしたりする場合には使わないソースコードは消すようにしましょう。例えば noto font をインストールしたあとに元の .zip
を消すとかそういうことです。
テストを書く
スゴーク簡単でも良いので、いつも使うライブラリを import + 定型的な処理を行なうコードを .py
で書いておいて build 時にテストすると精神的に良いです。
例えばですが Mecab を使って文字列パースするような処理を書いていると、pip でインストールされていることはもちろんですが apt の方でもちゃんと mecab が入っていることが確認できます。
def test_mecab(): import MeCab tokenizer = MeCab.Tagger('-Ochasen') tokenizer.parse('おはようせかい')
具体的にはこれを build したあとのイメージに入れて pytest
を実行してあげます。
docker run -v ./test.py:/analysis/ my-image pytest ./
僕の場合はイメージの build は gitlab CI 上で行っているので CI で build したあとに registry に登録する前にテストを動かすような構成にしています。(gitlab ci の yaml 形式ですが script にあるコマンドが順番に実行されるということだけわかれば雰囲気はわかるかと思います。)
.build-docker: stage: build before_script: - docker login -u gitlab-ci-token -p $CI_JOB_TOKEN $CI_REGISTRY script: - echo ${REGISTORY_URL} - export TAG=${REGISTORY_URL}/${IMAGE_NAME} - docker build -t ${TAG}:latest -f ${DOCKERFILE} ${CONTEXT} - docker tag ${TAG}:latest ${TAG}:${CI_COMMIT_SHA} - docker run -v ${PWD}:/analysis ${TAG}:latest pytest tests/${IMAGE_NAME} # < ここがそう - docker push ${TAG}:latest - docker push ${TAG}:${CI_COMMIT_SHA} variables: IMAGE_NAME: DOCKERFILE: # default values. if you use own values, override in variables scope. CONTEXT: . # ex.) registry.gitlab.com/atma_inc/analysis-template REGISTORY_URL: ${CI_REGISTRY}/${CI_PROJECT_NAMESPACE}/${CI_PROJECT_NAME}
こうしておくとテストに失敗すると登録自体が行われないので、登録されている == テストは通っているという安心感を得られます😃
例
自分がいつも使っているイメージ: https://gitlab.com/nyker510/analysis-template/-/blob/master/docker/cpu.Dockerfile
このリポジトリ analysis-template では dockerfile のコード管理 + master merge と tag を切った時に CI でビルド & registry に push ということをやっています。 内容に関しては昔詳しく書いたことがあるので、こちらも参考にしてみてください。
ドキュメントもあるよ!
こうしたほうがいい! 自分はこういうことをやっている とかあればコメント貰えると嬉しいですmm
JWT でのログインができるかどうかを確認するテスト
Django Rest Framework + allauth 使用時に JWT での Login ができているかどうかを念の為確認するテスト。
import pytest from django.urls import reverse from atma.shops.models import Shop from .factories import UserFactory @pytest.mark.django_db def test_authorized_by_jwt(client): """Json Web Token で login できる""" user = UserFactory() # type: User user.set_password('hogehoge') user.save() url = reverse('accounts_me') res = client.get(url) assert res.status_code == 401, res.json() auth_url = reverse('rest_login') res = client.post(auth_url, data={'email': user.email, 'password': 'hogehoge'}) assert res.status_code == 200, res.json() jwt_token = res.json().get('accessToken', None) assert jwt_token is not None header = {'HTTP_AUTHORIZATION': 'JWT {}'.format(jwt_token)} res = client.get(url, **header) assert res.status_code == 200, res.json()
- パスワードの設定は
user.set_password
で行なう - header を kwrgs で渡す
"HTTP_AUTHORIZATION"
を key にする
の3点がポイント
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 と事前に割り算をして計算量を減らせますが定式化通りやっています