10万件のセンサ計測、1000万枚の衛星画像ピクセル、1日数億トランザクションの株価ティック——現代の機械学習が向き合うデータ規模は、ガウス過程（GP）の本来の前提を完全に踏み越えています。GP回帰は「データ点が増えれば自動的に賢くなる」ノンパラメトリックな表現力と、予測の不確実性まで自然に出力する稀有な性質を持ちます。にもかかわらず、$N \times N$ のグラム行列を作り、その逆行列とコレスキー分解を要求するため、計算量は $O(N^3)$、メモリは $O(N^2)$。$N = 10000$ で既にメモリ800MB・学習時間は分オーダー、$N = 10^5$ ともなれば普通のワークステーションでは扱えません。
この壁を破るのが SVGP（Stochastic Variational Gaussian Process） です。SVGP は「全データ点ではなく、データの本質を $M \ll N$ 個の 誘導点（inducing points） で要約する」というスパース化と、「ELBO の確率的な勾配でミニバッチ最適化する」変分推論を組み合わせ、計算量を $O(NM^2)$、メモリを $O(M^2)$ に圧縮します。$M = 200$ なら $N = 10^7$ クラスまで現実的に学習でき、しかも GP 本来の不確実性定量化はそのまま保たれます。
この技術は実応用で広く使われています。たとえば衛星リモートセンシングでは1シーン数億ピクセルの土地被覆や気温分布の不確実性付き推定、宇宙ロボットの動力学学習では数万エピソード分の (state, action) ペアからの方策学習、気象予報では地表観測点のスパース化と長期トレンド推定、株価モデリングでは数年分の高頻度ティックデータからのボラティリティ推定など、いずれも「GP の表現力は欲しいが $N^3$ には耐えられない」典型例です。
本記事の内容
- GP 回帰の計算量問題と、スパースGP 系統の歴史的整理（SoR, FITC, VFE, SVGP）
- 誘導点 $\bm{u}$ の役割と「データの要約」としての解釈
- 変分推論によるELBOの導出（データ項 + KL項への分解）
- SVGP の最終形 $q(\bm{u}) = \mathcal{N}(\bm{m}, \bm{S})$ と確率勾配法での最適化
- Reparameterization トリック、自然勾配法、誘導点配置の戦略
- 多重出力GP、Deep GPs、Conditional Neural Processes など発展トピック
- PyTorch でのスクラッチ実装と、GPyTorch ベンチマーク、誘導点数 $M$ のトレードオフ
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
GP回帰の計算量問題 — なぜそのままではスケールしないのか
ガウス過程回帰の典型的な推論式を思い出します。観測データ $\mathcal{D} = \{(\bm{x}_i, y_i)\}_{i=1}^N$、ノイズ分散 $\sigma_n^2$、カーネル $k$、グラム行列 $\bm{K}_{NN} = [k(\bm{x}_i, \bm{x}_j)]$ に対して、新しい入力 $\bm{x}_*$ における予測分布は
$$ \begin{aligned} \mu_*(\bm{x}_*) &= \bm{k}_*^\top (\bm{K}_{NN} + \sigma_n^2 \bm{I})^{-1} \bm{y}, \\ \sigma_*^2(\bm{x}_*) &= k(\bm{x}_*, \bm{x}_*) – \bm{k}_*^\top (\bm{K}_{NN} + \sigma_n^2 \bm{I})^{-1} \bm{k}_* \end{aligned} $$
で与えられます。ここで $\bm{k}_* = [k(\bm{x}_*, \bm{x}_i)]_{i=1}^N$ は $\bm{x}_*$ と訓練点とのカーネル列ベクトルです。式の中央に居座る $(\bm{K}_{NN} + \sigma_n^2 \bm{I})^{-1}$ こそが計算量問題の主犯です。
この逆行列を毎回計算する代わりに、実装ではコレスキー分解 $\bm{K}_{NN} + \sigma_n^2 \bm{I} = \bm{L}\bm{L}^\top$ を取り、$\bm{L}^{-1}\bm{y}$ を前進代入で求めます。コレスキー分解自体が $O(N^3)$、メモリは $O(N^2)$。$N=10^4$ なら浮動小数点1要素8バイトで800MB、$N=10^5$ では80GB。PCのメモリに乗らないうえに、ハイパーパラメータ学習で何度もこの分解を回す必要があるので、実時間で破綻します。
直感的には、GP回帰は「全データ点を全データ点に対して比較する」処理であり、ペア数 $N^2$ が常に登場します。さらに連立方程式を解くために、その $N \times N$ 行列に対する「逆」を取らねばならない。連立方程式の解法は最良でも $O(N^{2.37})$ 程度（Strassen 系の理論的な下界に近い実用アルゴリズムは存在しない）で、実用上は $O(N^3)$ がボトルネックになります。
ところで、$N=10^4$ もの観測点は本当に「全部が独立に情報を持っている」のでしょうか。たとえば衛星画像で隣り合うピクセル100点は、ほぼ同じ情報を冗長に伝えているはずです。GP の予測も、似た入力 $\bm{x}_i \approx \bm{x}_j$ に対しては似た寄与しか与えません。だとすれば、「データ全体の情報を、もっと少数の代表点で要約できないか？」という発想が自然に湧きます。この要約点が、次節で導入する 誘導点 です。
スパースGPの歴史 — SoR から SVGP まで
スパースGPは2000年代から精力的に研究され、いくつかの近似が積み重なって SVGP に至っています。歴史を追うと、現代SVGP の設計判断がなぜ妥当なのかが見えてきます。
最初に登場したのは SoR (Subset of Regressors) で、GP の予測式の中で全データ $\bm{X}$ を $M$ 個の誘導点 $\bm{Z} \in \mathbb{R}^{M \times d}$ で置き換える素朴な近似でした。次にやってきたのが DTC (Deterministic Training Conditional) と FITC (Fully Independent Training Conditional) で、Snelson と Ghahramani（2006）の “Sparse Gaussian Processes using Pseudo-inputs” が決定的でした。FITC は「予測分布の対角成分を残し、非対角は誘導点経由でしか繋がらない」という対角補正を持ち、SoR の予測分散過小評価問題を改善しました。
Quiñonero-Candela と Rasmussen（2005）は “A Unifying View of Sparse Approximate Gaussian Process Regression” で、これらの近似がすべて 「事前分布側を改変する」 形に統一できることを示しました。一方、Titsias（2009）は完全に異なる道を切り開きます。”Variational Learning of Inducing Variables in Sparse Gaussian Processes” で提示された VFE (Variational Free Energy) は、「事前分布をいじるのではなく、事後分布 を変分近似する」立場を取り、誘導点 $\bm{Z}$ をハイパーパラメータではなく変分パラメータとして扱いました。これは過適合を理論的に抑える効果を持ち、現代SVGP の祖となります。
そして Hensman, Fusi, Lawrence（2013）の “Gaussian Processes for Big Data” が決定打を放ちます。VFE の ELBO（変分下界）を、データ点 $i$ ごとの和に分解できる形に書き直し、ミニバッチ確率勾配法で最適化できるようにしました。これが現代の SVGP（Stochastic Variational GP） です。$N = 10^9$ のデータでも、各ステップでミニバッチ $B \ll N$ 件しか触らないので、メモリは $O(M^2 + BM)$ で済みます。
歴史を一文で要約すれば、「事前分布を弄る FITC 系」から「事後分布を変分近似する VFE 系」へ、そして「変分下界をミニバッチ化した SVGP」へ、というのが流れです。次節からは、SVGP の心臓部である誘導点と変分推論をじっくり見ていきます。
誘導点（Inducing points）— データの「要約」とは何か
誘導点を導入する数学的なお膳立てを整えます。観測点での関数値を $\bm{f} = [f(\bm{x}_i)]_{i=1}^N \in \mathbb{R}^N$、誘導点 $\bm{Z} = \{\bm{z}_m\}_{m=1}^M$ での関数値を $\bm{u} = [f(\bm{z}_m)]_{m=1}^M \in \mathbb{R}^M$ とします。$\bm{f}$ と $\bm{u}$ は同じ GP から引いた値なので、両者の同時分布もガウス分布です。
$$ \begin{bmatrix} \bm{f} \\ \bm{u} \end{bmatrix} \sim \mathcal{N}\left(\bm{0}, \begin{bmatrix} \bm{K}_{NN} & \bm{K}_{NM} \\ \bm{K}_{MN} & \bm{K}_{MM} \end{bmatrix}\right) $$
ここで $\bm{K}_{NN} = [k(\bm{x}_i, \bm{x}_j)]$、$\bm{K}_{MM} = [k(\bm{z}_m, \bm{z}_{m’})]$、$\bm{K}_{NM} = [k(\bm{x}_i, \bm{z}_m)]$ です。ガウス分布の条件付き分布の公式から、$\bm{u}$ を所与としたときの $\bm{f}$ は
$$ p(\bm{f} | \bm{u}) = \mathcal{N}(\bm{f} \,|\, \bm{K}_{NM}\bm{K}_{MM}^{-1}\bm{u},\ \bm{K}_{NN} – \bm{K}_{NM}\bm{K}_{MM}^{-1}\bm{K}_{MN}) $$
で与えられます。この式の解釈が重要です。$\bm{K}_{NM}\bm{K}_{MM}^{-1}\bm{u}$ という項は、誘導点での関数値 $\bm{u}$ を「カーネル空間で補間」して観測点に持ってきた値と読めます。$\bm{K}_{NM}\bm{K}_{MM}^{-1}$ がまさに「誘導点から観測点への補間行列」の役割を果たしています。共分散 $\bm{K}_{NN} – \bm{K}_{NM}\bm{K}_{MM}^{-1}\bm{K}_{MN}$ は、補間誤差を表す残差共分散です。
ここでの鍵となる仮定が、SVGP（および VFE/FITC）の核です。「$\bm{u}$ が $\bm{f}$ の十分統計量である」つまり「$\bm{u}$ さえ知っていれば、データ $\bm{y}$ は $\bm{f}$ について何の追加情報も持たない」と仮定します。式で書くと、
$$ p(\bm{y} | \bm{u}) = \int p(\bm{y}|\bm{f}) p(\bm{f}|\bm{u})\,d\bm{f} $$
がデータの条件付き尤度の正しい姿だ、と認めるわけです。この仮定が成り立つのは厳密には $M = N$ かつ $\bm{Z} = \bm{X}$ のときだけですが、$M \ll N$ でも「データの本質的な変動が誘導点で十分捉えられている」ならば、近似は十分機能します。
直感的には、$M$ 個の誘導点はカーネル空間で「データの代表点」として配置され、補間行列 $\bm{K}_{NM}\bm{K}_{MM}^{-1}$ が観測点を再構成します。$M$ が多いほど補間誤差が小さくなり、$M \to N$ で厳密GP に収束します。$M$ をどう選ぶか、$\bm{Z}$ をどこに配置するかが、SVGP 設計の自由度になります。
ここまでで「誘導点が GP のスパース化を可能にする」イメージが固まりました。しかし、誘導点での関数値 $\bm{u}$ は私たちが直接観測できる量ではないので、これを推論する仕組みが必要です。そこで登場するのが変分推論です。
変分推論の枠組み — ELBO とその分解
ベイズ推論の目標は事後分布 $p(\bm{u}|\bm{y})$ を求めることですが、SVGP では真の事後分布の代わりに、扱いやすい変分分布 $q(\bm{u}) = \mathcal{N}(\bm{u}|\bm{m}, \bm{S})$ を用意し、これが真の事後分布に最も近くなるように $(\bm{m}, \bm{S})$ を学習します。「近さ」の指標が KL ダイバージェンスです。
$$ \mathrm{KL}\!\left[q(\bm{u})\,\|\,p(\bm{u}|\bm{y})\right] = \mathbb{E}_{q}\!\left[\log\frac{q(\bm{u})}{p(\bm{u}|\bm{y})}\right] $$
しかし真の事後 $p(\bm{u}|\bm{y})$ は計算困難（規格化定数 $p(\bm{y})$ が要る）です。そこで、対数周辺尤度 $\log p(\bm{y})$ を次のように分解します。
$$ \log p(\bm{y}) = \underbrace{\mathbb{E}_{q(\bm{u})}\!\left[\log\frac{p(\bm{y}, \bm{u})}{q(\bm{u})}\right]}_{\mathcal{L}\ (\text{ELBO})} + \mathrm{KL}\!\left[q(\bm{u})\,\|\,p(\bm{u}|\bm{y})\right] $$
左辺 $\log p(\bm{y})$ は $q$ に依存しない定数。右辺第2項の KL は非負。したがって ELBO（Evidence Lower BOund）$\mathcal{L}$ を $q$ について最大化することは、KL を最小化すること、すなわち $q(\bm{u})$ を真の事後分布に近づけることと等価です。これが変分推論の基本原理です。
ELBO の中身を展開します。同時分布 $p(\bm{y}, \bm{u}) = p(\bm{y}|\bm{u})p(\bm{u})$ なので、
$$ \mathcal{L} = \mathbb{E}_{q(\bm{u})}[\log p(\bm{y}|\bm{u})] – \mathrm{KL}[q(\bm{u})\,\|\,p(\bm{u})] $$
と分解できます。第1項を データ項（expected log-likelihood）、第2項を KL 項（事前分布への近さ） と呼びます。データ項は「変分事後分布の下でデータをどれだけうまく説明できるか」、KL 項は「事前分布から離れすぎないように」という正則化の役割を果たします。両者のバランスが取れる $q$ が ELBO の最大値を実現します。
ここで重要な観察があります。データ尤度 $p(\bm{y}|\bm{u}) = \int p(\bm{y}|\bm{f}) p(\bm{f}|\bm{u}) d\bm{f}$ は積分が残るので、扱いやすい形に変形が必要です。Jensen の不等式を使って下から押さえると、
$$ \log p(\bm{y}|\bm{u}) \geq \mathbb{E}_{p(\bm{f}|\bm{u})}[\log p(\bm{y}|\bm{f})] $$
となり、これを ELBO に代入することで、最終的に データ点ごとの和 に分解できます。
$$ \mathcal{L} = \sum_{i=1}^N \mathbb{E}_{q(f_i)}[\log p(y_i|f_i)] – \mathrm{KL}[q(\bm{u})\,\|\,p(\bm{u})] $$
ここで $q(f_i)$ は変分事後分布 $q(\bm{u})$ から誘導された $f_i$ の周辺分布です。データ項が「$N$ 個の独立項の和」になったことが、後でミニバッチ化を可能にする決定的な性質です。次節では、この $q(f_i)$ を具体的な式に落とします。
SVGPの完成形 — 数式と変分パラメータ
データ項 $\sum_i \mathbb{E}_{q(f_i)}[\log p(y_i|f_i)]$ の中身を具体化します。変分事後分布 $q(\bm{u}) = \mathcal{N}(\bm{u}|\bm{m}, \bm{S})$ と、条件付き分布 $p(\bm{f}|\bm{u}) = \mathcal{N}(\bm{f}|\bm{K}_{NM}\bm{K}_{MM}^{-1}\bm{u}, \bm{K}_{NN} – \bm{Q}_{NN})$（ただし $\bm{Q}_{NN} = \bm{K}_{NM}\bm{K}_{MM}^{-1}\bm{K}_{MN}$）を組み合わせると、$\bm{f}$ の周辺分布は
$$ q(\bm{f}) = \int p(\bm{f}|\bm{u})q(\bm{u})d\bm{u} = \mathcal{N}(\bm{f}|\bm{\mu},\bm{\Sigma}) $$
ただし
$$ \begin{aligned} \bm{\mu} &= \bm{K}_{NM}\bm{K}_{MM}^{-1}\bm{m}, \\ \bm{\Sigma} &= \bm{K}_{NN} – \bm{K}_{NM}\bm{K}_{MM}^{-1}(\bm{K}_{MM} – \bm{S})\bm{K}_{MM}^{-1}\bm{K}_{MN} \end{aligned} $$
となります。各データ点 $i$ について、$q(f_i)$ は $\bm{\mu}$ の $i$ 番目要素 $\mu_i$ と $\bm{\Sigma}$ の $(i,i)$ 対角要素 $\Sigma_{ii}$ から構成される1次元ガウス分布です。
$$ \mu_i = \bm{k}_{iM}^\top \bm{K}_{MM}^{-1}\bm{m}, \quad \Sigma_{ii} = k_{ii} – \bm{k}_{iM}^\top \bm{K}_{MM}^{-1}(\bm{K}_{MM} – \bm{S})\bm{K}_{MM}^{-1}\bm{k}_{iM} $$
ここで $\bm{k}_{iM} = [k(\bm{x}_i, \bm{z}_m)]_{m=1}^M$、$k_{ii} = k(\bm{x}_i, \bm{x}_i)$。重要なのは、$\bm{K}_{MM}^{-1}$ さえ前計算すれば、各 $i$ の計算が $O(M^2)$ で済むことです。$\bm{K}_{MM}^{-1}$ 自体は $O(M^3)$ ですが、これは $N$ に依存しません。
ガウス尤度 $p(y_i|f_i) = \mathcal{N}(y_i|f_i, \sigma_n^2)$ を仮定すると、データ項は閉形式で
$$ \mathbb{E}_{q(f_i)}[\log p(y_i|f_i)] = -\frac{1}{2}\log(2\pi\sigma_n^2) – \frac{1}{2\sigma_n^2}\big[(y_i – \mu_i)^2 + \Sigma_{ii}\big] $$
となります。右辺末尾の $\Sigma_{ii}$ が「変分事後分布の予測分散」を表し、これがあるおかげで不確実性が ELBO に自然に組み込まれます。
KL 項は両方ガウス分布なので閉形式です。
$$ \mathrm{KL}[\mathcal{N}(\bm{m},\bm{S}) \,\|\, \mathcal{N}(\bm{0}, \bm{K}_{MM})] = \frac{1}{2}\!\left[\mathrm{tr}(\bm{K}_{MM}^{-1}\bm{S}) + \bm{m}^\top \bm{K}_{MM}^{-1}\bm{m} – M + \log\frac{|\bm{K}_{MM}|}{|\bm{S}|}\right] $$
これで ELBO 全体が
$$ \mathcal{L} = -\frac{N}{2}\log(2\pi\sigma_n^2) – \frac{1}{2\sigma_n^2}\sum_{i=1}^N\!\left[(y_i – \mu_i)^2 + \Sigma_{ii}\right] – \mathrm{KL}[q(\bm{u})\,\|\,p(\bm{u})] $$
という、$N$ について和を取る形に書けました。ミニバッチ化はここから自然に出てきます。サイズ $B$ のミニバッチ $\mathcal{B} \subset \{1,\dots,N\}$ を取り、データ項を
$$ \hat{\mathcal{L}}_{\mathrm{data}} = \frac{N}{B}\sum_{i\in\mathcal{B}}\mathbb{E}_{q(f_i)}[\log p(y_i|f_i)] $$
と推定すれば、$\mathcal{L}$ の不偏な確率的推定量が得られ、確率勾配法（SGD, Adam）で最適化できます。学習パラメータは：
- 変分パラメータ: $\bm{m} \in \mathbb{R}^M$、$\bm{S} \in \mathbb{R}^{M \times M}$（正定値、コレスキー因子で表現するのが定石）
- 誘導点位置: $\bm{Z} \in \mathbb{R}^{M \times d}$
- カーネルハイパーパラメータ: $\sigma_f, \ell$ など
- 観測ノイズ: $\sigma_n^2$
これらをすべて勾配で同時最適化する——これが SVGP の完成形です。
Reparameterization と自然勾配法
SVGP は ELBO に解析形があり、ガウス尤度の場合は積分が閉形式で済むので、reparameterization トリックは厳密には不要です。しかし、非ガウス尤度（分類のベルヌーイ尤度、計数のポアソン尤度など）の場合、データ項 $\mathbb{E}_{q(f_i)}[\log p(y_i|f_i)]$ が閉形式で計算できなくなり、モンテカルロ近似が必要になります。このとき威力を発揮するのが reparameterization です。
$q(f_i) = \mathcal{N}(\mu_i, \Sigma_{ii})$ から $f_i$ をサンプルするのに、$\epsilon \sim \mathcal{N}(0, 1)$ を引いて
$$ f_i = \mu_i + \sqrt{\Sigma_{ii}}\,\epsilon $$
と変形します。こうすると $f_i$ は決定論的に $(\mu_i, \Sigma_{ii})$ の関数になり、ノイズ $\epsilon$ は分布パラメータと独立になります。勾配 $\nabla_{\bm{m}, \bm{S}, \bm{Z}, \theta}$ がサンプルを通して伝播でき、低分散な勾配推定量が得られます。VAE で広く知られるトリックですが、SVGP でも同じ原理で動作します。
もう一つ重要な工夫が 自然勾配法（natural gradient） です。Hensman らの SVGP は、変分パラメータ $(\bm{m}, \bm{S})$ の更新に通常の勾配ではなく、フィッシャー情報行列に基づく自然勾配を使うと、収束が劇的に速くなることを示しました。直感的には、$(\bm{m}, \bm{S})$ で張られるパラメータ空間は単純なユークリッド空間ではなく、確率分布のリーマン多様体なので、その幾何に沿った「真の最急方向」を辿るべきだ、という発想です。GPyTorch では NaturalVariationalDistribution として実装されており、典型的に普通のAdamより数倍速く収束します。
ハイパーパラメータ（カーネル長さスケール、ノイズ）と変分パラメータでステップサイズが大きく違うことが多いので、実装上は2つの最適化器を用意し、それぞれ別レートで更新する two-optimizer scheme が定番です。これを使うと、自然勾配で変分パラメータ、Adamでハイパーパラメータ、という分業が綺麗に決まります。
理論はここまでで完結です。次は、SVGP を実際に動かすうえで最も影響の大きい設計判断——誘導点をどこに配置するか——を見ていきます。
誘導点の選び方 — 配置戦略の比較
誘導点 $\bm{Z}$ は変分パラメータとして勾配で動かせますが、初期配置と数 $M$ の選び方が性能を大きく左右します。代表的な戦略を整理します。
ランダム選択: 訓練データから一様サブサンプリングする最単純な方法。実装は楽ですが、データ分布が不均一だと辺縁の予測が荒れます。k-meansクラスタリングは訓練データを $M$ クラスタに分け、各中心を $\bm{Z}$ にする方法。sklearn.cluster.KMeans を1行呼ぶだけで、ランダム選択より明らかに良い初期点を与えます。SVGPの初期化として広く採用されています。
勾配ベース最適化: $\bm{Z}$ を変分パラメータとして ELBO 勾配で動かす方法。Titsias の VFE で誘導点最適化が ELBO を理論的に厳密化することが保証されており、SVGP でも同様です。学習中に $\bm{Z}$ がデータの密集領域や関数の急峻な領域に集中する様子が観察できます。実用上は k-means 初期化 + 勾配微調整の併用がベストプラクティスです。
このほか、$\bm{K}_{MM}$ の条件数を保つため誘導点同士が近すぎないように制約する手法や、訓練データから ELBO 増加最大の点を逐次追加する Greedy 選択 もあります。$M$ の選び方は関数の滑らかさとデータ次元に依存し、経験則として1次元なら $M=20$〜$100$、$d=10$ 次元なら $M=200$〜$1000$、画像のような高次元では $M=1000$〜$5000$ が出発点。$M$ を増やしても ELBO が頭打ちになる点を実験的に探すのが定石です。
ここまでで SVGP の本体は出揃いました。実応用での拡張として、多重出力GP、Deep GPs、Conditional Neural Processes との関係を次に整理します。
拡張トピック — 多重出力GP、Deep GPs、CNPs
多重出力GP（MOGP）への拡張: 出力が $D$ 次元ベクトル $\bm{y}_i \in \mathbb{R}^D$ の場合、各出力を独立な GP でモデリングするのが最も素朴ですが、出力間の相関を捉えたいことが多くあります。Linear Model of Coregionalization (LMC) や Intrinsic Coregionalization Model (ICM) は、$D$ 個の潜在 GP を線形結合して各出力を作る枠組みで、SVGP との組み合わせが自然です。GPyTorch の MultitaskMultivariateNormal と LMCVariationalStrategy でこの構造を直接表現できます。
Deep Gaussian Processes (DGP): GP を層状に重ねたモデルで、$f_2(f_1(\bm{x}))$ のように関数合成で深い表現を作ります。各層の出力を次層の入力とすることで、ニューラルネットの「層」と「非線形性」を GP で実現します。学習は SVGP を各層に適用し、層ごとに変分分布を持たせます。Damianou と Lawrence（2013）が提案し、Salimbeni と Deisenroth（2017）の Doubly Stochastic Variational Inference でスケーラブルに学習可能になりました。1層 SVGP より表現力が高い反面、層数を増やすと最適化が難しくなる傾向があります。
Conditional Neural Processes (CNPs): Garnelo らが2018年に提案した、GP の「文脈点で関数を条件付けて新点を予測する」性質をニューラルネットで模倣する枠組み。文脈点をエンコーダで埋め込みベクトルに集約し、デコーダで予測分布を生成します。SVGP の誘導点 $\bm{Z}$ と「データ要約」という思想を共有しており、CNPs は実質的に「学習可能な誘導点要約 + ニューラル補間」と見ることもできます。Attention 機構を導入した Attentive NP や Latent NP など派生も多く、近年 GP との橋渡しが進んでいます。
これらの拡張に共通するのは、「SVGP の ELBO 構造を保ったまま、追加の構造（出力相関、層、潜在変数）を載せていく」という設計です。SVGP のスケーラビリティが、GP モデル全体を現代の機械学習エコシステムに乗せる土台になっています。応用先も急速に広がっており、衛星リモートセンシングでの土地被覆推定、ロボット動力学学習、気象予報のダウンスケーリング、株価のボラティリティ推定、大規模電子カルテからの病態推移モデリングなど、いずれも $N$ が数十万〜数億のスケールで SVGP が前提技術となっています。
理論と応用を見渡したところで、いよいよ手を動かす番です。次節以降、PyTorch でのスクラッチ実装、GPyTorch でのベンチマーク、誘導点数 $M$ のトレードオフ分析、多次元データへの拡張を順に試していきます。
Pythonでの実装 — PyTorchでスクラッチSVGP
まず最も基本となる SVGP の本体を PyTorch で実装します。理論の式と1対1で対応させ、ELBO の各項が何かを明示的にコメントします。
import torch, torch.nn as nn, math
class SVGP(nn.Module):
"""PyTorchによるSVGP（ガウス尤度・RBFカーネル）
変分: m (M,), L_chol (M,M)下三角で S = L L^T。誘導点Z, カーネル/ノイズも同時学習"""
def __init__(self, M, d, Z_init=None):
super().__init__()
self.M, self.d = M, d
if Z_init is None: Z_init = torch.randn(M, d)
self.Z = nn.Parameter(Z_init.clone()) # 誘導点 (M, d)
self.m = nn.Parameter(torch.zeros(M)) # 変分平均
self.L_chol = nn.Parameter(torch.eye(M) * 0.1) # 変分共分散コレスキー因子
self.log_sigma_f = nn.Parameter(torch.zeros(1)) # カーネル振幅 log
self.log_ell = nn.Parameter(torch.zeros(1)) # 長さスケール log
self.log_sigma_n = nn.Parameter(torch.zeros(1) - 1.0)# 観測ノイズ log
self.jitter = 1e-4
def kernel(self, A, B):
"""RBFカーネル k(a,b) = sigma_f^2 exp(-||a-b||^2 / (2 ell^2))"""
sf2 = self.log_sigma_f.exp().pow(2)
ell2 = self.log_ell.exp().pow(2)
d2 = (A.unsqueeze(1) - B.unsqueeze(0)).pow(2).sum(-1)
return sf2 * torch.exp(-0.5 * d2 / ell2)
ここまでがモデル本体の宣言です。L_chol を下三角行列に明示せずパラメータとして持っているのは、PyTorch の自動微分で扱いやすくするためです。実際に使うときに torch.tril() を取って下三角に落とします。jitter は $\bm{K}_{MM}$ にわずかに対角を加える数値安定化テクニックで、コレスキー分解の失敗を防ぎます。
次に ELBO の計算本体を実装します。データ項と KL 項を別々に出すことで、デバッグ時に各項がどう動いているか追えます。
def elbo(self, X, y, N_total):
"""ミニバッチ (X,y) のELBO不偏推定量。N_total: 全データ点数"""
B = X.size(0)
L = torch.tril(self.L_chol); S = L @ L.t()
sigma_n2 = self.log_sigma_n.exp().pow(2)
K_MM = self.kernel(self.Z, self.Z) + self.jitter * torch.eye(self.M)
K_NM = self.kernel(X, self.Z)
k_ii = self.log_sigma_f.exp().pow(2) * torch.ones(B) # RBFの対角
L_MM = torch.linalg.cholesky(K_MM)
A = torch.linalg.solve_triangular(L_MM, K_NM.t(), upper=False) # (M,B)
K_MM_inv_m = torch.cholesky_solve(self.m.unsqueeze(1), L_MM).squeeze(1)
mu = K_NM @ K_MM_inv_m # (B,)
# Sigma_ii = k_ii - k_iM^T K_MM^{-1} k_iM + k_iM^T K_MM^{-1} S K_MM^{-1} k_iM
Q_diag = (A * A).sum(0)
B_mat = torch.cholesky_solve(K_NM.t(), L_MM)
Sigma_ii = k_ii - Q_diag + (B_mat * (S @ B_mat)).sum(0)
# データ項（ガウス尤度の期待値・閉形式）
ll = -0.5*math.log(2*math.pi*sigma_n2.item()) \
- 0.5*((y - mu).pow(2) + Sigma_ii) / sigma_n2
data_term = (N_total / B) * ll.sum()
# KL項: KL(N(m,S) || N(0,K_MM))
tr_term = torch.diagonal(torch.cholesky_solve(S, L_MM)).sum()
quad = self.m @ torch.cholesky_solve(self.m.unsqueeze(1), L_MM).squeeze(1)
logdet_K = 2*torch.diagonal(L_MM).log().sum()
logdet_S = 2*torch.diagonal(L).log().sum()
kl = 0.5 * (tr_term + quad - self.M + logdet_K - logdet_S)
return data_term - kl
実装のポイントを3つ補足します。1つ目は torch.linalg.cholesky と cholesky_solve を使った数値安定化です。$\bm{K}_{MM}^{-1}$ を明示的に計算せず、コレスキー分解経由で連立方程式を解くことで条件数の悪さに耐えます。2つ目は (N_total / B) のスケーリングで、ミニバッチ ELBO を全データ ELBO の不偏推定量にしています。3つ目は Sigma_ii の計算で、$\bm{K}_{NN}$ の対角だけを使う $O(BM^2)$ 計算に抑えています。
予測関数も実装しておきます。$q(f_*) = \mathcal{N}(\mu_*, \sigma_*^2)$ をテスト点で評価します。
def predict(self, X_star):
"""テスト点 X_star (B*,d) での予測分布 (mu_*, sigma_*^2)"""
with torch.no_grad():
L = torch.tril(self.L_chol); S = L @ L.t()
K_MM = self.kernel(self.Z, self.Z) + self.jitter * torch.eye(self.M)
K_sM = self.kernel(X_star, self.Z)
k_diag = self.log_sigma_f.exp().pow(2) * torch.ones(X_star.size(0))
L_MM = torch.linalg.cholesky(K_MM)
mu_s = K_sM @ torch.cholesky_solve(self.m.unsqueeze(1), L_MM).squeeze(1)
B_mat = torch.cholesky_solve(K_sM.t(), L_MM)
A = torch.linalg.solve_triangular(L_MM, K_sM.t(), upper=False)
sigma_s2 = k_diag - (A*A).sum(0) + (B_mat * (S @ B_mat)).sum(0) \
+ self.log_sigma_n.exp().pow(2)
return mu_s, sigma_s2
予測分散には観測ノイズ $\sigma_n^2$ を最後に足しています。これは $y_* = f_* + \epsilon$ のノイズ込みの予測分布を返したいためです。関数 $f_*$ のみの不確実性が欲しい場合は最後の項を除外します。
N=10000データでの学習
合成データで SVGP の挙動を確かめます。1次元の非線形関数にノイズを乗せ、$N=10000$ 点で学習します。
import torch
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
torch.manual_seed(0); np.random.seed(0)
# 合成データ N=10000
N = 10000
X_np = np.random.uniform(-5, 5, (N, 1)).astype(np.float32)
y_np = (np.sin(X_np[:,0]) + 0.3*X_np[:,0] + 0.5*np.cos(2*X_np[:,0])
+ 0.2*np.random.randn(N)).astype(np.float32)
X, y = torch.from_numpy(X_np), torch.from_numpy(y_np)
# 誘導点を k-means で初期化
M = 64
Z_init = torch.from_numpy(
KMeans(n_clusters=M, n_init=10, random_state=0).fit(X_np)
.cluster_centers_.astype(np.float32))
model = SVGP(M=M, d=1, Z_init=Z_init)
opt = torch.optim.Adam(model.parameters(), lr=1e-2)
losses = []
for it in range(2000):
idx = torch.randint(0, N, (256,))
opt.zero_grad()
loss = -model.elbo(X[idx], y[idx], N_total=N)
loss.backward(); opt.step()
losses.append(loss.item())
if it % 200 == 0: print(f"iter {it:4d} -ELBO={loss.item():.2f}")
print(f"ell={model.log_ell.exp().item():.3f}, sigma_n={model.log_sigma_n.exp().item():.3f}")
このコードの肝は3点です。1つ目、k-means で誘導点を初期化することで、データの分布に沿った配置から始まります。2つ目、torch.randint での簡易ミニバッチサンプリング（厳密にやるなら DataLoader でシャッフル）。3つ目、-model.elbo(...) を損失とする（ELBO は最大化したいので符号反転）。$N=10000$ 点でも、$M=64$ 個の誘導点だけで内部行列を扱うので、1ステップあたり数ms程度で完了します。
学習が終わったら予測を可視化します。
# 予測と可視化
X_test = torch.linspace(-6, 6, 300).unsqueeze(1)
mu_test, var_test = model.predict(X_test)
mu_np, std_np = mu_test.numpy(), var_test.sqrt().numpy()
fig, axes = plt.subplots(1, 2, figsize=(11, 5))
axes[0].scatter(X_np[::20], y_np[::20], s=4, alpha=0.3, color='gray', label='Data')
axes[0].plot(X_test.numpy(), mu_np, 'b-', lw=2, label='SVGP mean')
axes[0].fill_between(X_test.numpy().ravel(), mu_np-2*std_np, mu_np+2*std_np,
alpha=0.3, color='blue', label='95% CI')
axes[0].scatter(model.Z.detach().numpy().ravel(), np.full(M, y_np.min()-0.5),
marker='|', s=80, color='red', label=f'Inducing (M={M})')
axes[0].set_title(f'SVGP (N={N}, M={M})'); axes[0].legend(fontsize=8); axes[0].grid(alpha=0.3)
axes[1].plot(losses, color='steelblue', alpha=0.7)
axes[1].set_yscale('log'); axes[1].set_title('Training loss (-ELBO)'); axes[1].grid(alpha=0.3)
plt.tight_layout(); plt.savefig('svgp_basic.png', dpi=140); plt.show()
このグラフから3点が読み取れます。左側のSVGP予測平均（青線）は元の非線形関数 $\sin x + 0.3x + 0.5\cos 2x$ をきれいに追えており、$M=64$ の誘導点だけで $N=10000$ のデータを十分要約できていることを示します。95%信頼区間（青の帯）はデータ密度が低い領域（$|x|>5$ あたり）で広がり、データ密度が高い領域で狭まる——GP本来の「データから離れるほど不確実」な性質が保たれています。誘導点（赤の縦線）はデータの定義域 $[-5, 5]$ に集中して配置され、k-means 初期化と勾配最適化で適切に位置調整された結果です。右側の学習曲線は対数スケールで滑らかに減少しており、確率勾配法でも安定して収束しています。
誘導点数Mのトレードオフ
次に、$M$ を変えたときの予測精度と計算時間のトレードオフを定量化します。
import time
def train_svgp(N, M, X, y, n_iter=1500, batch_size=256):
"""SVGPを学習しモデル・最終ELBO・所要時間を返す"""
kmeans = KMeans(n_clusters=M, n_init=5, random_state=0).fit(X.numpy())
Z = torch.from_numpy(kmeans.cluster_centers_.astype(np.float32))
m = SVGP(M=M, d=X.size(1), Z_init=Z)
opt = torch.optim.Adam(m.parameters(), lr=1e-2)
t0 = time.perf_counter()
for _ in range(n_iter):
idx = torch.randint(0, N, (batch_size,))
opt.zero_grad()
loss = -m.elbo(X[idx], y[idx], N_total=N); loss.backward(); opt.step()
return m, -loss.item(), time.perf_counter() - t0
Ms = [8, 16, 32, 64, 128, 256]
results = []
X_test = torch.linspace(-5, 5, 500).unsqueeze(1)
y_true = (np.sin(X_test.numpy()[:,0]) + 0.3*X_test.numpy()[:,0]
+ 0.5*np.cos(2*X_test.numpy()[:,0]))
for M in Ms:
m_, e_, t_ = train_svgp(N, M, X, y)
mu, _ = m_.predict(X_test)
mse = ((mu.numpy() - y_true) ** 2).mean()
results.append((M, e_, t_, mse))
print(f"M={M:3d} ELBO={e_:8.1f} time={t_:5.1f}s MSE={mse:.4f}")
Ms_arr = np.array([r[0] for r in results])
elbos = np.array([r[1] for r in results])
times = np.array([r[2] for r in results])
mses = np.array([r[3] for r in results])
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
for ax, ys, color, ylabel, title in [
(axes[0], elbos, 'navy', 'Final ELBO', 'Model quality vs M'),
(axes[1], mses, 'darkred', 'Test MSE', 'Prediction error vs M'),
(axes[2], times, 'darkgreen', 'Training time [s]', 'Computational cost vs M'),
]:
ax.plot(Ms_arr, ys, 'o-', color=color)
ax.set_xlabel('M (inducing points)'); ax.set_ylabel(ylabel); ax.set_title(title)
ax.set_xscale('log', base=2); ax.grid(alpha=0.3)
axes[1].set_yscale('log'); axes[2].set_yscale('log')
plt.tight_layout(); plt.savefig('svgp_M_tradeoff.png', dpi=140); plt.show()
3つのグラフからトレードオフが鮮明に読み取れます。左の ELBO は $M$ を増やすほど上昇しますが、$M=64$ あたりで頭打ちに近づきます。これは「64個の誘導点でこのデータの本質はほぼ表現できている」ことを意味し、それ以上は限界利益が小さい領域です。中央の MSE も同様の傾向で、$M=32$〜$64$ で急速に低下し、$M=128$ 以降は誤差軸が log なので僅かな改善しか見られません。右の学習時間は $M^2$ 程度に増加し、$M=256$ では $M=8$ の数十倍かかります。実用上は「ELBO が頭打ちになる手前」を選ぶのが定石で、このデータでは $M=64$ が最良コスパといえます。
多次元データへの拡張
最後に、入力が多次元の場合に SVGP がどう動くかを試します。$d=5$ 次元の非線形回帰問題を $N=20000$ 点で学習します。
# 5次元の非線形回帰問題 (N=20000, d=5)
torch.manual_seed(42)
N_md, d = 20000, 5
X_md = np.random.uniform(-2, 2, (N_md, d)).astype(np.float32)
w = np.array([1.0, -0.5, 0.8, 0.3, -0.6], dtype=np.float32)
y_md = (X_md @ w + np.sin(X_md[:,0])*np.cos(X_md[:,1]) + 0.3*X_md[:,2]**2
+ 0.2*np.random.randn(N_md).astype(np.float32))
X_md_t, y_md_t = torch.from_numpy(X_md), torch.from_numpy(y_md)
M_md = 200
Z_init = torch.from_numpy(
KMeans(n_clusters=M_md, n_init=5, random_state=42).fit(X_md)
.cluster_centers_.astype(np.float32))
model_md = SVGP(M=M_md, d=d, Z_init=Z_init)
opt_md = torch.optim.Adam(model_md.parameters(), lr=5e-3)
losses_md = []
for it in range(3000):
idx = torch.randint(0, N_md, (512,))
opt_md.zero_grad()
loss = -model_md.elbo(X_md_t[idx], y_md_t[idx], N_total=N_md)
loss.backward(); opt_md.step()
losses_md.append(loss.item())
if it % 300 == 0: print(f"iter {it:4d} -ELBO={loss.item():.2f}")
X_te = np.random.uniform(-2, 2, (2000, d)).astype(np.float32)
y_te = X_te @ w + np.sin(X_te[:,0])*np.cos(X_te[:,1]) + 0.3*X_te[:,2]**2
mu_te, _ = model_md.predict(torch.from_numpy(X_te))
mse_md = ((mu_te.numpy() - y_te)**2).mean(); r2 = 1 - mse_md/y_te.var()
print(f"5D: MSE={mse_md:.4f}, R^2={r2:.4f}")
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
axes[0].scatter(y_te, mu_te.numpy(), s=8, alpha=0.4)
mn, mx = y_te.min(), y_te.max()
axes[0].plot([mn, mx], [mn, mx], 'r--', label='y=x')
axes[0].set_xlabel('True y'); axes[0].set_ylabel('Predicted mu')
axes[0].set_title(f'5D SVGP: R^2={r2:.3f}'); axes[0].legend(); axes[0].grid(alpha=0.3)
axes[1].plot(losses_md, alpha=0.7, color='purple'); axes[1].set_yscale('log')
axes[1].set_xlabel('Iteration'); axes[1].set_ylabel('-ELBO'); axes[1].grid(alpha=0.3)
plt.tight_layout(); plt.savefig('svgp_5d.png', dpi=140); plt.show()
このグラフから2点が確認できます。左の散布図は真値 vs 予測で、点群が $y=x$ 線（赤破線）にきれいに沿っていることから、$R^2$ がほぼ1に近い良好な予測精度であることを示します。$N=20000$ という大規模データに対して、$M=200$ の誘導点のみで5次元の非線形関数を捉えられているわけです。右の学習曲線は対数スケールで4桁にわたって減少し、20000データでも安定して収束することを示しています。1次元の場合と同様、確率勾配法と SVGP の組み合わせが多次元でもそのまま機能する点が重要です。
GPyTorchによるSVGPベンチマーク
スクラッチ実装で原理を把握したところで、産業利用される実機ライブラリ GPyTorch で同じ問題を解き、性能を比較します。GPyTorch は GPU 対応・KeOps 統合・自然勾配などの高度な最適化を備えており、本番環境ではこちらが定番です。
import torch, gpytorch, time
class GPyTorchSVGP(gpytorch.models.ApproximateGP):
def __init__(self, Z):
vd = gpytorch.variational.CholeskyVariationalDistribution(Z.size(0))
vs = gpytorch.variational.VariationalStrategy(
self, Z, vd, learn_inducing_locations=True)
super().__init__(vs)
self.mean_module = gpytorch.means.ConstantMean()
self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())
def forward(self, x):
return gpytorch.distributions.MultivariateNormal(
self.mean_module(x), self.covar_module(x))
X = torch.from_numpy(X_np).contiguous(); y = torch.from_numpy(y_np).contiguous()
M_gp = 64
Z_gp = torch.from_numpy(
KMeans(n_clusters=M_gp, n_init=10, random_state=0).fit(X_np)
.cluster_centers_.astype(np.float32))
likelihood = gpytorch.likelihoods.GaussianLikelihood()
model_gp = GPyTorchSVGP(Z_gp); model_gp.train(); likelihood.train()
opt = torch.optim.Adam([{'params': model_gp.parameters()},
{'params': likelihood.parameters()}], lr=1e-2)
mll = gpytorch.mlls.VariationalELBO(likelihood, model_gp, num_data=X.size(0))
t0 = time.perf_counter()
for it in range(2000):
idx = torch.randint(0, X.size(0), (256,))
opt.zero_grad()
loss = -mll(model_gp(X[idx]), y[idx]); loss.backward(); opt.step()
if it % 200 == 0: print(f"iter {it:4d} -ELBO={loss.item():.3f}")
print(f"GPyTorch SVGP time: {time.perf_counter()-t0:.1f}s")
model_gp.eval(); likelihood.eval()
X_te = torch.linspace(-6, 6, 300).unsqueeze(1)
with torch.no_grad(), gpytorch.settings.fast_pred_var():
pred = likelihood(model_gp(X_te))
mu_gp, std_gp = pred.mean.numpy(), pred.stddev.numpy()
plt.figure(figsize=(10, 5))
plt.scatter(X_np[::20], y_np[::20], s=4, alpha=0.3, color='gray', label='Data')
plt.plot(X_te.numpy(), mu_gp, 'r-', lw=2, label='GPyTorch SVGP mean')
plt.fill_between(X_te.numpy().ravel(), mu_gp-2*std_gp, mu_gp+2*std_gp,
alpha=0.3, color='red')
Zl = model_gp.variational_strategy.inducing_points.detach().numpy().ravel()
plt.scatter(Zl, np.full(M_gp, y_np.min()-0.5), marker='|', s=80,
color='black', label=f'Inducing (M={M_gp})')
plt.title('GPyTorch SVGP on N=10000'); plt.legend(); plt.grid(alpha=0.3)
plt.tight_layout(); plt.savefig('gpytorch_svgp.png', dpi=140); plt.show()
このコードと出力から3点が確認できます。GPyTorch のクラス階層（ApproximateGP + VariationalStrategy + CholeskyVariationalDistribution）は SVGPの構成要素を綺麗に分離しており、自然勾配やマルチタスクへの拡張が VariationalStrategy の差し替えだけで済むモジュラ設計になっています。学習時間はスクラッチ実装と同程度（GPyTorch は内部最適化が効くため、大規模になるほど差が広がります）。予測曲線はスクラッチ実装とほぼ重なり、変分パラメータの最適化が両者で同じ最適解にたどり着いていることが分かります。本番環境では gpytorch.settings.fast_pred_var() のような最適化フラグや、KeOps バックエンドによるGPU高速化が威力を発揮します。
まとめ
本記事では、大規模データに対応するガウス過程の決定版近似である SVGP について、計算量問題の出発点から ELBO の導出、PyTorch スクラッチ実装、GPyTorch ベンチマーク、誘導点数のトレードオフまで一貫して解説しました。
- 計算量問題: 通常のGP回帰は $\bm{K}_{NN}$ のコレスキー分解で $O(N^3)$、メモリ $O(N^2)$。$N=10^4$ で実用上限、$N=10^5$ で破綻。
- 歴史の整理: SoR → DTC/FITC（Snelson-Ghahramani 2006）→ VFE（Titsias 2009）→ SVGP（Hensman et al. 2013）と発展。「事前分布をいじる」から「事後分布を変分近似」へ、そして「ミニバッチ化」へ。
- 誘導点の役割: $M$ 個の誘導点 $\bm{u} = f(\bm{Z})$ がデータの十分統計量として機能。補間行列 $\bm{K}_{NM}\bm{K}_{MM}^{-1}$ で観測点を再構成。
- ELBOの分解: $\mathcal{L} = \sum_i \mathbb{E}_{q(f_i)}[\log p(y_i|f_i)] – \mathrm{KL}[q(\bm{u})\|p(\bm{u})]$。データ項が $i$ ごとの和に分解できることがミニバッチ化を可能にする決定的性質。
- 変分パラメータ: $q(\bm{u}) = \mathcal{N}(\bm{m}, \bm{S})$ の $(\bm{m}, \bm{S})$、誘導点 $\bm{Z}$、カーネルハイパーをすべて勾配で同時最適化。
- 計算量: $O(NM^2)$、メモリ $O(M^2)$。$M=200$、$N=10^7$ クラスまで現実的に学習可能。
- 拡張: 多重出力GP（LMC/ICM）、Deep GPs、Conditional Neural Processes と自然につながる。
- 実装: PyTorch自前実装で原理を把握し、GPyTorch で本番運用。誘導点初期化は k-means、最適化は自然勾配 + Adam が定石。
GP の表現力と不確実性定量化を、現代の機械学習のスケールに乗せる技術——それが SVGP です。衛星リモートセンシング、ロボット動力学学習、気象予測、株価モデリングまで、応用範囲は今も急速に広がり続けています。
次のステップとして、以下の記事も参考にしてください。