NeurIPS2018読みメモ： CatBoost: unbiased boosting with categorical features（その-1：そもそも gradient boosting がわからない）

以下の論文を読みます。

Liudmila Prokhorenkova, Gleb Gusev, Aleksandr Vorobev, Anna Veronika Dorogush, Andrey Gulin. CatBoost: unbiased boosting with categorical features. In Proceedings of NIPS 2018. https://papers.nips.cc/paper/7898-catboost-unbiased-boosting-with-categorical-features

※ キャラクターは架空のものです。解釈誤りは筆者に帰属します。何もわかっていないのでおかしな点はご指摘ください。

f:id:cookie-box:20190101155733p:plain:w60

CatBoost…直訳すると「猫が励ます／後押しする」…なるほど、ツイッターなどをみると猫に励まされている方が大勢いますものね。

f:id:cookie-box:20190101160814p:plain:w60

違うよ？

はい。この論文は CatBoost という新しい「勾配ブースティング（gradient boosting）」のツールについて解説するものであるようです。こ、勾配ブースティングとは？勾配が励ますということですか？猫ちゃんが励ましてくれるのではなく？

いや、機械学習の文脈なんだから、「勾配ブースティング」というのは、モデルを理想のモデルに近づけるためにどんどん更新していくときに損失関数が小さくなる方向（勾配方向）に何か効果的に更新できる手法ということじゃないのかな…あ、ウィキペディアに「勾配ブースティング」という記事があるよ。

Gradient boosting - Wikipedia

これを読むと、勾配ブースティングは「予測モデルを複数個組み合わせて（アンサンブルさせて）最終的な予測を下す」ような学習をする場合の学習手法の一つらしい。このときの個々のモデルを「弱学習器」、組み合わせた全体のモデルを「強学習器」というみたいだね。

「弱学習器」を組み合わせて「強学習器」にする、ですか…でも、適当にその辺の弱い人々を100人集めてチームを組ませたところで藤井聡太七段に将棋で勝つことはできませんよね？

弱学習器に「訓練されていない」って意味はないからね！？ちゃんと訓練して！？というか弱学習器たちをどう訓練するべきかって手法が「勾配ブースティング」だからね。それでその勾配ブースティングでは、以下のように弱学習器 $h_m(x) \; (1 \leqq m \leqq M)$ を足し上げて強学習器 $F_M(x)$ にするらしい。各 $m$ ステップでは1時点前の強学習器 $F_{m-1}(x)$ の残差 $y - F_{m-1}(x)$ を埋めるような弱学習器 $h_m(x)$ を得ようとするイメージだと紹介されているね。

$\hat{y} = F_1 (x) \equiv h_1(x)$ で与えられる $\hat{y}$ がよい予測値になるような $h_1(x)$ を得る。
$\hat{y} = F_2 (x) \equiv F_1(x) + h_2(x)$ で与えられる $\hat{y}$ がよい予測値になるような $h_2(x)$ を得る。
$\hat{y} = F_3 (x) \equiv F_2(x) + h_3(x)$ で与えられる $\hat{y}$ がよい予測値になるような $h_3(x)$ を得る。
$\cdots$
$\hat{y} = F_M (x) \equiv F_{M-1}(x) + h_M(x)$ で与えられる $\hat{y}$ がよい予測値になるような $h_M(x)$ を得る。

