\newcommand{\argmax}{\mathop{\rm arg~max}\limits}
\newcommand{\argmin}{\mathop{\rm arg~min}\limits}

ガウス過程とベイズ最適化について

本記事では, optuna などのハイパーパラメータ最適化アルゴリズムで用いられる
ガウス過程とベイズ最適化 について説明します.

キーワード

確率過程
ガウス過程
ベイズ最適化
ハイパーパラメータ最適化

参考書

機械学習プロフェッショナルシリーズ
0 章に概要が書いてあり, 無料で読める. ページも 11 ページと少ない.

確率過程とは

時間などの条件によって変化する確率. 為替や株価など.
関数自体が確率によって変化 することを考える.
イメージ図

ガウス過程とは

各時刻の値の確率分布が, 正規分布で表される確率過程.
なぜ正規分布を使うのか?
→ どのような確率分布も, 平均をとってしまえば正規分布に近づくため (中心極限定理).
中心極限定理
ある母集団から $N$ 個のデータをサンプリングし, それらの平均値を確率変数 $Y$ とする.
$N \rightarrow \infty$ に近づくにつれ, 確率変数 $Y$ (を標準化したもの) は標準正規分布に近づく.

例: サイコロを $N \times M$ 回振り, $N$ 個を1グループとして平均値$Y$を求める. $Y \in \mathbb{R}^M$ をヒストグラム化する.
X = np.random.randint(1, 6 + 1, size=(n, M)) X = np.mean(X, axis=1) Y, _ = np.histogram(X, bins=6) Y = np.divide(Y, np.sum(Y)) # [0, 1] に標準化 # Yをプロット ...
はじめは一様分布の形だが, $n → \infty$ に近づくにつれ, だんだん正規分布の形に近づくことがわかる.
回帰問題
- 連続値を推定する問題.
- ↔ クラス分類 (離散値を推定する問題)
ガウス過程回帰
- ガウス過程を使って解く回帰問題

ガウス過程回帰とは

結論

カーネル関数 $k$, 既知の ${ (x_1, y_1), ..., (x_N, y_N) }$ が与えられたとき,
未知の $x^*$ についての $y^*$ の確率分布 $\mathcal{N}(\hat{\mu_{y^*}}, \hat{\Sigma_{y^*}})$ は

\begin{aligned}
     \\
    \hat{\mu_{y^*}}    &= k_*^T K^{-1} y \\
    \hat{\Sigma_{y^*}} &= k_{**} - k_*^T K^{-1} k_*
\end{aligned}

ただし $K$ は $N \times N$ の行列で, $K_{ij} = k(x_i, x_j)$

導出

訓練データ $\mathcal{D}$ とカーネル関数 $k(x_i, x_j)$ が与えられる. ただし, $y$ は簡単化のため $0$ とする.
```
\mathcal{D} = \{(x_1, y_1), (x_2, y_2), ..., (x_N, y_N)\}
```
※ $y$ の観測値にノイズ $\varepsilon \sim \mathcal{N}(0, \sigma^2)$ が含まれている場合を考慮して, カーネル関数 $k$ は観測ノイズを含めて考える.
```
k(x_i, x_j) = (x_i と x_j の距離) + \sigma^2 \delta(x_i, x_j)
```
```
\delta(x_i, x_j) =
\begin{cases}
1 &\text{if } x_i = x_j \\
0 &\text{if }  x_i \neq x_j
\end{cases}
```

このとき, ${y_1, ..., y_N}$ は, カーネル行列 $K$ を分散共分散行列にとる多変量正規分布から生成されたと考える.

\begin{aligned}
    y &= \mathcal{N}(0, K) \\
    K &=
        \begin{bmatrix}
            k(x_1, x_1) & \dots  & k(x_1, x_N) \\
            \vdots      & \ddots & \vdots \\
            k(x_N, x_1) & \dots  & k(x_N, x_N) \\
        \end{bmatrix}
\end{aligned}

未知のデータ $x^*, y^*$ もガウス過程に従うので, 以下のようにも表せる.

\begin{aligned}
    y' = [y_1 \dots y_N y^*]^T &= \mathcal{N}(0, K') \\
    K' &=
        \begin{bmatrix}
            \boxed{K}     & \boxed{k_{*}} \\
            \boxed{k_{*}^T} & \boxed{k_{**}}
        \end{bmatrix} \\
    k_{*} &= [k(x_1, x^*) \dots k(x_N, x^*)]^T \\
    k_{**} &= k(x^*, x^*) \\
