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に感動的にわかりやすく書かれているので、私は強くオススメをしています。

一致率90%以上はどれくらいあり得るのか計算してみる(ガバガバ設定)

将棋界で大騒ぎになっている疑惑について、連盟公式掲示板でのレスバトルをするための材料として手元で計算をしてみました。

計算の前提となっているパラメータがガバガバなので、白か黒かの材料にはなりようがないですが、とりあえず今出ている数字を使ってみたらどうなるかというメモです。

私は

  • 統計検定持ってません。変な計算してたら叱って下さい。
  • 調査設計については無知です。調査法がイケてない点については触れません。
  • 観る将なので、戦法や序盤・終盤でどれくらい手の広さが異なるか知りません。
  • 電王戦は見てるけど、ponaと技巧でどれくらい強さが違うかは知りません。

設定

千田先生のDropboxにある資料の数字から独自に作成します。 「千田棋譜 資料」内に当該ファイルがあります。 www.dropbox.com

  • ソフト指ししないプロ棋士の一手毎のソフト一致率は75%とする。
  • 一致するかどうかは系列相関なしで独立とする。

2点目の仮定について、一本道の変化だと一致率の高い局が多くなるはずなのでシロ派に不利な仮定と考えられます。一方、全部の手をソフト指しする必要はなく、ソフト指し側が一致率をコントロールできるという点ではシロ派に有利な仮定とも考えられます。 これは水掛け論になるので、計算を簡単にするために独立としておきます。

計算結果

一局で一致率が90%以上になる確率

序盤を通過した後で終局まで30手*1指すとします。27手以上一致すれば率が90%以上になります。

p   = 0.75
num = 30
lim = int(num * 0.9) - 1
q   = 1.0 - scipy.stats.binom.cdf(lim, num, p)

3.74%になりました。結構高いような気もします。

実績値と比べて検討もできますが、私にはヘビーなのでこの数字を使うことにします。

一致率の高い対局が複数出る確率

疑惑の7月以降12局*2で一致率90%以上が複数回現れる確率はどうなるでしょう。

r   = q
num = 12
s1  = 1.0 - scipy.stats.binom.cdf(0, num, r)
s2  = 1.0 - scipy.stats.binom.cdf(1, num, r)
s3  = 1.0 - scipy.stats.binom.cdf(2, num, r)

1局以上現れる確率は36.7%、2局以上現れる確率は7.2%、3局以上現れる確率は0.9%です。 三浦久保戦の一致率が高かったのも考慮して、2局あったとすると*3確率は7.2%です。

棋士を十数人集めれば、1人くらいは石に当たる奴が出てくる計算ですね。

90%という区切りが悪かった?

疑惑の4局は一致率が85%以上です。区切りを変えてみましょう。

p   = 0.75
num = 30
lim = int(num * 0.85)
q   = 1.0 - scipy.stats.binom.cdf(lim, num, p)

9.79%です。

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

12局中で4局以上現れる確率は2.4%です。

この設定だと5%は割ってきますが、12局がいいとこ取りなのをどう見るかでしょうか。

所感

最初から分かっていたことですが、一致率でクロ・シロの議論はできそうもありません。 しかし、当初想定していたよりはシロを支持できる数字が出たようにも思えます。

レスバトルで盛り上がる前に計算機を叩くよう心がけたいです。

*1:対局毎に終局までの手数は異なりますが、計算を簡単にするため30で固定しています。

*2:いいとこ取りになるので、この局数で良いのかという議論はあると思います。

*3:ある一局の一致率が高かったからといって、ヒット数を多くするのは適正ではないですが、サンプルも少ないのでとりあえず。

MXNetでmulti-input/multi-output

皆さんMXNet使っていますか? 年度初に著名データサイエンティストの記事が相次いで盛り上がった感がありましたが、もうChainerなりTensorFlowなりに移ってしまったのでしょうか…

MXNetはDeep Learningフレームワークの比較でドキュメントが弱いことをよく指摘されてるので、1ユーザとして草の根でお役立ち情報を発信していきたいです。

やりたいこと

形式の異なる複数のデータを入力として、複数の値を出力するモデルを学習したい。*1

  • Keras: Functional APIなるものを使って実現できるそうです。

Functional APIのガイド - Keras Documentation

  • Chainer: サポートされている模様。手続き的なフレームワークだと関数に通すだけなので難しいことは少なそう。

Google グループ

  • TensorFlow: サポートされていない模様。MXNetと同様に改造すればいけそう。

How to handle multi tensors input? · Issue #9 · tensorflow/serving · GitHub

参考にしたもの

MXNetで困ったら、まずはexampleを調べることをお勧めします。 英文を含めてもブログ記事の情報量は少なく、Issueには答えが書いていないものも多いです。

私もIssueを読み込みましたが、答えはexampleにありました。

mxnet/example/multi-task at master · dmlc/mxnet · GitHub

サンプルコード

github.com

  • 3ブロック目: 解像度の異なる画像2つと、ベクトル1つをインプットにしています。
