Ryzenマシン組んだのでXGBoostのGPU版やってみた

Ryzen 7 1700とGTX 1080 Tiでマシンを組んだので、動作確認がてらXGBoostGPU版を使ってみました。

タイムラインでそういう話題があったのでネタをパクったような形になってしまいましたが、私自身前からやってみたいと思っていたテーマであり、H2O.aiが最近担いでいる話でもあるので許して頂きたい。

構成

  • Ryzen 7 1700 (3.7GHz 8core/16thread)
  • GTX 1080 Ti FOUNDERS EDITION (1582MHz 11GB)

今のところ定格で動かしてます。 記事と前後してしまいますが、Deep Learning以外の機械学習アルゴリズムGPUで気軽に計算を高速化できるような予感があったのでCPUはケチりました。 将来GPUを追加する際の資金の足しにしたい。

インストー

公式のドキュメント

私はUbuntuの環境でmakeを使ってビルドしました。

CUDA

CUDAをインストールして

export PATH="$PATH:/usr/local/cuda-8.0/bin"
export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:/usr/local/cuda-8.0/lib64"

あたりのパスを通しておきます。nvccが動けば問題ないはずなので、cudnnは不要と思われます。

ビルド

make/config.mkファイルの該当行を

CUB_PATH = cub
PLUGIN_UPDATER_GPU = ON

と変更します。

後は通常と同様にmakeするだけです。

インストールについてもpython3 setup.py installとすれば問題ありませんでした。(Anaconda)

使い方

updaterパラメータを変更してGPUでの計算に切り替えることができます。 今までと同じパラメータを使っている限りではGPU版をビルドしても通常のCPUによる計算となります。

param['updater'] = 'grow_gpu'
  • GPU fast histgram
param['updater'] = 'grow_gpu_hist'
param['max_bin'] = 16

max_binLightGBMとかのアレなので、問題に合わせて設定してください。

サンプルをぼんやり眺めて、GPUが1つしかないのに

param['gpu_id'] = 1

とやるとエラーになって死ぬので注意です。

実験

CPU/GPUとexact/fast histgramで計算時間を比較しました。

データは以前の記事と同様にOttoのデータを利用しました。

github.com

f:id:marugari2:20170611104703p:plain

ntread CPU exact CPU hist GPU exact GPU hist
1 143.213904 43.093807 10.793109 7.385136
2 75.984772 24.485343 10.793109 7.385136
3 56.342093 18.662806 10.793109 7.385136
4 43.276191 14.980951 10.793109 7.385136
5 35.087812 12.747035 10.793109 7.385136
6 29.960601 11.476474 10.793109 7.385136
7 26.535089 11.396709 10.793109 7.385136
8 25.082611 9.916659 10.793109 7.385136
9 24.930102 11.761317 10.793109 7.385136
10 23.934886 11.543043 10.793109 7.385136
11 23.353146 11.724661 10.793109 7.385136
12 23.154665 11.528249 10.793109 7.385136
13 22.851917 11.614104 10.793109 7.385136
14 22.556055 11.890477 10.793109 7.385136
15 22.435674 11.309677 10.793109 7.385136
16 22.297691 11.825569 10.793109 7.385136

CPUよりGPU、exactよりもfast histgramが高速ですが、公表されているベンチマーク程は差が出ませんでした。

CPU fast histgramでも十分速い。X付きのにしときべきだったかなぁ…

nvidia-smiでモニタリングしたところGPUはまだまだ余裕ある感じだったので、データが巨大になった場合には差が出てくるのかもしれません。

また今更ですが、ntreadをコア数以上に設定しても速くなっていません。

LightGBMは?

LightGBMでもGPU対応されています。 公式のチュートリアル

インストールの仕方の説明でドライバをaptで入れてたので何となく回避しました。

OpenCLを使っているようなので、QuadroAMDのグラボを持ってる人には面白いかもしれません。

最後に