\end{aligned}

ここで, $y$ が与えられたときの $y*$ の条件付き確率 $p(y^*|y, x)$ を考えると,

\begin{aligned}
    p(y^*|y, x)  &= \mathcal{N}(\mu_{y^*}, \Sigma_{y^*}) \\
    \mu_{y^*}    &= k_*^T K^{-1} y \\
    \Sigma_{y^*} &= k_{**} - k_*^T K^{-1} k_*
\end{aligned}

が求まる.

プログラム

問題

$x_\text{train}$, $y_\text{train}$, $x_\text{test}$ が与えられる.
- $x_\text{train}$, $x_\text{test}$ はベクトルの配列
- $y_\text{train}$ はスカラーの配列
$x_\text{test}$ の各ベクトルに対応する出力 $y_\text{test}$ の確率分布を求めたい.
- それぞれの $x_i \in x_\text{test}$ において, $y(x_i) \sim \mathcal{N}(\square, \triangle)$ の $\square, \triangle$ はいくらと推定できるだろうか?

アルゴリズム

カーネル関数 $k(\cdot, \cdot)$ を設計
$x_\text{train}$ からカーネル行列$K$ を計算
$K[i,j] \leftarrow k(x_i, x_j)$
$x_\text{test}$ の各要素 $x_m$ について繰り返し:
1. $x_\text{train}$, $x_m$ から $k_*$ を計算
  $k_*[i] \leftarrow k(x_i, x_m)$
2. $x_m$ から $k_{**}$ を計算
  $k_{**} \leftarrow k(x_m, x_m)$
3. $K$, $k_*$, $k_{**}$ から $y_m$ の平均, 分散を計算
  - $\mu_{y_m} \leftarrow {k_*} {K^{-1}} {y_\text{train}}$
  - $\sigma^2_{y_m} \leftarrow {k_{**}}^2 - {k_*} {K^{-1}} {{k_*}^T}$

実行例

適当に設定した関数 $f(x)$ からトレーニングデータとしていくつかサンプリングし, ${x_\text{train}}$, ${y_\text{train}}$ とする.
$x^\ast \in [x_\text{min}, x_\text{max}]$ についての $y^\ast$ の確率分布をグラフにプロットする.

カーネル関数の設計

\text{k}(x_i, x_j|\theta_1,\theta_2,\theta_3) := {\theta_1} {\exp[-\frac{(x_i - x_j)^2}{\theta_2}]} + {\theta_3 \delta(x_i, x_j)}

\delta(x_i, x_j) =
    \begin{cases}
        1 &\text{if } x_i = x_j \\
        0 &\text{if } x_i \neq x_j
    \end{cases}

このカーネルは, 一般的に RBF カーネル と呼ばれる.
- SVM (サポートベクターマシン) 等の他の機械学習手法でも使われる.
他によく使われるカーネル関数として, Matern カーネルなどがある.
$\theta_1,\theta_2,\theta_3$ はハイパーパラメータ.
${\theta_3 \delta(x_i, x_j)}$ の項は "観測ノイズ" を表している.
(参考書では, カーネル関数に観測ノイズを含めている)

実行例 1

トレーニングデータの個数 $N$ をいろいろ変化させてみた結果

訓練データが密集している領域は分散が小さく, 密度が低い領域は分散が大きくなっている.
→ 訓練データに誤差(ノイズ) が含まれていないため.

実行例 2

カーネル関数のパラメータ $\theta$ をいろいろ変化させてみた結果

まとめ

カーネル関数のパラメータを変えると, 結果も変化する.
→ カーネル関数のパラメータも学習によって最適化する必要がある.

カーネル関数のパラメータθの最適化

最尤推定: 尤度関数を $\theta$ で微分して勾配を求め, 勾配法で最良の $\theta$ を求める.

訓練データ $\mathcal{D} = (X, y) = ([x_1, \dots, x_N], [y_1, \dots, y_N])$
- $x_i$ はベクトル, $y_i$ はスカラー

尤度関数は, 以下のように書き表せる.

