スパースガウス過程回帰（誘導点法）の理論とPython実装：Subset of RegressorsでO(n³)の壁を突破する

はじめに

https://yuhi-sa.github.io/posts/20260502_gaussian_process/1/ では、ガウス過程回帰（GPR）が観測データに対して閉形式の事後分布を持つこと、そしてコレスキー分解を使ってもなお学習に \(O(n^3)\) の計算量と \(O(n^2)\) のメモリを要することを見ました。同記事の「スパースGP概観」節では、\(m \ll n\) 個の**誘導点（inducing points）**を導入してカーネル行列を低ランク近似するという発想だけを紹介し、詳細は先送りにしていました。

本記事はその発展編です。誘導点法の中でも最も基本的な Subset of Regressors（SoR） を、ベイズ線形回帰との双対性から数式を導出し、numpyでフルスクラッチ実装します。そのうえで、scikit-learnのGaussianProcessRegressor（厳密GP）を正解として、

1. 誘導点数 \(M\) を増やしたときに予測が厳密GPへ収束していく様子
2. 学習・予測にかかる時間の高速化倍率
3. 誘導点の配置方法（ランダム／等間隔／k-means）が精度に与える影響
4. SoRに内在する予測分散の過小評価という既知の欠陥

を、実際にPythonを実行して得た数値だけで確認します。

誘導点法の基本アイデア

\(n\) 個の訓練点に対するカーネル行列 \(\mathbf{K} \in \mathbb{R}^{n \times n}\) の逆行列計算が \(O(n^3)\) で高コストになる根本原因は、\(\mathbf{K}\) がフルランク（階数 \(n\) ）であることです。誘導点法は、\(m \ll n\) 個の代表点 \(\mathbf{Z} = \{\mathbf{z}_1, \ldots, \mathbf{z}_m\}\) を選び、\(\mathbf{K}\) を階数 \(m\) の行列で近似することでこの問題を回避します：

ここで \(\mathbf{K}_{nm} \in \mathbb{R}^{n \times m}\) は訓練点と誘導点の間のカーネル行列、\(\mathbf{K}_{mm} \in \mathbb{R}^{m \times m}\) は誘導点同士のカーネル行列です。式(1)はNyström近似として知られ、\(\tilde{\mathbf{K}}\) を使う限り逆行列計算は \(m \times m\) 行列（Woodburyの恒等式により実質 \(O(nm^2)\) ）に縮小されます。

この節では、式(1)を天下り的に受け入れるのではなく、ベイズ線形回帰との双対性から同じ結果を導出します。これはhttps://yuhi-sa.github.io/posts/20260502_gaussian_process/1/のRKHS/Representer定理の節で見た「GP事後平均 = カーネルリッジ回帰の解」という双対関係の延長線上にある議論です。

導出：重み空間表現からSoRへ

誘導点 \(\mathbf{Z}\) を固定し、特徴写像

を定義します。これは \(m\) 個の誘導点との類似度を並べた行ベクトルです。ここで、関数を

というベイズ線形回帰としてモデル化します。重み \(\mathbf{w}\) の事前分散をわざわざ \(\mathbf{K}_{mm}^{-1}\) （誘導点間カーネル行列の逆行列）に選ぶのは、こうすると2点間の事前共分散が

となり、\(n\) 個の訓練点全体で評価すればちょうど式(1)の \(\tilde{\mathbf{K}}\) に一致するからです。つまり式(3)は、式(1)のNyström近似を再現する階数 \(m\) の退化GP事前分布を、\(m\) 個のパラメトリックな重みを持つ線形モデルとして書き下したものに他なりません。

あとは通常のベイズ線形回帰の公式をそのまま適用するだけです。ノイズモデル \(y_i = f(\mathbf{x}_i) + \varepsilon_i,\ \varepsilon_i \sim \mathcal{N}(0, \sigma_n^2)\) のもとで、計画行列 \(\Phi = \mathbf{K}_{nm}\) （\(\Phi_{ij} = k(\mathbf{x}_i, \mathbf{z}_j)\) ）、事前精度行列 \(\mathbf{K}_{mm}\) を使うと、事後精度行列は

となり、\(\mathbf{w}\) の事後平均は \(\bar{\mathbf{w}} = \frac{1}{\sigma_n^2}\mathbf{B}^{-1}\mathbf{K}_{mn}\mathbf{y}\) です。したがって新しい入力 \(\mathbf{x}_*\) における予測分布は、\(\mathbf{k}_{*m} := \phi(\mathbf{x}_*) = k(\mathbf{x}_*, \mathbf{Z})\) として、

