"Gradient Boosted Feature Selection" (Xu et al., KDD 2014) メモランダム - 渋谷駅前で働くデータサイエンティストのブログ

本日の輪読会で僕が担当した論文のメモランダムということで、置いときます。

概要

Gradient Boosted Feature Selection (Xu, Huang, Weinberger and Zheng, KDD 2014)

タイトルが示すように特徴量選択をやりたいというのが第一のモチベーションで、これをgradient boosted machineでやろうというお話。ちなみに昨今のKaggleやKDD cupではxgboostが大流行しているわけだが、元はと言えばこれは要はgradient boosted treeである。なので多分これはxgboostの枠組みでもやれるのではないかと勝手に期待してゐる。

1 INTRODUCTION / 2 RELATED WORK

特徴選択が機械学習にとって重要だというのは今更言わなくてもみんな分かっているであろうお話。最も広く用いられている特徴選択手法としてLasso（L1正則化項つき回帰）が挙げられるが、これは要は線形回帰の手法なので非線形では使えないとかそういう話がズラズラ続くので、この辺はばっさり割愛。

ということで、端的に言えば問題意識としては「線形でも非線形でもゴリゴリL1正則化みたいなことをやりたい」という一点に尽きていて、そのためには非線形でも強力な分類器を使うしかなくて、その点ではランダムフォレストかxgboostしかないよね、というのがここでの論点。そこだけ踏まえてこの次に行きましょう。

ちなみにここでは@myamada0さんのHSIC Lassoとかも紹介されているので、興味ある人はどうぞ。

（※L1正則化まわりは例えばこちらなど参考にしてみてください）

3 BACKGROUND

以下データセットは入力ベクトル $\{ \mathbf{x}_1, \cdots, \mathbf{x}_n \} \in \cal{R}^d$ と対応するラベル $\{ y_1, \cdots, y_n \} \in \cal{Y}$ から成り、ラベルの生成分布は未知とする。ラベルが二値・多クラス・連続値であるかどうかは問わない。ただしここでは基本的に二値 $\cal{Y} \in \{ -1, +1 \}$ に着目する。

3.1 Feature selection with the $l_1$ norm

Lassoの復習。L1正則化学習は以下のように書ける。

$min_{\mathbf{w}} \displaystyle \sum_{(\mathbf{x}_i, y_i)} \ell (\mathbf{x}_i, y_i, \mathbf{w}) + \lambda | \mathbf{w} |_1$ (1)

ここで $\ell (\cdot)$ は二乗損失で $\ell (\mathbf{x}_i, y_i, \mathbf{w}) = (\mathbf{w}^T \mathbf{x}_i - y_i )^2$ である。もちろん損失関数は他にも定義できて、例えば二値分類ならlog-lossつまり $\ell (\mathbf{x}_i, y_i, \mathbf{w}) = log( 1 + exp (y_i \mathbf{w}^T x_i) )$ となる。

3.2 The capped $l_1$ norm

L1正則化には目的が2つある。1つは汎化性能を上げること、そして特徴選択のスパース性を上げることである。ただしこの両者はある程度トレードオフの関係にあるため、最適化に当たっては何かしらの制約をかける必要がある。これを補正するためのアイデアとして、キャップつきL1ノルムという考え方がある。

$q_{\epsilon} (w_i) = min ( |w_i|, \epsilon)$ (2)

通常のL1ノルムに対してこの考え方を取り入れる利点は、ひとたび特徴量選択に用いられたらそれ以上ペナルティを与えない、即ち多くの特徴量に対してペナルティを与えたからといって小さな重みに不要なアドバンテージを与えない、という点である。これはL0ノルム（これは $\sum \delta(w_i)$ で表される：ただし $\delta(x) = 1 (x \neq 0), 0 (x = 0)$ ）の良い近似であるとも言える。つまり、特徴量を使うか否かに対してのみ罰則項として働き、重みの大きさには関与しないということでもある。(2)式の $\epsilon$ が十分に小さければ、即ち $\epsilon \leq min_i | w_i |$ であれば、選択された特徴量の次元数は $q_{\epsilon} (\mathbf{w}) / \epsilon$ で正しく表される。言い換えると、 $q_{\epsilon} (\mathbf{w})$ でペナルティを与えるということは、そもそも残す特徴次元数そのものに対してペナルティを与えるということに非常に近似的である。ただし、このキャップつきL1ノルムは非凸なので最適化は容易ではない。

