ガウス過程回帰（Gaussian Process Regression, GPR）は、関数の事前分布を直接扱う柔軟な非パラメトリック回帰手法です。しかし、その予測性能はハイパーパラメータの選び方で大きく変わります。長さスケールが小さすぎればデータのノイズに過剰に適合し、大きすぎればデータの構造を捉えられません。
では、最適なハイパーパラメータをどのように決めればよいのでしょうか。ガウス過程では、周辺尤度（marginal likelihood） を最大化するという原理的に美しい方法が使えます。これはベイズ的なモデル選択の考え方に基づいており、データへの適合と過剰な複雑さの間で自動的にバランスを取ってくれます。
本記事の内容
- ガウス過程回帰におけるハイパーパラメータの役割
- 周辺尤度 $p(\bm{y} | X, \bm{\theta})$ の導出
- 対数周辺尤度の3つの項の解釈（データ適合・複雑さペナルティ・定数）
- 対数周辺尤度の勾配の導出
- Python での実装（手動実装と scikit-learn の比較）
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ハイパーパラメータとは
ガウス過程回帰では、カーネル関数と観測ノイズの大きさがモデルの振る舞いを決定します。最もよく使われるカーネル関数であるRBFカーネル（ガウスカーネル、二乗指数カーネルとも呼ばれる） を例にとると、カーネル関数は次のように定義されます。
$$ k(\bm{x}, \bm{x’}) = \sigma_f^2 \exp\left( -\frac{\|\bm{x} – \bm{x’}\|^2}{2\ell^2} \right) $$
ここで、各パラメータの意味は以下のとおりです。
| パラメータ | 記号 | 意味 |
|---|---|---|
| 長さスケール | $\ell$ | 入力空間でどのくらいの距離のデータ点が互いに相関をもつかを制御する |
| 信号分散 | $\sigma_f^2$ | 関数値の変動幅（振幅）を制御する |
| ノイズ分散 | $\sigma_n^2$ | 観測に含まれるノイズの大きさを表す |
これら $\bm{\theta} = (\ell, \sigma_f^2, \sigma_n^2)$ がガウス過程回帰のハイパーパラメータです。通常の回帰モデルにおけるパラメータ（重み $\bm{w}$ など）とは異なり、データから直接学習される重みではなく、モデルの構造そのものを規定する量です。
観測ノイズを含めると、観測値 $\bm{y}$ の共分散行列は次のようになります。
$$ \bm{K}_y = \bm{K} + \sigma_n^2 \bm{I} $$
ここで、$\bm{K}$ はカーネル行列であり、$(i, j)$ 成分が $k(\bm{x}_i, \bm{x}_j)$ で与えられる $N \times N$ の行列です。$\bm{I}$ は $N \times N$ の単位行列です。
長さスケールの影響
長さスケール $\ell$ はガウス過程回帰の予測に最も大きな影響を与えるハイパーパラメータです。
- $\ell$ が小さい場合: カーネル関数の値が距離に対して急速に減衰するため、近くのデータ点のみが互いに影響します。結果として、予測関数は激しく振動し、データ点に過剰に適合します。
- $\ell$ が大きい場合: 遠くのデータ点まで影響が及ぶため、予測関数は滑らかになります。しかし、大きすぎるとデータの局所的な構造を捉えられなくなります。
このように、ハイパーパラメータの選択はガウス過程回帰の性能を左右するため、適切な最適化手法が必要です。
周辺尤度 $p(\bm{y}|X,\bm{\theta})$ の導出
ガウス過程回帰のハイパーパラメータ最適化の鍵となるのが、周辺尤度（marginal likelihood） です。「周辺」と呼ばれるのは、潜在関数値 $\bm{f}$ を積分消去（周辺化）して得られる尤度だからです。エビデンス（evidence）とも呼ばれます。
ベイズ推論の観点
ベイズ推論の枠組みでは、ハイパーパラメータの選択はモデル選択の問題として定式化されます。データ $\bm{y}$ が与えられたとき、ハイパーパラメータの事後分布はベイズの定理から次のように書けます。
$$ p(\bm{\theta} | \bm{y}, X) = \frac{p(\bm{y} | X, \bm{\theta}) \, p(\bm{\theta})}{p(\bm{y} | X)} $$
この式の分子に現れる $p(\bm{y} | X, \bm{\theta})$ が周辺尤度です。事前分布 $p(\bm{\theta})$ が一様分布であれば、事後分布を最大にすることは周辺尤度を最大にすることと等価です。
ガウス積分による解析的な導出
ガウス過程回帰では、関数 $f(\bm{x})$ にガウス過程の事前分布を置きます。
$$ f(\bm{x}) \sim \mathcal{GP}(0, k(\bm{x}, \bm{x’})) $$
$N$ 個のデータ点 $\bm{x}_1, \bm{x}_2, \dots, \bm{x}_N$ に対する関数値 $\bm{f} = (f(\bm{x}_1), \dots, f(\bm{x}_N))^T$ は、多変量ガウス分布に従います。
$$ \bm{f} | X, \bm{\theta} \sim \mathcal{N}(\bm{0}, \bm{K}) $$
観測にガウスノイズが加わるとすると、
$$ \bm{y} = \bm{f} + \bm{\varepsilon}, \quad \bm{\varepsilon} \sim \mathcal{N}(\bm{0}, \sigma_n^2 \bm{I}) $$
周辺尤度は、関数値 $\bm{f}$ を積分消去することで得られます。
$$ p(\bm{y} | X, \bm{\theta}) = \int p(\bm{y} | \bm{f}, X) \, p(\bm{f} | X, \bm{\theta}) \, d\bm{f} $$
被積分関数の各因子を明示的に書き下すと、
$$ p(\bm{y} | \bm{f}, X) = \mathcal{N}(\bm{y} | \bm{f}, \sigma_n^2 \bm{I}) = \frac{1}{(2\pi)^{N/2} |\sigma_n^2 \bm{I}|^{1/2}} \exp\left(-\frac{1}{2}(\bm{y} – \bm{f})^T (\sigma_n^2 \bm{I})^{-1} (\bm{y} – \bm{f})\right) $$
$$ p(\bm{f} | X, \bm{\theta}) = \mathcal{N}(\bm{f} | \bm{0}, \bm{K}) = \frac{1}{(2\pi)^{N/2} |\bm{K}|^{1/2}} \exp\left(-\frac{1}{2} \bm{f}^T \bm{K}^{-1} \bm{f}\right) $$
積分の中の指数部分を $\bm{f}$ について整理します。
$$ \begin{align} &-\frac{1}{2}\left[(\bm{y} – \bm{f})^T (\sigma_n^2 \bm{I})^{-1} (\bm{y} – \bm{f}) + \bm{f}^T \bm{K}^{-1} \bm{f}\right] \\ &= -\frac{1}{2}\left[\bm{f}^T \underbrace{\left(\sigma_n^{-2} \bm{I} + \bm{K}^{-1}\right)}_{\equiv \bm{A}} \bm{f} – 2\bm{f}^T \sigma_n^{-2} \bm{y} + \bm{y}^T \sigma_n^{-2} \bm{y}\right] \end{align} $$
$\bm{A} = \sigma_n^{-2} \bm{I} + \bm{K}^{-1}$ とおくと、$\bm{f}$ についての二次形式になっているので平方完成ができます。
$$ \begin{align} &= -\frac{1}{2}\left[(\bm{f} – \bm{A}^{-1}\sigma_n^{-2}\bm{y})^T \bm{A} (\bm{f} – \bm{A}^{-1}\sigma_n^{-2}\bm{y})\right] + \frac{1}{2}\bm{y}^T \sigma_n^{-2} \bm{A}^{-1} \sigma_n^{-2} \bm{y} – \frac{1}{2}\bm{y}^T \sigma_n^{-2} \bm{y} \end{align} $$
$\bm{f}$ に依存する部分はガウス分布の核（カーネル）の形をしているため、$\bm{f}$ についてガウス積分を実行すると、正規化定数として $(2\pi)^{N/2} |\bm{A}|^{-1/2}$ が出てきます。残りの $\bm{f}$ に依存しない項をまとめると、Woodbury の行列恒等式
$$ \bm{A}^{-1} = (\sigma_n^{-2}\bm{I} + \bm{K}^{-1})^{-1} = \bm{K} – \bm{K}(\bm{K} + \sigma_n^2 \bm{I})^{-1}\bm{K} $$
および行列式の関係
$$ |\sigma_n^2 \bm{I}| \cdot |\bm{K}| \cdot |\bm{A}|^{-1} = |\bm{K} + \sigma_n^2 \bm{I}| = |\bm{K}_y| $$
を用いることで、指数部分の $\bm{y}$ に関する二次形式が $\bm{y}^T \bm{K}_y^{-1} \bm{y}$ に帰着されます。最終的に、周辺尤度は次のように書けます。
$$ \begin{equation} p(\bm{y} | X, \bm{\theta}) = \frac{1}{(2\pi)^{N/2} |\bm{K}_y|^{1/2}} \exp\left( -\frac{1}{2} \bm{y}^T \bm{K}_y^{-1} \bm{y} \right) = \mathcal{N}(\bm{y} | \bm{0}, \bm{K}_y) \end{equation} $$
つまり、$\bm{f}$ と $\bm{\varepsilon}$ が独立なガウス分布に従うため、$\bm{y} = \bm{f} + \bm{\varepsilon}$ もガウス分布に従い、その共分散行列は $\bm{K} + \sigma_n^2 \bm{I}$ になるという結果です。ガウス過程の事前分布と尤度がどちらもガウス分布であるおかげで、潜在関数 $\bm{f}$ の積分を閉じた形で実行できたことがポイントです。
対数周辺尤度の式と3つの項の解釈
周辺尤度をそのまま最適化に使うと、指数関数の積により数値的に不安定になりやすいため、対数をとった対数周辺尤度（log marginal likelihood） を用います。対数は単調増加関数なので、最大値を与える $\bm{\theta}$ は変わりません。
$$ \begin{align} \log p(\bm{y} | X, \bm{\theta}) &= \log \left[ \frac{1}{(2\pi)^{N/2} |\bm{K}_y|^{1/2}} \exp\left( -\frac{1}{2} \bm{y}^T \bm{K}_y^{-1} \bm{y} \right) \right] \\ &= \underbrace{-\frac{1}{2} \bm{y}^T \bm{K}_y^{-1} \bm{y}}_{\text{データ適合項}} \underbrace{- \frac{1}{2} \log |\bm{K}_y|}_{\text{複雑さペナルティ項}} \underbrace{- \frac{N}{2} \log 2\pi}_{\text{定数項}} \end{align} $$
この式は3つの項に分解でき、それぞれに明確な意味があります。
第1項: データ適合項
$$ -\frac{1}{2} \bm{y}^T \bm{K}_y^{-1} \bm{y} $$
この項は、モデルがデータにどの程度適合しているかを測ります。$\bm{K}_y^{-1}$ はデータの相関構造を考慮した精度行列であり、この二次形式はマハラノビス距離の二乗に対応します。モデルがデータをよく説明できるとき（$\bm{y}$ が $\bm{K}_y$ で定まるガウス分布の中心付近にあるとき）、この項の絶対値は小さくなり、対数周辺尤度は大きくなります。
モデルを複雑にする方向（例えば $\bm{K}_y$ の固有値を大きくする方向）に動かすと、$\bm{K}_y^{-1}$ が小さくなるため、この項だけなら常に改善されます。つまり、この項だけを最大化するとモデルは際限なく複雑になります。
第2項: 複雑さペナルティ項
$$ -\frac{1}{2} \log |\bm{K}_y| $$
この項は、モデルの複雑さに対するペナルティとして機能します。$|\bm{K}_y|$ は共分散行列の行列式であり、モデルが許容する関数の「体積」を表しています。ハイパーパラメータがモデルを複雑にする方向に動くと（例えば $\ell$ が小さくなると）、$|\bm{K}_y|$ は大きくなり、このペナルティ項が対数周辺尤度を小さくする方向に作用します。
第1項のデータ適合を良くしようとすると第2項のペナルティが増え、逆にペナルティを減らそうとするとデータ適合が悪化する。この自動的なトレードオフが、周辺尤度最大化の本質です。
第3項: 定数項
$$ -\frac{N}{2} \log 2\pi $$
この項はハイパーパラメータに依存しないため、最適化の際には無視できます。多変量ガウス分布の正規化定数に由来するものです。
対数周辺尤度の勾配の導出
最適化アルゴリズム（L-BFGS法など）を適用するためには、対数周辺尤度のハイパーパラメータ $\theta_j$ に関する勾配が必要です。
行列微分の公式
まず、必要な行列微分の公式を確認しておきます。
逆行列の微分:
$$ \frac{\partial}{\partial \theta} \bm{K}_y^{-1} = -\bm{K}_y^{-1} \frac{\partial \bm{K}_y}{\partial \theta} \bm{K}_y^{-1} $$
行列式の対数の微分:
$$ \frac{\partial}{\partial \theta} \log |\bm{K}_y| = \mathrm{tr}\left( \bm{K}_y^{-1} \frac{\partial \bm{K}_y}{\partial \theta} \right) $$
各項の微分
対数周辺尤度 $\mathcal{L} = \log p(\bm{y} | X, \bm{\theta})$ の $\theta_j$ に関する偏微分を計算します。
第1項の微分:
$$ \begin{align} \frac{\partial}{\partial \theta_j} \left( -\frac{1}{2} \bm{y}^T \bm{K}_y^{-1} \bm{y} \right) &= -\frac{1}{2} \bm{y}^T \frac{\partial \bm{K}_y^{-1}}{\partial \theta_j} \bm{y} \\ &= -\frac{1}{2} \bm{y}^T \left( -\bm{K}_y^{-1} \frac{\partial \bm{K}_y}{\partial \theta_j} \bm{K}_y^{-1} \right) \bm{y} \\ &= \frac{1}{2} \bm{y}^T \bm{K}_y^{-1} \frac{\partial \bm{K}_y}{\partial \theta_j} \bm{K}_y^{-1} \bm{y} \end{align} $$
第2項の微分:
$$ \frac{\partial}{\partial \theta_j} \left( -\frac{1}{2} \log |\bm{K}_y| \right) = -\frac{1}{2} \mathrm{tr}\left( \bm{K}_y^{-1} \frac{\partial \bm{K}_y}{\partial \theta_j} \right) $$
第3項 はハイパーパラメータに依存しないので微分は0です。
勾配のまとめ
以上をまとめると、
$$ \begin{equation} \frac{\partial \mathcal{L}}{\partial \theta_j} = \frac{1}{2} \bm{y}^T \bm{K}_y^{-1} \frac{\partial \bm{K}_y}{\partial \theta_j} \bm{K}_y^{-1} \bm{y} – \frac{1}{2} \mathrm{tr}\left( \bm{K}_y^{-1} \frac{\partial \bm{K}_y}{\partial \theta_j} \right) \end{equation} $$
ここで、$\bm{\alpha} = \bm{K}_y^{-1} \bm{y}$ とおくと、第1項は $\frac{1}{2} \bm{\alpha}^T \frac{\partial \bm{K}_y}{\partial \theta_j} \bm{\alpha}$ と書けます。さらに、$\bm{\alpha} \bm{\alpha}^T$ は $N \times N$ 行列なので、トレースの性質 $\bm{a}^T \bm{B} \bm{a} = \mathrm{tr}(\bm{a}\bm{a}^T \bm{B})$ を用いて1つのトレースにまとめることができます。
$$ \begin{equation} \frac{\partial \mathcal{L}}{\partial \theta_j} = \frac{1}{2} \mathrm{tr}\left( \left( \bm{\alpha} \bm{\alpha}^T – \bm{K}_y^{-1} \right) \frac{\partial \bm{K}_y}{\partial \theta_j} \right) \end{equation} $$
この式が対数周辺尤度の勾配の最終形です。実装上は $\bm{\alpha} = \bm{K}_y^{-1} \bm{y}$ を一度計算しておけば、各ハイパーパラメータについてカーネル行列の偏微分 $\frac{\partial \bm{K}_y}{\partial \theta_j}$ さえ求めれば勾配が得られます。
RBFカーネルの偏微分
RBFカーネルの場合、各ハイパーパラメータに対する偏微分は以下のようになります。
長さスケール $\ell$ に関する偏微分:
$$ \frac{\partial k(\bm{x}_i, \bm{x}_j)}{\partial \ell} = \sigma_f^2 \exp\left( -\frac{\|\bm{x}_i – \bm{x}_j\|^2}{2\ell^2} \right) \cdot \frac{\|\bm{x}_i – \bm{x}_j\|^2}{\ell^3} = k(\bm{x}_i, \bm{x}_j) \cdot \frac{\|\bm{x}_i – \bm{x}_j\|^2}{\ell^3} $$
信号分散 $\sigma_f^2$ に関する偏微分:
$$ \frac{\partial k(\bm{x}_i, \bm{x}_j)}{\partial \sigma_f^2} = \exp\left( -\frac{\|\bm{x}_i – \bm{x}_j\|^2}{2\ell^2} \right) = \frac{k(\bm{x}_i, \bm{x}_j)}{\sigma_f^2} $$
ノイズ分散 $\sigma_n^2$ に関する偏微分:
$$ \frac{\partial (\bm{K}_y)_{ij}}{\partial \sigma_n^2} = \delta_{ij} $$
ここで、$\delta_{ij}$ はクロネッカーのデルタ（$i = j$ のとき1、それ以外は0）です。つまり、$\frac{\partial \bm{K}_y}{\partial \sigma_n^2} = \bm{I}$（単位行列）になります。
最適化の実装上のポイント
対数周辺尤度の勾配が得られたので、最適化アルゴリズムを用いてハイパーパラメータを求めます。実装上のポイントをいくつか紹介します。
コレスキー分解の利用
$\bm{K}_y^{-1}$ を直接計算するのではなく、コレスキー分解 $\bm{K}_y = \bm{L}\bm{L}^T$ を用いることで、数値的安定性と計算効率を向上させます。
- $\bm{\alpha} = \bm{K}_y^{-1}\bm{y}$ は、$\bm{L}\bm{z} = \bm{y}$ を前進代入で解き、$\bm{L}^T \bm{\alpha} = \bm{z}$ を後退代入で解くことで求まります。
- 行列式の対数は $\log|\bm{K}_y| = 2\sum_i \log L_{ii}$ で計算でき、桁あふれの心配がありません。
対数スケールでの最適化
ハイパーパラメータ $\ell, \sigma_f^2, \sigma_n^2$ はいずれも正の値をとる必要があります。$\log \ell, \log \sigma_f, \log \sigma_n$ に対して最適化を行うことで、制約なし最適化問題として扱えます。連鎖律により、対数スケールでの勾配は $\frac{\partial \mathcal{L}}{\partial \log \theta} = \theta \cdot \frac{\partial \mathcal{L}}{\partial \theta}$ で求まります。
多点初期化
対数周辺尤度は一般に多峰性（複数の極大を持つ）であるため、初期値によって異なる局所解に収束する可能性があります。複数のランダムな初期値から最適化を実行し、最も対数周辺尤度が大きい解を採用するのが実用的です。scikit-learn の GaussianProcessRegressor でも n_restarts_optimizer パラメータでこの多点初期化を設定できます。
最適化アルゴリズム
ガウス過程のハイパーパラメータ最適化では、L-BFGS-B法がよく用いられます。L-BFGS-B法は、ニュートン法を近似的に実行する準ニュートン法の一種で、勾配情報を用いて収束が速く、メモリ効率も良いという特徴があります。Pythonでは scipy.optimize.minimize の method='L-BFGS-B' で利用できます。
過学習と汎化のバランス — オッカムのかみそり的解釈
対数周辺尤度におけるデータ適合項と複雑さペナルティ項のトレードオフは、オッカムのかみそり（Occam’s razor） のベイズ的な実現と解釈できます。
オッカムのかみそりとは、「データを同程度にうまく説明できるモデルが複数あるならば、最も単純なモデルを選ぶべきだ」という原理です。周辺尤度最大化は、この原理を数学的に体現しています。
3つのモデルを考えてみましょう。
- 単純すぎるモデル（$\ell$ が非常に大きい）: ほぼ定数関数しか表現できません。複雑さペナルティは小さいですが、データへの適合が悪くなります。
- 適切な複雑さのモデル: データの構造を捉えつつ、不要な複雑さを持ちません。周辺尤度が最大になります。
- 複雑すぎるモデル（$\ell$ が非常に小さい）: 多様な関数を表現できますが、データへの適合は良くても、複雑さペナルティが大きくなります。
数学的に説明すると、周辺尤度 $p(\bm{y} | X, \bm{\theta})$ は確率密度関数であるため、データ空間全体で積分すると1になります。
$$ \int p(\bm{y} | X, \bm{\theta}) \, d\bm{y} = 1 $$
複雑なモデルは多様なデータパターンに確率を割り当てる必要があるため、各パターンに対する確率密度は薄くなります。逆に、単純なモデルは限られたパターンにしか確率を割り当てないため、合致するパターンには高い確率密度を集中できます。実際に観測されたデータに最も高い確率を割り当てるのが、「データを説明するのにちょうどよい複雑さ」のモデルです。
この性質により、周辺尤度最大化は交差検証などの外部的な評価手法を用いることなく、自動的に過学習を抑制する仕組みを内蔵しています。すべてのデータを学習に使いつつ、モデル選択も同時に行えるのがベイズ的アプローチの強みです。
Python実装
ここでは、手動でガウス過程回帰のハイパーパラメータ最適化を実装し、scikit-learn の結果と比較します。
スクラッチ実装
まず、対数周辺尤度とその勾配をスクラッチで実装し、scipy.optimize.minimize で最適化します。
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
# --- データ生成 ---
np.random.seed(42)
N = 30
X_train = np.sort(np.random.uniform(0, 5, N)).reshape(-1, 1)
y_train = np.sin(X_train).ravel() + 0.3 * np.random.randn(N)
X_test = np.linspace(-0.5, 5.5, 200).reshape(-1, 1)
# --- RBFカーネル行列の計算 ---
def rbf_kernel(X1, X2, length_scale, sigma_f):
"""RBFカーネル行列を計算する"""
dist_sq = np.sum(X1**2, axis=1).reshape(-1, 1) + \
np.sum(X2**2, axis=1).reshape(1, -1) - \
2 * X1 @ X2.T
return sigma_f**2 * np.exp(-0.5 * dist_sq / length_scale**2)
# --- 負の対数周辺尤度とその勾配を同時に計算 ---
def neg_log_marginal_likelihood(log_params, X, y):
"""
負の対数周辺尤度と勾配を返す（最小化用）
log_params: [log(ℓ), log(σ_f), log(σ_n)]
"""
ell = np.exp(log_params[0])
sigma_f = np.exp(log_params[1])
sigma_n = np.exp(log_params[2])
N = len(y)
# カーネル行列の構築
K = rbf_kernel(X, X, ell, sigma_f)
K_y = K + sigma_n**2 * np.eye(N) + 1e-8 * np.eye(N)
# コレスキー分解で安定に計算
try:
L = np.linalg.cholesky(K_y)
except np.linalg.LinAlgError:
return 1e10, np.zeros(3)
# alpha = K_y^{-1} y
alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))
# 対数周辺尤度の3つの項
data_fit = -0.5 * y @ alpha
complexity = -np.sum(np.log(np.diag(L))) # = -0.5 * log|K_y|
constant = -0.5 * N * np.log(2 * np.pi)
log_ml = data_fit + complexity + constant
# --- 勾配の計算 ---
K_y_inv = np.linalg.solve(L.T, np.linalg.solve(L, np.eye(N)))
M = np.outer(alpha, alpha) - K_y_inv # alpha alpha^T - K_y^{-1}
# 距離の二乗行列
dist_sq = np.sum(X**2, axis=1).reshape(-1, 1) + \
np.sum(X**2, axis=1).reshape(1, -1) - 2 * X @ X.T
# 各ハイパーパラメータに関するカーネル行列の偏微分
dK_dl = K * dist_sq / ell**3 # ∂K/∂ℓ
dK_dsf2 = K / sigma_f**2 # ∂K/∂(σ_f²)
dK_dsn2 = np.eye(N) # ∂K_y/∂(σ_n²) = I
# 勾配: 0.5 * tr(M @ ∂K_y/∂θ_j)
grad_l = 0.5 * np.trace(M @ dK_dl)
grad_sf2 = 0.5 * np.trace(M @ dK_dsf2)
grad_sn2 = 0.5 * np.trace(M @ dK_dsn2)
# 対数スケールでの勾配（連鎖律: ∂L/∂(log θ) = θ * ∂L/∂θ）
# log_params[1] = log(σ_f) なので ∂σ_f²/∂(log σ_f) = 2σ_f²
grad_log_ell = grad_l * ell
grad_log_sf = grad_sf2 * 2 * sigma_f**2
grad_log_sn = grad_sn2 * 2 * sigma_n**2
grads = np.array([grad_log_ell, grad_log_sf, grad_log_sn])
return -log_ml, -grads
# --- 最適化の実行 ---
log_params0 = np.array([np.log(1.0), np.log(1.0), np.log(0.1)])
result = minimize(
neg_log_marginal_likelihood,
log_params0,
args=(X_train, y_train),
method='L-BFGS-B',
jac=True # 関数が (値, 勾配) のタプルを返すことを指定
)
# 最適化されたハイパーパラメータ
ell_opt = np.exp(result.x[0])
sf_opt = np.exp(result.x[1])
sn_opt = np.exp(result.x[2])
print("=== スクラッチ実装の結果 ===")
print(f"長さスケール ℓ = {ell_opt:.4f}")
print(f"信号振幅 σ_f = {sf_opt:.4f}")
print(f"ノイズ σ_n = {sn_opt:.4f}")
print(f"対数周辺尤度 = {-result.fun:.4f}")
上記のコードのポイントを整理しておきましょう。
- コレスキー分解: $\bm{K}_y$ の逆行列を直接計算するのではなく、コレスキー分解 $\bm{K}_y = \bm{L}\bm{L}^T$ を用いて $\bm{\alpha} = \bm{K}_y^{-1} \bm{y}$ を安定に計算しています。行列式の対数も $\log |\bm{K}_y| = 2 \sum_i \log L_{ii}$ で計算でき、数値的に安定です。
- 対数スケール: ハイパーパラメータは対数スケールで最適化し、正の値の制約を自然に満たしています。
jac=True:minimizeに渡す関数が目的関数値と勾配を同時に返すことを指定しており、効率的に最適化できます。
scikit-learn との比較
次に、scikit-learn の GaussianProcessRegressor を使って同じデータで最適化を行い、結果を比較します。
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel, WhiteKernel
# カーネルの定義（ConstantKernel * RBF + WhiteKernel）
kernel = ConstantKernel(1.0, constant_value_bounds=(1e-3, 1e3)) \
* RBF(length_scale=1.0, length_scale_bounds=(1e-3, 1e3)) \
+ WhiteKernel(noise_level=0.01, noise_level_bounds=(1e-5, 1e1))
# ガウス過程回帰モデルの構築と学習
gpr = GaussianProcessRegressor(
kernel=kernel,
n_restarts_optimizer=10, # 多点初期化
random_state=42
)
gpr.fit(X_train, y_train)
print("\n=== scikit-learn の結果 ===")
print(f"最適化後のカーネル: {gpr.kernel_}")
print(f"対数周辺尤度 = {gpr.log_marginal_likelihood_value_:.4f}")
予測結果の可視化
最適化したハイパーパラメータを用いて、ガウス過程回帰の予測を可視化します。長さスケールが小さすぎる場合と大きすぎる場合も並べて比較します。
def gp_predict(X_train, y_train, X_test, ell, sigma_f, sigma_n):
"""ガウス過程回帰の予測分布を計算する"""
K = rbf_kernel(X_train, X_train, ell, sigma_f)
K_y = K + sigma_n**2 * np.eye(len(y_train)) + 1e-8 * np.eye(len(y_train))
K_s = rbf_kernel(X_train, X_test, ell, sigma_f)
K_ss = rbf_kernel(X_test, X_test, ell, sigma_f)
L = np.linalg.cholesky(K_y)
alpha = np.linalg.solve(L.T, np.linalg.solve(L, y_train))
# 予測平均
mu = K_s.T @ alpha
# 予測分散
v = np.linalg.solve(L, K_s)
cov = K_ss - v.T @ v
std = np.sqrt(np.diag(cov))
return mu, std
# --- 3つのハイパーパラメータ設定での予測 ---
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# (a) 最適化されたハイパーパラメータ
mu_opt, std_opt = gp_predict(X_train, y_train, X_test, ell_opt, sf_opt, sn_opt)
axes[0].fill_between(X_test.ravel(), mu_opt - 2*std_opt, mu_opt + 2*std_opt,
alpha=0.2, color='blue', label='95% confidence')
axes[0].plot(X_test, mu_opt, 'b-', label='Predictive mean')
axes[0].scatter(X_train, y_train, c='red', s=30, zorder=5, label='Training data')
axes[0].plot(X_test, np.sin(X_test), 'k--', alpha=0.5, label='True function')
axes[0].set_title(f'Optimized: l={ell_opt:.2f}, sf={sf_opt:.2f}, sn={sn_opt:.2f}')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].legend(fontsize=8)
# (b) 長さスケールが小さすぎる場合（過適合傾向）
mu_small, std_small = gp_predict(X_train, y_train, X_test, 0.1, sf_opt, sn_opt)
axes[1].fill_between(X_test.ravel(), mu_small - 2*std_small, mu_small + 2*std_small,
alpha=0.2, color='blue', label='95% confidence')
axes[1].plot(X_test, mu_small, 'b-', label='Predictive mean')
axes[1].scatter(X_train, y_train, c='red', s=30, zorder=5, label='Training data')
axes[1].plot(X_test, np.sin(X_test), 'k--', alpha=0.5, label='True function')
axes[1].set_title('Too small: l=0.1 (overfitting)')
axes[1].set_xlabel('x')
axes[1].legend(fontsize=8)
# (c) 長さスケールが大きすぎる場合（過度な平滑化）
mu_large, std_large = gp_predict(X_train, y_train, X_test, 10.0, sf_opt, sn_opt)
axes[2].fill_between(X_test.ravel(), mu_large - 2*std_large, mu_large + 2*std_large,
alpha=0.2, color='blue', label='95% confidence')
axes[2].plot(X_test, mu_large, 'b-', label='Predictive mean')
axes[2].scatter(X_train, y_train, c='red', s=30, zorder=5, label='Training data')
axes[2].plot(X_test, np.sin(X_test), 'k--', alpha=0.5, label='True function')
axes[2].set_title('Too large: l=10.0 (underfitting)')
axes[2].set_xlabel('x')
axes[2].legend(fontsize=8)
plt.tight_layout()
plt.show()
左のパネルでは最適化されたハイパーパラメータによる予測が真の関数をよく捉えていることがわかります。中央のパネルでは $\ell = 0.1$ と小さすぎるため、データ点の間で予測が激しく振動しています。右のパネルでは $\ell = 10.0$ と大きすぎるため、ほぼ定数に近い滑らかすぎる予測になってしまい、$\sin(x)$ の変動を捉えられていません。
対数周辺尤度のランドスケープ
ハイパーパラメータと対数周辺尤度の関係を可視化してみましょう。
# 長さスケールを変化させたときの対数周辺尤度
length_scales = np.logspace(-1, 1.5, 100)
log_mls = []
for ls in length_scales:
log_params = np.array([np.log(ls), np.log(sf_opt), np.log(sn_opt)])
neg_lml, _ = neg_log_marginal_likelihood(log_params, X_train, y_train)
log_mls.append(-neg_lml)
plt.figure(figsize=(8, 5))
plt.plot(length_scales, log_mls, 'b-', linewidth=2)
plt.axvline(ell_opt, color='r', linestyle='--',
label=f'Optimal l = {ell_opt:.3f}')
plt.xscale('log')
plt.xlabel('Length scale l')
plt.ylabel('Log marginal likelihood')
plt.title('Log Marginal Likelihood vs Length Scale')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このグラフでは、対数周辺尤度が長さスケール $\ell$ に対して明確なピークを持つことが確認できます。最適な $\ell$ の位置（赤い破線）で対数周辺尤度が最大となっており、これより小さくても大きくても対数周辺尤度は減少します。小さい側（過適合側）では複雑さペナルティが支配的となり、大きい側（過度な平滑化側）ではデータ適合項が支配的となって、対数周辺尤度が低下していることが読み取れます。
まとめ
本記事では、ガウス過程回帰のハイパーパラメータ最適化について、周辺尤度最大化の理論からPython実装までを解説しました。
- ガウス過程回帰のハイパーパラメータは、長さスケール $\ell$、信号分散 $\sigma_f^2$、ノイズ分散 $\sigma_n^2$ であり、予測性能を大きく左右する
- 周辺尤度 $p(\bm{y} | X, \bm{\theta})$ は、潜在関数 $\bm{f}$ をガウス積分で消去することで $\mathcal{N}(\bm{0}, \bm{K}_y)$ として解析的に得られる
- 対数周辺尤度は、データ適合項・複雑さペナルティ項・定数項の3つに分解でき、適合度と複雑さのトレードオフを自動的にバランスさせる（オッカムのかみそりの数学的実現）
- 対数周辺尤度の勾配は $\frac{\partial \mathcal{L}}{\partial \theta_j} = \frac{1}{2} \mathrm{tr}\left( (\bm{\alpha}\bm{\alpha}^T – \bm{K}_y^{-1}) \frac{\partial \bm{K}_y}{\partial \theta_j} \right)$ で解析的に得られ、L-BFGS-B法などで効率的に最適化できる
- コレスキー分解と対数スケールでの最適化により、数値的に安定な実装が可能
- スクラッチ実装と scikit-learn の
GaussianProcessRegressorで同等の結果が得られることを確認した
次のステップとして、以下のトピックも参考にしてください。
- ガウス過程分類: ガウス過程を分類問題に拡張する手法。尤度がガウス分布でなくなるため、ラプラス近似や変分推論が必要になります。
- ベイズ最適化: ガウス過程回帰を獲得関数と組み合わせて、ブラックボックス関数の大域的最適化を行う手法。ハイパーパラメータ最適化自体にも応用できます。