となります。逆行列は常に \(m \times m\) の \(\mathbf{B}\) に対してのみ計算すればよく、\(\mathbf{B}\) を構成するのに \(O(nm^2)\) 、コレスキー分解に \(O(m^3)\) ── 合わせて \(O(nm^2 + m^3)\) で、厳密GPの \(O(n^3)\) から大幅に削減されます（\(m \ll n\) のとき \(nm^2 \ll n^3\) ）。

式(6)・式(7)がまさにSoR（Subset of Regressors; Silverman, 1985; Smola & Bartlett, 2001）の予測式です。https://yuhi-sa.github.io/posts/20260704_gpr_practice/1/で紹介した実装コード（A = K_mm * sn2 + K_nm.T @ K_nm としてcho_solveする部分）は、\(\mathbf{A} = \sigma_n^2 \mathbf{B}\) の関係でこの式(5)・式(6)と完全に一致します。本記事ではこの式がどこから来るのかを、重み空間の視点から導出しました。

注意（既知の欠陥）：式(7)の分散は、\(\mathbf{x}_*\) がすべての誘導点から離れると \(\mathbf{k}_{*m} \to \mathbf{0}\) となり、\(\sigma^2_{\text{SoR}}(\mathbf{x}_*) \to 0\) に潰れます。厳密GPでは訓練データから離れるほど事前分散 \(\sigma_f^2\) に近づいて不確実性が増えるのに対し、SoRは逆に「誘導点から離れると自信満々になる」という直感に反する挙動を示します。これがFITC（Snelson & Ghahramani, 2006）が対角補正項を導入する動機です。後の実験でこの挙動を実際に可視化して確認します。

Python実装

式(5)〜(7)をそのままnumpyとscipyでフルスクラッチ実装します。数値安定性のため、\(\mathbf{B}\) の逆行列も明示的には計算せず、https://yuhi-sa.github.io/posts/20260502_gaussian_process/1/で解説したのと同じ理由でコレスキー分解（cho_factor/cho_solve）を使います。

import numpy as np
from scipy.linalg import cho_factor, cho_solve


class SparseGPRegressor:
    """Subset of Regressors (SoR) — 式(5)〜(7)のフルスクラッチ実装"""

    def __init__(self, length_scale=1.0, signal_var=1.0, noise_var=1e-2, jitter=1e-8):
        self.l = length_scale
        self.sf2 = signal_var
        self.sn2 = noise_var
        self.jitter = jitter

    def rbf(self, X1, X2):
        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, Z):
        """X: (n,1) 訓練入力, y: (n,), Z: (m,1) 誘導点"""
        self.Z = Z.copy()
        m = len(Z)
        K_mm = self.rbf(Z, Z) + self.jitter * np.eye(m)
        K_mn = self.rbf(Z, X)                       # (m, n)

        B = K_mm + (K_mn @ K_mn.T) / self.sn2        # 式(5): (m, m)
        self.B_chol_ = cho_factor(B, lower=True)

        Kmn_y = K_mn @ y                             # (m,)
        self.w_mean_ = cho_solve(self.B_chol_, Kmn_y) / self.sn2

    def predict(self, X_new):
        k_star_m = self.rbf(X_new, self.Z)           # (n_new, m)
        mu = k_star_m @ self.w_mean_                 # 式(6)
        v = cho_solve(self.B_chol_, k_star_m.T)       # (m, n_new)
        var = np.sum(k_star_m.T * v, axis=0)          # 式(7): k_*m^T B^-1 k_*m
        return mu, np.maximum(var, 1e-12)


def select_inducing_points(X, m, method="random", rng=None):
    rng = rng or np.random.default_rng(0)
    x = X.ravel()
    if method == "random":
        idx = rng.choice(len(x), size=m, replace=False)
        return X[idx].copy()
    if method == "uniform":
        return np.linspace(x.min(), x.max(), m).reshape(-1, 1)
    if method == "kmeans":
        from sklearn.cluster import KMeans
        km = KMeans(n_clusters=m, n_init=10, random_state=0).fit(X)
        return km.cluster_centers_
    raise ValueError(method)

fitが構築する行列は \(\mathbf{B}\) （\(m \times m\) ）だけで、\(n \times n\) の行列は一度も明示的に作りません。K_mn @ K_mn.Tの計算量は \(O(nm^2)\) ですが、これが誘導点法全体のボトルネックになります。