キャップつきL1ノルムは通常のL1 / L2ノルムと併用でき、適切な正則化係数 $\mu, \lambda \geq 0$ を付すことでトレードオフをコントロールできる。

$min_{\mathbf{w}} \displaystyle \sum_{(\mathbf{x}_i, y_i)} \ell (\mathbf{x}_i, y_i, \mathbf{w}) + \lambda |\mathbf{w}|_1 + \mu q_{\epsilon} (\mathbf{w})$ (3)

ここで $q_{\epsilon} (\mathbf{w})$ は $[ q_{\epsilon} (\mathbf{w}_1), \cdots, q_{\epsilon} (\mathbf{w}_d) ]$ である。

4 GRADIENT BOOSTED FEATURE SELECTION

(3)式それ自体はまだ線形のままなので、これを非線形にする。線形のものを非線形にするというとカーネル法やboostingが挙げられる。HSIC Lassoはカーネル法を用いているが、計算負荷が高い。そこでこの論文ではよりスケーラブルな方法ということでboostingを用いることにする。

Boostingそのものの説明は例えば朱鷺の森の解説記事でも読んでもらうとして、まず $\cal{H}$ を全ての可能な弱学習器としての回帰木の集合とする。もちろん精度にも木の数にも限りがあるので、 $|\cal{H}|$ は有限、ただしそれなりに大きな数とする。入力が $\phi (\mathbf{x}) = [ h_1 (\mathbf{x}), \cdots, h_{|\cal{H}|} (\mathbf{x}) ]^T$ なる $\cal{R}^{|\cal{H}|}$ にマッピングされているとして、ここでは(3)式の線形分類器を(4)式のように学習させる。

$min_{\mathbf{\beta}} \displaystyle \sum_{(\phi (\mathbf{x}_i), y_i)} \ell (\phi (\mathbf{x}_i), y_i, \mathbf{\beta}) + \lambda |\mathbf{\beta}|_1 + \mu q_{\epsilon}(\mathbf{\beta})$ (4)

ここで $\mathbf{\beta}$ は弱学習器の回帰木を選択するスパースなベクトルである。(4)式そのものはかなりの高次元になるが、 $\mathbf{\beta}$ がスパースなことを考えればその最適化は扱いやすい。一般性を失うことなく、集合 $\cal{H}$ における木はソートして最初の $\mathbf{\beta}$ に対する $T$ 個の入力が非ゼロだとすると、最終的な分類器は以下のようになる。

$H(\mathbf{x}) = \displaystyle \sum^T_{t=1} \beta_t h_t (\mathbf{x})$ (5)

Feature selection

(4)式は2つの罰則項を持つ。が、これはもはや重み付けではなく弱学習器としての木に対する正則化を行う式になっている。これがそのまま特徴量選択の意味合いを持つ。木のアンサンブルから選択される特徴量の個数をモデリングするために、ここで二値行列 $\mathbf{F} \in \{ 0, 1 \}^{d \times T}$ を定義する（ただし木 $h_t$ が特徴量 $f$ を用いる場合のみ $F_{ft} = 1$ とする）。これらを合わせると、特徴量 $f$ を選択するような木に対して与えられる全体の重み付けは(6)式のように表せる。

$\displaystyle \sum^T_{t=1} |F_{ft} \beta_t|$ (6)

ここで $q_{\epsilon}(\mathbf{\beta})$ を特徴量に紐付けられた実際の重みにペナルティを与えるように変形すると、最終的な最適化計画は以下のようになる。

$min_{\mathbf{\beta}} \ell (\mathbf{\beta}) + \lambda | \mathbf{\beta} |_1 + \mu \displaystyle \sum^d_{f=1} q_{\epsilon} (\displaystyle \sum^T_{t=1} |F_{ft} \beta_t| )$ (7)

