Optunaって便利だよね

ハイパーパラメータのチューニングの手法はいろいろとあるが、最近ではOptunaがよく使われる。
Optunaはベイズ最適化によって効率よくパラメータを探索してくれるため、グリッドサーチなどの他の手法に比べて高速だ。

と、ここまではよく聞く説明だけど
「ベイズ最適化って具体的に何をしてるの？」
と聞かれると詰まってしまう人は多いと思う（私です）。
この記事の目的はOptunaなどで使われるベイズ最適化の原理をちゃんと理解してすっきりすることだ。

ハイパーパラメータの探索

少し基本的なところからおさらいしよう。
機械学習の精度を上げるには適切な素性を加えて適切な学習をすることが大事だ。
そしてそれと並んでハイパーパラメータのチューニングも大事だ。
ハイパーパラメータの種類はモデルによって異なるが、学習する前に外部からあらかじめ設定する必要がある。
具体的には学習回数やモデルの複雑さを制御するパラメータなどだ。
適切なパラメータを見つけるころを「パラメータの探索」や「パラメータのチューニング」などという。
パラメータの探索には大きく分けて3種類の方法がある。

1. 手動で設定する

一番単純なのは手動でパラメータをいじって最適なパラメータを見つける方法だ。
簡単だが、組み合わせが増えると手動ですべてを探索するのは不可能なので、ある程度あたりをつけて調べる必要がある。
この方法で最適なパラメータを探索するにはモデルやパラメータに対する深い理解と経験（ヒューリスティクス）が必要なので、あまり実用的とは言えない。
Kaggleなどのデータコンペの実力者はこの方法でパラメータを設定する人もいるとか。

2. グリッドサーチなどの愚直な方法を使う

2つ目の方法はGrid Searchなどのナイーブな方法だ。
Grid Searchの仕組みは単純で、いわゆる総当たりでパラメータを探索する。
例えばパラメータA,Bがそれぞれ1-10の範囲を取りうるとすると、以下のようになる。

A,B = 
(1,1), (1,2) ... (1,10),
(2,1), (2,2) ... (2,10)
...
(10,1), (10,2)... (2,10)

こんな感じで全部のパターンにおいて最適化したいスコアを計算し、最もよかったときのパラメータを採用する。
格子(Grid)状に探索しているため　Grid Searchと呼ばれる。
弱点としては組み合わせが増えると時間がかかるという点だ。
n通りのパラメータとm通りのパラメータの組み合わせはn*m通りになる。
つまり、取りうるパターンの積になるので、変数が増えると組み合わせは爆発的に増えることになる。

他に乱数を用いて探索するRandomized Searchと呼ばれる方法もある。
乱数の分布をガウス分布などで指定できるので単純なGrid Searchよりも有効な場合はある。
しかし、結局場当たり的に探索していることには変わりない。
従って、パラメータが増えると探索に必要な回数は増えるし、逆に探索の回数が少なければ適切なパラメータにたどり着けない可能性がある。

3. Optunaなどのベイズ最適化を用いる方法

さて、ここで1と2のメリット、デメリットを今一度確認してみよう。
1のメリットは探索済みのパラメータを参照して、次のパラメータを決定できることだ。

例えば1-10の範囲を取るパラメータがあったとして、(2,4,6,8)を試したとしよう。
ここで2,4のスコアが高ければ
「お、最適はパラメータはこの辺にありそうだぞ」
とあたりをつけて3あたりを重点的に探索できる。
その反面、あたりの付け方が人間の判断に委ねられるので自動化はできない。

一方で2の手法は自動化できる反面探索範囲に融通が効かない。

ベイズ最適化を用いた方法はざっくり言えば両方の手法いいとこ取りをする。
つまり、探索範囲にあたりを付けつつ、自動で探索させることを目指す。

Expected Improvementの定式化

さて、なんとなくの方針がわかってきたところで、数式を定義していく。
inputとなるパラメータを$x$、outputとなるスコアを$y$とする。
$y$は小さいほど良いものとする。
僕たちは$y$を最小にするような$x$を見つけたい。
また、今までに$n$回探索を行ったとして、この結果をベースにしながら次の探索するパラメータを決定したい。
今までの$n$回分の結果を$D_n$と表すことにしよう。

ここで、事後確率を以下のように定義する。
$P\left(y \mid x, D_{n}\right)$
ここで、$P\left(y \mid x\right)$はxという条件のもとでのyが起きる確率密度を表す。
つまり、上の式は$D_{n}$　というデータがある時、あるxが与えられた時のyの確率密度を表している。

