未知の関数を観測データから推定し、しかも「どこが確からしくて、どこが当てにならないか」まで教えてほしい——こうした要求は、ベイズ最適化での次の実験点選び、高コストな数値シミュレーションの代理モデル、センサ較正の外挿など、いたるところで現れます。ガウス過程（GP）は、まさにこの「予測 + 不確実性」を同時に出せる強力な道具です。
ところが単層のガウス過程には、見過ごせない弱点があります。それは「関数の滑らかさがどこでも一様」だという仮定です。なめらかに変化する領域と、急激に立ち上がる段差が混在するような非定常な関数を、たった1つの長さスケールで表そうとすると、段差をなまらせるか、平坦部を暴れさせるか、どちらかになってしまう。深層学習が「層を重ねて表現力を上げる」ことで成功したのと同じ発想を、ガウス過程に持ち込めないか——この問いから生まれたのが 深層ガウス過程（Deep Gaussian Process, Deep GP） です。Damianou & Lawrence (2013) が提案したこのモデルは、ガウス過程を関数合成 $f = f_L \circ \cdots \circ f_1$ の形で多層に積み重ねます。
Deep GP が活きる場面は具体的です。たとえばベイズ最適化では、目的関数が非定常なとき単層GPの代理モデルが探索を誤誘導しますが、Deep GP なら局所的に滑らかさを変えて追従できます。また工学シミュレーションの代理モデルでは、有限要素解析やCFDのような高コスト計算を少数の評価点から学習し、不確実性つきで予測する用途で力を発揮します。
本記事の内容
- 単層GPの表現力の限界と、それを関数合成で乗り越える発想
- 各層をGPにする Deep GP の定式化と、なぜ周辺尤度が解析的に求まらないのか
- 周辺化不能性を回避する変分下界（ELBO）の導出と、層間の不確実性伝播
- doubly stochastic 変分推論（DSVI, Salimbeni & Deisenroth 2017）によるスケーラブルな学習
- 無限幅ニューラルネットワークとGPの対応、そして Deep GP との本質的な違い
- 過剰平滑化（over-smoothing）の問題と、平均関数による緩和
- GPyTorch を用いた Deep GP の実装と、回帰・ベイズ最適化・代理モデルへの応用
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ガウス過程回帰を初めから分かりやすく解説
- カーネル関数の種類と性質 — ガウス過程で使うカーネル
- 疎ガウス過程（Sparse Gaussian Processes）をわかりやすく解説
- Stochastic Variational Gaussian Process (SVGP) — 大規模データに対応するGP近似
単層ガウス過程の限界
なぜ単層では足りないのか
ガウス過程は、関数 $f(\bm{x})$ そのものに確率分布を置くモデルです。任意の有限個の入力点 $\bm{x}_1, \dots, \bm{x}_N$ における関数値 $\bm{f} = (f(\bm{x}_1), \dots, f(\bm{x}_N))^\top$ が、平均ベクトル $\bm{m}$ と共分散行列 $\bm{K}$ を持つ多変量正規分布に従う、と仮定します。
$$ \bm{f} \sim \mathcal{N}(\bm{m}, \bm{K}), \qquad K_{ij} = k(\bm{x}_i, \bm{x}_j) $$
ここでカーネル $k$ が「2点の関数値がどれだけ相関するか」を決めます。よく使うRBFカーネル $k(\bm{x}, \bm{x}’) = \sigma_f^2 \exp\!\big(-\|\bm{x}-\bm{x}’\|^2 / 2\ell^2\big)$ では、長さスケール $\ell$ が「どれくらい離れたら関数値が無相関になるか」を表します。問題は、この $\ell$ が入力空間全体で一定だという点です。つまり「関数の滑らかさはどこでも同じ」という定常性を暗黙に仮定しているのです。
現実の関数はそうとは限りません。たとえば材料の応力—ひずみ曲線は、弾性域では直線的でなめらかですが、降伏点を境に急に折れ曲がります。こうした関数を1つの $\ell$ で表そうとすると無理が生じます。
左図では長い長さスケール $\ell = 1.0$ を使ったため、$x=0$ にある段差がなめらかになまされてしまい、真の関数（赤）の鋭い立ち上がりを再現できていません。一方、右図は短い $\ell = 0.15$ を選んだので段差は捉えられるものの、本来平坦であるべき領域までギザギザに暴れ、信頼帯（薄い青）も広がっています。1つの長さスケールでは「段差を捉える」と「平坦部を落ち着かせる」を両立できない——これが単層GPの本質的な限界です。
深さで表現力を稼ぐという発想
ニューラルネットワークが成功した理由の一つは、単純な変換を何層も重ねることで、複雑な関数を「合成」によって表現した点にあります。同じ発想をGPに持ち込めないでしょうか。すなわち、入力をいったん別の空間へ写像し（第1層）、その写った先でさらにGP回帰する（第2層）。第1層が入力空間を伸び縮みさせて段差の手前で「座標を詰め」、段差を相対的にゆるやかにできれば、第2層は単純な定常GPでも全体として非定常な関数を表せるはずです。
この「GPを関数として合成する」アイデアこそが Deep GP です。次に、合成を数式できちんと定式化していきましょう。
ガウス過程の合成としての定式化
関数合成のモデル
Deep GP は、$L$ 個のガウス過程を入れ子に合成して作ります。入力 $\bm{X}$ から最終出力 $\bm{Y}$ までを、中間表現 $\bm{H}_1, \dots, \bm{H}_{L-1}$ を経由して写像します。
形式的には、各層の写像 $f_\ell$ が独立なガウス過程であるとして、次のように合成します。
$$ \begin{equation} \bm{Y} = f_L \circ f_{L-1} \circ \cdots \circ f_1(\bm{X}), \qquad f_\ell \sim \mathcal{GP}\big(m_\ell(\cdot),\, k_\ell(\cdot, \cdot)\big) \end{equation} $$
中間表現を明示すると、層ごとに次の確率的な遷移として書けます。
$$ \bm{H}_\ell = f_\ell(\bm{H}_{\ell-1}) + \bm{\epsilon}_\ell, \qquad \bm{H}_0 = \bm{X},\quad \bm{H}_L = \bm{Y} $$
ここで $\bm{\epsilon}_\ell$ は各層に入る微小なノイズです。重要なのは、$f_\ell$ が決定的な関数ではなくガウス過程からのサンプルだという点です。つまり中間層 $\bm{H}_\ell$ そのものが確率変数であり、潜在変数として扱われます。
各層が何をしているか
各層の役割を、入力空間がどう変形されるかという視点で見てみましょう。
左から、層1が入力 $x$ をなめらかに歪め（warp）、層2がその表現をさらに変形し、出力層が最終的な写像を作っています。各層単体は素直な滑らかな写像ですが、これらを合成すると、入力空間の一部を「引き伸ばし」別の部分を「圧縮」する効果が生まれます。たとえば段差の手前で座標を圧縮すれば、後段のGPから見れば段差はゆるやかな勾配に見え、定常カーネルでも扱えるようになる——これが Deep GP の表現力の源です。深さは、入力空間の局所的な伸縮（＝非定常性）を学習する自由度を与えるのです。
ここまでで「GPを合成すれば非定常性を扱える」という見通しが立ちました。しかし、ここに学習上の大きな壁が立ちはだかります。次にその壁——周辺尤度が解析的に計算できない問題——を見ていきます。
周辺化不能性という壁
なぜ尤度が解析的に求まらないのか
単層GP回帰の魅力の一つは、周辺尤度 $p(\bm{Y} \mid \bm{X})$ が解析的に求まることでした。観測 $\bm{Y} = \bm{f} + \bm{\epsilon}$ について、潜在関数値 $\bm{f}$ を積分消去すると、$\bm{Y}$ はやはりガウス分布になり、
$$ p(\bm{Y} \mid \bm{X}) = \mathcal{N}\big(\bm{Y} \mid \bm{m},\ \bm{K} + \sigma^2 \bm{I}\big) $$
と閉じた形で書けます。ガウス分布の線形変換と周辺化がまたガウス分布になる、という性質のおかげです。
ところが Deep GP では事情が一変します。出力の周辺尤度を求めるには、すべての中間層 $\bm{H}_1, \dots, \bm{H}_{L-1}$ を積分消去しなければなりません。
$$ p(\bm{Y} \mid \bm{X}) = \int p(\bm{Y} \mid \bm{H}_{L-1})\, p(\bm{H}_{L-1} \mid \bm{H}_{L-2}) \cdots p(\bm{H}_1 \mid \bm{X})\, d\bm{H}_1 \cdots d\bm{H}_{L-1} $$
ここで決定的に問題になるのが、各層の遷移確率 $p(\bm{H}_\ell \mid \bm{H}_{\ell-1})$ がガウス分布であっても、カーネル $k_\ell(\bm{H}_{\ell-1}, \bm{H}_{\ell-1})$ が前層の出力 $\bm{H}_{\ell-1}$ の非線形関数になっているという事実です。たとえばRBFカーネルなら、共分散行列の各要素に $\exp(-\|\bm{h}_i – \bm{h}_j\|^2 / 2\ell^2)$ が入り、$\bm{H}_{\ell-1}$ について非線形です。
つまり、積分の中で「ガウス分布の平均と共分散自体が、積分変数の非線形関数になっている」。これは線形ガウスの枠を完全に外れており、積分は解析的に実行できません。これを 周辺化不能性（intractability of marginalization） と呼びます。
モンテカルロでも素朴にはつらい
「解析的に無理なら数値積分すればよい」と思うかもしれません。しかし中間層 $\bm{H}_\ell$ はデータ点数 $N$ × 各層の次元 $D_\ell$ の高次元変数です。$N$ が数千、層が数層あれば、積分は数万次元になり、素朴なモンテカルロやグリッド積分では到底歯が立ちません。
そこで登場するのが、関数解析と最適化を橋渡しする変分推論です。「真の事後分布を直接求められないなら、扱いやすい分布族で近似し、その近似が真の分布にどれだけ近いかを最大化する」という戦略を取ります。次のセクションで、Deep GP 用の変分下界を導出します。
変分下界とベイズ的な層間伝播
誘導点による層の表現
変分推論に入る前に、各層のGPを誘導点（inducing points） で要約します。これは疎ガウス過程（疎ガウス過程）や SVGP（SVGP）で使われたのと同じ道具です。各層 $\ell$ に $M$ 個の誘導入力 $\bm{Z}_\ell$ と、そこでの誘導出力 $\bm{u}_\ell = f_\ell(\bm{Z}_\ell)$ を導入します。GPの事前分布のもとで $\bm{u}_\ell$ もガウス分布に従います。
$$ p(\bm{u}_\ell) = \mathcal{N}\big(\bm{u}_\ell \mid \bm{0},\ \bm{K}_{\bm{Z}_\ell \bm{Z}_\ell}\big) $$
誘導点を使う利点は2つあります。第一に、$N$ 個の関数値の代わりに $M \ll N$ 個の $\bm{u}_\ell$ で各層を要約でき、計算量が下がります。第二に、後で見るように層ごとの確率的な前向き計算を、$\bm{u}_\ell$ を条件づけたガウス分布として簡潔に書けます。
変分分布の設定
真の事後分布 $p(\{\bm{H}_\ell\}, \{\bm{u}_\ell\} \mid \bm{Y})$ は計算できません。そこで、扱いやすい変分分布 $q$ で近似します。各層の誘導出力に独立なガウス分布を置きます。
$$ q(\{\bm{u}_\ell\}) = \prod_{\ell=1}^{L} q(\bm{u}_\ell), \qquad q(\bm{u}_\ell) = \mathcal{N}\big(\bm{u}_\ell \mid \bm{m}_\ell,\ \bm{S}_\ell\big) $$
ここで $\bm{m}_\ell, \bm{S}_\ell$ が学習すべき変分パラメータです。中間層 $\bm{H}_\ell$ については、GPの条件付き分布をそのまま使い、$\bm{u}_\ell$ で条件づけた
$$ q(\bm{H}_\ell \mid \bm{H}_{\ell-1}, \bm{u}_\ell) = p(\bm{H}_\ell \mid \bm{H}_{\ell-1}, \bm{u}_\ell) $$
を採用します（変分分布の中で GP の事前条件付きをそのまま再利用する、という設計）。
ELBO の導出
変分推論の中核は、対数周辺尤度をイェンセンの不等式で下から押さえる証拠下界（Evidence Lower BOund, ELBO） です。対数尤度を変分分布の期待値として書き直し、イェンセンの不等式を適用します。
$$ \log p(\bm{Y}) = \log \int p(\bm{Y}, \{\bm{H}_\ell\}, \{\bm{u}_\ell\})\, d\bm{H}\, d\bm{u} $$
被積分関数に $q/q$ を掛けて期待値の形にし、$\log$ の凹性からイェンセンの不等式を使うと、
$$ \log p(\bm{Y}) \geq \mathbb{E}_{q}\!\left[ \log \frac{p(\bm{Y}, \{\bm{H}_\ell\}, \{\bm{u}_\ell\})}{q(\{\bm{H}_\ell\}, \{\bm{u}_\ell\})} \right] \equiv \mathcal{L}_{\mathrm{ELBO}} $$
が得られます。等号は $q$ が真の事後分布に一致するときに成り立ちます。ここで分子の同時分布を層ごとに分解し、$q$ の設定で中間層の条件付きが約分されることを使うと、ELBO は次の見通しのよい形にまとまります。
$$ \begin{equation} \mathcal{L}_{\mathrm{ELBO}} = \sum_{n=1}^{N} \mathbb{E}_{q(\bm{f}_n^L)}\!\big[\log p(y_n \mid \bm{f}_n^L)\big] – \sum_{\ell=1}^{L} \mathrm{KL}\big(q(\bm{u}_\ell)\,\|\,p(\bm{u}_\ell)\big) \end{equation} $$
この式の構造を絵で押さえておきましょう。
ELBO は2つの項に分かれます。第1項は再構成項（データ適合項） で、各データ点 $y_n$ を、最終層の関数値 $\bm{f}_n^L$ を通じてどれだけうまく説明できるかを測ります。第2項はKL正則化項で、各層の誘導点の変分事後 $q(\bm{u}_\ell)$ を事前 $p(\bm{u}_\ell)$ に近づけ、過剰適合を防ぎます。ELBO を最大化することは、対数周辺尤度を最大化することの近似であり、同時に変分分布を真の事後に近づけることに対応します。
不確実性が層を通って伝わる
ここで強調したいのは、再構成項の期待値 $\mathbb{E}_{q(\bm{f}_n^L)}[\cdot]$ が、すべての層を通った後の予測分布のもとで取られている点です。第1層の出力に不確実性があれば、それが第2層の入力分布になり、第2層の非線形写像で押し出され、最終層の予測分布に積み上がります。
左図は不確実性の「押し出し」を示します。層1の出力に分布 $q(h)$（紫）があり、それが層2の非線形関数 $f_2$ を通ると、出力側の分布（青）が非対称に変形されます。これは決定的なニューラルネットでは起こらない現象です。右図は、層を重ねるほど予測分散がどう変化するかを示しています。データが密集する中央では分散が小さく抑えられる一方、外挿域（灰色帯）では層を経るごとに不確実性が正直に積み上がります。「知らない領域では自信を持たない」という性質が、層を通しても保たれる——これが Deep GP の予測が単なる点予測より信頼できる理由です。
ELBO の形は得られましたが、再構成項の期待値はやはり層をまたぐ多重積分を含み、そのままでは計算できません。これをスケーラブルに評価する仕組みが、次に見る doubly stochastic 変分推論です。
Doubly Stochastic 変分推論（DSVI）
2つの確率性
Salimbeni & Deisenroth (2017) が提案した doubly stochastic 変分推論（DSVI） は、ELBO の評価を2種類のモンテカルロ近似で行います。「二重に確率的」という名前は、この2つの確率性に由来します。
第1の確率性はミニバッチによるデータのサブサンプリングです。ELBO の再構成項は $N$ 個のデータ点の和ですが、これをミニバッチ $\mathcal{B}$ で不偏推定します。
$$ \sum_{n=1}^{N} \mathbb{E}_{q}\big[\log p(y_n \mid \bm{f}_n^L)\big] \approx \frac{N}{|\mathcal{B}|} \sum_{n \in \mathcal{B}} \mathbb{E}_{q}\big[\log p(y_n \mid \bm{f}_n^L)\big] $$
これにより、$N$ が数十万を超えても1ステップの計算がミニバッチサイズに比例するだけで済みます。
第2の確率性は各層の前向きサンプリングです。最終層の予測分布 $q(\bm{f}_n^L)$ は、層をひとつずつサンプリングしながら前向きに伝播させることで近似します。各層で誘導点を条件づけた予測分布
$$ q(\bm{f}_n^\ell \mid \bm{f}_n^{\ell-1}) = \mathcal{N}\big(\bm{f}_n^\ell \mid \mu_\ell(\bm{f}_n^{\ell-1}),\ \Sigma_\ell(\bm{f}_n^{\ell-1})\big) $$
から1サンプルを引き、それを次の層の入力に渡す。これを最終層まで繰り返せば、$\bm{f}_n^L$ の1サンプルが得られます。期待値 $\mathbb{E}_q[\cdot]$ はこのサンプルで近似します。
再パラメータ化トリック
サンプリングを含む期待値を勾配で最適化するには、再パラメータ化トリック（reparameterization trick） が必要です。各層のガウスサンプルを
$$ \bm{f}_n^\ell = \mu_\ell(\bm{f}_n^{\ell-1}) + \bm{L}_\ell(\bm{f}_n^{\ell-1})\, \bm{\xi}, \qquad \bm{\xi} \sim \mathcal{N}(\bm{0}, \bm{I}) $$
と書きます。ここで $\bm{L}_\ell$ は $\Sigma_\ell$ のコレスキー因子です。乱数 $\bm{\xi}$ を分布の外に追い出すことで、サンプル $\bm{f}_n^\ell$ が変分パラメータの微分可能な関数になり、自動微分で勾配を通せるようになります。
全体のデータフロー
これらを組み合わせた DSVI の流れを整理します。
ミニバッチを抽出し（確率性①）、各層を再パラメータ化サンプリングで前向きに伝播させ（確率性②）、最終層で対数尤度を評価して不偏な ELBO 推定 $\hat{\mathcal{L}}$ を得ます。これを変分パラメータ $\phi$（誘導点 $\bm{Z}_\ell$、変分平均 $\bm{m}_\ell$、変分共分散 $\bm{S}_\ell$、カーネルのハイパーパラメータ）について自動微分し、確率的勾配上昇で更新します。
$$ \phi \leftarrow \phi + \eta\, \nabla_\phi \hat{\mathcal{L}} $$
この推定 $\hat{\mathcal{L}}$ は2つの確率性のもとで ELBO の不偏推定になっており、確率的勾配降下（SGD）の収束保証がそのまま適用できます。ミニバッチで $N$ への依存を断ち切り、層サンプリングで周辺化不能性を回避する——この2段構えが、Deep GP を大規模・多層でも学習可能にした鍵です。
DSVI によって Deep GP は実用的に学習できるようになりました。ここで一度立ち止まり、よく似た概念である「無限幅ニューラルネットワークとGPの対応」と Deep GP がどう違うのかを整理しておきましょう。
無限幅ニューラルネットワークとの関係
無限幅NNはGPになる
Deep GP を「GPのディープ版」と聞くと、「無限幅のディープニューラルネットがGPになる」という有名な結果（Neal 1996, Lee et al. 2018）を思い出すかもしれません。実際、隠れ層が1層で幅を無限大にしたニューラルネットは、出力が中心極限定理によってガウス過程に収束します。重みに独立な事前分布を置くと、出力は無数の独立な確率変数の和になり、その和がガウス分布に近づくためです。
左図のように、無限幅の1隠れ層NNは単層GPと等価になり、対応するカーネルは活性化関数から決まります。さらに深いNNを無限幅にすると、各層のカーネルが合成された NNGP（Neural Network Gaussian Process） が得られます。
Deep GP は NNGP とは別物
ここが肝心な点です。無限幅ディープNN（NNGP）の出力は、いくら層を重ねても周辺分布がガウス過程のままです。層の効果はカーネルの合成として吸収され、結局は1つの（複雑な）カーネルを持つ単層GPと等価になります。つまり「深さ」は新しいカーネルを作るだけで、本質的に単層GPの枠を出ません。
一方、Deep GP では各中間層 $\bm{H}_\ell$ が確率変数として残り、それを積分消去する前の段階で非線形性が効くため、出力の周辺分布は非ガウスになります。これこそが周辺化不能性の正体であると同時に、Deep GP が単層GPを超える表現力を持つ理由です。右図と表にまとめたとおり、両者は「深さの扱い方」が根本的に異なります。NNGP は深さをカーネルに吸収し、Deep GP は深さを潜在変数の階層として保持する。後者だからこそ、各層で完全にベイズ的な不確実性を扱えるのです。
この違いを押さえると、「Deep GP は単に表現力が高い」だけでなく、「深さに新たな自由度を持ち込んでいる」と理解できます。ただし、深ければ深いほどよいわけではありません。次に Deep GP に固有の落とし穴——過剰平滑化——を見ます。
過剰平滑化という落とし穴
深さが情報を消す
Deep GP を素朴に深くすると、思わぬ問題が起きます。平均関数をゼロに設定したGPを何層も合成すると、各層が入力の情報を少しずつ失い、深い層では出力が入力にほとんど依存しなくなってしまうのです。これを 過剰平滑化（over-smoothing） または表現の縮退と呼びます。
左図は、平均ゼロGPを素朴に合成したときの出力を、合成層数を変えて描いたものです。層数が増えるほど、出力は入力の細かな構造を失い、のっぺりした単調な形に縮退していきます。右図は、入力と出力の相関が深さとともにどう減衰するかを示します。平均ゼロの素朴な合成（赤）では相関がほぼ指数的に減衰し、深い層では入力情報が消えてしまいます。
平均関数で緩和する
Salimbeni & Deisenroth (2017) が提案した実用上の重要な工夫が、各層の平均関数を恒等写像に設定することです。中間層の次元が等しい場合、層の平均関数を $m_\ell(\bm{h}) = \bm{h}$ とすると、GPは「恒等写像からのずれ」をモデル化することになります。
$$ \bm{H}_\ell = \bm{H}_{\ell-1} + \big(\text{GPによる残差}\big) $$
これは深層学習における残差接続（ResNet）と同じ発想です。恒等平均があれば、たとえGPの残差がゼロでも入力情報がそのまま次層に流れるため、深くしても情報が消えません。右図の緑線が示すように、恒等平均関数を入れると相関が高い水準で保たれます。深さの恩恵を活かすには、平均関数の設計で情報経路を確保することが不可欠なのです。
過剰平滑化への対処を理解したところで、いよいよ Deep GP を実装し、その挙動を確かめましょう。
Pythonでの実装
GPyTorch による Deep GP の構築
ここでは、ガウス過程ライブラリ GPyTorch の Deep GP 機能を使って、2層の Deep GP を回帰問題に適用します。まずモデルを構築します。各層は誘導点つきの変分GP（SVGP相当）で、DSVI に基づいて学習されます。
import torch
import gpytorch
from gpytorch.models.deep_gps import DeepGPLayer, DeepGP
from gpytorch.means import ConstantMean, LinearMean
from gpytorch.kernels import ScaleKernel, RBFKernel
from gpytorch.variational import (
VariationalStrategy, CholeskyVariationalDistribution,
)
from gpytorch.distributions import MultivariateNormal
class ToyDeepGPLayer(DeepGPLayer):
"""1つのGP層。誘導点つき変分GPとして実装する。"""
def __init__(self, input_dims, output_dims, num_inducing=64, linear_mean=True):
# 誘導点を出力次元ごとに用意（バッチ次元 = 出力次元）
if output_dims is None:
inducing_points = torch.randn(num_inducing, input_dims)
batch_shape = torch.Size([])
else:
inducing_points = torch.randn(output_dims, num_inducing, input_dims)
batch_shape = torch.Size([output_dims])
variational_dist = CholeskyVariationalDistribution(
num_inducing, batch_shape=batch_shape)
variational_strategy = VariationalStrategy(
self, inducing_points, variational_dist, learn_inducing_locations=True)
super().__init__(variational_strategy, input_dims, output_dims)
# 恒等(線形)平均で過剰平滑化を緩和。出力層は定数平均。
self.mean_module = LinearMean(input_dims) if linear_mean else ConstantMean(batch_shape=batch_shape)
self.covar_module = ScaleKernel(
RBFKernel(batch_shape=batch_shape, ard_num_dims=input_dims),
batch_shape=batch_shape)
def forward(self, x):
mean_x = self.mean_module(x)
covar_x = self.covar_module(x)
return MultivariateNormal(mean_x, covar_x)
ここでのポイントは2つです。第一に、各層を CholeskyVariationalDistribution と VariationalStrategy で構成し、誘導点 $\bm{Z}_\ell$ の位置も学習対象（learn_inducing_locations=True）にしています。第二に、中間層には LinearMean（線形・恒等的な平均）を使い、前節で述べた過剰平滑化の緩和を実装している点です。出力層だけは定数平均にしています。
続いて、この層を2つ積み重ねて Deep GP 本体を組み立てます。
class TwoLayerDeepGP(DeepGP):
"""入力 -> 隠れ層(2次元) -> 出力(1次元) の2層 Deep GP。"""
def __init__(self, input_dim, hidden_dim=2):
hidden_layer = ToyDeepGPLayer(
input_dims=input_dim, output_dims=hidden_dim, linear_mean=True)
last_layer = ToyDeepGPLayer(
input_dims=hidden_dim, output_dims=None, linear_mean=False)
super().__init__()
self.hidden_layer = hidden_layer
self.last_layer = last_layer
self.likelihood = gpytorch.likelihoods.GaussianLikelihood()
def forward(self, x):
h = self.hidden_layer(x) # 第1層: 確率的な中間表現
y = self.last_layer(h) # 第2層: 中間表現から出力へ
return y
forward の中で self.hidden_layer(x) が確率分布を返し、それを self.last_layer(h) に渡しています。GPyTorch は内部でこの呼び出しを再パラメータ化サンプリングとして処理し、DSVI の前向き伝播を自動で行います。つまり前節で導出したサンプリングベースの層間伝播が、このコードの裏で動いています。
学習ループ（DSVIの ELBO 最大化）
学習では、DSVI の不偏 ELBO 推定をミニバッチで最大化します。GPyTorch では DeepApproximateMLL が層ごとの KL 項を含む ELBO を提供します。
import torch
from torch.utils.data import TensorDataset, DataLoader
from gpytorch.mlls import DeepApproximateMLL, VariationalELBO
def train_deep_gp(model, train_x, train_y, num_epochs=200, batch_size=128,
num_samples=8, lr=0.01):
model.train()
dataset = TensorDataset(train_x, train_y)
loader = DataLoader(dataset, batch_size=batch_size, shuffle=True)
optimizer = torch.optim.Adam(model.parameters(), lr=lr)
# ELBO: 再構成項 - 層ごとのKL（DSVIの目的関数）
mll = DeepApproximateMLL(
VariationalELBO(model.likelihood, model, num_data=train_x.size(0)))
for epoch in range(num_epochs):
epoch_loss = 0.0
for x_batch, y_batch in loader:
optimizer.zero_grad()
# num_samples 個の前向きサンプルで期待値をモンテカルロ近似
with gpytorch.settings.num_likelihood_samples(num_samples):
output = model(x_batch)
loss = -mll(output, y_batch) # ELBOの符号反転を最小化
loss.backward()
optimizer.step()
epoch_loss += loss.item()
if (epoch + 1) % 50 == 0:
print(f"epoch {epoch+1:3d} loss(-ELBO) = {epoch_loss/len(loader):.4f}")
return model
gpytorch.settings.num_likelihood_samples(num_samples) が、DSVI の第2の確率性（層サンプリング数）を指定しています。DataLoader のミニバッチが第1の確率性に対応します。損失 -mll(...) は ELBO の符号反転で、これを最小化することが ELBO の最大化に一致します。エポックごとに損失が減れば、ELBO が上がり、変分分布が真の事後に近づいている証拠です。
非定常関数での挙動確認
冒頭で見た「段差 + 滑らかな振動」という非定常関数に、単層GPと2層 Deep GP を当てはめて比較します。
import torch
import numpy as np
torch.manual_seed(0)
np.random.seed(0)
# 非定常な真の関数（段差 + 滑らかな振動）
def f_true(x):
return np.where(x < 0, 0.0, 1.0) + 0.15 * np.sin(3 * x)
X = np.linspace(-3, 3, 200).reshape(-1, 1)
y = f_true(X).ravel() + 0.03 * np.random.randn(200)
train_x = torch.tensor(X, dtype=torch.float32)
train_y = torch.tensor(y, dtype=torch.float32)
# 2層 Deep GP を学習
model = TwoLayerDeepGP(input_dim=1, hidden_dim=2)
model = train_deep_gp(model, train_x, train_y, num_epochs=200)
# 予測（サンプル平均と分散を取る）
model.eval()
test_x = torch.linspace(-3.2, 3.2, 300).reshape(-1, 1)
with torch.no_grad(), gpytorch.settings.num_likelihood_samples(50):
pred = model.likelihood(model(test_x))
mean = pred.mean.mean(0) # サンプル方向に平均
var = pred.variance.mean(0) # 予測分散
lower = mean - 2 * var.sqrt()
upper = mean + 2 * var.sqrt()
print("予測平均 shape:", mean.shape, " 段差付近の分散:", var[150].item())
このコードでは、num_likelihood_samples(50) で50回の前向きサンプルを引き、その平均と分散で予測分布を要約しています。Deep GP は段差付近（$x=0$）で予測平均が急峻に立ち上がり、平坦部では落ち着く——単層GPでは両立できなかった挙動を、層合成によって実現します。段差付近では中間層が座標を局所的に圧縮するため、後段GPから見れば段差がゆるやかな勾配に見え、定常カーネルでも鋭い遷移を表せるのです。
2層 Deep GP の事前サンプル
学習結果の理解を深めるために、Deep GP の事前分布からのサンプルパスを単層GPと比べてみましょう。これは、Deep GP がそもそもどんな関数の集合を「もっともらしい」と考えているかを可視化するものです。
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(5)
x = np.linspace(-3, 3, 250).reshape(-1, 1)
def rbf(Xa, Xb, ell=1.0, sf2=1.0):
d2 = (Xa - Xb.T) ** 2
return sf2 * np.exp(-0.5 * d2 / ell ** 2)
def sample_gp(X, ell, sf2, rng):
K = rbf(X, X, ell, sf2) + 1e-8 * np.eye(len(X))
L = np.linalg.cholesky(K)
return (L @ rng.standard_normal(len(X))).ravel()
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
# 単層GP（定常）
for _ in range(5):
axes[0].plot(x.ravel(), sample_gp(x, 0.6, 1.0, rng), lw=1.8, alpha=0.85)
axes[0].set_title("単層GPのサンプル（定常）")
# 2層 Deep GP（warp -> GP, 非定常）
for _ in range(5):
h = sample_gp(x, 0.8, 1.2, rng).reshape(-1, 1)
order = np.argsort(h.ravel())
s2 = sample_gp(h[order], 0.5, 1.0, rng)
y = np.empty_like(s2); y[order] = s2
axes[1].plot(x.ravel(), y, lw=1.8, alpha=0.85)
axes[1].set_title("2層 Deep GP のサンプル（非定常）")
for ax in axes:
ax.set_xlabel("$x$"); ax.set_ylabel("$f(x)$"); ax.grid(alpha=0.3)
plt.tight_layout(); plt.show()
左の単層GPのサンプルは、どの場所でも同じ滑らかさを持つ定常な曲線です。一方、右の2層 Deep GP のサンプルは、場所によって急に変化したりゆるやかになったりと、滑らかさが空間的に変わる非定常な曲線になっています。これが、Deep GP が単層GPより豊かな関数の集合を事前分布として持っていることの直接的な証拠です。学習では、この豊かな事前から観測データに合うものを変分推論で絞り込んでいるわけです。
実装で Deep GP の挙動を確認できたので、最後に実際の応用場面を整理しましょう。
応用
不確実性つき回帰・ベイズ最適化・代理モデル
Deep GP は、単なる予測精度だけでなく「不確実性の質」が問われる場面で真価を発揮します。代表的な3つの応用を見てみましょう。
左の (1) 不確実性つき回帰 では、非定常な関数を予測しつつ、データのない領域で信頼帯が広がり「ここは当てにならない」と正直に表現できています。これは安全性が重要な制御や医療の予測で重宝されます。
中央の (2) ベイズ最適化の代理モデル では、目的関数を少数の評価点から学習し、獲得関数（ここでは UCB）で次に試すべき点を選びます。Deep GP は目的関数が非定常なときに単層GPより正確な代理モデルとなり、無駄な評価を減らせます（ベイズ最適化の理論 も参照）。
右の (3) シミュレーション代理モデル では、有限要素解析やCFDのような高コストな数値計算を、少数の評価点から学習した Deep GP で代替します。非定常な応答（共振ピークなど）を持つシミュレーションでも、層合成によって精度よく近似でき、設計探索や感度解析を高速化できます。
これら3つに共通するのは、「予測の確からしさ」をモデルが自分で見積もれることです。決定的なニューラルネットは点予測しか出しませんが、Deep GP は各層の不確実性を伝播させた予測分布を返します。だからこそ「次にどこを調べるべきか」「どの予測を信じてよいか」という意思決定に直結するのです。
まとめ
本記事では、ガウス過程を関数合成で多層に積み重ねた深層ガウス過程（Deep GP）について、Damianou & Lawrence (2013) の定式化を軸に、理論の導出から実装・応用まで解説しました。
- 単層GPの限界: 単一の長さスケールは関数の滑らかさを一様と仮定するため、段差と平坦部が混在する非定常関数を表せない。
- 関数合成による定式化: Deep GP は $\bm{Y} = f_L \circ \cdots \circ f_1(\bm{X})$ と各層をGPにする。各層が入力空間を局所的に warp し、合成で非定常性を獲得する。
- 周辺化不能性: 各層のカーネルが前層出力の非線形関数になるため、中間層を積分消去できず周辺尤度が解析的に求まらない。これが Deep GP 学習の本質的な壁。
- 変分下界: ELBO は「再構成項 $-$ 層ごとの KL 項」に分解され、誘導点で各層を要約することで計算可能になる。期待値は全層を通った予測分布のもとで取られ、不確実性が層間を伝播する。
- DSVI: ミニバッチ（確率性①）と層の再パラメータ化サンプリング（確率性②）の二重のモンテカルロ近似で、ELBO の不偏勾配をスケーラブルに得る (Salimbeni & Deisenroth 2017)。
- NNとの関係: 無限幅NNはGP（NNGP）になるが周辺はガウスのまま。Deep GP は中間層を確率変数として保持するため周辺が非ガウスになり、深さに新たな自由度を持ち込む。
- 過剰平滑化: 素朴に深くすると入力情報が縮退する。各層の平均関数を恒等写像にする（残差接続的な発想）ことで緩和する。
Deep GP は、深層学習の表現力とベイズの不確実性定量化を一つのモデルに統合した、ベイズ深層学習の中核的なアプローチです。不確実性が意思決定に直結するベイズ最適化、シミュレーション代理モデル、安全性が問われる回帰問題で特に有効です。
次のステップとして、以下の記事も参考にしてください。