ガウス過程回帰（GPR）応用編：カーネル設計とハイパーパラメータ最適化の実務

ガウス過程回帰（GPR）のカーネル設計・ハイパーパラメータ最適化・Python実装を実務目線で解説。sklearn の GaussianProcessRegressor・Matern・ExpSineSquared・WhiteKernel・DotProduct によるカーネル合成（トレンド+周期+ノイズ分解）、log_marginal_likelihood 最大化と n_restarts_optimizer による局所解対策、予測分散と信頼区間・外挿の危険性、O(n³) 対策（Nyström/SoR/inducing points）までカバー。

6 min read

本記事の位置づけ

ガウス過程回帰の基礎では、GP の定義・事後分布の閉形式解・NumPy スクラッチ実装を扱いました。本記事はその実務続編です。基礎編で「RBF カーネル 1 本 + 手動ハイパーパラメータ」だった構成を、実データで通用するレベルへ引き上げます。

  • カーネル設計：RBF / Matérn / 周期 / 線形 / WhiteKernel の使い分けと、和・積による合成
  • ハイパーパラメータ最適化：対数周辺尤度最大化の実際と局所解対策（n_restarts_optimizer）
  • 予測分散の運用：信頼区間の読み方と外挿の危険性
  • スケーラビリティ：\(O(n^3)\) の壁と Nyström / inducing points による近似

コードはすべて scikit-learn 1.9 / NumPy 2.4 で実行検証済みです。

カーネル設計の実務