今日、我々はDeep Learningに留まらず勾配ブースティングにおいてもGPUモードをポチっとするだけで計算を高速にできるようになりました。 専門家以外にも扱い易い機械学習アルゴリズムGPUで使われるようになれば、ますますGPU需要が高まることが予想され、AI半導体メーカーの株価動向からは目が離せません。

勾配ブースティング落穂拾い - マルチクラス分類について

勾配ブースティングに特有の問題という訳ではないですが、それはマルチクラス分類をone-vs-allで解いていることに留意しましょうというのが今回の記事です。

2017/04/02追記 LightGBMではマルチクラス分類の際にsoftmaxとone-vs-allのどちらで解くか選択できるようです。 また、XGBoostではsoftmaxのみがサポートされています。 この記事では、マルチクラス分類のためにクラス数と同じ分類器を使っていることを指してone-vs-allとしてしまっています。ご注意ください。

マルチクラス分類について

ここではマルチクラス分類の問題を対象に、アンサンブル学習として代表的なランダムフォレストと勾配ブースティングとを比較します。

ランダムフォレスト

ランダムフォレストによるマルチクラス分類では学習器に分類木を用いています。 scikit-learnではDecisionTreeClassifierです。

コード

分類木はマルチクラスを自然に扱えるので、n_estimatorsと同じ数だけ木を作れば求めるモデルが構築できます。

勾配ブースティング

勾配ブースティングでは弱学習器に回帰木を用いているので事情が変わってきます。 マルチアウトプットの回帰木は一般的ではありませんので、問題をone-vs-allに変換して解くことになります。 そのため勾配ブースティングでは指定したn_estimatorsとクラス数の積の数だけ木を作る必要があります。

コード

勾配ブースティングは速度と精度を高いレベルで両立していますが、マルチクラス分類では

  • クラス数と比例して計算時間が増加する
  • Dark knowledge的なものを利用できない

という欠点が出てくるので、勾配ブースティングのメリットがスポイルされていないか留意する必要がありそうです。

数値例

勾配ブースティング(というかone-vs-all)が苦手そうな例を考えてみました。 セルの中は(クラスAの割合, クラスBの割合)を示しています。

sex/height
(5%, 10%) (15%, 10%) (20%, 20%)
(15%, 15%) (5%, 15%) (20%, 30%)
(20%, 25%) (20%, 25%)

この例にone-vs-allを使うと、クラスBについてsexで分割する以上のものが出てきません。 一方で分類木を使った場合は、クラスBについてsexで分割した後でクラスAの濃い部分を抽出するheightでの分割を発見できます。

マルチクラス分類をXGBoostのみで取組むのは微妙な感じがしますね。 あまり業務でスタッキングしたくないのですが、それも避けられないように思われます。 もしくはDeep Learningするか。

参考

ブースティングやDeep Learningが一般化した時代におけるマルチクラス分類の解き方についてはOtto Group Product Classification ChallengeのForumで盛んな議論があったようです。 実戦的なテクニックも数多いようなので、関心のあるかたはそちらで調査されるのが良いと思われます。

勾配ブースティング落穂拾い - 木の構築について

このシリーズについて

XGBoost芸人を自称してちょこちょこ活動をしてきたのですが、最近になって自分の理解の甘さを痛感するようになりました。

気になった箇所を散発的に復習しているのですが、その成果を備忘録として残しておこうと思います。 今のところ体系的にまとめるつもりはないので、これを読んでも勾配ブースティングの全体像はつかめませんので悪しからず。

今回のテーマ以外にはマルチクラス分類の際の挙動等に関心を持っています。

木の構築について

勾配ブースティングでは  M+1 回目のイテレーションで誤差

 \displaystyle
\sum_{i=1}^N L \left( y_i - \sum_{j=1}^M \eta_j f_j \left( x_i \right) \right)
= \sum_{i=1}^N L \left( y_i - f^M \left( x_i \right) \right)

の勾配を上手く表現した木を構築します。

この部分の処理についてscikit-learnとXGBoostでの違いを確認します。

scikit-learn

カステラ本に準拠した処理になっています。 勾配の計算は

 \displaystyle
\frac{\partial}{\partial f^M \left( x_i \right)} L \left( y_i - f^M \left( x_i \right) \right)

