「ベイズ最適化」を読む　～第2章ブラックボックス関数のベイズモデリング ②ガウス過程回帰モデル～

はじめに

データを使って仮説の生成と検証を行なうための方法であるベイズ最適化を学ぶために，今村秀明・松井孝太著「ベイズ最適化　ー適応的実験計画の基礎と実践ー」を読むことにした。

本記事は，「第2章　ブラックボックス関数のベイズモデリング」における，ガウス過程回帰モデルに関する読書メモである。

本書の紹介ページ

www.kindaikagaku.co.jp

2.2 ガウス過程回帰モデル

ベイズ最適化では，関数を近似する際にベイズモデリングを用いる。良く用いられるモデルとして，ベイズ線形回帰モデルとガウス過程回帰モデルが用いられる。

ガウス過程回帰は，森賀新・木田悠歩・須山敦志著「Pythonではじめるベイズ機械学習入門」を読んだ際にも出てきた話題であり，読書メモは以下にまとめた。
stern-bow.hatenablog.com

ベイズ線形回帰モデルからガウス過程回帰モデルへ

ベイズ線形回帰モデルでは，モデルパラメータ $\boldsymbol{w}$ と基底関数 $\boldsymbol{\phi}$ を用いて関数 $f$ をモデル化していた。しかしこの場合，入力 $x$ の次元や学習データの数が増えると，推定するべきモデルパラメータの数が指数関数的に増える(次元の呪い)という問題がある。

次元の呪いを回避するために，パラメータに依存しないノンパラメトリックなモデリング方法を用いることが考えられる。ガウス過程回帰は，ノンパラメトリックな方法に含まれる。

モデルパラメータ $\boldsymbol{w}$ が，多変量正規分布 $\mathcal{N}(\boldsymbol{0}, \Sigma)$ にしたがうとする。
このとき関数値ベクトル $\hat{\boldsymbol{f}} = (\boldsymbol{w}^T \boldsymbol{\phi}_1, \cdots, \boldsymbol{w}^T \boldsymbol{\phi}_n)^T = \boldsymbol{\Phi} \boldsymbol{w} \in \mathbb{R}^n$ がしたがう分布は，

$\begin{align} p(\hat{\boldsymbol{f}} | \boldsymbol{\Phi}) \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\Phi} \boldsymbol{\Sigma} \boldsymbol{\Phi}^T) \\ \end{align}$

になる。これは，関数値ベクトル $\hat{\boldsymbol{f}}$ が，モデルパラメータ $\boldsymbol{w}$ の線形変換になっているからである。
また関数値ベクトル $\hat{\boldsymbol{f}}$ がしたがう確率分布には，モデルパラメータ $\boldsymbol{w}$ が含まれていないので，ノンパラメトリックなモデルとなっている。

この式は， $n$ 次元の関数値ベクトル $\boldsymbol{f}$ が正規分布にしたがうという仮定を表している。また関数 $f(\boldsymbol{x})$ がガウス過程(Gaussian process)にしたがっていることを表し，

$\begin{align} \boldsymbol{f} \sim \mathcal{GP}(\mu, k) \\ \end{align}$

と表す。ただし $\mu$ は平均関数， $k$ は共分散関数またはカーネル関数と呼ぶ。