もう少し噛み砕いて説明する。
本来であればパラメータxが定まればスコアyも一意に定まるはずなので、確率分布を考えること自体がおかしなことにも見える。
しかし、ここではパラメータxが与えられた時にyがどうなるかは未知であるため、今までの試行結果$D_n$をベースとして「yはこの範囲の値をとりそう！」というのを確率密度で表してあげるわけだ。
yの確率分布なのでyで積分すれば1になる。

さて、これまでに出てきた式を使ってExpected Improvementと呼ばれる量を定義してみよう。

$E I_{D_{n}}(x)=\int_{-\infty}^{\infty} \max \left(y^{*}-y, 0\right) P\left(y \mid x, D_{n}\right) d y$

ここで$y^{*}$はベースラインとなる$y$の値とする。
具体的には得られた$n$個のyの上位$γ$%のスコアを$y^{*}$と定義する。
$γ$が上位５０パーセントであれば、今まで得られたデータに対して上位半分に入るかどうかの基準になるイメージだ。

$y$は小さい方がいいので$y^{*}-y$は大きいほど嬉しい。
この値が正になれば改善していることになるが、負になればむしろ悪化していることになる。
つまり、上の式の$\max(y^{*}-y, 0)$ はどれくらいスコアを改善できたかを表している。

さて、先ほども触れたように、パラメータxが与えられてもyの取りうる値は確率分布で表されるため、
改善の期待値を計算する。
期待値は確率密度をかけ算して積分すればいいので上の式のようになるというわけだ。

ベイズの定理を導入する

Expected Improvementを改善すればいいことは分かったが、
$P\left(y \mid x, D_{n}\right)$ をどうやって求めればいいのだろうか。

ここでベイズの定理を考える。
一般にAとBという事象が同時に起きる確率は$P(A)P(B)$と積で表すことができる。
これは (Aが起きる確率)*(Aが起こった時にBが起きる確率) = (Bが起きる確率)*(Bが起こった時にAが起きる確率)が等しいことを表しているので、以下のように書ける。

$P(A)P(B|A) = P(B)P(A|B)$
これを変形すると
$P(B|A) = \frac{P(B)P(A|B)}{P(A)}$
と書ける。これがベイズの定理だ。
(Aが起きたときにBが起こる確率)は(Bが起きたときにAが起きる確率)を用いて表すことができるということを表している。
ここで$P(B|A)$を事後確率、$P(A|B)$を事前確率と呼ぶ。

これを今回の式に当てはめてみよう。

$P\left(y \mid x, D_{n}\right)=\frac{P\left(x \mid y, D_{n}\right) P\left(y \mid D_{n}\right)}{P\left(x \mid D_{n}\right)}
$

今回求めたいのが左辺の事後確率だ。
ベイズの定理を使えばこれを事前確率$P(x|y)$から計算することができる。

ベイズの定理の解釈

ここで少し事前確率について少し考えてみる。
事後確率の式はパラメータ$x$が決まった時にスコアが$y$になる確率を示している。
これは本来の(パラメータを決める) → (スコアが決まる)という因果関係を表している。
一方で今回はn回の探索で得られた結果yから最適なxを推定したいので、向きが逆になる。

つまり、ベイズの定理を使うことでこれらを相互に変換することができるというわけだ。

Tree-Structured Parzen Estimator (TPE)

話を戻し、事前確率の計算方法を考える。
OptunaではTPEと呼ばれるベイズの定理を用いた方法を使って事前分布を計算している。
事前確率は以下の式だった。
$P\left(x \mid y, D_{n}\right)$

これはyを決めた時のxの確率分布を示している。
今回は得られたyが$y^{*}$より大きいかどうかがポイントとなるため、それぞれで場合分けした時のxの分布を考え、それぞれ$l(x),g(x)$ で表す。