基礎編で見たとおり、GP の性質はカーネル \(k(\mathbf{x}, \mathbf{x}')\) がすべて決めます。実務でのカーネル設計は「データに含まれる構造（滑らかさ・周期・トレンド・ノイズ）を列挙し、対応する部品を組み合わせる」作業です。

部品となる基本カーネル

sklearn.gaussian_process.kernels の主要カーネルを整理します。

カーネル式（概形）表現する構造
RBF(length_scale)\(\exp(-r^2 / 2\ell^2)\)無限回微分可能な滑らかな関数
Matern(length_scale, nu)式(1)滑らかさを \(\nu\) で制御
ExpSineSquared(ls, p)\(\exp\left(-\frac{2\sin^2(\pi r / p)}{\ell^2}\right)\)厳密な周期 \(p\) の繰り返し
DotProduct(sigma_0)\(\sigma_0^2 + \mathbf{x} \cdot \mathbf{x}'\)線形トレンド
WhiteKernel(noise_level)\(\sigma_n^2 \, \delta(\mathbf{x}, \mathbf{x}')\)観測ノイズ
ConstantKernel(constant)\(\sigma_f^2\)振幅スケール（積で使う）
RationalQuadratic(ls, alpha)\(\left(1 + \frac{r^2}{2\alpha\ell^2}\right)^{-\alpha}\)複数スケールの RBF 混合

ここで \(r = \|\mathbf{x} - \mathbf{x}'\|\) です。

Matérn カーネル：ν の選び方

Matérn カーネルは滑らかさパラメータ \(\nu\) を持ちます：

\[ k_{\nu}(r) = \sigma_f^2 \frac{2^{1-\nu}}{\Gamma(\nu)} \left(\frac{\sqrt{2\nu}\, r}{\ell}\right)^{\nu} K_{\nu}\left(\frac{\sqrt{2\nu}\, r}{\ell}\right) \tag{1} \]

\(\nu\) の実務的な選択指針は明快です。

  • \(\nu = 1/2\) （指数カーネル）：サンプルパスは連続だが微分不可能。ブラウン運動的なギザギザした信号（金融データ、粗い物理計測）向け
  • \(\nu = 3/2\) ：1 回微分可能。多くの実データに合う「デフォルト候補」
  • \(\nu = 5/2\) ：2 回微分可能。ベイズ最適化のサロゲートで最頻出（滑らかすぎず粗すぎない）
  • \(\nu \to \infty\) ：RBF に一致。無限回微分可能で、実データには「滑らかすぎる」ことが多い

\(\nu\) は勾配で最適化できない（閉形式の微分が \(\nu = 1/2, 3/2, 5/2\) 以外で高価）ため、候補を離散的に試して対数周辺尤度（LML）で比較するのが定石です。実行検証した比較コード：

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, WhiteKernel, ConstantKernel as C

rng = np.random.default_rng(42)
X = np.sort(rng.uniform(-4, 4, 40)).reshape(-1, 1)
y = np.sin(2 * X.ravel()) + 0.15 * np.abs(X.ravel()) * rng.standard_normal(40)

for nu in [0.5, 1.5, 2.5, np.inf]:
    kernel = C(1.0) * Matern(length_scale=1.0, nu=nu) + WhiteKernel(0.1)
    gpr = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=5,
                                   random_state=0)
    gpr.fit(X, y)
    print(f"nu={nu}: LML={gpr.log_marginal_likelihood_value_:.3f}")

実行結果：

nu=0.5: LML=-24.069
nu=1.5: LML=-20.297
nu=2.5: LML=-19.103
nu=inf: LML=-17.822

このデータ（滑らかな sin + ヘテロなノイズ）では \(\nu\) が大きいほど LML が改善しました。逆にギザギザした信号では小さい \(\nu\) が勝ちます。LML はカーネル構造そのものの比較にも使える——これが GP のモデル選択の基本作法です。

なお、RBF カーネルは SVM でも中心的な役割を果たしますが、SVM がマージン最大化の決定論的な道具としてカーネルを使うのに対し、GP はカーネルを共分散＝事前分布として使う点が本質的に異なります。

カーネルの和と積：合成の文法

カーネルは以下の閉包性を持ちます。

  • 和 \(k_1 + k_2\) ：独立な 2 つの成分の重ね合わせ（\(f = f_1 + f_2\) , \(f_i \sim \mathcal{GP}(0, k_i)\) ）
  • 積 \(k_1 \times k_2\) ：構造の変調（例：周期カーネル × RBF = 「周期パターンが徐々に変化する」局所周期性）

代表的な設計パターン：

データの構造カーネル合成
トレンド + 季節性 + ノイズDotProduct + ExpSineSquared + WhiteKernel
徐々に形が変わる周期（気温など）ExpSineSquared * RBF + WhiteKernel
長期トレンド + 短期変動RBF(long ls) + RBF(short ls) + WhiteKernel
振幅が入力に依存DotProduct * RBF

事前サンプルによる設計確認（sample_y）

合成したカーネルが意図した構造を表しているかは、フィット前に事前分布から関数をサンプリングして目視確認するのが最も確実です。GaussianProcessRegressor は fit せずに sample_y を呼ぶと事前分布からのサンプルを返します：

X_plot = np.linspace(0, 10, 200).reshape(-1, 1)
kernel = C(0.1) * DotProduct(sigma_0=1.0) + ExpSineSquared(1.0, periodicity=3.0)
gp_prior = GaussianProcessRegressor(kernel=kernel)  # fit しない = 事前分布
samples = gp_prior.sample_y(X_plot, n_samples=3, random_state=0)
# samples.shape == (200, 3)：意図どおり「傾き + 周期」の形か目視確認

「トレンドが乗った周期波形」が出てくれば設計は意図どおり、そうでなければ部品の組み合わせや初期値を見直します。数十行のコードでカーネルという抽象物を可視化できるため、合成カーネルを書くたびに実行する習慣をおすすめします。

多次元入力と ARD（自動関連度決定）

入力が多次元のときは、RBF(length_scale=[l1, l2, ...]) のように次元ごとに独立な length scale を持たせる ARD（Automatic Relevance Determination）が有効です。LML 最適化の結果、無関係な次元の length scale は自動的に大きくなり、その次元がモデルから実質的に除去されます。実行検証：

Xa = rng.uniform(-2, 2, (60, 2))
ya = np.sin(2 * Xa[:, 0]) + 0.05 * rng.standard_normal(60)  # x2 は無関係

kernel = C(1.0) * RBF(length_scale=[1.0, 1.0]) + WhiteKernel(0.1)
gpr = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10,
                               random_state=0)
gpr.fit(Xa, ya)
print(gpr.kernel_)
# 0.996**2 * RBF(length_scale=[0.925, 1e+05]) + WhiteKernel(noise_level=0.00233)

関係のある次元 \(x_1\) は \(\ell_1 = 0.925\) 、無関係な次元 \(x_2\) は \(\ell_2 = 10^5\) （上限に張り付き）——特徴量選択が周辺尤度最大化の副産物として自動で得られるわけです。学習後の length scale の逆数は、各特徴量の重要度指標として解釈できます。

実践：トレンド + 周期 + ノイズの分解

「線形トレンド + 周期 2.5 の正弦波 + ガウスノイズ」という合成データを、対応する 3 部品の和で分解してみます。

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import (
    ExpSineSquared, DotProduct, WhiteKernel, ConstantKernel as C,
)

rng = np.random.default_rng(42)
X = np.linspace(0, 10, 80).reshape(-1, 1)
y = (0.5 * X.ravel()                                  # 線形トレンド
     + 1.5 * np.sin(2 * np.pi * X.ravel() / 2.5)      # 周期 2.5
     + 0.3 * rng.standard_normal(80))                 # 観測ノイズ

kernel = (
    C(1.0) * DotProduct(sigma_0=1.0)                              # トレンド成分
    + C(1.0) * ExpSineSquared(length_scale=1.0, periodicity=3.0)  # 周期成分
    + WhiteKernel(noise_level=0.1)                                # ノイズ成分
)
gpr = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10,
                               normalize_y=True, random_state=0)
gpr.fit(X, y)
print(gpr.kernel_)   # 学習後のカーネル（ハイパーパラメータ最適化済み）
print(gpr.log_marginal_likelihood_value_)

実行結果：

0.299**2 * DotProduct(sigma_0=1.04)
  + 4.33**2 * ExpSineSquared(length_scale=9.08, periodicity=2.49)
  + WhiteKernel(noise_level=0.0201)
LML: 26.260

注目すべきは periodicity=2.49——初期値 3.0 から出発して真の周期 2.5 をデータだけからほぼ正確に復元しています。WhiteKernel の学習値はノイズ分散の推定値になっており、カーネル合成が「信号の加法分解 + ノイズ推定」を 1 回のフィットで同時に実現していることがわかります。

この「トレンド + 周期 + ノイズへの加法分解」という発想は、カルマンフィルタの状態空間モデル（レベル + 季節成分の状態分解）と双対的な関係にあります。実際、Matérn カーネルの GP は等価な状態空間表現を持ち、RTS スムーザで \(O(n)\) 推論できることが知られています。

ハイパーパラメータ最適化

対数周辺尤度：数式と直感

カーネルのハイパーパラメータ \(\boldsymbol{\theta}\) （length scale、振幅、周期、ノイズ分散など）は、対数周辺尤度（log marginal likelihood, LML）の最大化で決めます：

$$ \log p(\mathbf{y} \mid X, \boldsymbol{\theta}) = \underbrace{-\frac{1}{2}\mathbf{y}^T K_{\boldsymbol{\theta}}^{-1} \mathbf{y}}{\text{データ適合}} \underbrace{- \frac{1}{2} \log |K{\boldsymbol{\theta}}|}_{\text{複雑度ペナルティ}}

  • \frac{n}{2} \log 2\pi \tag{2} $$

ここで \(K_{\boldsymbol{\theta}} = K + \sigma_n^2 I\) です。第 1 項はデータへの当てはまり、第 2 項は「その \(\boldsymbol{\theta}\) が表現できる関数集合の大きさ」への罰則で、交差検証なしにオッカムの剃刀が働くのが GP の強みです。length scale を小さくすれば第 1 項は改善しますが \(|K_{\boldsymbol{\theta}}|\) が大きくなり第 2 項が悪化する——このトレードオフの釣り合い点が選ばれます。

勾配は閉形式で書けます：

\[ \frac{\partial}{\partial \theta_j} \log p(\mathbf{y} \mid X, \boldsymbol{\theta}) = \frac{1}{2} \mathrm{tr}\!\left[ \left( \boldsymbol{\alpha}\boldsymbol{\alpha}^T - K_{\boldsymbol{\theta}}^{-1} \right) \frac{\partial K_{\boldsymbol{\theta}}}{\partial \theta_j} \right], \quad \boldsymbol{\alpha} = K_{\boldsymbol{\theta}}^{-1} \mathbf{y} \tag{3} \]

scikit-learn はこの解析勾配を使い、\(\log \boldsymbol{\theta}\) 空間で L-BFGS-B を走らせます（スケールパラメータは対数空間の方が最適化が安定するため）。

局所解対策：n_restarts_optimizer

LML は \(\boldsymbol{\theta}\) について多峰です。典型的なのは「短い length scale + 小ノイズ（データを暗記）」と「長い length scale + 大ノイズ（信号をノイズ扱い）」という 2 つの解釈が別々の局所解になるケースで、周期カーネルでは周期の整数倍・整数分の 1 にも局所解が並びます。

対策は n_restarts_optimizer によるマルチスタートです。先ほどの合成カーネルで比較しました：

for nr in [0, 10]:
    gpr = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=nr,
                                   normalize_y=True, random_state=0)
    gpr.fit(X, y)
    print(f"n_restarts={nr}: LML={gpr.log_marginal_likelihood_value_:.3f}")

実行結果：

n_restarts= 0: LML=-85.099   # 初期値からの 1 回きり → 局所解に捕まる
n_restarts=10: LML= 26.260   # 10 回リスタート → 大域解（周期 2.49 を発見）

LML の差は 111 nat。リスタートなしでは周期構造を見つけられず、まったく別のモデルに落ちています。実務の指針：

  • 合成カーネル（特に周期を含む）では n_restarts_optimizer=10 以上を推奨
  • リスタートの初期値は各ハイパーパラメータの bounds から対数一様サンプリングされるため、bounds を物理的に妥当な範囲に絞ることが探索効率に直結する
  • gpr.log_marginal_likelihood(theta) で任意の \(\boldsymbol{\theta}\) の LML を手動評価できる。収束診断やプロファイル尤度の可視化に便利
theta_opt = gpr.kernel_.theta.copy()          # 最適化後（対数空間）
print(gpr.log_marginal_likelihood(theta_opt))  # 26.260
theta_bad = theta_opt.copy()
theta_bad[-1] = np.log(1.0)                    # ノイズ分散だけ 1.0 に固定
print(gpr.log_marginal_likelihood(theta_bad))  # -84.535 → ノイズ推定の重要性

その他の実務パラメータ

  • normalize_y=True：目的変数を平均 0・分散 1 に正規化。GP の事前平均は 0 なので、平均が大きくずれたデータでは必須級
  • alpha：対角に足す jitter（既定 \(10^{-10}\) ）。WhiteKernel でノイズを「学習」するか、alpha で「固定」するかは排他的に使い分ける（両方入れると同一の役割が重複する）
  • hyperparameter の固定：length_scale_bounds="fixed" で特定パラメータを最適化から除外できる

予測分散の使い方

信頼区間

predict(X, return_std=True) は予測平均 \(\mu_*\) と標準偏差 \(\sigma_*\) を返します。95% 信頼区間は \(\mu_* \pm 1.96\sigma_*\) です：

X_test = np.array([[5.0], [12.0], [20.0]])
mu, std = gpr.predict(X_test, return_std=True)
# x= 5.0 (内挿):  mean= 2.503  std=0.240  95% CI=[ 2.033,  2.974]
# x=12.0 (外挿):  mean= 4.582  std=0.247  95% CI=[ 4.098,  5.065]
# x=20.0 (外挿):  mean=10.199  std=0.277  95% CI=[ 9.656, 10.742]

return_cov=True を使えばテスト点間の完全な共分散行列が得られ、事後分布からの関数サンプリング（sample_y）や、複数点の同時信頼領域の構成に使えます。この予測分散は時系列異常検知にも直結します——観測値が予測区間 \(\mu_* \pm 3\sigma_*\) を外れたら異常、という動的しきい値がそのまま作れます。

外挿の危険性

上の結果で \(x = 20\) （学習域 \([0, 10]\) の遥か外）でも真値 10.0 に近い 10.199 を返しているのは、線形 + 周期という大域構造を持つカーネルだからです。同じ問題を RBF 単体で解くと様相が一変します：

kernel_rbf = C(1.0) * RBF(length_scale=1.0) + WhiteKernel(noise_level=0.1)
# 学習域 [0, 5] でフィット後：
# x= 2.5 (内挿): true= 2.188  mean= 2.102  std=0.083
# x= 6.0 (外挿): true= 2.249  mean= 1.096  std=1.285
# x=10.0 (外挿): true= 4.012  mean= 0.000  std=1.973

RBF は距離 \(\ell\) の数倍で相関がゼロになる局所カーネルなので、外挿では予測平均が事前平均 0 へ回帰し、標準偏差は事前の振幅 \(\sigma_f \approx 1.97\) まで戻ります。ここから 2 つの教訓が得られます。

  1. RBF/Matérn 系の外挿は「わかりません」という正直な答えを返す。std の膨張は故障ではなく仕様であり、この不確実性の申告こそ点推定モデルにない価値
  2. 外挿で構造を延長したいなら、その構造をカーネルに明示的に入れる（トレンドなら DotProduct、周期なら ExpSineSquared）。ただしそれは「構造が域外でも続く」という強い仮定を置くことと等価

なお LSTM や Transformer などのニューラル時系列モデルは外挿でも「自信ありげな」点予測を出し続けるため、予測区間が本質的に必要な用途（安全制約・在庫計画など）では GPR が今も第一候補です。モデル選択の全体像は時系列 ML ハブを参照してください。

ベイズ最適化との接続

GPR の応用先として最重要なのがベイズ最適化のサロゲートモデルです。評価コストの高いブラックボックス関数 \(f\) を GPR で近似し、予測平均 \(\mu_*\) （活用）と予測分散 \(\sigma_*\) （探索）を獲得関数で束ねて次の評価点を選びます：

\[ \mathbf{x}_{\text{next}} = \arg\max_{\mathbf{x}} \; \alpha_{\text{EI}}(\mathbf{x}), \quad \alpha_{\text{EI}}(\mathbf{x}) = \mathbb{E}\left[\max(f_{\text{best}} - f(\mathbf{x}), 0)\right] \tag{4} \]

本記事で扱った実務論点はすべてベイズ最適化の品質に直結します——カーネルは Matérn \(\nu = 5/2\) が事実上の標準、LML の局所解はサロゲートの劣化（＝無駄な評価）を招くためリスタート必須、そして獲得関数の「探索」項は予測分散そのものです。獲得関数（EI / UCB / PI）の導出と実装はベイズ最適化の基礎で詳説しています。

O(n³) の壁と近似手法

厳密な GPR は \(K_{\boldsymbol{\theta}} \in \mathbb{R}^{n \times n}\) の Cholesky 分解に \(O(n^3)\) 、保持に \(O(n^2)\) を要し、\(n \gtrsim 10^4\) で限界を迎えます。LML 最適化は分解を何十回も繰り返すため、体感的な上限はさらに低くなります。主要な近似ファミリーを整理します。

Nyström 近似 / SoR

\(m \ll n\) 個の誘導点（inducing points） \(Z = \{\mathbf{z}_1, \ldots, \mathbf{z}_m\}\) を選び、カーネル行列を低ランク近似します：

\[ K \approx \tilde{K} = K_{nm} K_{mm}^{-1} K_{mn} \tag{5} \]

SoR（Subset of Regressors）はこの \(\tilde{K}\) で GP を置き換え、Woodbury の恒等式により計算量を \(O(nm^2)\) に削減します。\(n = 2000, m = 100\) で実行検証したところ、予測平均の厳密 GP との差は最大 \(1.2 \times 10^{-7}\) ——滑らかな信号なら誘導点 5% で厳密解を実質再現できます：

from scipy.linalg import cho_factor, cho_solve

idx = rng.choice(n, m, replace=False)
Z = X[idx]                              # 誘導点（ランダム選択）
K_mm = kernel(Z, Z) + 1e-8 * np.eye(m)
K_nm = kernel(X, Z)
A = K_mm * sn2 + K_nm.T @ K_nm          # m×m に縮約
c, low = cho_factor(A)
alpha = cho_solve((c, low), K_nm.T @ y)
mu = kernel(X_test, Z) @ alpha          # O(nm^2) で予測平均

ただし SoR は予測分散を過小評価する（誘導点の外で不確実性がゼロに潰れる）欠点があります。

FITC / VFE / SVGP

  • FITC：SoR に対角補正を加え分散の潰れを緩和
  • VFE（変分自由エネルギー）：誘導点自体を変分下限の最大化で最適化。近似の「質」が LML 下限として定量化される
  • SVGP：VFE をミニバッチ化し \(n \sim 10^6\) 超へ。GPyTorch / GPflow で利用可能

実務の選択指針：\(n < 5{,}000\) なら厳密 GPR、\(n < 10^5\) なら VFE 系（誘導点 \(m = 100 \sim 1000\) ）、それ以上は SVGP かニューラルモデルへの乗り換えを検討、が目安です。

まとめ

論点実務の要点
カーネル選択Matérn \(\nu \in \{3/2, 5/2\}\) を既定に、LML で構造比較
カーネル合成和 = 成分の重ね合わせ、積 = 変調。トレンド+周期+WhiteKernel が定番
ハイパーパラメータLML 最大化 + n_restarts_optimizer>=10、bounds を物理的に絞る
予測分散95% CI は \(\mu \pm 1.96\sigma\) 。外挿での分散膨張は仕様（信頼せよ）
スケール\(n \gtrsim 10^4\) で誘導点近似（SoR → VFE → SVGP）

カーネルという言語でデータの構造を記述し、周辺尤度という基準でその記述を検証する——GPR の実務はこの往復運動です。

おすすめ書籍

ガウス過程と機械学習（持橋大地・大羽成征、講談社） Amazonで見る

ガウス過程の理論からベイズ最適化への応用まで、日本語で体系的に学べる決定版です。

ベイズ推論による機械学習入門（須山敦志、講談社） Amazonで見る

ベイズ推論の枠組みをゼロから積み上げる入門書で、本記事の背後にある確率モデリングの考え方を補強できます。

※ 上記は Amazon アソシエイトのリンクです。

関連記事

参考文献

  • Rasmussen, C. E., & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press. （特に Ch. 4 カーネル、Ch. 5 モデル選択、Ch. 8 近似手法）
  • Duvenaud, D. (2014). Automatic Model Construction with Gaussian Processes. PhD thesis, University of Cambridge. （カーネル合成の文法）
  • Titsias, M. (2009). “Variational Learning of Inducing Variables in Sparse Gaussian Processes.” AISTATS 2009. （VFE）
  • Hensman, J., Fusi, N., & Lawrence, N. D. (2013). “Gaussian Processes for Big Data.” UAI 2013. （SVGP）
  • scikit-learn documentation: Gaussian Processes

この記事に関連する無料ツール

DevToolBox Moving Average Calculator SMA/EMA/WMA をブラウザで即座に計算 DevToolBox JSON Formatter JSONを整形・圧縮・バリデーション DevToolBox Regex Tester 正規表現をリアルタイムテスト

関連記事