ベイズ最適化（Bayesian Optimization）は、評価コストが高いブラックボックス関数の最適化を効率的に行う手法です。例えば、ニューラルネットワークの学習率やバッチサイズといったハイパーパラメータを探索する場面を考えましょう。1回の学習に数時間かかるモデルに対して、グリッドサーチやランダムサーチで何百通りも試すのは現実的ではありません。ベイズ最適化は、「これまでの評価結果から、次にどの点を試せば最も有益か」をガウス過程の予測分布に基づいて賢く判断することで、少ない評価回数で良い解を見つけます。
ベイズ最適化の応用先はハイパーパラメータチューニングだけではありません。風洞実験や化学実験の条件最適化、CFD や有限要素法を用いたシミュレーションベースの設計最適化など、「関数の形は分からないが、1回の評価にコストがかかる」あらゆる問題に適用できます。
本記事の内容
- ベイズ最適化のフレームワーク（代理モデル + 獲得関数のループ）
- 獲得関数 PI, EI, UCB の数学的な導出
- 探索と活用のトレードオフ
- Python でのスクラッチ実装と各ステップの可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ガウス過程回帰では、観測データが与えられたとき、未知の入力 $\bm{x}_*$ に対する予測分布が平均 $\mu(\bm{x}_*)$ と分散 $\sigma^2(\bm{x}_*)$ を持つガウス分布として得られます。本記事の獲得関数はすべてこの予測平均と予測分散を使って構成されるため、ガウス過程回帰の仕組みを理解しておくことが前提となります。
ベイズ最適化のフレームワーク
ベイズ最適化は、次の2つの要素を軸に構成されます。
- 代理モデル（Surrogate Model）: 目的関数 $f(\bm{x})$ を近似する確率的な回帰モデル。通常はガウス過程回帰を使う
- 獲得関数（Acquisition Function）: 代理モデルの予測に基づいて「次にどの点を評価すべきか」を定量化する関数
最適化したい目的関数 $f(\bm{x})$ の解析的な式は不明であり、勾配も使えない状況を想定します。目標は次の最小化問題を解くことです。
$$ \bm{x}^{*} = \arg\min_{\bm{x} \in \mathcal{X}} f(\bm{x}) $$
ベイズ最適化のアルゴリズムは、以下のループを繰り返します。
- 初期点をいくつか選び、目的関数 $f(\bm{x})$ を評価する
- 評価済みデータ $\mathcal{D}_n = \{(\bm{x}_i, y_i)\}_{i=1}^{n}$ を使ってガウス過程回帰モデルを学習する
- 獲得関数 $\alpha(\bm{x})$ を最大化する点 $\bm{x}_{n+1} = \arg\max_{\bm{x} \in \mathcal{X}} \alpha(\bm{x})$ を求める
- 目的関数を評価する: $y_{n+1} = f(\bm{x}_{n+1})$
- データを更新 $\mathcal{D}_{n+1} = \mathcal{D}_n \cup \{(\bm{x}_{n+1}, y_{n+1})\}$ し、ステップ 2 に戻る
このループを評価回数の予算に達するまで繰り返します。
なぜガウス過程を使うのか
代理モデルにはさまざまな手法が考えられますが、ベイズ最適化ではガウス過程回帰が最も広く採用されています。その理由は、ガウス過程回帰が予測値だけでなく予測の不確実性（分散）も同時に出力してくれるからです。
ガウス過程回帰では、未知の入力 $\bm{x}$ に対する予測分布が次のガウス分布として得られます。
$$ f(\bm{x}) \mid \mathcal{D}_n \sim \mathcal{N}(\mu_n(\bm{x}),\; \sigma_n^2(\bm{x})) $$
データが密に観測されている領域では $\sigma_n^2(\bm{x})$ が小さくなり（自信がある）、データが少ない領域では $\sigma_n^2(\bm{x})$ が大きくなります（不確実性が高い）。この予測平均と予測分散の両方を使って、次にどこを評価すべきかを判断するのが獲得関数の仕事です。
獲得関数1: PI（Probability of Improvement）
最もシンプルな獲得関数が PI（Probability of Improvement）です。「ある候補点 $\bm{x}$ を評価したとき、現在の最良値を改善する確率はどれくらいか？」を測ります。
PI の定式化
これまでに観測された目的関数の最小値を $f^+ = \min_{i=1,\dots,n} y_i$ とします（最小化問題を考えます）。候補点 $\bm{x}$ における改善確率は、
$$ \text{PI}(\bm{x}) = P(f(\bm{x}) < f^+ \mid \bm{x}, \mathcal{D}_n) $$
と定義されます。
PI の導出
ガウス過程の予測分布より $f(\bm{x}) \sim \mathcal{N}(\mu(\bm{x}), \sigma^2(\bm{x}))$ です。確率変数 $f(\bm{x})$ を標準化します。
$$ Z = \frac{f(\bm{x}) – \mu(\bm{x})}{\sigma(\bm{x})} \sim \mathcal{N}(0, 1) $$
すると、改善確率は以下のように計算できます。
$$ \begin{align} \text{PI}(\bm{x}) &= P(f(\bm{x}) < f^+) \\ &= P\!\left(\frac{f(\bm{x}) - \mu(\bm{x})}{\sigma(\bm{x})} < \frac{f^+ - \mu(\bm{x})}{\sigma(\bm{x})}\right) \\ &= \Phi\!\left(\frac{f^+ - \mu(\bm{x})}{\sigma(\bm{x})}\right) \end{align} $$
ここで $\Phi(\cdot)$ は標準正規分布の累積分布関数です。よって、PI の式は次のようになります。
$$ \boxed{\text{PI}(\bm{x}) = \Phi(\gamma), \quad \gamma = \frac{f^+ – \mu(\bm{x})}{\sigma(\bm{x})}} $$
PI は直感的に分かりやすい指標ですが、「どれだけ改善するか」ではなく「改善するか否か」のみを見ています。そのため、わずかな改善であっても確率が高ければ優先してしまい、現在の最良値付近ばかりを調べる活用に偏った振る舞いをする傾向があります。
獲得関数2: EI（Expected Improvement）
EI（Expected Improvement）は PI の弱点を克服した獲得関数で、「改善量の期待値」を最大化します。ベイズ最適化において最も広く使われている獲得関数です。
EI の定式化
改善量（Improvement）を次のように定義します。
$$ I(\bm{x}) = \max(f^+ – f(\bm{x}), \; 0) $$
$f(\bm{x})$ が $f^+$ より小さければ改善量は $f^+ – f(\bm{x}) > 0$、そうでなければ改善量は 0 です。EI はこの改善量の期待値として定義されます。
$$ \text{EI}(\bm{x}) = \mathbb{E}\bigl[\max(f^+ – f(\bm{x}), \; 0)\bigr] $$
EI の解析解の導出
EI の大きな利点は、ガウス過程の予測分布のもとで解析的に閉じた形が得られることです。以下、丁寧に導出します。
$f(\bm{x}) \sim \mathcal{N}(\mu, \sigma^2)$ と略記します（$\mu = \mu(\bm{x})$, $\sigma = \sigma(\bm{x})$）。また、$\phi(\cdot)$ を標準正規分布の確率密度関数、$\Phi(\cdot)$ を累積分布関数とします。
期待値の定義に従い、$f(\bm{x}) < f^+$ の範囲のみが寄与するため、
$$ \text{EI}(\bm{x}) = \int_{-\infty}^{f^+} (f^+ – y) \cdot \frac{1}{\sigma\sqrt{2\pi}} \exp\!\left(-\frac{(y – \mu)^2}{2\sigma^2}\right) dy $$
ここで $z = \frac{y – \mu}{\sigma}$ と標準化します。$y = \mu + \sigma z$, $dy = \sigma \, dz$ であり、$y = f^+$ のとき $z = \frac{f^+ – \mu}{\sigma}$ です。この値を $\gamma = \frac{f^+ – \mu}{\sigma}$ とおきます。
$$ \begin{align} \text{EI}(\bm{x}) &= \int_{-\infty}^{\gamma} (f^+ – \mu – \sigma z) \cdot \frac{1}{\sqrt{2\pi}} e^{-z^2/2} \, dz \\ &= \int_{-\infty}^{\gamma} (f^+ – \mu – \sigma z) \, \phi(z) \, dz \end{align} $$
ここで $f^+ – \mu = \sigma \gamma$ であることに注意して、2つの積分に分解します。
$$ \text{EI}(\bm{x}) = \sigma\gamma \int_{-\infty}^{\gamma} \phi(z) \, dz – \sigma \int_{-\infty}^{\gamma} z \, \phi(z) \, dz $$
第1項 について、$\int_{-\infty}^{\gamma} \phi(z) \, dz = \Phi(\gamma)$ なので、
$$ \text{第1項} = \sigma\gamma \, \Phi(\gamma) = (f^+ – \mu) \, \Phi(\gamma) $$
第2項 を計算します。$\phi(z) = \frac{1}{\sqrt{2\pi}} e^{-z^2/2}$ の微分が次の関係を満たすことを利用します。
$$ \frac{d\phi}{dz} = -z \, \phi(z) \quad \Longleftrightarrow \quad z \, \phi(z) = -\frac{d\phi}{dz} $$
これを代入すると、
$$ \begin{align} \int_{-\infty}^{\gamma} z \, \phi(z) \, dz &= \int_{-\infty}^{\gamma} \left(-\frac{d\phi}{dz}\right) dz \\ &= -\bigl[\phi(z)\bigr]_{-\infty}^{\gamma} \\ &= -\phi(\gamma) + \underbrace{\phi(-\infty)}_{= 0} \\ &= -\phi(\gamma) \end{align} $$
よって、第2項は、
$$ \text{第2項} = -\sigma \cdot (-\phi(\gamma)) = \sigma \, \phi(\gamma) $$
以上を合わせると、EI の解析解が得られます。
$$ \boxed{\text{EI}(\bm{x}) = (f^+ – \mu(\bm{x})) \, \Phi(\gamma) + \sigma(\bm{x}) \, \phi(\gamma), \quad \gamma = \frac{f^+ – \mu(\bm{x})}{\sigma(\bm{x})}} $$
ここで、$\Phi(\cdot)$ は標準正規分布の累積分布関数、$\phi(\cdot)$ は標準正規分布の確率密度関数です。$\sigma(\bm{x}) = 0$ の場合（既に観測済みの点）は $\text{EI}(\bm{x}) = 0$ と定義します。
EI の直感的な解釈
EI の式は2つの項からなり、それぞれ異なる役割を持っています。
- 第1項 $(f^+ – \mu) \, \Phi(\gamma)$: 予測平均 $\mu$ が現在の最良値 $f^+$ より小さいとき（つまり改善が見込めるとき）に大きくなります。これは活用（Exploitation）に対応する項です。
- 第2項 $\sigma \, \phi(\gamma)$: 予測分散 $\sigma$ が大きい場所で大きくなります。まだデータが少ない領域を積極的に調べる探索（Exploration）に対応する項です。
このように、EI は探索と活用のバランスを自動的かつ自然に取ってくれます。
獲得関数3: UCB（Upper Confidence Bound）
UCB（Upper Confidence Bound）は、探索と活用のトレードオフをハイパーパラメータ $\kappa$ で明示的に制御する獲得関数です。GP-UCB（Gaussian Process Upper Confidence Bound）とも呼ばれます。
UCB の定式化
最小化問題では、LCB（Lower Confidence Bound）を用います。信頼区間の下限を獲得関数として使い、これを最小化する点を次の評価点に選びます。
$$ \text{LCB}(\bm{x}) = \mu(\bm{x}) – \kappa \, \sigma(\bm{x}) $$
最大化問題の場合は、信頼区間の上限を獲得関数として使い、これを最大化します。
$$ \boxed{\text{UCB}(\bm{x}) = \mu(\bm{x}) + \kappa \, \sigma(\bm{x})} $$
ここで $\kappa > 0$ は探索と活用のバランスを制御するハイパーパラメータです。
- $\kappa$ が大きい: $\sigma(\bm{x})$ の寄与が大きくなり、不確実性の高い領域を優先する（探索寄り）
- $\kappa$ が小さい: $\mu(\bm{x})$ が支配的になり、現時点で最良と予測される領域を優先する（活用寄り）
典型的には $\kappa = 2.0$ 前後が使われます。理論的には、$\kappa$ を反復回数 $t$ に応じて増加させるスケジュールを用いると、累積リグレットに関する理論保証が得られることが知られています（Srinivas et al., 2010）。
UCB の利点
UCB の式はシンプルであり、EI のような積分の導出を必要としません。また、$\kappa$ を直接調整できるため、「この問題ではもっと探索を重視したい」「活用を強めたい」といった意図を反映しやすいという利点があります。
探索と活用のトレードオフ
ベイズ最適化において最も本質的な概念が、探索（Exploration）と活用（Exploitation）のトレードオフです。
- 活用: 現時点で最良と予測される領域を集中的に評価する。局所的な最適解を精密に求めることができるが、大域的により良い解を見落とす可能性がある
- 探索: まだデータが少なく不確実性が高い領域を評価する。大域的な最適解を発見できる可能性があるが、すでに良いと分かっている領域を十分に活用できない
獲得関数を設計するとは、すなわちこのトレードオフをどう制御するかを決めることに他なりません。各獲得関数の特性を改めてまとめます。
| 獲得関数 | 式 | 探索と活用のバランス |
|---|---|---|
| PI | $\Phi(\gamma)$ | 活用に偏りやすい。改善確率のみを見るため、わずかな改善でも優先する |
| EI | $(f^+ – \mu)\Phi(\gamma) + \sigma\phi(\gamma)$ | 自動的にバランスを取る。第1項が活用、第2項が探索に対応 |
| UCB | $\mu + \kappa\sigma$ | $\kappa$ で明示的に制御可能。$\kappa$ が大きいほど探索寄り |
実用上は、まず EI を試し、必要に応じて UCB で $\kappa$ を調整するアプローチが一般的です。
アルゴリズム全体のステップ
ここで、ベイズ最適化の全体像をアルゴリズムとして整理しておきましょう。
入力: 目的関数 $f$, 探索空間 $\mathcal{X}$, 初期点数 $n_0$, 最大評価回数 $N$
- 初期化: $n_0$ 個の点 $\{\bm{x}_1, \dots, \bm{x}_{n_0}\}$ をランダムに（またはラテン超方格法で）選び、$y_i = f(\bm{x}_i)$ を評価する。データ集合 $\mathcal{D}_0 = \{(\bm{x}_i, y_i)\}_{i=1}^{n_0}$ を構築する
- ループ（$t = n_0, n_0 + 1, \dots, N – 1$）:
- データ $\mathcal{D}_t$ を用いてガウス過程回帰を学習する
- 探索空間上の各点で獲得関数 $\alpha(\bm{x})$ を計算する
- $\bm{x}_{t+1} = \arg\max_{\bm{x} \in \mathcal{X}} \alpha(\bm{x})$ を求める
- 目的関数を評価する: $y_{t+1} = f(\bm{x}_{t+1})$
- データを更新する: $\mathcal{D}_{t+1} = \mathcal{D}_t \cup \{(\bm{x}_{t+1}, y_{t+1})\}$
- 出力: 観測された最良の点 $\bm{x}^* = \arg\min_i y_i$ を返す
ステップ 2-3 の獲得関数の最大化は、代理モデル（ガウス過程）の評価が $f$ の評価に比べて圧倒的に安価なため、グリッドサーチや L-BFGS-B などの通常の最適化手法で高速に行えます。
Python 実装
ここからは、ガウス過程回帰と EI 獲得関数を使ったベイズ最適化を Python でスクラッチ実装します。1次元のテスト関数を対象に、各ステップで代理モデルと獲得関数を可視化し、最適解に収束していく様子を確認しましょう。
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
# --- RBFカーネル ---
def rbf_kernel(X1, X2, length_scale=1.0, variance=1.0):
"""RBF（ガウス）カーネルを計算する"""
sqdist = (np.sum(X1**2, axis=1).reshape(-1, 1)
+ np.sum(X2**2, axis=1).reshape(1, -1)
- 2 * X1 @ X2.T)
return variance * np.exp(-0.5 * sqdist / length_scale**2)
# --- ガウス過程回帰 ---
def gpr_predict(X_train, y_train, X_test,
length_scale=1.0, variance=1.0, noise=1e-6):
"""ガウス過程回帰の予測平均と予測標準偏差を返す"""
K = rbf_kernel(X_train, X_train, length_scale, variance) \
+ noise * np.eye(len(X_train))
K_s = rbf_kernel(X_train, X_test, length_scale, variance)
K_ss = rbf_kernel(X_test, X_test, length_scale, variance)
K_inv = np.linalg.inv(K)
mu = K_s.T @ K_inv @ y_train # 予測平均
cov = K_ss - K_s.T @ K_inv @ K_s # 予測共分散
sigma = np.sqrt(np.maximum(np.diag(cov), 0))
return mu.ravel(), sigma
# --- EI獲得関数 ---
def expected_improvement(mu, sigma, f_best):
"""Expected Improvement を計算する（最小化問題）"""
with np.errstate(divide='ignore', invalid='ignore'):
gamma = (f_best - mu) / sigma
ei = (f_best - mu) * norm.cdf(gamma) + sigma * norm.pdf(gamma)
ei[sigma < 1e-10] = 0.0
return ei
# --- 目的関数（最小化対象のテスト関数） ---
def objective(x):
"""複数の極値を持つ関数"""
return np.sin(3 * x) + x**2 * 0.1 - 0.5 * np.sin(7 * x)
ベイズ最適化の実行と可視化
上記の関数を使って、ベイズ最適化を6ステップ実行し、各ステップの代理モデルと EI を可視化します。
np.random.seed(42)
# 探索範囲
x_min, x_max = -3.0, 3.0
X_grid = np.linspace(x_min, x_max, 500).reshape(-1, 1)
y_true = objective(X_grid).ravel()
# 初期点（2点からスタート）
X_train = np.array([[-2.0], [2.0]])
y_train = objective(X_train).ravel()
# GPR のハイパーパラメータ
ls, var = 0.8, 2.0
n_iterations = 6
fig, axes = plt.subplots(n_iterations, 2, figsize=(14, 4 * n_iterations))
for i in range(n_iterations):
# ガウス過程回帰で予測
mu, sigma = gpr_predict(X_train, y_train, X_grid,
length_scale=ls, variance=var)
f_best = np.min(y_train)
# EI を計算
ei = expected_improvement(mu, sigma, f_best)
# 次の評価点を決定
x_next_idx = np.argmax(ei)
x_next = X_grid[x_next_idx].reshape(1, -1)
y_next = objective(x_next).ravel()
# --- 左: 代理モデル ---
ax1 = axes[i, 0]
ax1.plot(X_grid, y_true, 'k--', alpha=0.4, label='True function')
ax1.plot(X_grid, mu, 'b-', label='GP mean')
ax1.fill_between(X_grid.ravel(), mu - 2 * sigma, mu + 2 * sigma,
alpha=0.2, color='blue', label='GP ±2σ')
ax1.scatter(X_train, y_train, c='red', s=60, zorder=5,
label='Observations')
ax1.axvline(x=x_next.item(), color='green', ls=':', lw=1.5,
label=f'Next: x={x_next.item():.2f}')
ax1.set_title(f'Step {i + 1}: Surrogate Model')
ax1.set_xlabel('x')
ax1.set_ylabel('f(x)')
ax1.legend(loc='upper left', fontsize=8)
ax1.set_xlim(x_min, x_max)
# --- 右: EI獲得関数 ---
ax2 = axes[i, 1]
ax2.plot(X_grid, ei, 'g-', linewidth=2)
ax2.fill_between(X_grid.ravel(), 0, ei, alpha=0.3, color='green')
ax2.axvline(x=x_next.item(), color='green', ls=':', lw=1.5,
label=f'Next: x={x_next.item():.2f}')
ax2.set_title(f'Step {i + 1}: Expected Improvement')
ax2.set_xlabel('x')
ax2.set_ylabel('EI(x)')
ax2.legend(loc='upper right', fontsize=8)
ax2.set_xlim(x_min, x_max)
# データに追加
X_train = np.vstack([X_train, x_next])
y_train = np.append(y_train, y_next)
plt.tight_layout()
plt.savefig('bayesian_optimization_steps.png', dpi=150,
bbox_inches='tight')
plt.show()
# 最終結果
best_idx = np.argmin(y_train)
print(f"最良の観測点: x = {X_train[best_idx].item():.4f}")
print(f"最良の目的関数値: f(x) = {y_train[best_idx]:.4f}")
このコードを実行すると、6ステップの可視化が得られます。各行の左が代理モデル（ガウス過程の予測平均と信頼区間）、右が EI 獲得関数です。緑色の破線が次に評価される点を示しています。
ステップごとの挙動を観察すると、ベイズ最適化の特徴が見て取れます。
- 初期ステップ: 不確実性が高い領域（信頼区間が広い場所）を優先的に探索する。EI の第2項 $\sigma \phi(\gamma)$ が支配的に効いている
- 中盤ステップ: 良い値が見つかった付近と、まだ未知の領域の両方にEIのピークが現れ、探索と活用がバランスする
- 後半ステップ: 代理モデルが真の関数をよく近似し、最小値付近を精密に探索する。EI の第1項 $(f^+ – \mu)\Phi(\gamma)$ が支配的になる
まとめ
本記事では、ベイズ最適化の理論と Python 実装について解説しました。
- ベイズ最適化は代理モデル（ガウス過程回帰）と獲得関数の2つの要素からなるフレームワークである
- 獲得関数 PI は改善確率を測る最もシンプルな指標であるが、活用に偏りやすい
- 獲得関数 EI は改善量の期待値であり、$\text{EI} = (f^+ – \mu)\Phi(\gamma) + \sigma\phi(\gamma)$ という解析解を持つ。第1項が活用、第2項が探索に対応し、探索と活用のバランスを自然に取る
- 獲得関数 UCB は $\mu + \kappa\sigma$ の形で、ハイパーパラメータ $\kappa$ により探索と活用のバランスを明示的に制御できる
- Python でのスクラッチ実装により、ベイズ最適化が少ない評価回数で効率的に最適解を見つける様子を確認した
実務でベイズ最適化を利用する場合は、scikit-optimize の gp_minimize や Optuna、BoTorch などのライブラリが便利です。これらのライブラリでは多次元探索空間、カーネルの自動選択、並列評価などが実装されており、大規模な問題にも対応できます。本記事で学んだ理論的基盤があれば、ライブラリのパラメータ設定やトラブルシューティングにおいても正しく判断できるようになるでしょう。
次のステップとして、以下の記事も参考にしてください。