Blog

Research

シンプルなベイズ最適化について

このエントリーをはてなブックマークに追加

Black-box関数最適化

ベイズ最適化

みなさんこんにちは！ 6月からAI Labに異動してきた野村（@nmasahiro_）と申します．

先月まではSparkで集計基盤を作ったりといったデータエンジニアリング業務に従事していましたが，現在は数理最適化の研究を行っています．

今回は機械学習手法のハイパーパラメータ最適化によく利用されているベイズ最適化について勉強してみたので，自分なりの理解をまとめておこうと思います．

まだ勉強中の身ですので，もし気になる点などあればご指摘いただけると嬉しいです！

 

はじめに

Black-box関数最適化問題というのは，目的関数の代数的な表現，つまり，数式自体が与えられず，解に対する目的関数値のみを利用することができる問題です．

機械学習手法のハイパーパラメータ最適化などがBlack-box関数最適化問題の枠組みに属します．例としてDeep Learningのハイパーパラメータ最適化問題を考えた場合，バッチサイズ，ネットワーク構造（中間層の数，各層のユニット数，活性化関数の種類など），ドロップアウト率，最適化手法の種類，学習率などが最適化すべきハイパーパラメータに対応します．

Black-box関数最適化では解の1回の評価の計算コストが大きいことを前提としていることが特徴です．前述のDeep Learningの例を考えると，1回の評価に数時間かかることが容易に考えられると思います．よって，探索にかかる評価回数をいかに少なくできるかが重要となります．

なお，本記事では常に最小化問題を仮定します．

 

ベイズ最適化

ベイズ最適化（Bayesian Optimization）は有力なBlack-box関数最適化手法の1つであり，近年では機械学習手法のハイパーパラメータ最適化によく利用されています．ベイズ最適化は，Grid Search，Random Searchと比較して効率よく優れた解を求められることが報告されています[Snoek, 2015]． Googleはベイズ最適化でクッキーを作ったりしてます．楽しそう．

ベイズ最適化では評価に大きな計算コストのかかる目的関数の代わりにAcquisition functionと呼ばれる代理関数を最適化することで，次の探索点を選択します．このAcquisition functionを構築するために，Gaussian Processが使われることが多いです．ベイズ最適化には様々なvariantが存在しますが，本記事では，文献[Brochu, 2010]で説明されているシンプルなベイズ最適化のモデルを説明しようと思います．

 

全体のアルゴリズム

ベイズ最適化は，以下の操作を繰り返すことで探索を行います．

  • Gaussian Processによって予測した平均，分散を使って計算したAcquisition functionを最適化して，次の探索点を得る
  • 得られた探索点を評価し，評価値を得る
  • 探索点と評価値の組を保存し，Gaussian Processを更新する

つまり，ベイズ最適化では，ある解の”良さ”をAcquisition functionを通して評価し，そのAcquisition functionをGaussian Processによって更新するという訳です．それでは，このGaussian ProcessとAcquisition functionが何者なのかを見ていきます．

 

Gaussian Process

Gaussian Processは関数の空間で議論をすることから，正規分布の無限次元への拡張と言われることがあります．正規分布が確率変数上の分布であり，平均ベクトルと共分散行列により一意に定められるのに対し，Gaussian Processは関数上の分布であり，平均関数と共分散関数により定められます．

$$f({\bf x}) \sim \mathcal{GP}(m({\bf x}), k({\bf x}, {\bf x}^{\prime}))$$

直感的には，Gaussian Processは，ある\({\bf x}\)に対してスカラー\(f({\bf x})\)を返すのではなく，\({\bf x}\)における\(f\)の値\(f({\bf x})\)が従う正規分布（つまり，平均と分散）を返すイメージです．

共分散行列（カーネル関数）\(k\)としてよく使われるのがsquared exponential functionであり，以下の式で表されます．\({\bf x},{\bf x}^{\prime}\)が似ていると1,似てないと0に近くなるので，類似度を表現するようなものと考えることができます．カーネル関数にハイパーパラメータを導入したり，squared exponential function以外のカーネル関数を使ったりといったことをよくするのですが，ここでは省略します．ハイパーパラメータ最適化手法のハイパーパラメータを決定するのが難しい問題は根強いですね．．．

$$k({\bf x}, {\bf x}^{\prime}) = \exp \Bigl( -\frac{1}{2}\| {\bf x} – {\bf x}^{\prime} \|^2 \Bigr).$$