となり、これを各サンプルのラベル扱いにして DecisionTreeRegressor に投げます。

コード

このとき当然 DecisionTreeRegressor は2乗誤差*1を最小にするよう葉に属するサンプルにスコアを付けます。つまり平均値となっています。

ロス関数がRMSEのケースではそれで問題ないのですが、他の関数ではズレが生じます。 そのため、木の分岐は DecisionTreeRegressor の結果を使うのですが、最終的なスコアには補正をかけて対応します。

例えば、ロス関数がMAEの場合は、葉に属するサンプルの中央値をスコアとして採用するのが、木の構造を所与としたときの最適となります。

コード

XGBoost

XGBoostでは正則化項が追加されているので最適化の目的関数が異なります。式の導出は別に譲りますが以下のようになります。

 \displaystyle
\sum_{j=1}^M \left( G_j w_j + \alpha \left| w_j \right| + \frac{1}{2} \left( H_j + \lambda \right) w_j^2 \right) + \gamma T

木の構造を所与としたとき最適な  w_j^*

 \displaystyle
w_j^* = \frac{l \left( G_j , \alpha \right)}{H_j + \lambda} \\
\mathrm{where} \quad l(x, y) =
\begin{array}{l}
x - y \quad \mathrm{if} \ x \geq y \\
x + y \quad \mathrm{else}
\end{array}

となります。これは最適化の1階条件

 \displaystyle
G_j + \alpha \mathrm{sign} \left( w_j^* \right) + \left( H_j + \lambda \right) w_j^* = 0 \\
\displaystyle
w_j^* = -\frac{G_j + \alpha \mathrm{sign} \left( w_j^* \right)}{H_j + \lambda}

によります。ここで  w_j^* \leq 0 となるためには  G_j \geq \alpha でなければなりません。

これを代入すれば目的関数から  w_j が消去できます。 したがってXGBoostでは、ある分割についてはじめから最適なスコアを付ける前提で、最適な分岐を探索していることになります。

ポエム

XGBoostは通常の勾配ブースティングより精度が良いとされており、その原因は主としてL1/L2正則化にあるとされています。 確かに通常の勾配ブースティングでもインサンプルの誤差を0にするのは簡単なので、汎化性能の高さは正則化によるものと考えるのが自然です。

しかしながら、今回見てきたようにXGBoostは通常のものより精度の高い近似計算をしています。 個人的には正則化がそれほど効いていると思えないので、最適化の工夫によって高い汎化性能が実現しているのではないかと疑っています。 XGBoostの収束が速い分だけ、同水準のインサンプル誤差で比較したときに森の複雑性がより小さいと期待できるのではないでしょうか。

参考

*1:オプションで絶対誤差も指定可能

4教科で珍しく良い点取った生徒のカンニングを疑ってはいけない、「竜王を弁護する」について

はじめに

将棋ソフト不正疑惑に関する事件について新たな記事が出ていたので、また計算機を叩いてみました。

三浦九段不正疑惑について、渡辺明竜王を弁護する

2017/02/27追記 当該記事は削除されてしまっています。一時、ハフィントンポストが全文を掲載してくれていたのですが、そちらも削除されています。ご覧なりたい方は自身で魚拓を漁ってください。

三浦九段については晴れて「シロ」判定がされて、これから名誉回復をどうしていくかという段で、このような記事が出てくるのには違和感があります。

ネット上でもあまり評判良くないようなので、今更内容に噛み付く必要もないのですが、データジャーナリズムの手習いとして公開してみます。

前提

ここでは、記事中にある一致率の推定値を真の値として扱います。

一致率については、将棋ソフトや局面の違いによって値に大きな差異が出る可能性が指摘されています。 そのような数字を基に組み立てたロジックは無意味だと一蹴することもできるのでしょうが、一致率を受け入れた上でも反論が可能と私は考えました。

設定

記事中では、ある対局で高い一致率の示される確率について2つの数字が提示されています。

  • 10%: 「当方の計算では、出現率でいえば、良くて10%程度」
  • 25%: 「第三者委員会の報告は、(中略)上位4分の1に入るような出来」