もう少しきちんと説明すると、いま、考えている損失関数 $L$ の平均を最小化するような関数 $\hat{F}$ がほしい。 $\hat{F} = \underset{F}{\rm arg \, min} \, E_{x,y} \bigl[L \bigl( y, F(x) \bigr) \bigr]$ この $\hat{F}$ をどう手に入れようかということになるけど、ここではある関数クラス $\mathcal{H}$ に属する関数 $h_m \in \mathcal{H}$ を $M$ 個もってきてその重み付き和（＋定数）として得ることにする。 $\displaystyle \hat{F}(x) = \sum_{m=1}^{M} \gamma_m h_m(x) + {\rm const.}$ それでこの $\{h_m\}$ を得るのに、ここでは以下のように訓練データ上での損失を最小化することにしたい。 $\displaystyle \begin{split} F_0(x) &= \underset{\gamma}{\rm arg \, min} \sum_{n=1}^{N} L(y_n, \gamma) \\ F_m(x) &= F_{m-1}(x) + \underset{h_m \in \mathcal{H}}{\rm arg \, min} \left[ \sum_{n=1}^{N} L\bigl( y_n, F_{m-1}(x_n) + h_m(x_n) \bigr) \right] \end{split}$ （☆）けど、これで任意の損失関数 $L$ について最適な $\{h_m\}$ を得るのは一般に不可能なので、 $F_m(x)$ の更新式として上式の代わりに以下の式を用いることにする。これが勾配ブースティングだって。 $\displaystyle \begin{split} F_m(x) &= F_{m-1}(x) - \gamma_m \sum_{n=1}^N \nabla_{F_m-1} L \bigl(y_n, F_{m-1}(x_n) \bigr) \\ \gamma_m &= \underset{\gamma}{\rm arg \, min} \sum_{n=1}^{N} L \bigl(y_n, F_{m-1}(x_n) - \gamma \nabla_{F_{m-1}} L \bigl(y_n, F_{m-1}(x_n) \bigr) \end{split}$ …思うに、本当は各ステップで最も上手く残差を埋める $h_m$ を得たいけど（☆）、たぶん弱学習器がとりうる関数空間 $\mathcal{H}$ 内でベストな $h_m$ を見つけるのって難しい。そもそも ☆ は $h_m$ を陽に示してないし。だから代替案として、「1時点前の強学習器 $F_{m-1}$ での $L$ の勾配方向 $\nabla_{F_m-1} L \bigl(y_n, F_{m-1}(x_n) \bigr)$ 」を正解ラベルとして、 $\mathcal{H}$ 内でこれを最もよく当てる $h_m$ を見つけることにする。これだって厳密な最適解を見つけるのは難しそうだけど、いま「 $\mathcal{H}$ から最もよく正解ラベルを当てる弱学習器を選ぶ」枠組みというのはあるはずだから（学習できる弱学習器を選んだはずだから）、それにあてはめればいい。問題は選んだ弱学習器 $h_m$ を強学習器にどれだけ反映するかの係数 $\gamma_m$ だけど、これは「1次元だから適宜損失が最小になるように探索して」って感じだね。そしたら強学習器を $F_m(x) = F_{m-1}(x) - \gamma_m h_m$ と更新すればいい。

え、えっと？要するに、2人目以降は前の人までの判断のミスを埋めることを目指して学習するんですね？そうですね…例えばスパムメール分類器を学習するとして、5つの訓練データの正解ラベルを $(1,1,0,0,0)$ とでもしましょう。1人目が訓練データの特徴からこの正解ラベルを目指して学習した結果が $(1,0,0,0,1)$ であったとします。この場合、2人目は1人目のミスを埋めるように学習しなければならないので、2人目にとっての正解ラベルは $(0,1,0,0,−1)$ ですね、おそらく（符号はなんとでもなるのでどちら向きでもいいですが）。さらに3人目は1人目と2人目の判断を最適に合体させたような強学習器のミスを埋めることを目指すことになります…でもこれ、ミスを埋めていくといわれると確かに強学習器がどんどん改善されていきそうですが、1時点前までの強学習器の損失っぷりが「弱学習器で埋められるようなもの」でなければなりませんよね…？例えば個々の弱学習器が線形回帰モデルで損失関数が2乗誤差の和だったら、1人目ができる限り頑張って学習した時点で残差は平均ゼロのノイズであるはずです。平均ゼロのノイズに線形回帰モデルをフィッティングしても傾きもバイアスもゼロですから、この残差を2人目が埋め合わせようというのはできないのでは…？

1人目と2人目が参照する訓練データが全く同じなら損失が2乗誤差の線形回帰モデルではそうなりそうだね。でもロジスティック回帰による線形分類なら違いそうかな。ロジスティック回帰結果の真のラベルとの誤差を再びロジスティック回帰することはできそうだから。もっとも、こういう組み合わせて強くする学習に向いてるモデルの最たる例が「決定木」らしい。ウィキペディアの「勾配ブースティング」にも「勾配ブースティングは典型的には固定サイズの決定木を弱学習器として使用される」とあるね。

決定木？何ですかそれは？ご神木から得た神託を予測値とするのですか？

違うかな。ウィキペディアに Decision tree learning という項目があるね。

Decision tree learning - Wikipedia

上の記事の最初の図に決定木の例があるよ。"is sex male?" を根ノードにした木構造をしているね。これは「タイタニックの乗客が生存したか否か」の2値分類だ。各乗客の特徴として「性別」「年齢」「ともに乗船した配偶者または兄弟の人数」が与えられているらしい（この木で使用されていないだけで他の特徴もあるかもしれないけど）。各乗客は根ノードの質問から出発して、各ノードの質問の答えに応じていずれかの終端ノードに振り分けられる。緑色の終端ノードに振り分けられたら「生存」、赤色の終端ノードに振り分けられたら「死亡」と判定する。例えば女性の乗客なら真っ先に一番右の終端ノードに振り分けられて「生存」と判定されるけど、実際女性の乗客の生存率は 73% みたいだね。女性ならば「生存」といっておけば 73% は当たる。27% は外してしまうけど。というか図のたった3つの内部ノードだけで全乗客に対して 80% 以上の正解率が出るんだね（ 0.36 * 0.73 + 0.02 * 0.89 + 0.02 * (1 - 0.05) + 0.61 * (1 - 0.17) = 0.8059 ）。タイタニックの乗客で生存したのは雑にいうと「女性か、または9歳以下でかつともに乗船した兄弟が2人以下の男児」だったらしい。

同じ9歳以下の男児でも兄弟の人数でこうも生存率が変わりますか…3人以上の兄弟と乗船した9歳男児は下の兄弟の避難を優先させたのか、あるいは三等船室の乗客であった確率が高かったのでしょうか…。しかし、「タイタニック」はもう21年前の映画なんですね。今時の高校生はタイタニックと聞いて「あー沈没した船ね」となるのでしょうか。

…それで、この図のような決定木をつくる方法だけど、まず根ノードでは最も説明力の高い特徴の閾値によって分岐させて、そこから先は分岐させたデータの中で最も説明力の高い特徴の閾値によって分岐させて…といった手続きによるみたいだね。そうすればなるべくバランスのよい木はつくれるけど、データを決定付けるよい閾値を拾いきれるとは限らない。だから、いくつかの異なる木を組み合わせたくなってくる。組み合わせる方法の1つは訓練データの1部ずつをつかっていくつも木をつくって多数決させるというもので、この手法は「バギング」といってランダムフォレストもこの一種だ。組み合わせる別の方法は、バギングのように並列に木を組み合わせるのではなくて、直列に組み合わせる。それぞれの木はそれまでの木のミスを補うように学習する。これが「ブースティング」で、その具体的な強力なアルゴリズムの一つが「勾配ブースティング」だね。

なるほど。先ほどの「勾配ブースティング」の手続きを、決定木を弱学習器として実行するわけですか。決定木であったとしても、各弱学習器 $h_m$ はそれまでの強学習器の残差を正解ラベルとして学習し、得られた $h_m$ を損失が最小になるような係数 $\gamma_m$ で強学習器に足し合わせるというのは変わりませんね。

なんだけど、決定木の場合はもう少し改良できるってかいてある。決定木って特徴空間を直方体の領域に分割するよね。この直方体ごとに $\gamma_m$ を変えても構わない。というか、 $h_m$ の領域分割さえわかれば、 $h_m$ が各領域に割り当てた値に興味もなくて、各領域にどれだけの定数を足し合わせれば強学習器の損失が最小になるかを探索するだけだね。

…確かに全ての領域について同じ係数で足してやらなければならないということもありませんね。新しい弱学習器は、特徴空間内の一部の領域の判定精度を向上させるが、その領域の外はむしろそれまでの強学習器のままでよく余計に何か足さないでほしい、ということもあるでしょうし。…しかし、こうウィキペディアだけ読んでいると雲をつかむようです。実際に何か決定木を学習してみることはできないでしょうか。

決定木というかむしろ勾配ブースティング決定木のパッケージはいくつもあるから、何か動かしてみればいいんじゃないかな。

そうですね、では xgboost というパッケージのR版をつかってみましょう（スクリプトは記事の下部）。データは…タイタニックのデータを拾ってくるのも面倒なので iris でいいです。木の数は指定できるのですね。しかし木を2つに指定しても6つできてしまうようですが…ああ、3クラス分類だから弱学習器1つにつき木が3つなのですね。面倒なので2クラスに絞りましょう。しかし、2つ目の弱学習器によって1つ目の弱学習器単体よりも分類精度がよくなる特徴と木の深さの組み合わせがなかなかないですね…特徴をがく片の幅と長さだけに絞り、木の深さを4にすると申し訳程度に訓練誤差が減少しますね。

[1]	train-error:0.260000 
[2]	train-error:0.240000 
   pred
y    0  1
  0 34 16
  1  8 42
        Feature       Gain     Cover Frequency
1: Sepal.Length 0.94330762 0.7858799 0.6666667
2:  Sepal.Width 0.05669238 0.2141201 0.3333333

f:id:cookie-box:20190201004758p:plain:w640

しかし木まで図示してくれるのですか。便利なものですね。

こうなると特徴空間の分割もみたいよね。

特徴空間の分割ですか…これでどうでしょう。白丸が "versicolor"、黒丸が "virginica" のデータで、破線が各学習器の分割境界、青が濃いほど "versicolor" と判定、赤が濃いほど "virginica" と判定という意味です。赤と青の濃さは上のグラフの Value に比例しています。一番左の図が1つ目の弱学習器、真ん中の図が2つ目の弱学習器で、一番右の図は結果の強学習器です。2つの弱学習器の図を単純に重ね合わせたら赤と青が相殺しないのがいけてないですね…。

それどうやって描いたの？？

上のグラフを見て自力で長方形を描きましたね…。

そっか…。

（次回があれば）つづく

スクリプト

library(xgboost)
max.depth <- 4 # 木の最大深さ
nrounds <- 2 # 木の数

### iris データの準備と特徴カラムと損失関数の定義
df <- data.frame(iris)
# いきなり3クラス分類だとわかりづらいので setosa を除外して2クラス分類にする
target_labels <- c("versicolor", "virginica") # "setosa"
num_class <- 1 # 3 # 3クラス以上でない場合は1を指定する
objective <- "binary:logistic" # "multi:softmax" 
df <- df[df$Species %in% target_labels, ]
df$Species <- factor(df$Species, levels=target_labels) # クラスを絞ったので factor 型を修正
# 特徴も2つのみに絞る
feats <- c("Sepal.Length", "Sepal.Width") # "Petal.Length", "Petal.Width"


# 説明変数と目標ラベル
x <- as.matrix(df[,feats])
y <- as.integer(df[,"Species"]) - 1

# 訓練する
model <- xgboost(data=x, label=y, max.depth=max.depth, nrounds=nrounds, 
                 objective=objective, num_class=num_class, verbose=1)

# 訓練したモデルの判定結果（訓練データ上）
pred <- predict(model, x)
if (length(target_labels) == 2) { # 2クラス分類にした場合は予測結果が 0～1 なので閾値で 0 or 1 に
    pred <- ifelse(pred < 0.5, 0, 1)
}
print(table(y, pred)) # 混同行列

# 各特徴の重要度を取得する
imp <- xgb.importance(feats, model)
print(imp) # 重要度

# モデルを図示
library(DiagrammeR)
print(xgb.plot.tree(feature_names=feats, model=model))
print(xgb.dump(model))