ガウス過程（Gaussian Process）は柔軟なノンパラメトリックモデルですが、データ数 $N$ に対して予測分布の計算量が $O(N^3)$ になるという欠点があります。$N$ が数千を超えると、実用的な計算時間で推論を行うことが困難になります。
疎ガウス過程（Sparse Gaussian Processes）は、この計算量の問題を解決するための手法です。少数の「誘導点（inducing points）」を導入し、計算量を $O(NM^2)$（$M \ll N$）に削減します。
本記事では、疎ガウス過程の基本的なアイデアから数学的な導出、Pythonでの実装までを解説します。
本記事の内容
- 標準ガウス過程の計算量の問題
- 誘導点による近似の考え方
- 変分疎ガウス過程（SVGP）の導出
- Pythonでの実装と比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
標準ガウス過程の計算量
まず、標準的なガウス過程回帰の計算量を確認しましょう。訓練データ $\bm{X} = \{x_1, \dots, x_N\}$ と観測値 $\bm{y} = \{y_1, \dots, y_N\}$ が与えられたとき、新しい入力 $x_*$ に対する予測分布は次のようになります。
$$ \begin{align} \mu_* &= \bm{k}_*^T (\bm{K} + \sigma_n^2 \bm{I})^{-1} \bm{y} \\ \sigma_*^2 &= k_{**} – \bm{k}_*^T (\bm{K} + \sigma_n^2 \bm{I})^{-1} \bm{k}_* \end{align} $$
ここで $\bm{K}$ は $N \times N$ のグラム行列です。この逆行列の計算に $O(N^3)$ の計算量がかかります。
| データ数 $N$ | 計算量のオーダー | 実用性 |
|---|---|---|
| 100 | $10^6$ | 問題なし |
| 1,000 | $10^9$ | やや重い |
| 10,000 | $10^{12}$ | 非常に重い |
| 100,000 | $10^{15}$ | 事実上不可能 |
誘導点による近似
疎ガウス過程の基本的なアイデアは、$N$ 個の訓練データ全てを使う代わりに、$M$ 個（$M \ll N$）の誘導点（inducing points） $\bm{Z} = \{z_1, \dots, z_M\}$ を導入することです。
イメージとしては、大量のデータの特徴を少数の代表的な点で要約するということです。誘導点での関数値 $\bm{u} = f(\bm{Z})$ を介して予測を行うことで、計算量を $O(NM^2)$ に削減します。
誘導点の位置 $\bm{Z}$ はハイパーパラメータとして最適化の対象になります。
変分疎ガウス過程（SVGP）の導出
変分疎ガウス過程（Stochastic Variational Gaussian Process, SVGP）は、変分推論の枠組みで疎ガウス過程を定式化する手法です。
同時分布
ガウス過程の事前分布のもとで、訓練データでの関数値 $\bm{f}$ と誘導点での関数値 $\bm{u}$ の同時分布は次のようになります。
$$ \begin{equation} p(\bm{f}, \bm{u}) = p(\bm{f}|\bm{u}) \, p(\bm{u}) \end{equation} $$
ここで、
$$ \begin{align} p(\bm{u}) &= \mathcal{N}(\bm{u} | \bm{0}, \bm{K}_{MM}) \\ p(\bm{f}|\bm{u}) &= \mathcal{N}(\bm{f} | \bm{K}_{NM}\bm{K}_{MM}^{-1}\bm{u}, \, \tilde{\bm{K}}) \end{align} $$
ただし $\bm{K}_{MM}$ は誘導点間のカーネル行列（$M \times M$）、$\bm{K}_{NM}$ は訓練データと誘導点間のカーネル行列（$N \times M$）、$\tilde{\bm{K}} = \bm{K}_{NN} – \bm{K}_{NM}\bm{K}_{MM}^{-1}\bm{K}_{MN}$ です。
変分下界（ELBO）
誘導点の変分分布 $q(\bm{u}) = \mathcal{N}(\bm{u}|\bm{m}, \bm{S})$ を導入します。ここで $\bm{m}$ は $M$ 次元の平均ベクトル、$\bm{S}$ は $M \times M$ の共分散行列です。
ELBOは次のように導出されます。
$$ \begin{equation} \ln p(\bm{y}) \geq \text{ELBO} = \sum_{i=1}^{N} E_{q(f_i)}\left[\ln p(y_i|f_i)\right] – D_{\text{KL}}(q(\bm{u}) \| p(\bm{u})) \end{equation} $$
ここで $q(f_i)$ は次のように計算されます。
$$ \begin{equation} q(f_i) = \int p(f_i|\bm{u}) \, q(\bm{u}) \, d\bm{u} = \mathcal{N}(f_i | \mu_i, \sigma_i^2) \end{equation} $$
$$ \begin{align} \mu_i &= \bm{k}_i^T \bm{K}_{MM}^{-1} \bm{m} \\ \sigma_i^2 &= k_{ii} – \bm{k}_i^T \bm{K}_{MM}^{-1}(\bm{K}_{MM} – \bm{S})\bm{K}_{MM}^{-1}\bm{k}_i \end{align} $$
ただし $\bm{k}_i = \bm{K}_{Mi}$ は誘導点と $x_i$ 間のカーネルベクトルです。
ガウス尤度の場合
尤度が $p(y_i|f_i) = \mathcal{N}(y_i|f_i, \sigma_n^2)$ の場合、期待値が解析的に計算できます。
$$ \begin{equation} E_{q(f_i)}[\ln p(y_i|f_i)] = -\frac{1}{2}\ln(2\pi\sigma_n^2) – \frac{1}{2\sigma_n^2}\left((y_i – \mu_i)^2 + \sigma_i^2\right) \end{equation} $$
予測分布
学習後の予測分布は、新しい入力 $x_*$ に対して次のようになります。
$$ \begin{align} \mu_* &= \bm{k}_*^T \bm{K}_{MM}^{-1} \bm{m} \\ \sigma_*^2 &= k_{**} – \bm{k}_*^T \bm{K}_{MM}^{-1}(\bm{K}_{MM} – \bm{S})\bm{K}_{MM}^{-1}\bm{k}_* \end{align} $$
計算量は $O(M^2)$ であり、$M$ が小さければ高速に予測できます。
Pythonでの実装
疎ガウス過程をスクラッチで実装し、標準ガウス過程と比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
np.random.seed(42)
# RBFカーネル
def rbf_kernel(X1, X2, length_scale=1.0, variance=1.0):
"""RBFカーネル行列を計算"""
sq_dist = np.sum(X1**2, axis=1, keepdims=True) - 2 * X1 @ X2.T + np.sum(X2**2, axis=1)
return variance * np.exp(-0.5 * sq_dist / length_scale**2)
# 標準ガウス過程回帰
def gp_predict(X_train, y_train, X_test, noise=0.1, length_scale=1.0):
"""標準GPの予測"""
N = len(X_train)
K = rbf_kernel(X_train, X_train, length_scale) + noise**2 * np.eye(N)
K_s = rbf_kernel(X_train, X_test, length_scale)
K_ss = rbf_kernel(X_test, X_test, length_scale)
L = np.linalg.cholesky(K)
alpha = np.linalg.solve(L.T, np.linalg.solve(L, y_train))
mu = K_s.T @ alpha
v = np.linalg.solve(L, K_s)
var = np.diag(K_ss) - np.sum(v**2, axis=0)
return mu, var
# 疎ガウス過程（SVGP）
class SparseGP:
def __init__(self, X_inducing, noise=0.1, length_scale=1.0):
self.Z = X_inducing # 誘導点
self.M = len(X_inducing)
self.noise = noise
self.length_scale = length_scale
# 変分パラメータ
self.m = np.zeros(self.M)
self.S = np.eye(self.M)
def fit(self, X_train, y_train, n_iter=100, lr=0.01):
"""ELBO最大化による学習"""
N = len(X_train)
M = self.M
elbo_history = []
for it in range(n_iter):
# カーネル行列の計算
Kmm = rbf_kernel(self.Z, self.Z, self.length_scale) + 1e-6 * np.eye(M)
Knm = rbf_kernel(X_train, self.Z, self.length_scale)
L_mm = np.linalg.cholesky(Kmm)
Kmm_inv = np.linalg.solve(L_mm.T, np.linalg.solve(L_mm, np.eye(M)))
# q(f_i)の平均と分散
A = Knm @ Kmm_inv # N x M
mu_f = A @ self.m # N
# 各データ点の分散
Knn_diag = np.ones(N) # variance=1のRBFカーネルの対角要素
AS = A @ self.S
var_f = Knn_diag - np.sum(A * (Knm @ Kmm_inv), axis=1) + np.sum(AS * A, axis=1)
var_f = np.maximum(var_f, 1e-6)
# ELBO計算
# 尤度項
like_term = -0.5 * N * np.log(2 * np.pi * self.noise**2) \
- 0.5 / self.noise**2 * (np.sum((y_train - mu_f)**2) + np.sum(var_f))
# KLダイバージェンス
L_S = np.linalg.cholesky(self.S + 1e-6 * np.eye(M))
kl = 0.5 * (np.trace(Kmm_inv @ self.S)
+ self.m @ Kmm_inv @ self.m
- M
+ np.sum(np.log(np.diag(L_mm))) * 2
- np.sum(np.log(np.diag(L_S))) * 2)
elbo = like_term - kl
elbo_history.append(elbo)
# 自然勾配に基づく更新
Lambda = A.T @ A / self.noise**2 + Kmm_inv
self.S = np.linalg.inv(Lambda + 1e-6 * np.eye(M))
self.m = self.S @ (A.T @ y_train / self.noise**2)
return elbo_history
def predict(self, X_test):
"""予測分布の計算"""
M = self.M
Kmm = rbf_kernel(self.Z, self.Z, self.length_scale) + 1e-6 * np.eye(M)
Ksm = rbf_kernel(X_test, self.Z, self.length_scale)
Kmm_inv = np.linalg.inv(Kmm)
A = Ksm @ Kmm_inv
mu = A @ self.m
Kss_diag = np.ones(len(X_test))
var = Kss_diag - np.sum(A * (Ksm @ Kmm_inv), axis=1) + np.sum((A @ self.S) * A, axis=1)
var = np.maximum(var, 1e-6)
return mu, var
# データ生成
N = 500
X_train = np.sort(np.random.uniform(-5, 5, N)).reshape(-1, 1)
y_true = np.sin(X_train.ravel()) + 0.5 * np.cos(2 * X_train.ravel())
y_train = y_true + np.random.normal(0, 0.3, N)
X_test = np.linspace(-6, 6, 200).reshape(-1, 1)
# 標準GP
mu_full, var_full = gp_predict(X_train, y_train, X_test, noise=0.3)
# 疎GP（誘導点数 M=20）
M = 20
Z = np.linspace(-5, 5, M).reshape(-1, 1)
sgp = SparseGP(Z, noise=0.3, length_scale=1.0)
elbo_hist = sgp.fit(X_train, y_train, n_iter=50)
mu_sparse, var_sparse = sgp.predict(X_test)
# 可視化
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# 左: 標準GP
axes[0].scatter(X_train, y_train, s=3, alpha=0.3, color='gray', label='Training data')
axes[0].plot(X_test, mu_full, 'b-', linewidth=2, label='Full GP mean')
axes[0].fill_between(X_test.ravel(),
mu_full - 2*np.sqrt(var_full),
mu_full + 2*np.sqrt(var_full),
alpha=0.2, color='blue')
axes[0].set_title(f'Full GP (N={N})')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].legend()
# 中央: 疎GP
axes[1].scatter(X_train, y_train, s=3, alpha=0.3, color='gray', label='Training data')
axes[1].plot(X_test, mu_sparse, 'r-', linewidth=2, label='Sparse GP mean')
axes[1].fill_between(X_test.ravel(),
mu_sparse - 2*np.sqrt(var_sparse),
mu_sparse + 2*np.sqrt(var_sparse),
alpha=0.2, color='red')
axes[1].scatter(Z, np.zeros(M), marker='^', s=100, color='green',
zorder=5, label=f'Inducing points (M={M})')
axes[1].set_title(f'Sparse GP (M={M})')
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].legend()
# 右: ELBO推移
axes[2].plot(elbo_hist, 'o-', markersize=3, color='purple')
axes[2].set_xlabel('Iteration')
axes[2].set_ylabel('ELBO')
axes[2].set_title('ELBO Convergence')
plt.tight_layout()
plt.show()
左のグラフは標準ガウス過程（$N=500$ のデータ全てを使用）、中央は疎ガウス過程（$M=20$ の誘導点のみ使用）の結果です。緑の三角が誘導点の位置を示しています。少数の誘導点でも標準GPに近い予測が得られることがわかります。
まとめ
本記事では、疎ガウス過程（Sparse Gaussian Processes）について解説しました。
- 標準ガウス過程は計算量 $O(N^3)$ が問題となり、大規模データに適用できない
- 疎ガウス過程は $M$ 個の誘導点を導入し、計算量を $O(NM^2)$ に削減する
- 変分疎ガウス過程（SVGP）はELBOの最大化として定式化され、誘導点の位置も最適化できる
- 予測分布は $O(M^2)$ で計算でき、誘導点数が少なくても標準GPに近い精度が得られる
次のステップとして、以下の記事も参考にしてください。