はっきりと「ある対局で高い一致率の示される確率」と書かれている訳ではないですが、文脈からそう判断して問題ないと思われます。

私も以前の記事で計算をしたのですが、ある対局で一致率85%以上となる確率を9.79%と概算しています。 不安定な一致率の推定値から計算したものなので、これが正しいと考えてもいませんが、ここでは三浦九段に不利な10%を採用してみます。

確率計算

弁護記事では「そんな出来の良い将棋が、疑いをもたれた短期間に、4局も」あったと指摘しています。「短期間」が何を指すか明らかではありませんが、私の前回記事では某まとめサイト情報を参考に12局としていました。

12局中4局以上で高い一致率が示される確率

r   = 0.1
num = 12
lim = 3
s   = 1.0 - scipy.stats.binom.cdf(lim, num, r)

2.56%となります。 教科書にあるような有意水準5%だと、通常では起こりにくいことが起こっていて疑わしいとなります。

以前の記事では、12局というのが恣意的に切り取った一部なのでシロともクロとも言えないとコメントしていましたが、よくよく考えれば恣意的な切り取りについても数値的に評価できそうなので、その計算を今回追加します。

一致率の高くなる期間を切り取れる確率

三浦九段の2015年度対局数は42局です。 42局の中に一致率の高い4局を含む12局が出現するか、という問題に変わると、この確率は当然ながら2.56%より高くなります。

解析的に解く方法が思いつかないので、モンテカルロ法を使ってみます。エレガントな解法があれば教えてください。

確率10%のクジを42回引いて、結果のrolling_sumの最大値が3を超えるかで判定しています。

num_ext  = 42
num_iter = 10000
res = np.zeros(num_iter)
for ii in range(num_iter):
    samp = np.random.choice([0, 1], num_ext, replace=True, p=[1.0-r, r])
    seq  = pd.Series(samp).rolling(window=num).sum()
    res[ii] = seq.max()
s = (res > lim).mean()

モンテカルロ法なので結果にブレはありますが、17%程度という結果になりました。

トップ棋士は、ベースの一致率が高くなりがちというだけでなく、対局数つまり試行回数が増えてしまうため、一致率による疑いをかけられやすい構造にあるようです。

計算は以下にまとめています。 github.com

主張

サンプルを恣意的に切り取ることで、17%程度発生するありふれた事象を2%程度に見せることができます。 我々が統計の結果を作る・見るときには、このようなテクニックが使われていないか常に注意する必要があります。

いつもソフトを使ってる訳でないので、一致率の高い4局が問題で、12局だとか42局は関係ないという論を張っている人もいるようですが、これは誤りです。 無実の人が同じ試行をしたときに起こりそうもない結果が出たときのみ疑いをかけるのが正しい姿勢です。 羽生さんの対局から一致率の高い1局を引っ張ってきて、いつもはソフト指ししていないけれど、当該局ではソフト指ししていたという主張は馬鹿げているでしょう。

また、弁護記事では試験のカンニングについて言及されていますが、もし教師の方がいらっしゃたら記事のような手続きで生徒を疑うのはやめてください。個人的に疑いを持って物的証拠を探すのはセーフかもしれませんが、疑いを同僚教師に話したりするのはアウトです。 過去の偏差値からは取れそうもない点数を取る程度の事象であれば、「卒業までに受ける試験の数×カンニングが疑われるような素行の生徒数」の回数の試行があれば、簡単に発生してしまうのが確率というものです。

最期に

最近、話題のBIGについて計算したのですが、そちらでは直感通りのつまらない結果になりました。

MXNetを中心としたCustom Loss functionの話

これはDeepLearning Advent Calendar 2016の12月24日分の記事です。

Deep Learning一般の内容ではないので申し訳ないですが、来年はMXNet Advent Calendarやりたいですねという願いもあり、空き枠に滑り込みました。