既に述べたように、もし $\epsilon$ が十分に小さければ（即ち $\epsilon \leq min_f |\sum^T_{t=1} F_{ft} \beta_t|$ ）、 $\mu = 1 / \epsilon$ と置けて、なおかつ特徴選択のためのペナルティがそのまま用いられた特徴次元数と一致する。

4.1 Optimization

(7)式の最適化計画は、非凸な上に微分不可能である。けれどもここではgradient boostingを使って最適化してみる。まず $\cal{L}(\mathbf{\beta})$ を損失関数とし、 $\nabla \cal{L}(\mathbf{\beta})_t$ を $\beta_t$ に対する勾配とする。Gradient boostingはcoordinate descentと同様に捉えることができる。即ち毎ステップごとに最急勾配を持つ次元を更新していく、というやり方である。まず全ての弱学習器回帰木の集合 $\cal{H}$ をnegation closed、即ち各々の $h \in \cal{H}$ に対して $-h \in \cal{H}$ も存在するとする。これにより常に負の勾配だけを追いかけることができ、常に $\mathbf{\beta}$ は増大し続ける。これで必ず非負の最適な $\beta$ が存在することになる。次元 $t^{\ast}$ において最急負勾配を探索するという操作は(8)式のように表せる。

$t^{\ast} = argmin_t \nabla \cal{L}(\beta)_t$ (8)

以下、全ての可能な弱学習器回帰木に対してイテレートせずとも近似的に最小化できる戦略について示す。

$l_1$ -regularization

最適化における各ステップは固定されたステップサイズ $\alpha > 0$ の分だけ $\mathbf{\beta}$ の1つの次元に対して増えるため、 $\mathbf{\beta}$ のL1ノルムは $T$ 回のイテレーション後 $| \mathbf{\beta} |_1 = \alpha T$ となる。これは言い換えると、 $\mathbf{\beta}$ のL1ノルムにペナルティを加えるということは $T$ 回イテレーションでearly stoppingをかけるのと等価だということを意味する。そこで、以下 $\lambda | \mathbf{\beta} |_1$ の項は落とし、代わりに等価なハイパーパラメータとして $T$ を導入する。

Gradient Decomposition

あるイテレーション $T' + 1$ において最急勾配を見つけるためには、勾配を2つの要素に分解する必要がある。1つは損失関数 $\ell()$ に対するもの、もう1つはキャップつきL1ノルム罰則項に対するものである。