検証：誘導点数を訓練点数と一致させると厳密GPに一致するか

式(6)・式(7)の導出が正しければ、誘導点をすべての訓練点にする（\(\mathbf{Z} = \mathbf{X}\) 、\(m = n\) ）と、SoRは近似ではなく厳密GPと完全に一致するはずです（\(\mathbf{K}_{mm} = \mathbf{K}_{nn} = \mathbf{K}\) となり低ランク近似が恒等写像に退化するため）。これを実際に確認しました。

rng = np.random.default_rng(0)
n = 40
X = rng.uniform(-3, 3, size=(n, 1))
y = np.sin(X.ravel()) + 0.1 * rng.standard_normal(n)
Xs = np.linspace(-3, 3, 25).reshape(-1, 1)

exact = ExactGP(length_scale=1.0, signal_var=1.0, noise_var=0.01)   # 通常のGP（cho_factorベース）
exact.fit(X, y)
mu_exact, var_exact = exact.predict(Xs)

sor = SparseGPRegressor(1.0, 1.0, 0.01, jitter=1e-10)
sor.fit(X, y, Z=X)   # m == n
mu_sor, var_sor = sor.predict(Xs)

print("max |mu_exact - mu_sor|   =", np.max(np.abs(mu_exact - mu_sor)))
print("max |var_exact - var_sor| =", np.max(np.abs(var_exact - var_sor)))

実行結果：

max |mu_exact - mu_sor|   = 2.529596088152175e-09
max |var_exact - var_sor| = 4.067293168930064e-10

差はどちらも \(10^{-9}\) 〜\(10^{-10}\) オーダーで、これは浮動小数点誤差の範囲です。式(6)・式(7)の導出が正しいことがコードレベルでも裏付けられました。

実験1：誘導点数MによるRMSE収束と高速化の実測

いよいよ本題です。\(n=5000\) の1次元回帰データで、scikit-learnのGaussianProcessRegressor（厳密GP、optimizer=Noneでハイパーパラメータは固定）を正解として、誘導点数 \(M \in \{10, 25, 50, 100, 200\}\) （配置はk-means）でSoRを学習・予測し、(a) 真の関数へのRMSE、(b) 厳密GP予測へのRMSE、(c) 学習+予測にかかる時間、を比較しました。

import time
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel, WhiteKernel

N, LENGTH_SCALE, SIGNAL_VAR, NOISE_STD = 5000, 0.7, 1.0, 0.2
NOISE_VAR, DOMAIN, N_TEST = NOISE_STD**2, (-10.0, 10.0), 500

rng = np.random.default_rng(42)

def true_f(x):
    return np.sin(x) + 0.5 * np.sin(3.0 * x) + 0.1 * x

def rmse(a, b):
    return float(np.sqrt(np.mean((a - b) ** 2)))

X_train = rng.uniform(*DOMAIN, size=(N, 1))
y_train = true_f(X_train.ravel()) + NOISE_STD * rng.standard_normal(N)
X_test = np.linspace(*DOMAIN, N_TEST).reshape(-1, 1)
f_test_true = true_f(X_test.ravel())

kernel = ConstantKernel(SIGNAL_VAR, "fixed") * RBF(length_scale=LENGTH_SCALE, length_scale_bounds="fixed") \
         + WhiteKernel(noise_level=NOISE_VAR, noise_level_bounds="fixed")

gp = GaussianProcessRegressor(kernel=kernel, optimizer=None)
t0 = time.perf_counter(); gp.fit(X_train, y_train); t_fit = time.perf_counter() - t0
t0 = time.perf_counter(); mu_exact, _ = gp.predict(X_test, return_std=True); t_pred = time.perf_counter() - t0
# ... 7回反復し中央値を採用（timing_probeと同一手法）

for M in [10, 25, 50, 100, 200]:
    Z = select_inducing_points(X_train, M, method="kmeans", rng=np.random.default_rng(0))
    sgp = SparseGPRegressor(LENGTH_SCALE, SIGNAL_VAR, NOISE_VAR, jitter=1e-8)
    sgp.fit(X_train, y_train, Z)
    mu_sor, _ = sgp.predict(X_test)
    # RMSE, 実行時間を記録（詳細は本文の表を参照）

実行結果（学習+予測時間は7回実行の中央値）：