今，観測データ\(\mathcal{D}_{1:t} = \{ {\bf x}_{1:t}, y_{1:t} \}\)があるとします．求めたいのは\(f_{t+1} = f({\bf x}_{t+1})\)の予測分布\(P(f_{t+1}|\mathcal{D}_{1:t}, {\bf x}_{t+1})\)です．

まず，ガウス過程の性質より，\({\bf f}_{1:t}\)と\(f_{t+1}\)の同時分布は以下で与えられます．

$$P({\bf f}_{1:t+1}) = P({\bf f}_{1:t}, f_{t+1}) = \mathcal{N}\Biggl( {\bf 0}, \left[
\begin{array}{cc}
{\bf K} & {\bf k} \\
{\bf k}^{\rm T} & k({\bf x}_{t+1}, {\bf x}_{t+1})
\end{array}
\right] \Biggr)$$

ここで，
$${\bf k} = [k({\bf x}_{t+1}, {\bf x}_1), k({\bf x}_{t+1}, {\bf x}_2), \cdots, k({\bf x}_{t+1}, {\bf x}_t)]$$
です．

それでは，予測分布\(P(f_{t+1}|\mathcal{D}_{1:t}, {\bf x}_{t+1}) = \mathcal{N}(\mu_{{\bf x}_{t+1}|{\bf x}_{1:t}}, \sigma_{{\bf x}_{t+1}|{\bf x}_{1:t}}^2)\)の\(\mu, \sigma\)を求めていこうと思います．

ここで，以下の条件つき正規分布についての補題を証明なしで使います（詳細はPRML上巻のpp.82-85参照）

$$\mu_{{\bf a}|{\bf b}} = \mu_{{\bf a}} + \Sigma_{{\bf a} {\bf b}} \Sigma_{{\bf b} {\bf b}}^{-1}({\bf x}_{{\bf b}} – \mu_{{\bf b}}),\\
\Sigma_{{\bf a}|{\bf b}} = \Sigma_{{\bf a} {\bf a}} – \Sigma_{{\bf a} {\bf b}} \Sigma_{{\bf b} {\bf b}}^{-1} \Sigma_{{\bf b} {\bf a}}.$$

すると，\(\mu_{{\bf x}_{t+1}|{\bf x}_{1:t}}, \sigma_{{\bf x}_{t+1}|{\bf x}_{1:t}}^2\)は以下の式で求められます．

$$\mu_{{\bf x}_{t+1}|{\bf x}_{1:t}} = {\bf k}^{\rm T}{\bf K}^{-1}{\bf f}_{1:t}, \\
\sigma_{{\bf x}_{t+1}|{\bf x}_{1:t}}^2 = k({\bf x}_{t+1}, {\bf x}_{t+1}) – {\bf k}^{\rm T}{\bf K}^{-1}{\bf k}.$$

結果をまとめると，以下のようになります．

$$\begin{align}
P(f_{t+1}|\mathcal{D}_{1:t}, {\bf x}_{t+1}) &= \mathcal{N}(\mu_{{\bf x}_{t+1}|{\bf x}_{1:t}}, \sigma_{{\bf x}_{t+1}|{\bf x}_{1:t}}^2) \\
&= \mathcal{N}({\bf k}^{\rm T}{\bf K}^{-1}{\bf f}_{1:t}, k({\bf x}_{t+1}, {\bf x}_{t+1}) – {\bf k}^{\rm T}{\bf K}^{-1}{\bf k}).
\end{align}$$

 

サンプル数2のとき，3つ目の予測値はどうなるか

上の式だけ見ても，イマイチGaussian Processが予測しようとしているものがわからなかったので，サンプル数が2，つまり\({\bf x}_{1:2}\)が得られているとき，3つ目の解の評価値\(f({\bf x}_3)\)の平均\(\mu_{{\bf x}_{3} | {\bf x}_{1:2}}\)を具体的に計算してみます．

まず，\({\bf K}\)について，\(k({\bf x}_1, {\bf x}_1) = k({\bf x}_2, {\bf x}_2) = 1\)であり，
$$ {\bf K} = \left[
\begin{array}{cc}
1 & k({\bf x}_1, {\bf x}_2) \\
k({\bf x}_1, {\bf x}_2) & 1
\end{array}
\right]$$となるので，

$$\begin{align}
{\bf K}^{-1} = \frac{1}{1 – k({\bf x}_1, {\bf x}_2)^2} \left[
\begin{array}{cc}
1 & – k({\bf x}_1, {\bf x}_2) \\
– k({\bf x}_1, {\bf x}_2) & 1
\end{array}
\right].
\end{align}$$

よって，