TensorFlowとChainerが強くてMXNetは肩身が狭いですが、AWSという錦の御旗も得たのでユーザ増加に期待したいです。 MXNetはドキュメントが弱いとよく言われますが、最近はドキュメントサンプルコードも充実してきています。 9月頃はドキュメントがリンク切れしておりサードパーティーのドキュメントを読んでいたくらいなので随分良くなりましたよ。

さて、今回はCustom Loss functionについて紹介しようと思います。機能自体は前からあったようですが今月になって簡単な資料が公開され、一般ユーザにも存在が知られるようになったものです。

ロス関数について

ここでは、モデル出力値とラベルの差を表現する関数で、勾配法によって最小化する対象をロス関数と呼びます。 どのライブラリでも回帰なら2乗誤差、分類ならクロスエントロピーが使われているはずです。

注意が必要なのは、モデル精度の最終的な評価をするための指標とロス関数は別になるケースもあるということです。 分類問題ではクロスエントロピーよりもAccuracyやAUCで精度評価することが多いと思いますが、そういった場合でもロス関数にはクロスエントロピーが選ばれてクロスエントロピーを最小化するのが通常です。

評価指標をAUCに変えたら、それを最適化してくれると思っているユーザは結構多く、私もそういった方への案内誘導をボランティアでやっています。*1

一方で、ロス関数に拘ったからといって報われるとは限らず、絶対誤差で評価されるタスク等でも2乗誤差を最小化する方が上手くいくケースも多いのが実際と思われます。

ごく限定的な状況ですが、標準のロス関数では対応できないような問題のためにロス関数をカスタムするインターフェースが各ライブラリから提供されています。

MXNetでの事情

MXNetのLinearRegressionOutputSoftmaxOutputといった出力層向けのシンボルは、出力とロス関数が一体化しておりパラメータでロス関数を変更できません。 ロス関数を変えたい場合はMakelossというシンボルを利用する必要があります。

回帰の例

LinearRegressionOutputを使った通常の方法は以下のようになります。これはMSEを最小化します。

data  = mx.sym.Variable('data')
label = mx.sym.Variable('label')
fc    = mx.sym.FullyConnected(data, num_hidden=3)
act   = mx.sym.Activation(fc, act_type='relu')
fco   = mx.sym.FullyConnected(act, num_hidden=1)
net1  = mx.sym.LinearRegressionOutput(fco, label)

同じようにMSEの最小化をCustom Loss functionでやると以下のようになります。

loss = mx.sym.square(
    mx.sym.Reshape(fco, shape=(-1,)) - label
)
net2 = mx.sym.MakeLoss(loss)

mx.sym.squareやelement-wise minusの導関数はMXNetが知っているので自動微分を使って上手く最小化してくれる訳です。 numpy等を使って好き勝手にロス関数を作ることはできません。

注意 MakeLossにmx.sym.squareした結果を渡すとpredictした際にも2乗した結果が出力されるそうです。

custom loss symbol in R/Python · Issue #3368 · dmlc/mxnet · GitHub

確かにどれが出力として求められているかはMXNetで判断しにくいかもしれません。これについては使い易いように修正がされるそうです。

分類の例

MXNetのLogisticRegressionOutputは名前のような挙動をしないことが知られています。

RPubs - MXnet LogisticRegressionOutput test

SoftmaxOutputを使うことが推奨されていますが、これは手前のレイヤーで隠れユニット数をカテゴリ数と同じだけ確保する必要があります。 2値分類のときに隠れユニットを2つ持っておくのは冗長な気がするので、Custom Loss functionでロジスティック回帰を組んでみようというのがこの試みです。

SoftmaxOutputを使った通常の方法は以下のようになります。これはLog Lossを最小化します。

data  = mx.sym.Variable('data')
label = mx.sym.Variable('sm_label')
fc    = mx.sym.FullyConnected(data, num_hidden=3)
act   = mx.sym.Activation(fc, act_type='relu')
fco   = mx.sym.FullyConnected(act, num_hidden=2)
net1  = mx.sym.SoftmaxOutput(fco, name='sm')

ロジスティック回帰らしく書いた場合は以下のようになります。