\(n=5000\) に対して誘導点をわずか2%（\(M=100\) ）にしただけで、厳密GPとの予測平均の差は \(4\times10^{-7}\) まで縮み、実用上は厳密GPと区別がつかなくなります。一方 \(M=10\) （0.2%）では真の関数へのRMSEが厳密GP比で約23倍悪化しており、明らかに近似が粗すぎます。速度面では、\(M=100\) でも厳密GPの337倍高速、\(M=10\) の極端な設定では2042倍高速でした。\(M\) を増やすほど精度は上がりますが、\(\mathbf{K}_{mn}\mathbf{K}_{nm}\) の構築コスト（\(O(nm^2)\) ）のため実行時間も \(M\) にほぼ2乗で効いてきて、高速化倍率は徐々に下がっていきます——これは理論通りの \(O(nm^2 + m^3)\) vs \(O(n^3)\) のトレードオフです。

左図は縦軸が対数スケールで、SoRの厳密GPへのRMSE（赤）が \(M\) の増加とともに指数的に縮小し、真の関数へのRMSE（青）は \(M=50\) 付近で厳密GPの水準（黒破線）に達して頭打ちになる様子を示しています。これは、\(M=50\) 以降の誤差はもはや近似誤差ではなく、データに含まれるノイズ由来の本質的な下限であることを意味します。右図は同じ \(M\) に対する学習+予測の実行時間で、青の折れ線に付した倍率が厳密GP（黒破線）に対する高速化率です。

実験2：誘導点の配置方法（ランダム／等間隔／k-means）の比較

誘導点をどこに置くかも近似精度を左右します。あえて誘導点数を絞った \(M=25\) （実験1でまだ厳密GPに完全収束していない設定）で、ランダム選択・等間隔（np.linspace）・k-meansクラスタリング中心の3方式を比較しました。ランダム選択は結果が乱数シードに依存するため、20個のシードで平均・標準偏差を取っています。

M = 25
Z_uniform = select_inducing_points(X_train, M, method="uniform")
Z_kmeans  = select_inducing_points(X_train, M, method="kmeans")

rmse_random = []
for seed in range(20):
    Z_random = select_inducing_points(X_train, M, method="random", rng=np.random.default_rng(seed))
    sgp = SparseGPRegressor(LENGTH_SCALE, SIGNAL_VAR, NOISE_VAR, jitter=1e-8).fit(X_train, y_train, Z_random)
    mu, _ = sgp.predict(X_test)
    rmse_random.append(rmse(mu, f_test_true))

実行結果：

exact GP RMSE_vs_true = 0.01765
M=25
  uniform : RMSE_vs_true=0.03633  RMSE_vs_exact=0.03886
  kmeans  : RMSE_vs_true=0.02879  RMSE_vs_exact=0.02683
  random  : RMSE_vs_true=0.25582 +/- 0.10313 (best=0.07622, worst=0.52803, n=20)

k-means（RMSE 0.0288）が等間隔配置（0.0363）よりわずかに優れ、ランダム選択（平均0.2558、標準偏差0.1031）は両者よりはるかに悪く、しかもシードによるばらつきが非常に大きいことがわかります。今回の実験データは入力を一様分布からランダムに生成しているため、\(M=25\) 個をランダムに間引くと、20区間中いくつかの領域に誘導点がまったく置かれない「穴」が生じやすくなります。等間隔配置とk-meansはどちらもこの穴を避けるように働くため安定して低いRMSEを達成しますが、ランダム選択は運が悪いと（今回の実験のworst=0.528のように）厳密GPの30倍近いRMSEになることもあります。

右図のヒストグラムは訓練入力の密度（一様分布なのでほぼ平坦）、その上の3本の記号列が各手法の誘導点位置です。ランダム選択（オレンジ）は \(x \approx -1\) 付近など明らかに間隔が空いている箇所があるのに対し、等間隔（緑）とk-means（紫）はどちらも定義域全体にほぼ均等に配置されていることが視覚的にも確認できます。

SoRの既知の欠陥：誘導点から離れた場所での分散の過小評価

導出の節で触れた通り、式(7)の予測分散は \(\mathbf{k}_{*m} \to \mathbf{0}\) のとき \(0\) に潰れます。これを実際に確認するため、誘導点をわざと定義域の左半分 \([-10, -3]\) にしか置かず（\(M=15\) 、等間隔）、訓練データ自体は定義域全体 \([-10, 10]\) に一様に存在する（つまり厳密GPは右半分でも十分な情報を持つ）状況で、右半分 \(x > 2\) での予測標準偏差を厳密GPとSoRで比較しました。