$$\begin{align}
\label{eq:sample2:mu}
{\bf k}^{\rm T} {\bf K}^{-1} {\bf f}_{1:2} = \frac{1}{1 – k({\bf x}_1, {\bf x}_2)^2} & \Bigl\{ (k({\bf x}_3, {\bf x}_1) – k({\bf x}_3, {\bf x}_2) k({\bf x}_1, {\bf x}_2)) f({\bf x}_1) \nonumber \\
&+ (k({\bf x}_3, {\bf x}_2) – k({\bf x}_3, {\bf x}_1) k({\bf x}_1, {\bf x}_2)) f({\bf x}_2) \Bigr\}
\end{align}$$

上式から，以下のように考えることができます．

  • \(f({\bf x}_3)\)の平均である\(\mu_{{\bf x}_{3} | {\bf x}_{1:2}}\)は\(f({\bf x}_1), f({\bf x}_2)\)の線形和になっている．
  • カーネル関数について，近いほど大きい値を取ったことを思い出すと，\({\bf x}_3\)が\({\bf x}_1\)に近いほど\(f({\bf x}_1)\)を考慮し，\({\bf x}_3\)が\({\bf x}_2\)に近いほど\(f({\bf x}_2)\)を考慮しているように見える．

 

Acquisition function

ついにGaussian Processによって予測された評価値の平均，分散を使うときがやってきました．目的関数の代わりになるようなAcquisition functionを定義する必要があるのですが，それには以下のような性質を満たしていることが理想であると考えられます．

  • 良いかもしれないけどまだ探していないところに評価を与える（探索）
  • ある程度良いということがわかってるところに評価を与える（活用）

よく知られているAcquisition functionとして以下の3つがあるので，それぞれについて説明していきます．

  • Probability of Improvement (PI)
  • Expected Improvement (EI)
  • Lower Confidence Bound (LCB)

 

Probability of Improvement (PI)
PIでは改善確率を指標として用います． \(1:t\)でのベストな探索点を\({\bf x}^{-} = {\rm arg} \min_{{\bf x}_i \in {\bf x}_{1:t}}f({\bf x}_i)\)とすると，
新たな探索点\({\bf x}\)が\({\bf x}^{-}\)を改善する確率PI\(({\bf x})\)は以下で計算されます．

$$\begin{align}
{\rm PI}({\bf x}) &= P(f({\bf x}) \leq f({\bf x}^{-})) \\
&= \int_{-\infty}^{f({\bf x}^{-})} \mathcal{N}(f({\bf x})|\mu_{{\bf x}}, \sigma_{\bf x}^2) df \\
&= \Phi \Bigl( \frac{f({\bf x}^{-}) – \mu({\bf x})}{\sigma({\bf x})} \Bigr)
\end{align}$$

ここで\(\Phi\)は標準正規分布の累積分布関数(cdf)です．

しかし，この方法には注意点があります．それは，改善さえすればPIが良くなってしまうという点です．
つまり，無限小の改善であってもPIは良くなってしまうため，不確実だがもっといいかもしれないところより，
無限小の改善だけど改善はしそうというところに探索が集まってしまうのです．つまり，活用に振りすぎているということになります．
これを改善しようとして，trade-off parameter \(\xi \geq 0\)を用いる方法が提案されています．

$$\begin{align}
{\rm PI}({\bf x}) &= P(f({\bf x}) \leq f({\bf x}^{-} – \xi)) \\
&= \int_{-\infty}^{f({\bf x}^{-}) – \xi} \mathcal{N}(f({\bf x})|\mu_{\bf x}, \sigma_{\bf x}^2) \\
&= \Phi \Bigl( \frac{f({\bf x}^{-}) – \xi – \mu({\bf x})}{\sigma({\bf x})} \Bigr)
\end{align}$$

\(\xi\)はユーザパラメータなので決めるのは（個人的には）なかなか難しいと思ってます．先ほども少し書かせていただいた「ハイパーパラメータ最適化手法のハイパーパラメータ設定難しい説」ですね．

最初は\(\xi\)を大きくして探索をさせるようにし，だんだん\(\xi\)を小さくして活用に振らせるスケジューリングが経験的に良いと報告されているようです．

 

Expected Improvement (EI)
PIは改善確率だけを見ていたことにより，無限小の改善だけど改善はしそうというところが評価されやすくなるという問題点が存在しました．

EIは改善確率に加えて改善量（の期待値）も見るということがPIと異なります．以下で具体的に計算をしていきます．

 

改善量\({\rm I}({\bf x})\)を以下の式で定義します．

$$\begin{align}
{\rm I}({\bf x}) = \max \{ 0, f_{{\bf x}^{-}} – f_{t+1}({\bf x}) \}
\end{align}$$

