Optuna(TPE)のアルゴリズム理解 ー Part 3 ー

はじめに

part1,part2からの続編になります．

2.獲得関数

獲得関数とはある点に対して，探索時に良い評価値を見つける可能性に応じて値をつける関数のこと．GPではこの獲得関数の設定によって探索性能が大きく変わる．

TPEでは論文にあるとおり，獲得関数を以下のように定義している．

EI_{y^*}(x) = \int_{-\infty}^{y^*} (y^* - y)p(y|x) dy = \int_{-\infty}^{y^*} (y^* - y)\frac{p(x|y)p(y)}{p(x)} dy

最後の式変形はベイズの定理．意味としてはある $x$ を与えたときに評価値 $y$ の値が 閾値 $y^{*}$ よりもどの程度良くなるか，という期待値．獲得関数に関してはこちらのページの34~39あたりが参考になります．

これをさらに式変形を行うことで最終的に $EI_{y^{*}}(x)$ の $x$ のみに依存する項は $\frac{l(x)}{g(x)}$ に比例することがわかる．そのため，TPEでは上で求めた確率密度関数を用いて $\frac{l(x)}{g(x)}$ を最大化するように探索を行う．

3.サンプリング

最後にサンプリング．GPのように目的関数を直接モデル化し，空間上のあらゆる点の予測値がわかる方法は最強だと思われがちだが，次元が大きくなると予測しなければいけない点の数自体が膨大になり，目的関数をモデル化してもそのモデル内でどこが最高の獲得関数を取るかわからない(正確には全てを予想するだけの時間がないだけであり，時間を無限にかけることで実行可能なのですかね？)．そのため，GPではモデル内の探索自体にCMA-ESを用いている．

TPEでは以下のようにサンプリングを行っている．
参考はoptuna/samplers/random.py と optuna/samplers/tpe/sampler.py です．

以下では，正規分布関数の累積分布関数を