x1 = np.zeros((num_train, 1, 8, 8))
x2 = np.zeros((num_train, 1, 16, 32))
x3 = np.zeros((num_train, 10))
y = np.zeros((num_train, num_cls))
z = [y[:,ii] for ii in range(y.shape[1])]
  • 6ブロック目: 公式のDataIterを使うとシンボルの名前が嫌になります。インプットのデフォルト名はdataで、アウトプットのデフォルト名はsoftmax_labelになります。出力層の名前を自分で勝手に決めるとエラーが出るのはこいつが原因です。
print(mx_dat0.provide_data)
print(mx_dat0.provide_label)
  • 8ブロック目: DataIterを改造して、名前を付けやすくしています。

  • 9ブロック目: 適当なネットワークを組みました。

data0 = mx.sym.Variable('input0')
data1 = mx.sym.Variable('input1')
data2 = mx.sym.Variable('input2')
def get_symbol(sym, prefix=''):
    net = mx.sym.BatchNorm(sym,
                           name=prefix+'_bn')
    net = mx.sym.FullyConnected(net,
                                name=prefix+'_fc',
                                num_hidden=3)
    return net
net0 = mx.sym.Flatten(data0)
net1 = mx.sym.Flatten(data1)
net2 = get_symbol(data2, 'i2')
netc  = mx.sym.Concat(*[net0, net1, net2])
fc   = []
out  = []
for ii in range(num_cls):
    fc.append(mx.sym.FullyConnected(netc,
                                    name='fc'+str(ii),
                                    num_hidden=2))
    out.append(mx.sym.SoftmaxOutput(fc[ii],
                                    name='clf'+str(ii)))
net = mx.sym.Group(out)

10ブロック目でプロットしたネットワーク図でinput1input2が見えませんが、インプットのデータを弄るとエラーが出るので、これで大丈夫のはずです。

注意

公式でもmulti-inputやmulti-outputのためのIF改善がTODOに挙げられています。 また、MXNet自体がNNVMとの機能分割で大工事中ですし、将来のバージョンアップで諸々変わってしまう可能性があります。

v1.0 Stable Release TODO List · Issue #2944 · dmlc/mxnet · GitHub

その他

MXNetのLogisticRegressionOutputは名前から想像するような動きをしてくれないので、SoftmaxOutputを使った方が良さそうです。

*1:Deep Learningは人間の脳を模倣した仕組みなので、複数種類のインプット(視覚・聴覚・記憶等)から複数のアウトプット(話す内容・身振り・手振り等)を得るためのルールを学習できます。(嘘)

分類性能とMAE

名刺お疲れ様でした。 終盤はGPUメモリ足りないとイライラしていましたが、足りないのは創意工夫でした。

私はMXNetを使って取り組んでました。multi-inputとかmulti-outputのやり方を泣きながら調べたのも良い思い出です。

その話はまた時間があればまとめます。今日は軽いネタを一つ。

本題

今回のコンペはマルチラベル問題で、評価指標はMAEでした。 サンプル数が  N で、そのインデックスを  i とします。同様にクラス数が  M で、そのインデックスを  j とします。  i 番目のサンプルの クラス  j に対するフラグを  y_{i,j}、予測値を  f_{i,j} と表します。

MAE(Mean Absolute Error)は

{\displaystyle
\mathrm{MAE} =
\frac{1}{NM} \sum_{1 \leq i \leq N} \sum_{1 \leq j \leq M}
\left| y_{i,j} - f_{i,j} \right|
}

で計算します。

このとき、 f は確率ではなくフラグにしないと殆どの場合で損するはずです。 MAEにはmedianを返すのがお約束なので、これ以上は不要な説明かもしれませんが、泥臭くやってみます。

何らかの機械学習で得た予測値  g_{i,j} があるとして、これを  \Delta_{i,j} で補正した  f_{i,j} = g_{i,j} + \Delta_{i,j} を予測値にします。

補正したことによるMAEの改善幅は

 {\displaystyle
q_{i,j} \Delta_{i,j} - \left( 1 - q_{i,j} \right) \Delta_{i,j}
= \Delta_{i,j} \left( 2 q_{i,j} - 1 \right)
}

となります。  q_{i,j} y_{i,j} = 1 となる確率とか確信度とかそういった数字です。 先の機械学習モデルがまともなものであれば、 q_{i,j} g_{i,j} で代用できるでしょう。 そうすると、 g_{i,j} \geq 0.5 では  f_{i,j} = 1 にすべきですし、 g_{i,j} \leq 0.5 では  f_{i,j} = 0 にすべきです。

したがって、評価指標は実際にはAccuracyになってしまっていたんじゃないか、というお話です。

以下のチュートリアルでも、特に説明はありませんが確率でなくフラグを出力しています。

コンテストチュートリアル(Sansan)

余談

分類器の性能をMAEで評価するのは意味がないと個人的には考えるのですが、そのようなペーパーも少ないながら発見できます。 以下のような表で比較をしていたりします。

classifier Accuracy MAE RMSE
clf1
clf2

著者がインド系の方が多いので、インドの偉い先生の好みだったりするんでしょうか。 日本人がMAEで分類性能を評価しているものを見つけてオッと思ったら指導教官がインド系の先生だったりもしました。