ここで，maxオペレータにより\({\rm I}\)の取り得る値の範囲は\(0\)以上となるため，最終的に改善量\({\rm I}\)の期待値を計算する際には\(0\)以上の範囲のみを考えれば良いことになります．よって，\({\rm I}\)が従う分布についても\(0\)以上の範囲のみについて考えることにします．

すると，\(f_{t+1}({\bf x}) \sim \mathcal{N}(\mu({\bf x}), \sigma^2({\bf x}))\)より，
\({\rm I}({\bf x}) \sim \mathcal{N}(f({\bf x}^{-}) – \mu({\bf x}), \sigma^2({\bf x}))\)．
よって，\({\rm I}\)の尤度は，

$$\begin{align}
P(I) = \frac{1}{\sqrt{2\pi}\sigma({\bf x})} \exp \Bigl( -\frac{1}{2} \frac{({\rm I} – f({\bf x}^{-}) + \mu({\bf x}))^2}{\sigma^2({\bf x})} \Bigr)
\end{align}$$
となる．

よって，\({\rm I}\)の期待値は，
\begin{align}
\label{eq:ei:1}
\mathbb{E}(I) = \int_{{\rm I}=0}^{{\rm I}=\infty} {\rm I} \frac{1}{\sqrt{2\pi}\sigma({\bf x})} \exp \Bigl( -\frac{1}{2} \frac{({\rm I} – f({\bf x}^{-}) + \mu({\bf x}))^2}{\sigma^2({\bf x})} \Bigr) d{\rm I}
\end{align}

ここで，\(t=\frac{{\rm I} – f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})}\)とおくと，
\({\rm I}=0\)で\(t=\frac{- f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})}\)，\({\rm I}=\infty\)で\(t=\infty\)．
また，\(dt = \frac{1}{\sigma({\bf x})} d{\rm I}, {\rm I} = \sigma({\bf x})t + f({\bf x}^{-}) – \mu({\bf x})\)より，

$$\begin{align}
\mathbb{E}(I) &= \int_{t=\frac{- f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})}}^{t=\infty} (\sigma({\bf x})t + f({\bf x}^{-}) – \mu({\bf x})) \frac{1}{\sqrt{2\pi}\sigma({\bf x})} \exp \Bigl( -\frac{1}{2} t^2 \Bigr) \sigma({\bf x}) dt \\
\label{eq:ei:2}
&= \sigma({\bf x}) \int_{t=\frac{- f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})}}^{t=\infty} \frac{t}{\sqrt{2\pi}} \exp \Bigl( -\frac{1}{2} t^2 \Bigr) dt \\
\label{eq:ei:3}
& \ \ \ \ + (f({\bf x}^{-}) – \mu({\bf x})) \int_{t=\frac{- f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})}}^{t=\infty} \frac{1}{\sqrt{2\pi}} \exp \Bigl( -\frac{1}{2} t^2 \Bigr) dt
\end{align}$$

ここで，上式について，\(t^2 = u\)とおくと\(t dt = \frac{1}{2} du\)
よって，

$$\begin{align}
\sigma({\bf x}) \int_{t=\frac{- f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})}}^{t=\infty} \frac{t}{\sqrt{2\pi}} \exp \Bigl( -\frac{1}{2} t^2 \Bigr) dt &= \sigma({\bf x}) \int_{u=t^2}^{u=\infty} \frac{1}{\sqrt{2\pi}} \exp (-\frac{1}{2}u) \cdot \frac{1}{2} du \\
&= \frac{\sigma({\bf x})}{\sqrt{2\pi}} \Bigl[- \exp(-\frac{1}{2}u) \Bigr]_{u=t^2}^{u=\infty}\\
&= \frac{\sigma({\bf x})}{\sqrt{2\pi}} \Bigl[- \exp(-\frac{1}{2}\Bigl[ \frac{{\rm I} – f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})} \Bigr]^2 ) \Bigr]_{{\rm I}=0}^{{\rm I}=\infty} \\
&= \frac{\sigma({\bf x})}{\sqrt{2\pi}} \Bigl(0 + \exp(-\frac{1}{2} \Bigl[\frac{- f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})} \Bigr]^2) \Bigr) \\
&= \sigma({\bf x}) \cdot \frac{1}{\sqrt{2\pi}} \exp(-\frac{1}{2} \Bigl[\frac{- f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})} \Bigr]^2) \\
&= \sigma({\bf x}) \cdot \mathcal{N}(\frac{-f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})} | 0, 1) \\
&= \sigma({\bf x}) \cdot \phi \Bigl( \frac{-f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})} \Bigr)
= \sigma({\bf x}) \cdot \phi \Bigl( \frac{f({\bf x}^{-}) – \mu({\bf x})}{\sigma({\bf x})} \Bigr)
\end{align}$$

