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

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

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

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

## キーワード

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

## 参考書

![reference_book_thumbnail.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191898/bb618eb1-fa82-2515-8e96-175acc0127a4.jpeg)

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

## 確率過程とは

- 時間などの条件によって変化する確率. 為替や株価など.
- **関数自体が確率によって変化** することを考える.
  **イメージ図**
  ![stochastic_process_diagram.drawio.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191898/abc39ddd-f120-dd1e-b301-e3ebe2c4cd7e.png)


## ガウス過程とは

- 各時刻の値の確率分布が, 正規分布で表される確率過程.
- なぜ正規分布を使うのか?
  → どのような確率分布も, 平均をとってしまえば正規分布に近づくため (**中心極限定理**).

  :::note info
  **中心極限定理**
  ある母集団から $N$ 個のデータをサンプリングし, それらの平均値を確率変数 $Y$ とする.
  $N \rightarrow \infty$ に近づくにつれ, 確率変数 $Y$ (を標準化したもの) は標準正規分布に近づく.

  例: サイコロを $N \times M$ 回振り, $N$ 個を1グループとして平均値$Y$を求める. $Y \in \mathbb{R}^M$ をヒストグラム化する.

    ```python
    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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191898/335181e5-7f0d-f2da-7c76-cbaded9dc209.png)
    はじめは一様分布の形だが, $n → \infty$ に近づくにつれ, だんだん正規分布の形に近づくことがわかる.
  :::

- **回帰問題**

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

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

## ガウス過程回帰とは

### 結論

![gpr_question_diagram.drawio.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191898/0aa8f409-7e09-6739-186c-435e94ac20a7.png)

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

```math
\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)$

### 導出

0. 訓練データ $\mathcal{D}$ とカーネル関数 $k(x_i, x_j)$ が与えられる. ただし, $y$ は簡単化のため $0$ とする.

   ```math
   \mathcal{D} = \{(x_1, y_1), (x_2, y_2), ..., (x_N, y_N)\}
   ```

   ※ $y$ の観測値にノイズ $\varepsilon \sim \mathcal{N}(0, \sigma^2)$ が含まれている場合を考慮して, カーネル関数 $k$ は観測ノイズを含めて考える.

   ```math
   k(x_i, x_j) = (x_i と x_j の距離) + \sigma^2 \delta(x_i, x_j)
   ```

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

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

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

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

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

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

   ```math
   \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$ はいくらと推定できるだろうか?

#### アルゴリズム

0. カーネル関数 $k(\cdot, \cdot)$ を設計
1. $x_\text{train}$ から カーネル行列$K$ を計算
   $K[i,j] \leftarrow k(x_i, x_j)$
2. $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$ の確率分布をグラフにプロットする.

##### カーネル関数の設計

```math
\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)}
```

```math
\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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191898/edd61980-5f55-8809-d2a3-0ad8146c5876.png)
![gpr_prediction_example-20-1.0-0.4-0.1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191898/8dc06dce-03d7-3b0b-fc88-4803e4009951.png)
![gpr_prediction_example-40-1.0-0.4-0.1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191898/8ddd509f-e420-cf6d-9c50-5ad7b2497d12.png)

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

##### 実行例 2

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

![gpr_prediction_example-10-0.5-0.4-0.1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191898/7aef2b1e-0188-c4e2-97ad-40bc8c2f1e94.png)
![gpr_prediction_example-10-1.0-0.8-0.1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191898/90d2403b-f56f-f6f6-7d8e-560809205661.png)
![gpr_prediction_example-10-1.0-0.4-0.2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191898/1fc9ca8b-00b3-e54f-689f-6f94655e227c.png)

###### まとめ

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

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

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

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

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

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

   ```math
   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. 尤度関数の対数を取り, パラメータごとに微分する (省略).

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

:::note info
**多変量正規分布の密度関数 [^1]**

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

```math
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191898/9b49fa40-48b9-a7a6-4121-83437a3d890d.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^?)$ を追加

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

### 獲得関数 (acquisition function)

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

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

- Probability of Improvement (PI)
- Expected Improvement (EI)
- GP Upper Confidence Band (GP-UCB)
  ```math
  a_\text{UCB}(x) = \mu(x) + k_N \sigma(x)
  ```

など

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

```python
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191898/d5af639b-4345-f530-6378-09c3b6b27183.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