data  = mx.sym.Variable('data')
label = mx.sym.Variable('label')
fc    = mx.sym.FullyConnected(data, num_hidden=3)
act   = mx.sym.Activation(fc, act_type='relu')
fco   = mx.sym.FullyConnected(act, num_hidden=1)
p     = mx.sym.Activation(fco, act_type='sigmoid')
eps   = 1e-6
p     = mx.sym.minimum(mx.sym.maximum(p, eps), 1.0-eps)
q     = 1.0 - p
lp    = mx.sym.Reshape(mx.sym.log(p), shape=(-1,))
lq    = mx.sym.Reshape(mx.sym.log(q), shape=(-1,))
loss  = label * lp + (1.0 - label) * lq
net2  = mx.sym.MakeLoss(loss)

出力を確率に直すやり方がわからないので、ちゃんと動いているかわからないのですが、 方針自体はこれで良いはず…

中途半端になってしまいましたが、これ以上掘るのは効率が悪いので、本日はここまでとします。 手元のJupyter Notebookはこちらです。

github.com

やっぱり

以前に比べればドキュメントが強化されているとはいえ、最近どういう機能が追加されたか等を個人で調べるのは大変です。 ユーザが増えて互助していければ素敵ですね。

*1:こんなことは働き盛りの人間がやる仕事ではないので、ボット化したいなぁ

LightGBMのPythonパッケージ触ってみた

DMLCひとりアドベントカレンダー0日目の記事です。 強い競合が現れたということで、DMLCとは直接関係ないですがLightGBMについて紹介します。

LightGBMとは

勾配ブースティング木の高速な実装としてXGBoostが有名ですが、Microsoftの開発した更に高速な実装がLightGBMです。 実験によるとXGBoostの数倍高速だということです。

パフォーマンス向上の大きな要因は、決定木の分割をexactに探索するのではなく、ヒストグラムのビン単位にして計算を間引いていることにあるようです。 ビン単位にすると精度が落ちそうな気もしますが、汎化性能を考えたら特徴量が粗いくらいでちょうど良いのだろうと腹落ちしています。

同じアプローチでFastBDTが先に発表されましたが開発が滞ってしまっている一方、LightGBMでは最近Pythonパッケージがベータリリースされました。 今日、我々がPythonで勾配ブースティングをする際にはXGBoostかLightGBMの2択*1となります。

導入

ビルドはこちら、インストールはこちらを参考に。

XGBoostとそれほど変わりません。

使い方

Exampleより

xgboost.DMatrix的なものを作って…

# create dataset for lightgbm
lgb_train = lgb.Dataset(X_train, y_train)
lgb_eval  = lgb.Dataset(X_test, y_test, reference=lgb_train)

パラメータをディクショナリに格納して…

# specify your configurations as a dict
params = {
    'task' : 'train',
    'boosting_type' : 'gbdt',
    'objective' : 'regression',
    'metric' : {'l2', 'auc'},
    'num_leaves' : 31,
    'learning_rate' : 0.05,
    'feature_fraction' : 0.9,
    'bagging_fraction' : 0.8,
    'bagging_freq': 5,
    'verbose' : 0
}

trainするだけ。XGBoostと同じ流れですね。

# train
gbm = lgb.train(params,
                lgb_train,
                num_boost_round=20,
                valid_sets=lgb_eval,
                early_stopping_rounds=5)

predictpandas.DataFrameを直接受け入れてくれるのには注意が必要です。

# predict
y_pred = gbm.predict(X_test, num_iteration=gbm.best_iteration)

XGBoostのパラメータとの対応は以下のようになります。

2016/12/18追記 tksさんに指摘を頂いたので記述を修正しました。 github.com

LightGBM XGBoost 候補 備考
objective objective regression, binary, multiclass, ...
boosting_type booster gbdt, dart 線形モデルはなし
max_depth max_depth num_leavesの制約を受けるので注意
learning_rate learning_rate
feature_fraction colsample_bytree
bagging_fraction subsample bagging_freqを1以上にしないと有効にならない
min_gain_to_split gamma
lamnda_l1 alpha
lambda_l2 lambda
metric eval_metrix l1, l2, auc, ...

