ガウス過程回帰の基礎とPython実装

ガウス過程とは

ガウス過程（Gaussian Process, GP） は、関数そのものを確率変数として扱う確率モデルです。直感的には、有限次元の多変量正規分布を「無限次元」へ拡張したものと見なせます。

通常の正規分布が「ベクトル \(\mathbf{y} \in \mathbb{R}^n\) の分布」を定めるのに対し、GP は「関数 \(f: \mathcal{X} \rightarrow \mathbb{R}\) の分布」を定めます。GP の定義は次の通りです：任意の有限個の入力点 \(\{\mathbf{x}_1, \ldots, \mathbf{x}_n\}\) に対して、対応する関数値ベクトル \([f(\mathbf{x}_1), \ldots, f(\mathbf{x}_n)]^T\) が多変量正規分布に従うとき、\(f\) はガウス過程に従うといいます。

GP は平均関数 \(m(\mathbf{x})\) とカーネル関数（共分散関数）\(k(\mathbf{x}, \mathbf{x}')\) の2つで完全に特徴付けられます：

実用上は \(m(\mathbf{x}) = 0\) とすることが多く、関数の形状はカーネルによって決まります。

カーネル関数（事前分布の設計）

カーネル関数 \(k(\mathbf{x}, \mathbf{x}')\) は2点間の関数値の共分散を与え、関数の滑らかさ・周期性などの事前知識を埋め込みます。

RBF（二乗指数）カーネル

最も広く使われるのが RBF カーネル（または二乗指数カーネル）です：

• \(\sigma_f^2\) : 出力の分散（関数の振幅）
• \(\ell\) : length scale。距離 \(\ell\) だけ離れた点間の相関がどの程度残るかを決める

\(\ell\) が大きいほど関数は滑らかになり、小さいほど局所的に激しく変動します。RBF カーネルから生成される関数は無限回微分可能です。

Matern カーネル

実問題ではしばしば「滑らかすぎ」が問題になります。Matern カーネルは滑らかさパラメータ \(\nu\) で微分可能性を制御できます：

ここで \(K_{\nu}\) は第二種変形ベッセル関数です。\(\nu = 1/2\) で指数カーネル、\(\nu \to \infty\) で RBF に収束します。実用では \(\nu = 3/2, 5/2\) がよく用いられます。

観測データへの条件付け（事後分布）

GP の真価は、観測データを与えると 解析的に閉形式で事後分布が得られる 点にあります。

\(n\) 個の観測 \(\mathcal{D} = \{(\mathbf{x}_i, y_i)\}_{i=1}^n\) にノイズを仮定します：

新しい入力 \(\mathbf{x}_*\) における関数値 \(f_*\) と観測ベクトル \(\mathbf{y}\) の同時分布は次式で書けます：

ここで \(K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)\) 、\(\mathbf{k}_* = [k(\mathbf{x}_1, \mathbf{x}_*), \ldots, k(\mathbf{x}_n, \mathbf{x}_*)]^T\) です。多変量正規分布の条件付けの公式から、事後予測分布が閉形式で求まります：

式(6)が予測平均、式(7)が予測分散です。データが密な領域では \(\sigma_*^2\) が小さくなり、データから離れた領域では \(\sigma_f^2\) に近づく 「データに引き寄せられる挙動」 が自動的に得られます。これがニューラルネットなどの点推定モデルにはない、GP の本質的な強みです。

周辺尤度とハイパーパラメータ

カーネルのハイパーパラメータ \(\boldsymbol{\theta} = (\sigma_f^2, \ell, \sigma_n^2)\) は、対数周辺尤度（log marginal likelihood） の最大化により学習できます：

第1項がデータ適合度、第2項が モデル複雑度のペナルティ（オッカムの剃刀）として働きます。GP がベイズ的に過適合を抑制するのはこの項のおかげです。

Python実装：1次元 GP 回帰

NumPy と SciPy だけで GP 回帰をスクラッチ実装します。Cholesky 分解で逆行列計算を安定化します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import cho_factor, cho_solve

class GaussianProcessRegressor:
    def __init__(self, length_scale=1.0, signal_var=1.0, noise_var=1e-4):
        self.l = length_scale
        self.sf2 = signal_var
        self.sn2 = noise_var

    def rbf(self, X1, X2):
        """式(2): RBFカーネル"""
        d2 = (np.sum(X1**2, axis=1, keepdims=True)
              - 2.0 * X1 @ X2.T
              + np.sum(X2**2, axis=1))
        return self.sf2 * np.exp(-0.5 * d2 / self.l**2)

    def fit(self, X, y):
        self.X_train = X.copy()
        self.y_train = y.copy()
        K = self.rbf(X, X) + self.sn2 * np.eye(len(X))
        # Cholesky で安定化
        self.L_, self.lower_ = cho_factor(K, lower=True)
        self.alpha_ = cho_solve((self.L_, self.lower_), y)

    def predict(self, X_new):
        """式(6), (7): 事後予測平均と分散"""
        k_star = self.rbf(self.X_train, X_new)
        mu = k_star.T @ self.alpha_
        v = cho_solve((self.L_, self.lower_), k_star)
        var = self.sf2 - np.sum(k_star * v, axis=0)
        var = np.maximum(var, 1e-10)
        return mu, var

    def log_marginal_likelihood(self):
        """式(8): 対数周辺尤度"""
        n = len(self.y_train)
        log_det = 2.0 * np.sum(np.log(np.diag(self.L_)))
        return (-0.5 * self.y_train @ self.alpha_
                - 0.5 * log_det
                - 0.5 * n * np.log(2 * np.pi))

サンプリングと予測

ノイズあり観測から事後分布を構築し、95%信頼区間を可視化します。

# 真の関数（ノイズなし）
def true_f(x):
    return np.sin(2.0 * x) + 0.3 * x

# 観測データ
np.random.seed(0)
X_train = np.array([[-3.0], [-2.0], [-0.5], [1.0], [2.5], [3.5]])
y_train = true_f(X_train.ravel()) + 0.1 * np.random.randn(len(X_train))

gp = GaussianProcessRegressor(length_scale=1.0, signal_var=1.0, noise_var=0.01)
gp.fit(X_train, y_train)

X_test = np.linspace(-5.0, 5.0, 300).reshape(-1, 1)
mu, var = gp.predict(X_test)
sigma = np.sqrt(var)

plt.figure(figsize=(10, 5))
plt.plot(X_test, true_f(X_test.ravel()), "k--", label="True function")
plt.plot(X_test, mu, "b-", label="GP mean")
plt.fill_between(X_test.ravel(), mu - 1.96 * sigma, mu + 1.96 * sigma,
                 alpha=0.2, color="blue", label="95% CI")
plt.scatter(X_train, y_train, c="red", s=60, zorder=5, label="Observations")
plt.xlabel("x"); plt.ylabel("f(x)")
plt.title("Gaussian Process Regression (RBF kernel)")
plt.legend(); plt.grid(True); plt.tight_layout()
plt.savefig("gp_regression.png", dpi=150)

print(f"log marginal likelihood: {gp.log_marginal_likelihood():.4f}")

観測点付近では信頼区間が狭まり、データから離れるほど不確実性（青い帯）が広がる挙動が確認できます。これは式(7)の第1項 \(k(\mathbf{x}_*, \mathbf{x}_*) = \sigma_f^2\) が支配的になるためです。

ハイパーパラメータの影響

length scale \(\ell\) とノイズ分散 \(\sigma_n^2\) は GP の挙動を大きく変えます。

length scale \(\ell\)

• \(\ell\) が小さい：局所的に激しく変動、データ点間の補間も急峻。過適合のリスク。
• \(\ell\) が大きい：滑らかな近似、データ点を貫通せず平均化される。過小適合のリスク。

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for ax, l in zip(axes, [0.3, 1.0, 3.0]):
    gp = GaussianProcessRegressor(length_scale=l, signal_var=1.0, noise_var=0.01)
    gp.fit(X_train, y_train)
    mu, var = gp.predict(X_test)
    sigma = np.sqrt(var)
    ax.plot(X_test, true_f(X_test.ravel()), "k--")
    ax.plot(X_test, mu, "b-")
    ax.fill_between(X_test.ravel(), mu - 1.96*sigma, mu + 1.96*sigma,
                    alpha=0.2, color="blue")
    ax.scatter(X_train, y_train, c="red", s=40, zorder=5)
    ax.set_title(f"length_scale = {l}")
    ax.grid(True)
plt.tight_layout()
plt.savefig("gp_length_scale.png", dpi=150)

ノイズ分散 \(\sigma_n^2\)

\(\sigma_n^2\) が小さいと GP は観測点を厳密に通り、大きいと観測値からの乖離を許容して滑らかな曲線になります。式(6)の \((\mathbf{K} + \sigma_n^2 \mathbf{I})^{-1}\) の数値安定化 jitter としても重要です（典型値 \(10^{-6} \sim 10^{-4}\) ）。

実用上は、式(8)の対数周辺尤度を scipy.optimize.minimize で最大化してハイパーパラメータを自動調整します。

計算コストと近似手法

GP の最大の弱点は \(O(n^3)\) の計算量（カーネル行列の Cholesky 分解）と \(O(n^2)\) のメモリ です。\(n \approx 10^4\) あたりが純粋な GP の実用上限とされます。

大規模データには以下の近似が用いられます：

• Inducing point methods（FITC, VFE, SVGP）：\(m \ll n\) 個の代表点で近似し \(O(nm^2)\) に削減
• Random Fourier Features：カーネルを有限次元のランダム特徴で近似
• Local approximation：近傍データのみで予測

ベイズ最適化への接続

GP の予測平均（式6）と予測分散（式7）は、ベイズ最適化のサロゲートモデルとしてそのまま再利用されます。獲得関数（Expected Improvement, UCB など）は \(\mu_*\) と \(\sigma_*\) の関数として定義され、不確実性が高くかつ改善が期待される点を次の評価点として選びます。

たとえば UCB 獲得関数は次のようにシンプルに書けます：

GP がベイズ最適化の根幹を支えていることがわかります。詳細は次の記事を参照してください。

関連記事

• ベイズ最適化の基礎とPython実装 - 本記事のGPをサロゲートモデルとして用い、獲得関数（EI/UCB/PI）で効率的に大域最適解を探索する手法を解説しています。
• MCMC入門：メトロポリス・ヘイスティングス法とギブスサンプリング - GPの背景にあるベイズ推論の計算手法。
• 確率的勾配降下法からAdamまで - GPで学習されたハイパーパラメータをチューニング対象とする勾配ベース最適化。
• サポートベクターマシン（SVM）の理論とPython実装 - 同じくRBFカーネルを用いる代表的な手法。GPの確率的予測とSVMの決定論的マージン最大化を対比できます。
• 粒子群最適化（PSO）の理論とPython実装 - 大量評価で探索する集団最適化。GPサロゲートを用いるベイズ最適化との評価コスト観点での使い分けが理解できます。

参考文献

• Rasmussen, C. E., & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press.
• Murphy, K. P. (2023). Probabilistic Machine Learning: Advanced Topics. MIT Press, Ch. 18.
• Snelson, E., & Ghahramani, Z. (2006). “Sparse Gaussian Processes using Pseudo-inputs.” NeurIPS 2006.
• Hensman, J., Fusi, N., & Lawrence, N. D. (2013). “Gaussian Processes for Big Data.” UAI 2013.