機械学習な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 と事前に割り算をして計算量を減らせますが定式化通りやっています
python の DateTime・Timezone と Django での取扱い
django でアプリケーションを作っていて timezone 周りで困ったことがあったので、初めて django および python での時間の取り扱いについて真面目に調べてみた記録です。
はじめに・ われわれが使っている時計の時間とは何か
まず時計が示している時間とは何かということを考えてみましょう。
世界中に時計がひとつしかなければ(すべての人間がひとつの時計を正として生活していれば)特に問題はなくその時計の時刻自体を YYYY-m-dd HH:mm
などで保存したものが常に正しいので何も考えることはありません。 皆同じものを参照しているので不整合はありません。
しかし実際には世界には時差があります。時差があるというのはどういうことかというと時計を文字通り読み取った値が同じでも、物理的には違う時間なことがあるということです。
例えば東京の 4/23 18:00 とニューヨークの 4/23 5:00 は物理的には同じ時間を表していますが、時計は違う時刻を指しています。それは東京・ニューヨークが別の Timezone に属していているためです。
Timezone ごとに、標準時 (イギリスのあたり) の時計に対して何時間ずらすのかが決まっています。例えば東京では標準時から +9H ・ ニューヨークなら -4H ときまっているので、先の 18:00 と 5:00 が同じ物理時間を指していることが確認できます (差分が 9 + 4 = 13 時間であるからです)。
python で時刻を扱う
これらの問題は python で時刻を取り扱う場合でも同様なので、 単に数値としての時刻と Timezone の情報を持った時刻を区別して取り扱うように実装されています。 前者を naive time, 後者を aware time と呼びます。
naive time は単に datetime として作成したときの状態です。標準ライブラリの datetime で現在時刻をとったような状態では naive time です。
from datetime import datetime naive = datetime.now() # datetime.datetime(2020, 4, 23, 10, 15, 27, 886221)
aware な時刻かどうかは datetime object が持っている tz
によって判定されています。naive な datetime は tz を持っていません。
naive.tzname() is None # True
tzを設定するひとつの方法が naive.astimezone
です。この関数は tzinfo
instance を引数にとってその timezone の aware な時刻のオブジェクトを返します。
tzinfo の設定は pytz
を使うのが楽ちんです。例えば UTC 時刻として登録する場合は以下のようになります。
import pytz naive.astimezone(pytz.UTC) # datetime.datetime(2020, 4, 23, 10, 15, 27, 886221, tzinfo=<UTC>)
一度 aware になってしまうと naive な時刻とは引き算をすることができません。
くどいですが区別するのは Timezone 込の情報と Timezone 無しのものは比較することができないからです。 これら2つを同様に扱ってしまってナイーブに引き算などすると例えば東京とニューヨークを同様に扱って引き算できてしまうので、先の例だと時間差がゼロになってしまいます。
実際 naive time と aware time を引き算すると TypeError
が発生します。親切ですね。
naive - naive.astimezone(pytz.UTC) --------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-24-98bec8f1194b> in <module> ----> 1 naive - naive.astimezone(pytz.UTC) TypeError: can't subtract offset-naive and offset-aware datetimes
astimezone
で変換さた時刻は、元の naive な時刻を標準時で持つような tzinfo の時刻となります。
そのためここで tokyo や newyork を設定しても引き算すると時刻の差はゼロです。
tokyo_tz = pytz.timezone('Asia/Tokyo') nyc_tz = pytz.timezone('America/New_York') naive.astimezone(tokyo_tz) - naive.astimezone(nyc_tz) # datetime.timedelta(0)
もし「東京の時計で naive
な時刻」を持つような時間を作りたい場合には timezone.locale を使います。
例えば東京とニューヨークそれぞれで同じ時計が表示されている状態の差分を見る場合は以下のようになります。期待通り時差 +9H から -4H を引くので 13時間のズレになっていることがわかります。
# 東京とニューヨークで文字盤上同じ時刻を引き算すると時差と一致する diff = tokyo_tz.localize(naive) - nyc_tz.localize(naive) diff.total_seconds() / 60 ** 2 # -13.0
Django での時間の取り扱い
Django では USE_TZ=True
のとき時差を考慮した処理を行なうようになります。すなわち django が python object として扱う時刻はすべて aware な datetime になるということです。反対に USE_TZ=False
にしていると tz を設定しない naive な datetime を使います。この場合現在時刻を知りたければ datetime
を使うことができます。
from datetime import datetime now = datetime.now()
一方で USE_TZ=True
の場合は aware な方法で now を作成しなくてはいけません。このような場合のために django では django.utils.timezone
というモジュールがありその内部に now
が実装されています。
from django.utils import timezone aware_now = tiemzone.now()
この内部では datetime.now を呼び出したあとに UTC を tz に設定しています。 このように django では時差考慮するとなった時には基本的に UTC として取り扱います。これによってデータベースに保存する際にすべての時刻が UTC であることが担保されますので、引き算等の演算での不整合が発生することを防ぐことができます。
# django/utils/timezone.py # UTC time zone as a tzinfo instance. utc = pytz.utc # 中略 def now(): """ Return an aware or naive datetime.datetime, depending on settings.USE_TZ. """ if settings.USE_TZ: # timeit shows that datetime.now(tz=utc) is 24% slower return datetime.utcnow().replace(tzinfo=utc) else: return datetime.now()
webapplication ではリクエスト元のユーザーが違う timezone に属していることがあります。このような場合にはリクエストごとに Timezone を変更したくなりますが、この用途のために django には activate という関数があり、これを呼び出すとデフォルトの TimeZone を上書きすることができます。
from django.utils.timezone import activate activate('Your-Favorite-Timezone')
テンプレートエンジンや DRF のシリアライザなどユーザーに情報を返す際にはここで設定された TimeZone での locale datetime への変換がされるので、ユーザーの TimeZone の時計で見た時の時刻をユーザーごとに設定して返すことができます。
Example: time から datetime への変換
「TimeField を今日の時間とみなして、その時間と DateTimeField の値を比較したい」といいう要望があったとしましょう。 Time を date に設定する関数 combine を使うと以下のように書くことができるように思えますが USE_TZ
の場合エラーになります。それは to_today_datetime
が返している datetime
が naive だからです。
from django.utils import timezone def to_today_datetime(time, now): today = now.today() dt = timezone.datetime.combine(today, time) return dt time = # TimeField の値. Time 型 now = timezone.now() # Log という model があったとする created_at = Log.objects.first().created_at # created_at は DateTimeField created_at - to_today_datetime(time, now) # TypeError: can't subtract offset-naive and offset-aware datetimes
では utc
にすれば OK かというとそうでもありません。 なぜなら今日と認識しているのはリクエストを送ったユーザーですから、ユーザーにとっての時間として表現する必要があるからです。 今設定されている timezone は get_current_timezone
で取得することができます。 得られた tzinfo を使って localize することで時間をユーザーのtimezoneでtimeを時計で見たときの時刻に変更することができます。
def to_today_datetime(time, now): today = now.today() dt = timezone.datetime.combine(today, time) # これでは UTC の人にとって `dt` が時計に表示されているときの時刻になる。 # dt = dt.astimezone(timezone.utc) # 今のリクエスト context での timezone を取得して、その timezone で `dt` が表示された時の時間に直す tz = timezone.get_current_timezone() dt = tz.localize(dt) return dt