P\left(x \mid y, D_{n}\right)=\left\{\begin{array}{ll}
l\left(x \mid D_{n}\right) & \text { if } y<y^{*} \\
g\left(x \mid D_{n}\right) & \text { if } y \geq y^{*}
\end{array}\right.

以下の図は[1]から引用した概念図だ。
n回のサンプルで得られた結果を$y^{*}$以上と以下で分割する。

$l(x), g(x)$はそれぞれの分布にあたるため、今回の例だと以下のようになる[1]。
なお、実際には得られるのはデータ点の集合であり連続的な関数ではないため、データ点から関数を推定する必要がある。
ここでは詳細は省くが、Parzen Estimator (カーネル密度推定法)という方法を用いることで推定することができる。

ここで先程の事前分布の式を$l(x), g(x)$を使って書き下してみる。


P\left(y \mid x, D_{n}\right)=
\frac{P\left(x \mid y, D_{n}\right) P\left(y \mid D_{n}\right)}{P\left(x \mid D_{n}\right)}　\\


=\left\{\begin{array}{cl}
\frac{l\left(x \mid D_{n}\right) \cdot P\left(y \mid D_{n}\right)}{\gamma l\left(x \mid D_{n}\right)+(1-\gamma) g\left(x \mid D_{n}\right)} & \text { if } y<y^{*} \\
\frac{g\left(x \mid D_{n}\right) \cdot P\left(y \mid D_{n}\right)}{\gamma l\left(x \mid D_{n}\right)+(1-\gamma) g\left(x \mid D_{n}\right)} & \text { if } y \geq y^{*}
\end{array}\right.

事前分布$P(x|y)$の部分がそれぞれ$l(x), g(x)$で置き換えられ、分母もγを使って$l(x), g(x)$の和で表せる。

さて、これらを使って僕たちが最大化したいExpected Improvementを計算してみよう。

E I_{D_{n}}(x)=
\int_{-\infty}^{\infty} \max \left(y^{*}-y, 0\right) P\left(y \mid x, D_{n}\right) d y \\
=\int_{-\infty}^{y^{*}} \left(y^{*}-y\right) P\left(y \mid x, D_{n}\right) d y \\

=\int_{-\infty}^{y^{*}} \left(y^{*}-y\right) \frac{l\left(x \mid D_{n}\right) \cdot P\left(y \mid D_{n}\right)}{\gamma l\left(x \mid D_{n}\right)+(1-\gamma) g\left(x \mid D_{n}\right)}  d y \\

=\left(\gamma+\frac{g\left(x \mid D_{n}\right)}{l\left(x \mid D_{n}\right)}(1-\gamma)\right)^{-1}
\int_{-\infty}^{y^{*}} \left(y^{*}-y\right) P\left(y \mid D_{n}\right)  d y \\


=\left(\gamma+\frac{g\left(x \mid D_{n}\right)}{l\left(x \mid D_{n}\right)}(1-\gamma)\right)^{-1}\left\{\gamma y^{*}-\int_{-\infty}^{y^{*}} y P\left(y \mid D_{n}\right) d y\right\}

やや複雑な式変形なので丁寧に書いた[2]。
1行目から2行目にかけてはmax関数が0でないのは$y<y^{*}$のときのみなのを考慮して変形。
2行目から3行目はベイズの定理を用いた先程の式を代入している。
3行目から4行目でyの積分の定数部分をくくっている。
4行目から5行目では$P(y)$の$y<y^{*}$での積分が$\gamma$になることを使っている。

ベイズ最適化の方法

さて、いろいろと式変形をしてきたわけだが、この式で重要なのが$\frac{g(x)}{l(x)}$の部分だ。
式全体を最大化するためには$\frac{g(x)}{l(x)}$が小さいほうが良い。

つまり、それまでに得られた結果をもとに$\frac{g(x)}{l(x)}$が小さくなるような点を次の探索候補し、それを繰り返すことで最適なパラメータを見つけることができるというわけだ。

以下の図のようなイメージ[1]。

そもそもyが基準よりも良いスコアの時の分布が$l(x)$で悪い時の分布が$g(x)$なので、$l(x)$が$g(x)$に比べて大きくなるように探索していくとイメージすることもできる。
直感的にも正そうだ。

Optunaの場合

Optunaの場合はObjectiveクラスの中で最適化したい関数の定義やパラメータの範囲を指定する。
これはこの記事におけるxやyの定義にあたる。
そして、study = optuna.create_study()でそれまでの結果を参照しながら最適なパラメータを探索していく。
study.best_paramsで最適なパラメータを得ることができる。
詳しい実装方法はこの記事の範疇を超えるのでこれくらいで。

まとめ

Optunaの使い方の記事はいっぱいあるが、原理の詳細に突っ込んだ記事はあまりなかったのでまとめてみた
これでOptunaについてより深い理解ができるかも？

参考文献

[1] Kaggleで勝つデータ分析の技術
説明及び図を参考にした
[2] https://www.slideshare.net/hskksk/hyperopt
式変形の参考にしたスライド