F(x)=\int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}\sigma}\exp\Biggl(-\frac{(x'-\mu)^2}{2\sigma^2}\Biggl) dx' 
= \frac{1}{2}\Biggl(1+{\rm erf }\ \frac{x-\mu}{\sqrt{2}\sigma} \Biggl)

という形で導出し，表記する．また，erfはガウスの誤差関数である．(対数正規分布も似たような形になる．)

def F(x, mu, sigma):
  import scipy.special

  idx = (x - mu) / np.max(np.sqrt(2) * sigma, 1e-12)
  return  0.5 * ( 1 + scipy.special.erf(idx) )

####[1]連続値と離散値のハイパパラメータ
#####1.カーネル密度推定

集合 $L$ とハイパパラメータの定義域，集合 $G$ とハイパパラメータの定義域を事前に $ParzenEstimator()$ に代入することでカーネル密度推定を利用した $L$, $G$ それぞれの確率密度関数のパラメータ
$w^L_i=\{w^l_{i,1},w^l_{i,2},...\},\ w^G_i=\{w^g_{i,1},w^g_{i,1},...\}$
$\sigma^L_i=\{\sigma^l_{i,1},\sigma^l_{i,2},...\},\ \sigma^G_i=\{\sigma^g_{i,1},\sigma^g_{i,2},...\}$
を推定(というよりは経験則に基づき決定)．

ただし，集合 $L, G$ の $i$ 番目のハイパパラメータ要素をそれぞれ $x^L_i=\{x^l_{i,1},x^l_{i,2},...\},\ x^G_i=\{x^g_{i,1},x^g_{i,2},...\}$ とする．
このとき
$l_i(x_i) = \sum_{j} w^l_{i,j} K(x_i|x^l_{i,j},\sigma^l_{i,j})$
$g_i(x_i) = \sum_{j} w^g_{i,j} K(x_i|x^g_{i,j},\sigma^g_{i,j})$

ただし，カーネル関数は
$K(x|\mu,\sigma) = \frac{1}{\sqrt{2\pi}\sigma}\exp(-\frac{(x-\mu)^2}{2\sigma^2})$
対数スケールに変換している場合は
$K(x|\mu,\sigma) = \frac{1}{\sqrt{2\pi}\sigma x}\exp(-\frac{(\ln x-\mu)^2}{2\sigma^2})$
としている．

#####2.サンプリング
KDEの結果を $\_sample\_from\_gmm()$ に代入．

$\pi^L_i=\{\pi^l_{i,1},\pi^l_{i,2},...\}$
$\pi^l_{i,j} = \frac{w^l_{i,j}}{\sum_{k} w^l_{i,k}}\ ,\ \sum_{j} \pi^l_{i,j} = 1$
を生成する．
上の $\pi^l_{i,j}$ という確率で $j$ が生成される乱数 $r$ を生成し，その乱数 $r$ に対応する観測値 $x^l_{i,r}$ と バンド幅$\sigma^l_{i,r}$ をハイパパラメータとして持つ正規分布に従う乱数 $r_{norm}$ を1つ生成する．つまり，$r_{norm} 〜 N(x^l_{i,r},\bigl(\sigma^l_{i,r})^2\bigl)$なる値を生成する(対数スケールなら対数を取ったものに対して正規乱数を発生させる)．$r_{norm}$が定義域内に収まる場合サンプリング集合 $S_i$ に追加する．

この操作を24個(デフォルト値)のサンプル値を得るまで繰り返し，サンプリング集合 $S_i=\{s_{i,1}, s_{i,2},...,s_{i,24}\}$ を得る．

#####3.サンプリングに基づいた対数尤度算出
######A.連続値の場合
ここでは直接サンプル値を引数として取るときの確率密度関数値を求めに行っている．また，$i$ はハイパパラメータ番号，$j$ はサンプル番号，$k$ は観測変数に関する番号を指す．

規格化定数

a_i^l=\ln \Biggl( \int_{x_{i,min}}^{x_{i,max}} \sum_{j}  w^l_{i,j} K(x_i|,x^l_{i,j},\sigma^l_{i,j}) dx_i \Biggl)

a_i^g=\ln \Biggl( \int_{x_{i,min}}^{x_{i,max}} \sum_{j}  w^g_{i,j} K(x_i|,x^g_{i,j},\sigma^g_{i,j}) dx_i \Biggl)

mahalanobis距離
$D^l_{i,j,k}=\Bigl( \frac{s_{i,j}-x^l_{i,k}}{\sigma^l_{i,k}}\Bigl)^2$
$D^g_{i,j,k}=\Bigl( \frac{s_{i,j}-x^g_{i,k}}{\sigma^g_{i,k}}\Bigl)^2$
jacobian行列(ここでは対数スケール変換なので単純にsample値で割るだけ．ちなみにこの計算はHyperoptでは存在しない．)
通常スケールの場合
$J_{i,j,k}=1$
対数スケールの場合
$J_{i,j,k}=\frac{1}{s_{i,j}}$

対数確率密度関数値

\ln l_{i}(s_{i,j})

=\ln \Biggl( \sum_{k}\exp\Bigl(-\frac{D^l_{i,j,k}}{2} + \ln \frac{w^l_{i,k} J_{i,j,k}}{\sqrt{2\pi}a_i^l\sigma^l_{i,k}}\Bigl) \Biggl)

=\ln \Biggl( \sum_{k}\frac{w^l_{i,k} J_{i,j,k}}{\sqrt{2\pi}a_i^l\sigma^l_{i,k}}\exp\Bigl(-\frac{D^l_{i,j,k}}{2} \Bigl) \Biggl)

\ln g_{i}(s_{i,j})

=\ln \Biggl( \sum_{k}\exp\Bigl(-\frac{D^g_{i,j,k}}{2} + \ln \frac{w^g_{i,k} J_{i,j,k}}{\sqrt{2\pi}a_i^g\sigma^g_{i,k}}\Bigl) \Biggl)

=\ln \Biggl( \sum_{k}\frac{w^g_{i,k} J_{i,j,k}}{\sqrt{2\pi}a_i^g\sigma^g_{i,k}}\exp\Bigl(-\frac{D^g_{i,j,k}}{2} \Bigl) \Biggl)

上の式を求めるために以下のclassmethodを利用するが，関数内に定義されているmは実際に計算を行うと相殺されるので計算上は必要なし．つまり，"return np.log(np.exp(x).sum(axis=1))"を計算しているのと同値．

@classmethod
def _logsum_rows(cls, x):
  # type: (np.ndarray) -> np.ndarray

  x = np.asarray(x)
  m = x.max(axis=1)
  return np.log(np.exp(x - m[:, None]).sum(axis=1)) + m

######B.離散値の場合
確率密度関数の定義域内の積分値は1ではないため，実際にはKDEによる確率密度モデルを定義域全域で積分した値によって割ることによって確率密度関数を規格化している．また，離散値では離散値自体の増加方向，減少方向に対して意味を持つので，カーネルの重なり部分を考慮するために一つの離散値間の距離(ここでは $q$ と仮定している．整数値の場合は $q=1$ となる．)を用いて，実際に離散値が取る値 $\pm \frac{q}{2}$の区間で積分をし，均している(?)．

ここで求めている値は結果的に以下の値となる．

\ln l_{i}(s_{i,j})

=\ln \Biggl( \int_{s_{i,j}-\frac{q}{2}}^{s_{i,j}+\frac{q}{2}} \sum_{j}  w^l_{i,j} K(x_i|,x^l_{i,j},\sigma^l_{i,j}) dx_i \Biggl)
- \ln \Biggl( \int_{x_{i,min}}^{x_{i,max}} \sum_{j}  w^l_{i,j} K(x_i|,x^l_{i,j},\sigma^l_{i,j}) dx_i \Biggl)

\ln g_{i}(s_{i,j})

=\ln \Biggl( \int_{s_{i,j}-\frac{q}{2}}^{s_{i,j}+\frac{q}{2}} \sum_{j}  w^l_{i,j} K(x_i|,x^g_{i,j},\sigma^g_{i,j}) dx_i \Biggl)
- \ln \Biggl( \int_{x_{i,min}}^{x_{i,max}} \sum_{j}  w^g_{i,j} K(x_i|,x^g_{i,j},\sigma^g_{i,j}) dx_i \Biggl)

\ln EI_{y^*,i} \propto \ln \Biggl(\frac{l_{i}(x_i)}{g_{i}(x_i)} \Biggl) = \ln l_{i}(x_i) - \ln g_{i}(x_i)

####[2]カテゴリカル変数
#####1.確率密度関数推定
カテゴリカル変数は数字の増減に対して幾何的な意味を持たないため，単純な重み付きヒストグラムによって確率密度関数を推定する．

集合 $L, G$ はそれぞれ $x^L_i=\{x^l_{i,1},x^l_{i,2},...\},\ x^G_i=\{x^g_{i,1},x^g_{i,2},...\}$ とする．また，各観測点に対する重みを$w^L_i=\{w^l_{i,1},w^l_{i,2},...\},\ w^G_i=\{w^g_{i,1},w^g_{i,1},...\}$
とする．

最初に集合 $L, G$ で各カテゴリーが何回ずつ出たか重み付きで数え上げる．(以下に例を示す．)

例：$x_i=\{1,6,5,3,3,5,2,2,3,3\},\ w_i=\{1,1,1,1,1,1,1,1,1,1\}$
$カテゴリkの重み付き出現回数=\sum_{j \ for \ all \ x_{i,j}=k} w_{i,j}$

カテゴリ種類	カテゴリ1	カテゴリ2	カテゴリ3	カテゴリ4	カテゴリ5	カテゴリ6	試行回数
出現回数	1	2	4	0	2	1	10

次に各カテゴリーの出現確率が0になることを避けるため，各カテゴリーの出現回数に1を加える．

カテゴリ種類	カテゴリ1	カテゴリ2	カテゴリ3	カテゴリ4	カテゴリ5	カテゴリ6	合計
出現回数 + 1	2	3	5	1	3	2	16

そして，全てのカテゴリーの総出現回数で割ることによって，規格化する．

カテゴリ種類	カテゴリ1	カテゴリ2	カテゴリ3	カテゴリ4	カテゴリ5	カテゴリ6	合計
出現確率	$\frac{1}{8}$	$\frac{3}{16}$	$\frac{5}{16}$	$\frac{1}{16}$	$\frac{3}{16}$	$\frac{1}{8}$	1

#####2.サンプリング
上の操作によって得た重み付きヒストグラム(一番最後の表のように各カテゴリの出現確率を導出したもの)に従った疑似乱数を既定の個数生成し，サンプル点とする．

#####3.サンプリングに基づいた対数尤度算出
サンプリングによって得た点に対応する重み付きヒストグラムの高さ(つまり，最後の表の出現確率)がこのサンプル点が正しいクラス($L$ or $G$)に分類される尤度であるから，これらに対して対数を取ることによって対数尤度を得る．

###サンプル点の採択
最終的に上で得た各サンプル点の対数獲得関数値の最良値を取ったサンプル点を採択し，これを次回の評価点とする．

以下のコードに関しても結果的に引数を与えた状態で呼び出し元に返しているので，実は"return np.asarray([samples[best]] * samples.size)"の"* samples.size"は要らないのではないかと思っています．

@classmethod
def _compare(cls, samples, log_l, log_g):
# type: (np.ndarray, np.ndarray, np.ndarray) -> np.ndarray

  samples, log_l, log_g = map(np.asarray, (samples, log_l, log_g))
  if samples.size:
    score = log_l - log_g
    if samples.size != score.size:
      raise ValueError("The size of the 'samples' and that of the 'score should be same. But (samples.size, score.size) = ({}, {})".format(samples.size, score.size))

    best = np.argmax(score)
    return np.asarray([samples[best]] * samples.size)
  else:
    return np.asarray([])