LightGBMにもmax_depthパラメータが用意されていますが、オフィシャルにはnum_leavesで葉の数を制御します。

また、特徴量をまとめるビンの数をmax_binで設定できます。

詳細は公式のWikiを調べて下さい。

カテゴリ変数のサポート

LightGBMではカテゴリ変数がサポートされています。 One-hot encodingが不要になったので、メモリ的に有利であったり、カラム方向のサンプリングやfeature importanceが直感的になる等のメリットがありそうです。

下の実験からのコード片ですが、traincategorical_featureパラメータにカテゴリ変数として扱うカラムのインデックスをリストで渡してやるだけです。 カテゴリ変数はintのみ対応しているので、文字だったりする場合は0始まりのコードに置き換える必要があります。

obj_lgb = lgb.train(
    prm_lgb, dt_lgb_c, num_boost_round=num_round,
    valid_sets=dv_lgb_c,
    categorical_feature=list(range(len(x.columns.values))))

処理を完全には追っていないのですが、rpartのようなカテゴリ併合をしてくれる訳ではないようで、数値変数であればx <= thresholdで分割されるところが、x is categoryで分割されるようになるだけのようです。

また以下の実験の出力からなのですが、数値変数として扱うと以下のように小数点の閾値で分割されます。

threshold=1.5 3.5 2.5 0.5 0.5 0.5 0.5 0.5 4.5 14.5 0.5 1.5 0.5 0.5 1.5 2.5 0.5 0.5 0.5 0.5 0.5 4.5 0.5 0.5 0.5 0.5 1.5 0.5 2.5 2.5

一方、カテゴリ変数として扱うとデータにある値そのもので分割されます。

threshold=0 0 0 0 1 1 3 2 0 1 0 0 0 1 1 0 0 3 1 1 0 3 0 0 0 0 1 0 1 0

実験

手元でもパフォーマンスを確認してみようと、Ottoのデータで実験してみました。

github.com

モデル MLogLoss 計算時間
XGBoost 0.570366 65.20
LightGBM 0.522378 16.02
LightGBM Cat 0.671286 15.78

この設定ではLightGBMが4倍以上速い結果となりました。精度もLightGBMの方が良好です。

全変数をカテゴリ変数として扱ったLightGBM Catの有り難みがないように見えますが、One-hot encodingによってカラム数が膨らんでいる場合には計算時間の短縮が実感できるはずです。 精度が悪いのは、今回のデータが順序性を持っているためと思われます。

最後に

XGBoost4jのような使い方をすると話は別かもしれませんが、クライアントマシンでせこせこ計算する分にはLightGBMで良さそうです。 開発も活発なので、XGBoostにしかないような機能も続々と取り込まれている上に、コードもスッキリしており改造が比較的容易になっています。

CNTKなんかもパフォーマンス良いようですし、巻き返しに力の入っているMicrosoftさんのLightGBMには将来性も期待できるのではないでしょうか。

次は23日にMXNetのネタを考えています。

*1:scikit-learnのGradientBoostingもMAEで分割ができる等、面白そうではあるんですが…

ピュアオーディオ的な解像度で語る乱数のコク

周回遅れも甚だしいですが、「乱数のコク」について思う所があったので計算をしてみました。

オリジナル

一様乱数と比べて、5つの一様乱数の平均の方がコクがあるとのことです。 シミュレーションで得たヒストグラムは以下のようになりました。

def rand_koku(n, seed):
    np.random.seed(seed)
    u1 = np.random.rand(n)
    u2 = np.random.rand(n)
    u3 = np.random.rand(n)
    u4 = np.random.rand(n)
    u5 = np.random.rand(n)
    return (u1, (u1+u2+u3+u4+u5)/5.0)

f:id:marugari2:20161108002149p:plain

平均を取った側は正規分布のような釣り鐘型になっています。 果たしてこれをコクと呼べるでしょうか。2次モーメントも異なっていますし、分布の形も別物です。 オーディオで言ったらSONYBOSEくらい差があります。