また，

$$\begin{align}
(f({\bf x}^{-}) – \mu({\bf x})) \int_{t=\frac{- f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})}}^{t=\infty} \frac{1}{\sqrt{2\pi}} \exp \Bigl( -\frac{1}{2} t^2 \Bigr) dt &= (f({\bf x}^{-}) – \mu({\bf x})) (1 – \int_{t=-\infty}^{t=\frac{- f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})}} \frac{1}{\sqrt{2\pi}} \exp \Bigl( -\frac{1}{2} t^2 \Bigr) dt) \\
&= (f({\bf x}^{-}) – \mu({\bf x})) \Bigl(1 – \Phi \Bigl( \frac{-f({\bf x}^{-}) + \mu({\bf x})}{\sigma({\bf x})} \Bigr) \Bigr) \\
&= (f({\bf x}^{-}) – \mu({\bf x})) \cdot \Phi \Bigl( \frac{f({\bf x}^{-}) – \mu({\bf x})}{\sigma({\bf x})} \Bigr)
\end{align}$$

まとめると，

$$\begin{align}
\mathbb{E}({\rm I}) &= \begin{cases}
\sigma({\bf x}) \cdot \phi (Z) + (f({\bf x}^{-}) – \mu({\bf x})) \cdot \Phi (Z) & {\rm if}\ \sigma({\bf x}) > 0 \\
0 & {\rm if}\ \sigma({\bf x}) = 0
\end{cases}\\
Z &= \frac{f({\bf x}^{-}) – \mu({\bf x})}{\sigma({\bf x})}
\end{align}$$

 

Upper Confidence Bound (UCB)
BanditでおなじみUCBです．ここで，\(\kappa ( \geq 0)\)はユーザパラメータです．

\(\mu({\bf x})\)が活用，\(\sigma({\bf x})\)が探索に対応する項であり，\(\kappa\)によって探索と活用のトレードオフを取っています．

 

\begin{align}
{\rm UCB}({\bf x}) = \mu({\bf x}) + \kappa \sigma({\bf x})
\end{align}

 

おわりに

本記事ではシンプルなベイズ最適化のモデルについて説明しました．個人的に，今回取り上げたベイズ最適化において改善できそう（あるいは改善したい）と感じた部分は以下になります．

  • サンプル数\(n\)に対して時間計算量が\(O(n^3)\)のため，収束性が重要な問題に対しても多数のサンプル生成が難しい．
  • 設定すべきパラメータが多い．例えば，カーネル関数（とそれに用いるハイパーパラメータ），Acquisition functionの最適化手法どれにするかなど，職人芸が必要な部分がある．
    • Acquisition functionはnon-convexであるため，局所探索に重きを置いている手法を複数回繰り返す，あるいは，大域的探索が可能な手法を使う必要があると考えられる．
  • 離散変数をうまく扱う枠組みがない．

次回は以下のようなことを行いたいと考えています．

  • 上記の課題点を解決できているベイズ最適化のvariantがあるかどうか調査
    • ベイズ最適化のサーベイ論文[Shahriari, 2016]読んでみる
    • GPをDNNに置き換えた[Snoek, 2015]は気になる
  • ハイパーパラメータに対してrobust, あるいはsensitiveな機械学習手法の洗い出し
  • 他のBlack-box関数最適化手法との比較実験

読んで下さり，ありがとうございました．

 

引用文献

  • E.Brochu, V.M.Cora, N.de Freitas : A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning, \(arXiv:1012.2599\) (2010).
  • B. Shahriari, K. Swersky, Z. Wang, R. Adams, and N. de Freitas : Taking the human out of the loop: A Review of Bayesian Optimization. \({\it Proc. \ of \ the \ IEEE, (1)}\), 12/2015(2016)
  • J.Snoek, H.Larochelle, R.P.Adams : Practical Bayesian Optimization of Machine Learning Algorithms, \({\it Advances in Neural Information Processing System 25}\), pp.2951-2959(2012).
  • J.Snoek, O.Rippel, K.Swersky, R.Kiros, N.Satish, N.Sundaram, M.M.A.Patwary, Prabhat, and R.P.Adams : Scalable Bayesian optimization using Deep Neural Networks, \({\it In \ Proc. \ of \ ICML’15}\) (2015)

Author

アバター
nomura