p(y|X,\theta)
= \mathcal{N}(y|0, K_\theta)
= \frac{ \exp{[-\frac{1}{2} y^T {K_\theta}^{-1} y ]} }{ \sqrt{(2 \pi)^N |K_\theta|} }

尤度関数の対数を取り, パラメータごとに微分する (省略).

\frac{\partial L}{\partial \theta} =
\begin{cases}
    \frac{\partial L}{\partial \theta_1} = ... \\
    \frac{\partial L}{\partial \theta_2} = ... \\
    \frac{\partial L}{\partial \theta_3} = ...
\end{cases}

多変量正規分布の密度関数 ¹

$X \sim \mathcal{N}(\mu, \Sigma)$ のとき, 確率密度関数は

f_X(x_1, ..., x_k) = \frac{ \exp{[-\frac{1}{2} (x-\mu)^T \Sigma^{-1} (x-\mu) ]} }{ \sqrt{(2 \pi)^k |\Sigma|} }

実行例 3
カーネル関数のパラメータを最適化した結果

左上: 最適化前 / 右上: 最適化後
左下: 対数尤度(訓練データ) / 右下: 対数尤度(全データ)

⇒ 最適化後はより真値にフィットし, 尤度が高くなっている.

ガウス過程回帰の問題点

ガウス過程回帰ではカーネル行列の逆行列 $K^{-1}$ を求める必要があり, 計算コストが大きい.
(訓練データ $\mathcal{D}$ のサイズを $N$ とすると, およそ $O(N^3)$ の計算量)

ベイズ最適化

アルゴリズム

評価関数 $f(x)$ によって計算される評価値がなるべく大きくなるような $x$ を求めたい.
($x^* = \argmax_{x \in \mathcal{X}} [{f(x)}]$ を求めたい.)

ただし, $f(x)$ の評価には大きなコストがかかるため, なるべく不要な評価をすることなく 探索したい.

★ ガウス過程回帰を用いることで, 未知の入力 $x^?$ に対する評価値 $f(x^?)$ の確率分布を推測でき, どの入力が最適値の探索に有用か目処が立てられる.

具体的には, 次の 1 ~ 3 を繰り返すことで, 最適値に近づける.

次に探索する入力 $x^?$ を策定
選んだ入力 $x^?$ を評価 (評価値 $y^?$ を計算)
ガウス過程回帰モデルの訓練データ $\mathcal{D}$ に $(x^?, y^?)$ を追加

ベイズ最適化以外のハイパーパラメータ探索法
ランダムサーチ, グリッドサーチ (空間を格子状に区切って探索) など

獲得関数 (acquisition function)

次に探索する入力 $x^?$ の策定において, 具体的には
$\mu$ と $\sigma$ を適当に組み合わせて定義した 獲得関数 が最大になるような $x^?$ を選定する.

獲得関数は強化学習の報酬関数に似ており, 探索 (Exploration) と 活用 (Exploitation) のトレードオフをうまく汲み取った関数を作る.
(なるべく期待値 $\mu$ が高く, 分散 $\sigma$ が大きい (=未探索) $x$ についての値が大きくなるように設計する.)

Probability of Improvement (PI)
Expected Improvement (EI)

GP Upper Confidence Band (GP-UCB)

a_\text{UCB}(x) = \mu(x) + k_N \sigma(x)

など

ベイズ最適化のプログラム例 (optuna)

import optuna

# ハイパーパラメータの評価関数を定義
def f(x: float):
    ...

def criteria(trial: optuna.Trial):
    x = trial.suggest_float("x", -4, 4)  # `float` のハイパーパラメータを定義
    # 他にも `trial.suggest_int()`, `trial.suggest_categorical()` など
    # 様々なパラメータタイプが実装されている.
    return f(x)

study = optuna.create_study(direction="maximize")
study.optimize(criteria, n_trials=N)

# 各トライアルで評価したパラメータ $x$ をプロット (`trials`)
# 各トライアルの段階で発見した最良のパラメータ ($x^*$) をプロット (`best_trial`)
...

💭 前半は未探索の領域の探索を重視し, 後半は最適値が存在すると思われる領域の探索を重視している (ように見える).
→ 強化学習の ε-Greedy 法と同じ考え方

https://ja.wikipedia.org/wiki/%E5%A4%9A%E5%A4%89%E9%87%8F%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83 ↩