$\nabla \cal{L} (\mathbf{\beta})_t = \frac{\partial \ell}{\partial \beta_t} + \mu \displaystyle \sum^d_{f=1} \nabla q_{\epsilon} ( \displaystyle \sum^{\mathit{T}'}_{t=1} \mathit{F}_{ft} \beta_t )$ (9)

（ここで元の(7)式からは第2項後ろの絶対値記号を外してある：理由は簡単で $F_{ft}, \beta_t$ とも非負であるため）

$q_{\epsilon} ( \displaystyle \sum^{T'}_{t=1} F_{ft} \beta_t )$ の勾配は $\sum_t F_{ft} \beta_t = \epsilon$ なる点ではwell-definedではない。だが $\beta_t$ が減少することはないため、(10)式に見えるように右側極限値を取ることができる。

$\nabla q_{\epsilon} ( \displaystyle \sum^{T'}_{t=1} F_{ft} \beta_t ) = \begin{cases} F_{ft}, ~ if ~ \sum_t F_{ft} \beta_t < \epsilon \\ 0, ~ if ~ \sum_t F_{ft} \beta_t \geq \epsilon \end{cases}$ (10)

仮に $\epsilon = \alpha$ そして $\alpha > 0$ がステップサイズとすれば、 $\sum_t F_{ft} \beta_t > \epsilon$ は特徴量 $f$ が過去のイテレーションで既にある木の中で使われている場合にのみ成り立つ。ここで特徴量 $f$ がまだ使われていない場合に $\phi_f = 1$ 、そうでなければ $\phi_f = 0$ で表すとする。これを用いれば、(10)式の $\nabla q_{\epsilon} ( \sum^{T'}_{t=1} F_{ft} \beta_t )$ を場合分けせずに1つの勾配で表せる。

$\nabla \cal{L} ( \mathbf{\beta} )_t = \frac{\partial \ell}{\partial \beta_t} + \mu \displaystyle \sum^d_{f=1} \phi_f \mathit{F}_{ft}$ (11)

ただし $\phi_f F_{ft} = 1$ となるのは木 $t$ から最初に特徴量 $f$ が選択された時に限る。言い換えると、第2項は多くの新しい（過去に選択されていない）特徴量を使ってしまうような木に対して効果的にペナルティを与えるということでもある。

Greedy Tree construction

(11)式により、ようやくどの木に対しても勾配を計算できるようになった。ただし最適な $t^{\ast}$ を見つけるには依然として全ての木を探索するしかない。そこで、ここではその問題を既に定式化した損失関数を最小化させるベストの木 $h_t$ を探索する問題への転換を目指す。これは古典的な決定木のCARTアルゴリズムで近似できる。

そのためには(11)式の第1項をchain ruleでその手の問題になるように別々の微分の積同士の関係にバラす必要がある。その結果が(12)式である。ただし $H ( \mathbf{x}_i )$ はi番目の入力である。

$\nabla \cal{L} ( \mathbf{\beta} )_t = \displaystyle \sum^n_{i=1} \frac{\partial \ell}{\partial \mathit{H} ( \mathbf{x}_i )} \frac{\partial \mathit{H} ( \mathbf{x}_i )}{\partial \beta_t} + \mu \displaystyle \sum^d_{f=1} \phi_f \mathit{F}_{ft}$ (12)

ここで $H ( \mathbf{x}_i ) = \mathbf{\beta}^T \mathbf{h} ( \mathbf{x}_i )$ は全ての $h_t ( \mathbf{x}_i )$ のただの線形和で、全ての学習データに対する予測値であることに注意。よって $\frac{\partial \mathit{H} ( \mathbf{x}_i )}{\partial \beta_t} = h_t ( \mathbf{x}_i )$ である。仮に負の勾配 $g_i = - \frac{\partial \ell}{\partial \mathit{H} ( \mathbf{x}_i )}$ を導入すれば、(12)式はさらに(13)式のように最定式化できる。

$h_t = argmin_{h_t \in \cal{H}} \displaystyle \sum^n_{i=1} -g_i h_t ( \mathbf{x}_i ) + \mu \displaystyle \sum^d_{f=1} \phi_f \mathit{F}_{ft}$ (13)

なおここでは $\cal{H}$ は正規化された木に限った（即ち $\sum_i h^2_t ( \mathbf{x}_i ) = 1$ ）。この(13)式に2つの定数項、即ち $\frac{1}{2} \sum_i h^2_t ( \mathbf{x}_i ), \frac{1}{2} \sum_i g^2_i$ を加えれば、以下の式が完成する。

$h_t = argmin_{h_t \in \cal{H}} \frac{1}{2} \displaystyle \sum^n_{i=1} (g_i - h_t ( \mathbf{x}_i ))^2 + \mu \displaystyle \sum^d_{f=1} \phi_f \mathit{F}_{ft}$ (14)

これはまさにペナルティつき二乗損失そのもの、"impurity function"である。そしてその解は回帰木に対するgreedy CARTアルゴリズム（つまりxgboostの原型に当たるアルゴリズム）で効率的に導くことができる。(14)式の第1項は損失関数の負の勾配に最も合うように特徴量を分割し、第2項は過去のイテレーションで用いられた特徴量を切り分ける働きを担っている（筆者注：boostingとしての要素は第2項の $\phi_f$ が担っていて、過去に一度学習に使われた場合は以後使われないという操作を反映している）。以上のアルゴリズムを擬似コードに改めたものがAlgorithm 1である。

4.2 Structured Feature Selection

ここでは事前知識としてある程度のスパース性が特徴量に存在している場合にどうするかという話題がなされている（らしい）が、GBFSではその点にも対応できるという話がされている。この後のRESULTSパートに絡む話で言うと、事前に特徴量をある個数のbagにまとめられるケースに、できるだけ少ないbag（特徴量の特定のかたまり）のみで精度を出せるようにしたい、みたいな話への対応という話題。ということで詳細は面倒なので割愛。

5 RESULTS

検証に用いた環境はLinux 2.6.32.x86_64のオンプレのデスクトップで6コアIntel i7 (2.66 GHz)、96 GB RAMだそうです。

5.1 Synthetic data

これはもうFig.1を見たまんまで、こんな非線形のものを綺麗に分離するなんて非線形カーネル SVMでもないと無理では。。。というところをGBFSならできまぁす！というお話。

5.2 Structured feature selection

Princetonの"Colon"というデータセット（結腸がんのデータセット）に対して適用してみたもの。がん化した結腸が40サンプル、健常な結腸が22サンプルあって、6500のヒト遺伝子について測定がなされている。これは分子生物学的に9つのクラスタ(bag)に分けられるので、そのうちのどれが結腸がんにより強く関与するか？を特定したいという課題。そこでFig.2を見るとこれまた一目瞭然で、他の既存手法に比べてGBFSがずば抜けて綺麗に1つのbagに原因遺伝子を絞り込めているのが分かる。

5.3 Benchmark data sets

もう面倒なので図だけ貼っていく。詳細が気になる人は本文読んで下さい。

上記各データセットのholdoutに対するclassification error rates一覧。GBFS強いですねー。

Fig.4はKDD cup 1999でのデータに適応したもの。特徴量を絞り込めた上で尚且つerror rateを抑えられているのが分かる。Fig. 5はGBFSによって選んだ特長量でガウシアンカーネル SVMを学習させたらテストデータへのerror rateが一番低くなったよ的な主張。Fig.6はSMK-CAN-187という187行しかないのに19,993次元もあるデータに対するerror rateを見たもので、これもGBFSが良いスコアを叩き出しているのが見て取れる。

6,7,8 DISCUSSION etc.

果てしなく面倒なので全て割愛。興味のある人は本文が公開されているのでご自身でどぞー。ただ、ちょっと思ったのが「GBFSすげーって言ってもランダムフォレストで特徴量選択してもスコア大して変わらなくね？」。。。感想はご随意に。

9 Misc

DMLCのxgboostのGitHub見に行ったら、こんなissueが立ってました。

tomtung commented on 19 Sep

as described in http://alicezheng.org/papers/gbfs.pdf

Basically it seems like: when measuring the goodness of a tree node split (for learning a regression tree), add a tunable constant penalty if the feature used for splitting hasn't been used by existing trees.

Seems quite simple, and can be a nice alternative for people who use feature importance to select features.

お前がやれよ、とか言われそうｗ

読後感など

正則化についてはちろっと以前もブログで書いたことがありますが、要はカーネル法で非線形化すると重くて嫌だからboostingでやりました（baggingでやるというとランダムフォレストが先にあるので）という話なのかなと。むしろそこが出発点で、途中の誤差関数の変形のところはそこを遡ってやってるみたいな印象を受けました。

ちなみに輪読会の席で出た疑問としては「結局 $T$ をハイパーパラメータにしたのは筋悪では？」というもの。確かに、可能な全ての弱学習機回帰木 $\cal{H}$ のどこまでを網羅するかという範囲だけでコントロールするのは、やや面白くないかもという印象はあります。でも、かといってこれをどうしたらもっと綺麗になるかというとアイデアもなく。。。

そして、RESULTSのところを見ると予測精度面では結構ランダムフォレストの特徴量選択による結果とあんまり変わらない部分もあるというのがｗ　とは言えColonデータでの特徴量選択のシャープさには見事なものがあると思ったので、やっぱりxgboostに実装してもらって普及させる方が狙いとしては重要になるのかも？ってところですかねぇ。。。

概要

1 INTRODUCTION / 2 RELATED WORK

3 BACKGROUND

3.1 Feature selection with the norm

3.2 The capped norm

4 GRADIENT BOOSTED FEATURE SELECTION

Feature selection

4.1 Optimization

-regularization

Gradient Decomposition

Greedy Tree construction

4.2 Structured Feature Selection

5 RESULTS

5.1 Synthetic data

5.2 Structured feature selection

5.3 Benchmark data sets

6,7,8 DISCUSSION etc.

9 Misc

読後感など

3.1 Feature selection with the $l_1$ norm

3.2 The capped $l_1$ norm

$l_1$ -regularization