LoginSignup
0
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@s059ff

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

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

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

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

キーワード

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

参考書

reference_book_thumbnail.jpg

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

確率過程とは

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

ガウス過程とは

  • 各時刻の値の確率分布が, 正規分布で表される確率過程.

  • なぜ正規分布を使うのか?
    → どのような確率分布も, 平均をとってしまえば正規分布に近づくため (中心極限定理).

    中心極限定理
    ある母集団から $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をプロット
...

    central_limit_theorem_example.png
    はじめは一様分布の形だが, $n → \infty$ に近づくにつれ, だんだん正規分布の形に近づくことがわかる.

  • 回帰問題

    • 連続値を推定する問題.
    • ↔ クラス分類 (離散値を推定する問題)

  • ガウス過程回帰

    • ガウス過程を使って解く回帰問題

ガウス過程回帰とは

結論

gpr_question_diagram.drawio.png

カーネル関数 $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)$

導出

  1. 訓練データ $\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}

  2. このとき, ${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}

  3. 未知のデータ $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}

  4. ここで, $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$ はいくらと推定できるだろうか?

アルゴリズム

  1. カーネル関数 $k(\cdot, \cdot)$ を設計
  2. $x_\text{train}$ から カーネル行列$K$ を計算
    $K[i,j] \leftarrow k(x_i, x_j)$
  3. $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$ をいろいろ変化させてみた結果

gpr_prediction_example-10-1.0-0.4-0.1.png
gpr_prediction_example-20-1.0-0.4-0.1.png
gpr_prediction_example-40-1.0-0.4-0.1.png

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

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

gpr_prediction_example-10-0.5-0.4-0.1.png
gpr_prediction_example-10-1.0-0.8-0.1.png
gpr_prediction_example-10-1.0-0.4-0.2.png

まとめ

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

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

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

  1. 訓練データ $\mathcal{D} = (X, y) = ([x_1, \dots, x_N], [y_1, \dots, y_N])$

    • $x_i$ はベクトル, $y_i$ はスカラー

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

    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|} }

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

    \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}

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

$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
カーネル関数のパラメータを最適化した結果

gpr_training_example.png

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

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

ガウス過程回帰の問題点

ガウス過程回帰ではカーネル行列の逆行列 $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 を繰り返すことで, 最適値に近づける.

  1. 次に探索する入力 $x^?$ を策定
  2. 選んだ入力 $x^?$ を評価 (評価値 $y^?$ を計算)
  3. ガウス過程回帰モデルの訓練データ $\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`)
...

optuna_example.png

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

  1. 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 ↩

0
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?