乱数のコクが感じられる例

今回は2種類の正規乱数を用意しました。

1つはBox-Muller法による正規乱数です。 これはnumpyの正規乱数生成にも用いられているベーシックなものです。

def normal_box_muller(n, seed):
    np.random.seed(seed)
    nm = int(n / 2)
    u  = np.random.rand(n)
    u1 = u[:nm]
    u2 = u[nm:n]
    r  = -2.0 * np.log(u1)
    v  = 2.0 * np.pi * u2
    z1 = np.sqrt(r) * np.cos(v)
    z2 = np.sqrt(r) * np.sin(v)
    z  = np.concatenate([z1, z2, -z1, -z2])
    return z

もう1つはMoroの逆関数法による正規乱数です。 正規分布の累積分布関数の逆関数を利用して、一様乱数から正規乱数を得ることができます。

def normal_moro(n, seed):
    np.random.seed(seed)
    a0  = 2.50662823884
    a1  = -18.61500062529
    a2  = 41.39119773534
    a3  = -25.44106049637
    b0  = -8.47351093090
    b1  = 23.08336743743
    b2  = -21.06224101826
    b3  = 3.13082909833
    c0  = 0.3374754822726147
    c1  = 0.9761690190917186
    c2  = 0.1607979714918209
    c3  = 0.0276438810333863
    c4  = 0.0038405729373609
    c5  = 0.0003951896511919
    c6  = 0.0000321767881768
    c7  = 0.0000002888167364
    c8  = 0.0000003960315187
    u   = np.random.rand(n)
    y   = u - 0.5
    f1  = (np.abs(y) < 0.42)
    y1  = y[f1]
    r   = y1 ** 2
    num = ((a3*r+a2)*r+a1)*r+a0
    den = (((b3*r+b2)*r+b1)*r+b0)*r+1.0
    x1  = y1 * num / den
    f2  = np.logical_not(f1)
    u2  = u[f2]
    y2  = y[f2]
    r   = np.where(y2 > 0.0, 1.0 - u2, u2)
    r   = np.log(-np.log(r))
    x2  = c0+r*(c1+r*(c2+r*(c3+r*(c4+r*(c5+r*(c6+r*(c7+r*c8)))))))
    x2  = np.where(y2 < 0.0, -x2, x2)
    z1  = np.concatenate([x1, x2])
    z   = np.concatenate([z1, -z1])
    return z

アルゴリズムはどちらも Paul Glasserman, Monte Carlo Methods in Financial Engineering を参考にしています。

乱数を200,000,000個発生させて得たヒストグラムは以下のようになりました。

  • 両側0.5シグマ f:id:marugari2:20161108003322p:plain 両側0.5シグマの範囲では差を見て取ることができません。

  • 片側4シグマ f:id:marugari2:20161108003533p:plain 同じ一様乱数系列を用いていますが、正規分布への変換法によって片側4シグマの裾の分布が異なっています。

英語版WikipediaではBox-Muller法の弱点について記載があります。 裾の分布に拘る場合はMoroの逆関数のようなものを使うことが推奨されているようです。

何気なくrandnで正規乱数を得ることも多いですが、そのアルゴリズムの違いによって得られる結果が異なることを見てきました。 こういったところにこそ「乱数のコク」があるのではないでしょうか。

最後に

正規分布は非常に基礎的な確率分布ですが、それに従う乱数を得るためのアルゴリズムは(私にとっては)そんなに簡単ではありません。 正規乱数でさえそうなのです。

ベイズ統計では、複雑な事前分布と複雑な尤度関数の積からなる事後分布に従うサンプルを生成したいケースがあるのですが、そのような複雑な分布関数にはBox-Mullerのような解析的なアプローチもMoroのような精度の良い近似関数も存在しないことが常です。 そのようなマトモでない分布に従う乱数を生成するためにMarcov Chain Monte Carloが必要となっている訳です。

こういった話はRobert and Casella, Monte Carlo Statistical Methodsに感動的にわかりやすく書かれているので、私は強くオススメをしています。