Z = np.linspace(-10, -3, 15).reshape(-1, 1)   # 誘導点は左半分だけに配置
sgp = SparseGPRegressor(LENGTH_SCALE, SIGNAL_VAR, NOISE_VAR, jitter=1e-8).fit(X_train, y_train, Z)
mu_sor, var_sor = sgp.predict(X_test)
std_sor = np.sqrt(var_sor)

mask_right = X_test.ravel() > 2.0
print("mean std (exact GP), x>2  :", std_exact[mask_right].mean())
print("mean std (SoR),      x>2  :", std_sor[mask_right].mean())
print("mean std (SoR),      x<-8 :", std_sor[X_test.ravel() < -8].mean())

実行結果：

mean std (exact GP), x>2   (no inducing pts there): 0.2009
mean std (SoR),      x>2   (no inducing pts there): 0.0000
mean std (SoR),      x<-8  (near inducing pts):      0.0190

誘導点の近く（\(x < -8\) ）ではSoRの標準偏差0.0190は厳密GPと同程度の水準ですが、誘導点がない右半分（\(x > 2\) ）では厳密GPが0.2009と相応の不確実性を保持しているのに対し、SoRは0.0000まで潰れています。訓練データ自体は右半分にも豊富にあるにもかかわらず、SoRが「誘導点から離れている」というだけの理由で過信するのは危険な挙動です。

図の右半分でSoRの平均（赤線）が0に張り付き、信頼区間（赤帯）がほぼ見えなくなっているのに対し、厳密GP（黒線・灰色帯）は真の関数の変動を捉え続けています。この欠陥を対角補正で緩和するのがFITC（Snelson & Ghahramani, 2006）で、さらに誘導点自体を変分下限の最大化で最適配置するのがVFE（Titsias, 2009）、それをミニバッチ化して \(n \sim 10^6\) 超に拡張するのがSVGP（Hensman, Fusi & Lawrence, 2013）です。これらはhttps://yuhi-sa.github.io/posts/20260502_gaussian_process/1/の「スパースGP概観」節で名前だけ触れた手法ですが、本記事で見た「SoRの分散はなぜ潰れるのか」を理解していれば、FITC以降の改良がどの問題を解いているのかが自然につながります。

まとめ：実務での使い分け

SoRは実装がもっとも単純で、\(n \times n\) 行列を一度も作らずに厳密GPを高い精度で近似できることが、本記事の実測から確認できました。ただし予測分散の較正には弱点があるため、不確実性の値そのものを意思決定に使う場面（ベイズ最適化の獲得関数など）では、FITCやVFE系の手法への切り替えを検討する価値があります。

関連記事

• ガウス過程回帰の基礎とPython実装 - 本記事の前提となる、厳密GPの理論・導出・数値安定性を解説した基礎編です。
• ガウス過程回帰の実践：カーネル設計と応用 - カーネル設計・ハイパーパラメータ最適化を含む実務的な観点からSoRを含む近似手法を概観しています。本記事はそこで簡潔に触れられたSoRの導出と実験を深掘りする内容です。
• ベイズ最適化の基礎とPython実装 - GPをサロゲートモデルとして使うベイズ最適化では、評価回数が数百〜数千を超えると本記事のスケーラビリティの議論が直接関係してきます。
• SVMのカーネル設計：Mercerの定理・グラム行列の正定値性とRandom Fourier Features - Random Fourier Featuresはカーネル行列自体を低ランクな明示的特徴写像で近似する手法で、誘導点法とは異なる角度からカーネル法の\(O(n^2)\) 〜\(O(n^3)\) の壁に挑む「親戚」にあたります。

参考文献

• Rasmussen, C. E., & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press.
• Silverman, B. W. (1985). “Some Aspects of the Spline Smoothing Approach to Non-Parametric Regression Curve Fitting.” Journal of the Royal Statistical Society: Series B, 47(1), 1-21.
• Smola, A. J., & Bartlett, P. (2001). “Sparse Greedy Gaussian Process Regression.” NeurIPS 2000.
• Snelson, E., & Ghahramani, Z. (2006). “Sparse Gaussian Processes using Pseudo-inputs.” NeurIPS 2006. （FITC）
• 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）
• Quiñonero-Candela, J., & Rasmussen, C. E. (2005). “A Unifying View of Sparse Approximate Gaussian Process Regression.” Journal of Machine Learning Research, 6, 1939-1959.