カーネル関数を $k(\boldsymbol{x}, \boldsymbol{x}' ) = \boldsymbol{\phi} ( \boldsymbol{x} )^T \Sigma \boldsymbol{\phi} ( \boldsymbol{x} )$ と表現すると，基底関数 $\boldsymbol{\phi}$ を明示的に表さなくてもよいことが分かる。

ガウス過程の基本的な性質

カーネル関数として，ガウスカーネル $k(x, x') = \exp(-\theta | x- x'|^2)$ を用いると，ブラックボックス関数 $f$ に対して，

入力[tex; x, x']が近ければ， $f(x), f(x')$ も近い値を取る。
また入力[tex; x, x']が離れていれば， $f(x), f(x')$ はほとんど独立であるように振舞う。

という性質をモデルに反映することが分かる。

数値実験

本書P36 図2.6を参考に，ガウス過程事前分布モデルのサンプルパスと，観測データで条件付けたガウス過程事後分布モデルのサンプルパスを描画した。
データが得られることにより，そのデータ周りでの誤差が低減されることが確認できた。

上記の分析用pythonコードはこちら。

クリックで展開

import numpy as np
from scipy.stats import multivariate_normal

np.random.seed(0)

# カーネル関数
def kernel(x, x_prime):
  return np.exp(-0.5 * np.abs(x - x_prime)**2)

# 入力値
x = np.linspace(0, 10, 100)

# 観測データがない場合 ##########
mu = np.zeros(len(x))
cov = np.zeros((len(x), len(x)))
for i in range(len(x)):
  for j in range(len(x)):
    cov[i, j] = kernel(x[i], x[j])


# 観測データがある場合 #####
x_obs = np.array([4.0])
y_obs = np.array([1.0])

# 観測データを含むカーネル行列
cov_obs = np.zeros((len(x) + 1, len(x) + 1))
for i in range(len(x) + 1):
  for j in range(len(x) + 1):
    if i < len(x) and j < len(x):
      cov_obs[i, j] = kernel(x[i], x[j])
    elif i == len(x):
      cov_obs[i,j] = kernel(x_obs[0], x[j]) if j < len(x) else kernel(x_obs[0],x_obs[0])
    elif j == len(x):
        cov_obs[i,j] = kernel(x[i], x_obs[0])

# 条件付きガウス分布を計算
cov_x_obs = cov_obs[:len(x),len(x):]
cov_obs_x = cov_obs[len(x):,:len(x)]
cov_obs_obs = cov_obs[len(x):,len(x):]
mu_cond = mu + cov_x_obs @ np.linalg.inv(cov_obs_obs) @ (y_obs - mu[-1])
cov_cond = cov - cov_x_obs @ np.linalg.inv(cov_obs_obs) @ cov_obs_x
std_cond = np.sqrt(np.diag(cov_cond))
upper_cond = mu_cond + 1.96 * std_cond
lower_cond = mu_cond - 1.96 * std_cond


# グラフ描画 #####
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()  # seabornスタイルを設定
fig, axes = plt.subplots(1, 2, figsize=(20, 8))

# 観測データがある場合 ########
# サンプルパスを描画
np.random.seed(0)
for i in range(3):
  f = np.random.multivariate_normal(mu, cov)
  axes[0].plot(x, f)

# 信用区間を計算
std = np.sqrt(np.diag(cov))
upper = mu + 1.96 * std
lower = mu - 1.96 * std
axes[0].fill_between(x, upper, lower, alpha=0.2)

axes[0].set_title("Gaussian Process (No Observation Data)")
axes[0].set_xlabel("x")
axes[0].set_ylabel("f(x)")
axes[0].set_xlim(0, 10)
axes[0].set_ylim(-4, 4)


# 観測データがない場合 ########
# サンプルパスを描画
np.random.seed(0)
for i in range(3):
  f = np.random.multivariate_normal(mu_cond, cov_cond)
  axes[1].plot(x, f)

# 信用区間を計算
axes[1].fill_between(x, upper_cond, lower_cond, alpha=0.2)
axes[1].plot(x_obs,y_obs,'ro')

axes[1].set_title("Gaussian Process (With Observation Data)")
axes[1].set_xlabel("x")
axes[1].set_ylabel("f(x)")
axes[1].set_xlim(0, 10)
axes[1].set_ylim(-4, 4)

plt.tight_layout()
plt.show()

ガウス過程モデルの推論

ベイズ線形回帰モデルの時と同様に，ガウス過程回帰モデルでも，新しい入力 $\boldsymbol{x}_*$ が得られたときの出力 $y_*$ を推論することに興味がある。

関数 $f$ に対して，平均関数が0のガウス過程事前分布を仮定する。すなわち， $f \sim \mathcal{GP}(0, k(\boldsymbol{x}, \boldsymbol{x}'))$ である。

観測誤差がない場合

$n$ 個の入力データ $\boldsymbol{x}_i, i=1, \cdots, n$ を用いて，関数値ベクトル $\boldsymbol{f} = (f(\boldsymbol{x}_1), \cdots, f(\boldsymbol{x}_n))^T$ を算出する。
このとき，新しい入力 $\boldsymbol{x}_*$ と，これに対する関数値 $f_*$ について，

$\begin{align} \begin{bmatrix} \boldsymbol{f} \\ f_* \end{bmatrix} \sim \mathcal{N} \left( \boldsymbol{0}, \begin{bmatrix} \boldsymbol{K}_n & \boldsymbol{k}_n(\boldsymbol{x}_*) \\ \boldsymbol{k}_n(\boldsymbol{x}_*)^T & k(\boldsymbol{x}_*, \boldsymbol{x}_*) \end{bmatrix} \right) \end{align}$

となる。ただし， $\boldsymbol{K}_n$ は $K_{ij} = k(\boldsymbol{x}_i, \boldsymbol{x}_j)$ となるカーネル行列であり，また

$\begin{align} \boldsymbol{k}_n(\boldsymbol{x}_*) = (k(\boldsymbol{x}_*, \boldsymbol{x}_1), \cdots, k(\boldsymbol{x}_*, \boldsymbol{x}_n)^T) \\ \end{align}$

である。

これは $\boldsymbol{f}$ と $f_*$ の同時分布になっている。いま興味がある $f_*$ は，条件付き確率分布 $p(f_* | \boldsymbol{f})$ を求めればよいが，これは多変量正規分布の条件付き確率の性質を用いれば求めることができる。

参考 :
- 「Pythonではじめるベイズ機械学習入門」を読む～第3章回帰モデル⑤(ガウス過程回帰モデル：ガウス尤度)～ - jiku log
- 多変量正規分布の条件付き分布の思い出し方(増補改訂版) #統計検定 - jiku log

観測誤差がある場合

観測誤差がある場合，すなわち $y_i = f(\boldsymbol{x}_i) + \varepsilon_i, \varepsilon_i \sim \mathcal{N}(0, \sigma^2)$ の場合を考える。
このとき， $\boldsymbol{y} = (y_1, \cdots, y_n)^T$ について，

$\begin{align} \boldsymbol{y} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{K}_n + \sigma^2 I_n) \\ \end{align}$

となるので，

$\begin{align} \begin{bmatrix} \boldsymbol{y} \\ f_* \end{bmatrix} \sim \mathcal{N} \left( \boldsymbol{0}, \begin{bmatrix} \boldsymbol{K}_n + \sigma^2 I_n & \boldsymbol{k}_n(\boldsymbol{x}_*) \\ \boldsymbol{k}_n(\boldsymbol{x}_*)^T & k(\boldsymbol{x}_*, \boldsymbol{x}_*) \end{bmatrix} \right) \end{align}$

となる。

ガウス過程の平均関数

ベイズ最適化の文脈では，ガウス過程の平均関数 $\mu$ の設定方法はあまり大きな問題にはならない。
これは，

ベイズ最適化では，分散の方に興味がある。
0でない平均関数 $\mu$ を設定して $f_*$ の事後分布を計算しても，分散には平均関数が含まれない。

ためである。
そのためベイズ最適化では，平均関数には定数0が置かれることが多い。

カーネル関数のハイパーパラメータ

ガウスカーネル $k(x, x') = \exp(-\theta | x- x'|^2)$ に含まれるハイパーパラメータ $\theta$ を最適化するためには，周辺尤度最大化法を用いる。
すなわちブラックボックス関数 $f$ を周辺化によって積分消去し，この周辺尤度を $\theta$ について最大化する。
詳細は，本書の第3章に紹介されている。