機械学習のモデルを構築するとき、学習率やドロップアウト率、正則化の強さなど、調整すべきハイパーパラメータが数多くあります。ニューラルネットワークの学習には数時間から数日かかることもあり、あらゆる組み合わせを総当たりで試すグリッドサーチは現実的ではありません。たとえば、5つのハイパーパラメータについてそれぞれ10個の候補値を試すだけで $10^5 = 100{,}000$ 回の評価が必要になり、1回に1時間かかるとすると約11年を要します。
一方、ランダムサーチは候補を無作為に選ぶため効率は改善しますが、「過去の試行結果を参考にして次の試行を決める」という知的な戦略がありません。人間の熟練エンジニアは、過去の実験結果を踏まえて「このあたりのパラメータが良さそうだ」と推測し、同時に「まだ試していない領域にもっと良い値があるかもしれない」と探索的な判断もします。
ベイズ最適化（Bayesian Optimization）は、この人間の直感的な最適化プロセスを数学的に定式化したものです。ガウス過程（Gaussian Process）で目的関数のサロゲートモデル（代理モデル）を構築し、獲得関数（acquisition function）で「次にどこを評価すべきか」を決定します。少ない評価回数で効率的に最適解に近づけるのが最大の特長です。
ベイズ最適化を理解すると、以下のような応用が開けます。
- ハイパーパラメータ最適化: ニューラルネットワークの学習率、層数、ユニット数などの自動調整。Optuna や HyperOpt の背後にある理論
- 実験計画: 材料探索や化学合成の条件最適化。1回の実験が高コストな場面で威力を発揮
- ロボティクス: ロボットの制御パラメータのチューニング。実機実験は壊れるリスクもあるため、少ない試行で最適化したい
- シミュレーション最適化: CFD（計算流体力学）やFEM（有限要素法）のパラメータ最適化
本記事の内容
- ベイズ最適化の全体像と動機
- ガウス過程回帰の復習
- 獲得関数（PI, EI, UCB）の定義と導出
- EIの解析解の導出
- Pythonでのスクラッチ実装（1次元関数の最適化）
- グリッドサーチ・ランダムサーチとの比較
- 多次元への拡張と実用上の注意点
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ガウス過程回帰を初めから分かりやすく解説 — ガウス過程の基礎理論
- 正規分布とは？定義やグラフの見方を徹底解説 — 正規分布の性質
- ベイズ推定とは？仕組みについてわかりやすく解説 — ベイズ推論の基本
ベイズ最適化の全体像
フレームワーク
ベイズ最適化は以下の2つの構成要素から成ります。
- サロゲートモデル（代理モデル）: 少ない評価点から目的関数の全体像を推定する確率的モデル。一般にガウス過程（GP）が使われます
- 獲得関数（acquisition function）: サロゲートモデルの予測を基に、「次にどこを評価すべきか」を決める基準
ベイズ最適化のアルゴリズムは非常にシンプルです。
- 初期点をいくつか評価する
- 以下を繰り返す: – (a) これまでの評価結果でサロゲートモデル（GP）を更新する – (b) 獲得関数を最大化する点 $\bm{x}_{\text{next}}$ を求める – (c) $\bm{x}_{\text{next}}$ で目的関数を評価し、データに追加する
- 最良の評価点を返す
このシンプルなループの中に、探索と活用のバランス（exploration-exploitation trade-off）という深い構造が隠されています。「活用」とは、現在最も良さそうな領域をさらに詳しく調べること。「探索」とは、まだ調べていない領域を試して、もっと良い解が隠れていないか確認することです。良い獲得関数は、この2つを自動的にバランスします。
なぜガウス過程が最適なのか
サロゲートモデルにはさまざまな選択肢がありますが、ガウス過程が最も広く使われる理由は以下の3つです。
- 予測の不確実性が得られる: GPは予測平均だけでなく予測分散（不確実性）も提供します。獲得関数はこの不確実性情報を必要とします
- 少数データでも良い推定ができる: GPはノンパラメトリックモデルであり、データ数に応じて柔軟に関数の複雑さを調整します
- 解析的に扱いやすい: GPの予測分布は正規分布であり、獲得関数の期待値を解析的に計算できます
では、ガウス過程回帰の基本を簡潔に復習しましょう。
ガウス過程回帰の復習
基本定義
ガウス過程（Gaussian Process, GP）は関数の確率分布であり、平均関数 $m(\bm{x})$ とカーネル関数 $k(\bm{x}, \bm{x}’)$ で特徴づけられます。
$$ f(\bm{x}) \sim \mathcal{GP}(m(\bm{x}), k(\bm{x}, \bm{x}’)) $$
任意の有限個の入力点 $\bm{x}_1, \dots, \bm{x}_N$ に対して、対応する関数値 $f(\bm{x}_1), \dots, f(\bm{x}_N)$ は多変量正規分布に従います。
$$ \begin{pmatrix} f(\bm{x}_1) \\ \vdots \\ f(\bm{x}_N) \end{pmatrix} \sim \mathcal{N}\left(\begin{pmatrix} m(\bm{x}_1) \\ \vdots \\ m(\bm{x}_N) \end{pmatrix}, \begin{pmatrix} k(\bm{x}_1, \bm{x}_1) & \cdots & k(\bm{x}_1, \bm{x}_N) \\ \vdots & \ddots & \vdots \\ k(\bm{x}_N, \bm{x}_1) & \cdots & k(\bm{x}_N, \bm{x}_N) \end{pmatrix}\right) $$
簡単のため平均関数 $m(\bm{x}) = 0$ とし、カーネル関数にはRBF（ガウス）カーネルを使います。
$$ \begin{equation} k(\bm{x}, \bm{x}’) = \sigma_f^2 \exp\left(-\frac{|\bm{x} – \bm{x}’|^2}{2\ell^2}\right) \end{equation} $$
ここで $\sigma_f^2$ は出力の分散スケール、$\ell$ は長さスケール（関数の滑らかさを制御するパラメータ）です。
予測分布
$N$ 個の観測データ $\mathcal{D} = \{(\bm{x}_i, y_i)\}_{i=1}^{N}$（$y_i = f(\bm{x}_i) + \epsilon_i$, $\epsilon_i \sim \mathcal{N}(0, \sigma_n^2)$）が与えられたとき、新しい入力 $\bm{x}_*$ での予測分布は、
$$ f(\bm{x}_*) | \mathcal{D} \sim \mathcal{N}(\mu(\bm{x}_*), \sigma^2(\bm{x}_*)) $$
ここで、
$$ \mu(\bm{x}_*) = \bm{k}_*^T (\bm{K} + \sigma_n^2 \bm{I})^{-1} \bm{y} $$
$$ \sigma^2(\bm{x}_*) = k(\bm{x}_*, \bm{x}_*) – \bm{k}_*^T (\bm{K} + \sigma_n^2 \bm{I})^{-1} \bm{k}_* $$
$\bm{K}$ は訓練データ間のカーネル行列、$\bm{k}_*$ は新しい点と訓練データ間のカーネルベクトルです。
予測平均 $\mu(\bm{x}_*)$ は訓練データの近くでは観測値に近い値をとり、予測分散 $\sigma^2(\bm{x}_*)$ は訓練データから離れるほど大きくなります。この「データがある場所では自信を持ち、データがない場所では不確実性を認める」という性質が、ベイズ最適化の探索-活用バランスの基盤となります。
ガウス過程の予測分布を使って、次に評価すべき点をどのように決めるのか。その役割を担うのが獲得関数です。
獲得関数の理論
獲得関数の役割
獲得関数 $\alpha(\bm{x})$ は、ガウス過程の予測平均 $\mu(\bm{x})$ と予測標準偏差 $\sigma(\bm{x})$ を入力として、各点 $\bm{x}$ の「評価する価値」をスカラー値で表します。次に評価する点は $\bm{x}_{\text{next}} = \arg\max_{\bm{x}} \alpha(\bm{x})$ として決定されます。
良い獲得関数は以下の2つの要素をバランスします。
- 活用（exploitation）: $\mu(\bm{x})$ が大きい点を選ぶ（現在の最良推定の近くを深掘り）
- 探索（exploration）: $\sigma(\bm{x})$ が大きい点を選ぶ（不確実な領域を調査）
代表的な獲得関数として、PI（Probability of Improvement）、EI（Expected Improvement）、UCB（Upper Confidence Bound）の3つを解説します。
PI（Probability of Improvement）
直感
PIは「現在の最良値よりも良い値が得られる確率」を獲得関数とする、最もシンプルな方法です。
数学的定義
現在の最良観測値を $f^+ = \max_{i} y_i$ とします。点 $\bm{x}$ での改善確率は、
$$ \begin{equation} \text{PI}(\bm{x}) = P(f(\bm{x}) > f^+) = \Phi\left(\frac{\mu(\bm{x}) – f^+}{\sigma(\bm{x})}\right) \end{equation} $$
ここで $\Phi(\cdot)$ は標準正規分布の累積分布関数です。
この式の導出は素直です。GPの予測分布の下で $f(\bm{x}) \sim \mathcal{N}(\mu(\bm{x}), \sigma^2(\bm{x}))$ なので、$f(\bm{x}) > f^+$ となる確率は、標準化変量 $Z = (f(\bm{x}) – \mu(\bm{x}))/\sigma(\bm{x})$ が $(f^+ – \mu(\bm{x}))/\sigma(\bm{x})$ より大きい確率、すなわち $\Phi((\mu(\bm{x}) – f^+)/\sigma(\bm{x}))$ です。
PIの問題点
PIは「改善の確率」を最大化しますが、「改善の大きさ」を考慮しません。たとえば、ある点で改善確率が80%だが改善量は微小な場合と、別の点で改善確率が30%だが大幅な改善が期待できる場合、PIは前者を選びます。これは多くの場面で非効率です。PIは活用に偏りすぎる傾向があるのです。
この問題を解決するのがEI（Expected Improvement）です。
EI（Expected Improvement）
直感
EIは「改善量の期待値」を獲得関数とします。改善確率だけでなく、「どれだけ改善するか」も考慮するため、PIよりもバランスの取れた探索が可能です。
数学的定義と導出
改善量（Improvement）を $I(\bm{x}) = \max(f(\bm{x}) – f^+, 0)$ と定義します。EIはこの改善量の期待値です。
$$ \text{EI}(\bm{x}) = \mathbb{E}[\max(f(\bm{x}) – f^+, 0)] $$
$f(\bm{x}) \sim \mathcal{N}(\mu, \sigma^2)$（ここでは表記を簡略化して $\mu = \mu(\bm{x})$, $\sigma = \sigma(\bm{x})$ とします）の下でこの期待値を計算します。
$$ \text{EI}(\bm{x}) = \int_{f^+}^{\infty} (f – f^+) \cdot \frac{1}{\sigma}\phi\left(\frac{f – \mu}{\sigma}\right) df $$
ここで $\phi(\cdot)$ は標準正規分布の確率密度関数です。$z = (f – \mu)/\sigma$ と変数変換すると、$f = \mu + \sigma z$, $df = \sigma \, dz$ であり、$f = f^+$ のとき $z = (f^+ – \mu)/\sigma$ です。$Z = (\mu – f^+)/\sigma$ と置くと、
$$ \text{EI}(\bm{x}) = \int_{-Z}^{\infty} (\mu + \sigma z – f^+) \phi(z) \, dz $$
この積分を2つに分解します。
$$ \text{EI}(\bm{x}) = (\mu – f^+) \int_{-Z}^{\infty} \phi(z) \, dz + \sigma \int_{-Z}^{\infty} z \phi(z) \, dz $$
第1項の積分は $\Phi(Z)$ です。
第2項の積分を計算しましょう。$z\phi(z) = -\phi'(z)$（$\phi(z) = \frac{1}{\sqrt{2\pi}}e^{-z^2/2}$ の微分）を使うと、
$$ \int_{-Z}^{\infty} z \phi(z) \, dz = \left[-\phi(z)\right]_{-Z}^{\infty} = \phi(-Z) = \phi(Z) $$
最後の等号は $\phi$ の偶関数性（$\phi(-z) = \phi(z)$）によります。
以上をまとめると、EIの解析解が得られます。
$$ \begin{equation} \text{EI}(\bm{x}) = (\mu(\bm{x}) – f^+)\Phi(Z) + \sigma(\bm{x})\phi(Z), \quad Z = \frac{\mu(\bm{x}) – f^+}{\sigma(\bm{x})} \end{equation} $$
EIの解析解の解釈
この式の2つの項には明確な意味があります。
- 第1項 $(\mu – f^+)\Phi(Z)$: 活用（exploitation）を表します。予測平均が現在の最良値を上回る度合い $(\mu – f^+)$ に、改善が起こる確率 $\Phi(Z)$ を掛けたもの
- 第2項 $\sigma \phi(Z)$: 探索（exploration）を表します。予測の不確実性 $\sigma$ が大きいほど、偶然良い値が得られる可能性が高まることを反映しています
$\sigma(\bm{x}) = 0$ のとき（すでに観測済みの点）、EI = 0 になることも確認できます。また、$\sigma(\bm{x}) \to \infty$ のとき EI $\to \infty$ となるので、完全に未知の領域は必ず探索されます。
UCB（Upper Confidence Bound）
直感
UCBは予測の「楽観的な上限」を獲得関数とします。これは「最も楽観的に考えたときに最良の結果が期待できる点を選ぶ」という戦略です。
数学的定義
$$ \begin{equation} \text{UCB}(\bm{x}) = \mu(\bm{x}) + \kappa \sigma(\bm{x}) \end{equation} $$
ここで $\kappa > 0$ は探索と活用のバランスを制御するパラメータです。$\kappa$ を大きくすると探索重視、小さくすると活用重視になります。
UCBの利点は、$\kappa$ のスケジューリングに関する理論的な結果が得られていることです。GP-UCBアルゴリズムでは、$\kappa$ を反復回数 $t$ に応じて適切に増やすことで、累積後悔（regret）の理論的な上界が保証されます。
3つの獲得関数の比較
| 獲得関数 | 長所 | 短所 | 調整パラメータ |
|---|---|---|---|
| PI | 計算が簡単 | 活用に偏りすぎる | なし（$\xi$ で調整可） |
| EI | バランスが良い、解析解あり | 探索がやや弱い場合がある | なし（$\xi$ で調整可） |
| UCB | 理論的保証がある | $\kappa$ の設定が必要 | $\kappa$ |
実用上はEIが最も広く使われています。バランスの良さに加え、解析解があるため計算効率も高いです。
理論が整ったところで、Pythonでベイズ最適化をスクラッチ実装してみましょう。
Pythonでのスクラッチ実装
ガウス過程と獲得関数の実装
1次元の目的関数を最小化する問題（符号を反転して最大化問題として扱う）をベイズ最適化で解きます。
import numpy as np
from scipy.stats import norm
from scipy.optimize import minimize_scalar
class GaussianProcess:
"""ガウス過程回帰"""
def __init__(self, length_scale=1.0, signal_var=1.0, noise_var=1e-6):
self.l = length_scale
self.sf2 = signal_var
self.sn2 = noise_var
self.X = None
self.y = None
self.K_inv = None
def kernel(self, x1, x2):
"""RBFカーネル"""
sq_dist = np.subtract.outer(x1, x2)**2
return self.sf2 * np.exp(-sq_dist / (2 * self.l**2))
def fit(self, X, y):
"""訓練データでGPを学習"""
self.X = np.array(X)
self.y = np.array(y)
K = self.kernel(self.X, self.X) + self.sn2 * np.eye(len(self.X))
self.K_inv = np.linalg.inv(K)
def predict(self, X_new):
"""予測平均と予測標準偏差を返す"""
X_new = np.array(X_new)
k_star = self.kernel(self.X, X_new)
k_ss = self.kernel(X_new, X_new)
mu = k_star.T @ self.K_inv @ self.y
cov = k_ss - k_star.T @ self.K_inv @ k_star
std = np.sqrt(np.maximum(np.diag(cov), 1e-10))
return mu, std
上のコードでは、RBFカーネルを用いたガウス過程回帰のクラスを定義しました。predict メソッドは予測平均と予測標準偏差を返し、これらが獲得関数の入力となります。
次に、3つの獲得関数を実装します。
def acquisition_ei(mu, std, f_best, xi=0.01):
"""Expected Improvement"""
with np.errstate(divide='warn'):
Z = (mu - f_best - xi) / std
ei = (mu - f_best - xi) * norm.cdf(Z) + std * norm.pdf(Z)
ei[std < 1e-10] = 0.0
return ei
def acquisition_pi(mu, std, f_best, xi=0.01):
"""Probability of Improvement"""
Z = (mu - f_best - xi) / std
pi = norm.cdf(Z)
pi[std < 1e-10] = 0.0
return pi
def acquisition_ucb(mu, std, kappa=2.0):
"""Upper Confidence Bound"""
return mu + kappa * std
獲得関数の実装は EI の解析解をそのまま反映しています。$\xi$ は微小な定数で、すでに最良値と同程度の改善を避ける効果があり、探索をわずかに促進します。
ここまでの部品を組み合わせて、ベイズ最適化のメインループを実装しましょう。
ベイズ最適化ループの実装と可視化
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
class GaussianProcess:
def __init__(self, length_scale=1.0, signal_var=1.0, noise_var=1e-6):
self.l = length_scale
self.sf2 = signal_var
self.sn2 = noise_var
def kernel(self, x1, x2):
sq_dist = np.subtract.outer(x1, x2)**2
return self.sf2 * np.exp(-sq_dist / (2 * self.l**2))
def fit(self, X, y):
self.X = np.array(X)
self.y = np.array(y)
K = self.kernel(self.X, self.X) + self.sn2 * np.eye(len(self.X))
self.K_inv = np.linalg.inv(K)
def predict(self, X_new):
X_new = np.array(X_new)
k_star = self.kernel(self.X, X_new)
k_ss = self.kernel(X_new, X_new)
mu = k_star.T @ self.K_inv @ self.y
cov = k_ss - k_star.T @ self.K_inv @ k_star
std = np.sqrt(np.maximum(np.diag(cov), 1e-10))
return mu, std
def acquisition_ei(mu, std, f_best, xi=0.01):
Z = (mu - f_best - xi) / std
ei = (mu - f_best - xi) * norm.cdf(Z) + std * norm.pdf(Z)
ei[std < 1e-10] = 0.0
return ei
# 目的関数（最大化する）
def objective(x):
return -(x - 2)**2 * np.sin(3 * x + 1) + 2
# 最適化ループ
np.random.seed(42)
bounds = (0, 5)
x_grid = np.linspace(bounds[0], bounds[1], 500)
# 初期点
X_observed = np.array([0.5, 2.5, 4.5])
y_observed = objective(X_observed)
gp = GaussianProcess(length_scale=0.5, signal_var=2.0, noise_var=1e-4)
n_iterations = 8
fig, axes = plt.subplots(4, 2, figsize=(14, 20))
for iteration in range(n_iterations):
ax = axes[iteration // 2, iteration % 2]
# GPフィッティング
gp.fit(X_observed, y_observed)
mu, std = gp.predict(x_grid)
# EI計算
f_best = np.max(y_observed)
ei = acquisition_ei(mu, std, f_best)
# 次の評価点
x_next_idx = np.argmax(ei)
x_next = x_grid[x_next_idx]
# プロット: 目的関数とGP
ax.plot(x_grid, objective(x_grid), 'k--', linewidth=1, alpha=0.5, label='目的関数')
ax.plot(x_grid, mu, 'b-', linewidth=2, label='GP予測平均')
ax.fill_between(x_grid, mu - 1.96*std, mu + 1.96*std,
alpha=0.15, color='blue')
ax.scatter(X_observed, y_observed, c='red', s=60, zorder=5, label='観測点')
ax.axvline(x_next, color='green', linestyle=':', linewidth=2, alpha=0.7)
# EIのスケーリング表示
ei_scaled = ei / (np.max(ei) + 1e-10) * (np.max(mu + 1.96*std) - np.min(mu - 1.96*std)) * 0.3
ax.fill_between(x_grid, np.min(mu - 1.96*std), np.min(mu - 1.96*std) + ei_scaled,
alpha=0.3, color='green', label='EI')
ax.set_title(f'反復 {iteration + 1}: 次の点 x={x_next:.2f}', fontsize=12)
ax.set_xlabel('x')
ax.set_ylabel('f(x)')
if iteration == 0:
ax.legend(fontsize=8, loc='upper right')
# 新しい点を追加
y_next = objective(x_next)
X_observed = np.append(X_observed, x_next)
y_observed = np.append(y_observed, y_next)
plt.tight_layout()
plt.savefig("bayesian_optimization_iterations.png", dpi=150, bbox_inches="tight")
plt.show()
print(f"\n最終結果:")
print(f" 最良の x = {X_observed[np.argmax(y_observed)]:.4f}")
print(f" 最良の f(x) = {np.max(y_observed):.4f}")
8回の反復を通じて、ベイズ最適化の探索過程を視覚的に追跡できます。
-
初期段階（反復1-2）: GPは3つの初期点から目的関数の大まかな形を推定しています。予測の不確実性（青い帯）がデータのない領域で大きく、EI（緑の領域）もこの不確実な領域で高い値を示しています。最適化は探索を優先しています
-
中期段階（反復3-5）: 探索により得られた情報でGPの推定が改善され、目的関数の形状がより正確に捉えられています。EIは最適値付近（活用）と未探索領域（探索）にピークを持ち、自動的にバランスを取っています
-
後期段階（反復6-8）: GPの推定は目的関数にかなり近づいており、EIは最適値付近に集中しています。予測の不確実性が全体的に小さくなり、最適解の周辺を精密に探索する段階に移行しています
特筆すべきは、ベイズ最適化がたった8回の評価で目的関数の最大値をかなり正確に見つけていることです。目的関数は複雑な非線形関数ですが、GPによるサロゲートモデルとEIによる賢い探索点選択により、効率的な最適化が実現されています。
次に、ベイズ最適化の効率性をグリッドサーチやランダムサーチと定量的に比較してみましょう。
グリッドサーチ・ランダムサーチとの比較
3つの最適化手法の比較実験
同じ目的関数に対して、ベイズ最適化、グリッドサーチ、ランダムサーチの性能を比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
class GaussianProcess:
def __init__(self, length_scale=0.5, signal_var=2.0, noise_var=1e-4):
self.l = length_scale
self.sf2 = signal_var
self.sn2 = noise_var
def kernel(self, x1, x2):
sq_dist = np.subtract.outer(x1, x2)**2
return self.sf2 * np.exp(-sq_dist / (2 * self.l**2))
def fit(self, X, y):
self.X = np.array(X)
self.y = np.array(y)
K = self.kernel(self.X, self.X) + self.sn2 * np.eye(len(self.X))
self.K_inv = np.linalg.inv(K)
def predict(self, X_new):
X_new = np.array(X_new)
k_star = self.kernel(self.X, X_new)
k_ss = self.kernel(X_new, X_new)
mu = k_star.T @ self.K_inv @ self.y
cov = k_ss - k_star.T @ self.K_inv @ k_star
std = np.sqrt(np.maximum(np.diag(cov), 1e-10))
return mu, std
def acquisition_ei(mu, std, f_best, xi=0.01):
Z = (mu - f_best - xi) / std
ei = (mu - f_best - xi) * norm.cdf(Z) + std * norm.pdf(Z)
ei[std < 1e-10] = 0.0
return ei
def objective(x):
return -(x - 2)**2 * np.sin(3 * x + 1) + 2
# 真の最大値を見つける
x_fine = np.linspace(0, 5, 10000)
true_max = np.max(objective(x_fine))
n_trials = 50
n_evals = 20
bounds = (0, 5)
results_bo = np.zeros((n_trials, n_evals))
results_random = np.zeros((n_trials, n_evals))
results_grid = np.zeros((n_trials, n_evals))
for trial in range(n_trials):
np.random.seed(trial)
# --- ベイズ最適化 ---
X_bo = np.random.uniform(bounds[0], bounds[1], 3)
y_bo = objective(X_bo)
gp = GaussianProcess()
x_grid_bo = np.linspace(bounds[0], bounds[1], 500)
for i in range(n_evals):
if i < 3:
results_bo[trial, i] = np.max(y_bo[:i+1])
else:
gp.fit(X_bo[:i], y_bo[:i])
mu, std = gp.predict(x_grid_bo)
f_best = np.max(y_bo[:i])
ei = acquisition_ei(mu, std, f_best)
x_next = x_grid_bo[np.argmax(ei)]
X_bo = np.append(X_bo, x_next) if len(X_bo) <= i else X_bo
y_bo = np.append(y_bo, objective(x_next)) if len(y_bo) <= i else y_bo
results_bo[trial, i] = np.max(y_bo[:i+1])
# --- ランダムサーチ ---
X_rand = np.random.uniform(bounds[0], bounds[1], n_evals)
y_rand = objective(X_rand)
for i in range(n_evals):
results_random[trial, i] = np.max(y_rand[:i+1])
# --- グリッドサーチ ---
X_grid = np.linspace(bounds[0], bounds[1], n_evals)
np.random.shuffle(X_grid) # 順序をランダムに
y_grid = objective(X_grid)
for i in range(n_evals):
results_grid[trial, i] = np.max(y_grid[:i+1])
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: 最良値の推移
ax = axes[0]
evals_range = np.arange(1, n_evals + 1)
for results, label, color in [
(results_bo, 'ベイズ最適化 (EI)', 'blue'),
(results_random, 'ランダムサーチ', 'orange'),
(results_grid, 'グリッドサーチ', 'green'),
]:
mean = np.mean(results, axis=0)
std = np.std(results, axis=0)
ax.plot(evals_range, mean, color=color, linewidth=2, label=label)
ax.fill_between(evals_range, mean - std, mean + std, color=color, alpha=0.15)
ax.axhline(true_max, color='black', linestyle=':', linewidth=1, label=f'真の最大値 ({true_max:.2f})')
ax.set_xlabel('評価回数')
ax.set_ylabel('これまでの最良値')
ax.set_title('最適化の収束比較 (50回平均)')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# 右: 最終結果の箱ひげ図
ax = axes[1]
final_results = [results_bo[:, -1], results_random[:, -1], results_grid[:, -1]]
bp = ax.boxplot(final_results, labels=['ベイズ最適化', 'ランダム', 'グリッド'],
patch_artist=True)
colors_bp = ['lightblue', 'lightyellow', 'lightgreen']
for patch, color in zip(bp['boxes'], colors_bp):
patch.set_facecolor(color)
ax.axhline(true_max, color='black', linestyle=':', linewidth=1, label='真の最大値')
ax.set_ylabel('最終最良値')
ax.set_title(f'{n_evals}回評価後の最良値分布')
ax.legend()
plt.tight_layout()
plt.savefig("bo_vs_random_vs_grid.png", dpi=150, bbox_inches="tight")
plt.show()
比較実験の結果から、以下のことが明らかになります。
-
収束速度: ベイズ最適化は他の2つの手法よりも明確に速く最適値に収束しています。特に5-10回の評価で大きな差が生まれており、評価コストが高い場面でのベイズ最適化の優位性が確認できます
-
最終的な精度: 20回の評価後の最良値の分布（箱ひげ図）を見ると、ベイズ最適化は中央値が高く、ばらつきも小さいことがわかります。これは、ベイズ最適化が安定して良い解を見つけられることを意味しています
-
ランダムサーチ vs グリッドサーチ: 1次元の問題ではグリッドサーチが比較的良い性能を示しますが、多次元になるとグリッドの密度が指数的に薄くなるため、ランダムサーチの方が有利になることが知られています
次に、多次元への拡張と実用上の注意点をまとめましょう。
多次元への拡張と実用上の注意点
多次元GPの扱い
1次元から多次元への拡張は概念的には単純で、カーネル関数の引数がスカラーからベクトルに変わるだけです。ARD（Automatic Relevance Determination）カーネルを使うと、各次元ごとに異なる長さスケールを学習できます。
$$ \begin{equation} k(\bm{x}, \bm{x}’) = \sigma_f^2 \exp\left(-\sum_{d=1}^{D} \frac{(x_d – x_d’)^2}{2\ell_d^2}\right) \end{equation} $$
ここで $\ell_d$ は第 $d$ 次元の長さスケールです。$\ell_d$ が大きい次元は目的関数への影響が小さく、$\ell_d$ が小さい次元は影響が大きいことを意味します。この性質は自動的な特徴選択として機能します。
高次元での課題
ベイズ最適化には次元の呪いがあり、一般に $D \leq 20$ 程度の問題に適しています。その理由は以下の通りです。
- GPの計算コスト: $N$ 個のデータに対してカーネル行列の逆行列を計算するため $O(N^3)$ の計算量がかかる
- 獲得関数の最大化: 高次元空間での最大化は困難。勾配法やマルチスタートの最適化が必要
- サロゲートモデルの精度: 高次元空間を正確にモデル化するには多くのデータが必要
高次元問題への対処としては、ランダム埋め込み（REMBO）、部分空間への射影、加法的GPモデルなどの手法が研究されています。
カーネルハイパーパラメータの最適化
実用上最も重要なのが、GPのカーネルハイパーパラメータ（$\sigma_f^2$, $\ell$, $\sigma_n^2$）の設定です。適切なハイパーパラメータを設定しないと、GPのサロゲートモデルが目的関数を正しく近似できず、最適化の効率が大幅に低下します。
最も一般的な方法は、対数周辺尤度の最大化（Type-II最尤推定）です。
$$ \log p(\bm{y} | \bm{X}, \theta) = -\frac{1}{2}\bm{y}^T(\bm{K} + \sigma_n^2\bm{I})^{-1}\bm{y} – \frac{1}{2}\log|\bm{K} + \sigma_n^2\bm{I}| – \frac{N}{2}\log(2\pi) $$
この最大化はGPを更新するたびに（ベイズ最適化の各反復で）実行するのが理想的ですが、計算コストとのトレードオフがあります。
実用上のヒント
- 初期点: ラテン超方格サンプリング（LHS）で均等に分散した初期点を選ぶ
- 正規化: 入力と出力を適切にスケーリングする（各次元を [0, 1] に正規化するなど）
- ノイズへの対応: 目的関数にノイズがある場合、$\sigma_n^2 > 0$ を設定する
- 獲得関数の選択: 一般にはEIで十分。探索が不足する場合はUCBに切り替える
- 終了条件: EIの最大値が十分小さくなったら終了
まとめ
本記事では、ベイズ最適化の理論をガウス過程と獲得関数の観点から解説しました。
- ベイズ最適化の構造: サロゲートモデル（ガウス過程）で目的関数を近似し、獲得関数で次に評価する点を賢く選ぶ
- ガウス過程: 関数の確率分布を提供し、予測平均（最良の推定）と予測分散（不確実性）の両方を得られる
- 獲得関数: PI（改善確率）、EI（期待改善量）、UCB（楽観的上限）の3つを導出した。EIの解析解は活用項と探索項に明確に分解される
- 比較実験: ベイズ最適化はグリッドサーチやランダムサーチよりも収束が速く、少ない評価回数で良い解に到達する
- 多次元への拡張: ARDカーネルによる自動的な特徴選択が可能だが、20次元程度が実用的な上限
ベイズ最適化は「評価コストの高い関数の最適化」という実用的な問題を、ベイズ統計の美しい理論で解決する手法です。ガウス過程による不確実性の定量化と、獲得関数による探索-活用バランスの自動調整が、その効率性の鍵となっています。
次のステップとして、以下の記事も参考にしてください。
- ガウス過程回帰を初めから分かりやすく解説 — GPの理論をより深く理解する
- ベイズ予測分布とは？事後分布から未知データを予測する理論と実装 — 予測分布の一般理論
- PyMCで学ぶ確率的プログラミング — ベイズモデリングの実践入門 — ベイズ推論を手軽に実装する