XGBoostにDart boosterを追加しました

はじめに

XGBoostにBoosterを追加しました。 以下のようなIssueを見つけ、興味があったので実装してみたものです。 github.com

今のところ、Dart boosterの仕様について私以外の誰も把握していないはずなので、皆さんに使って頂きたく解説記事を書きます。*1

モチベーション

論文の

Boosted Treesでは誤差を潰すために回帰木を大量に作ります。

木の数が多くなったときに残っている誤差は小さいですが、以降に構築される回帰木はその些末な誤差にフィッティングされることになります。

これは効率悪いように思われ、イテレーションの終盤においても一定の影響力を持った回帰木を作りたいということで、NN系でよく用いられる(ものとは趣きが異なるように私には思える)DropoutをBoosted Treesに転用しています。

自分の

  • XGBoostにはいつもお世話になっているのでcall-for-contributionに反応したかった
  • ライブラリ使ってるだけでしょ、という煽りが癪
  • Gitできないので5日でクビになってしまう
  • XGBoostエバンジェリストとしてSIerさんに拾ってもらいたい

実装について

DARTにおける回帰木の作り方は通常と同じです。 構築した回帰木に対するDropoutやスケーリングの処理を追加すれば良いので、回帰木をコントロールしているGradientBooster(今回は特にGBTree)を拡張するのみです。

ドロップ

Dart::PredictDropTreesをコールしています。 ここでntree_limitに0がセットされていると、既に構築した回帰木の幾つかをドロップしたleaf scoreを計算します。

したがって、学習結果を用いて予測値を得る場合にはntree_limitに正の値(num_roundと同じ数)を設定する必要があります。 ここはAPIを弄りたくなかったので複雑な仕様になっています。

スケーリング

ドロップしたleaf scoreに対して通常と同じように回帰木を構築します。

ドロップした回帰木と追加の回帰木で平均を取るためにDart::CommitModel内でNormalizeTreesをコールしています。

パラメータ

booster

dartを指定します。 内部的にはGBTreeを継承していますので、パラメータもgbtreeと同じものが設定できます。 Predメソッドでのバッファを無効化しているので、gbtreeと比べて学習が遅くなっています。

Dart boosterに固有のパラメータを以下で解説します。

sample_type

ドロップする回帰木のサンプリング方法を指定するオプションです。 デフォルトはuniformで、このとき一様にドロップする回帰木が選ばれます。 weightedに指定すると、影響力の大きな回帰木が選ばれやすくなります。

normalize_type

DARTではドロップする回帰木と新たに学習した回帰木の平均を取ります。 このときの平均の取り方を指定します。 デフォルトはtreeで、ドロップした回帰木のそれぞれ1本が新たに学習した回帰木と同じ重みを持つとして計算します。

 {\displaystyle
D = \sum_{i=1}^k F_i \sim \tilde{F}
}

 {\displaystyle
\alpha \left( \sum_{i=1}^k F_i + \frac{\lambda}{k} \tilde{F} \right)
\sim \alpha \left( 1 + \frac{\lambda}{k} \right) D = D, \quad
\alpha = \frac{k}{k + \lambda}
}

learning_rateを小さく、num_roundを大きくした場合に影響力が小さくなりすぎるので、forestというオプションを用意しています。 ドロップした回帰木の合計と新たに学習した回帰木が同じ重みを持つとして計算します。

 {\displaystyle
\alpha \left( \sum_{i=1}^k F_i + \lambda \tilde{F} \right) \sim \alpha \left( 1 +\lambda \right) D = D, \quad
\alpha = \frac{1}{1 + \lambda}
}

rate_drop

既存の回帰木がドロップされる確率を設定します。

skip_drop

論文ではlearning_rateが1で固定されているので問題となりませんが、learning_rateを小さな値に設定した場合には足踏みを繰り返して学習が進みません。 あるイテレーションではドロップせずに、通常のGBTreeのようにboostingできるような仕様にしています。 イテレーション毎にドロップをしない確率をskip_dropで設定します。

数値実験

こちら

パラメータによってはgbtreeより良い精度を出せているようです。

今後について

Boosted Treesでは通常learning_rateを小さく設定しているため、既にbaggingが効いてしまっています。 論文のアルゴリズムのように、ドロップした回帰木と同じような値を新たな回帰木に学習させて、それらの平均を取っているだけでは劇的な性能改善は望めないでしょう。

過学習は怖いですが、バリデーションのスコアを参照してセレクションを効かせるような仕組みに発展させたいです。*2

その他

  • NormalizeTreesの計算が間違っていたらご指摘ください*3
  • Boosted Trees周辺で数学的に難しくない論文があれば教えて下さい、XGBoostベースで組んでみます

*1:私の英語力の問題で、日本語記事が世界最速です。DARTでライバルに差をつけるなら今ですよ。

*2:卒業論文くらいにはちょうど良いテーマだと思いますが大学生の皆さん如何でしょう?

*3:既に、Pull requestが通ってから1日で修正をしてしまっている…