18.3. ガウス過程推論

この節では、前節で導入した GP 事前分布を用いて、事後推論を行い予測を作成する方法を示す。まずは回帰から始める。回帰では、推論を 閉形式 で行うことができる。これは、実際にガウス過程をすぐ使い始めるための「GP の要点」セクションである。まず、基本的な操作をすべてゼロから実装し、その後で GPyTorch を導入する。これにより、最先端のガウス過程を扱ったり、深層ニューラルネットワークと統合したりすることがずっと便利になる。これらのより高度な話題については、次の節で詳しく扱う。その節では、近似推論が必要となる設定――分類、点過程、あるいは非ガウス尤度全般――についても考える。

18.3.1. 回帰における事後推論

観測 モデルは、学習したい関数 \(f(x)\) と、観測値 \(y(x)\) を結びつける。どちらもある入力 \(x\) によって添字付けられる。分類では、\(x\) は画像の画素であり、\(y\) は対応するクラスラベルである。回帰では、\(y\) は通常、地表温度、海面水位、\(CO_2\) 濃度などの連続値出力を表す。

回帰では、出力は潜在的なノイズのない関数 \(f(x)\) に、独立同分布なガウス雑音 \(\epsilon(x)\) が加わったものだと仮定することがよくある。

(18.3.1)\[y(x) = f(x) + \epsilon(x),\]

ここで \(\epsilon(x) \sim \mathcal{N}(0,\sigma^2)\) である。\(\mathbf{y} = y(X) = (y(x_1),\dots,y(x_n))^{\top}\) を訓練観測のベクトル、\(\textbf{f} = (f(x_1),\dots,f(x_n))^{\top}\) を、訓練入力 \(X = {x_1, \dots, x_n}\) で問い合わせた潜在的なノイズのない関数値のベクトルとする。

ここでは \(f(x) \sim \mathcal{GP}(m,k)\) を仮定する。これは、任意の関数値の集まり \(\textbf{f}\) が、平均ベクトル \(\mu_i = m(x_i)\) と共分散行列 \(K_{ij} = k(x_i,x_j)\) をもつ同時多変量ガウス分布に従うことを意味する。RBF カーネル \(k(x_i,x_j) = a^2 \exp\left(-\frac{1}{2\ell^2}||x_i-x_j||^2\right)\) は、標準的な共分散関数の選択肢である。記法を簡単にするため、平均関数は \(m(x)=0\) と仮定する。導出は後で容易に一般化できる。

入力の集合

(18.3.2)\[X_* = x_{*1},x_{*2},\dots,x_{*m}.\]

で予測したいとする。すると、\(x^2\) と \(p(\mathbf{f}_* | \mathbf{y}, X)\) を求めたいことになる。回帰設定では、\(\mathbf{f}_* = f(X_*)\) と \(\mathbf{y}\) の同時分布を求めた後、ガウスの恒等式を使ってこの分布を便利に求めることができる。

式 (18.3.1) を訓練入力 \(X\) で評価すると、\(\mathbf{y} = \mathbf{f} + \mathbf{\epsilon}\) である。ガウス過程の定義（前節参照）より、\(\mathbf{f} \sim \mathcal{N}(0,K(X,X))\) である。ここで \(K(X,X)\) は、あり得るすべての入力対 \(x_i, x_j \in X\) に対して共分散関数（別名 カーネル）を評価して得られる \(n \times n\) 行列である。\(\mathbf{\epsilon}\) は単に \(\mathcal{N}(0,\sigma^2)\) からの iid サンプルからなるベクトルなので、分布は \(\mathcal{N}(0,\sigma^2I)\) である。したがって \(\mathbf{y}\) は、2つの独立な多変量ガウス変数の和であり、分布は \(\mathcal{N}(0, K(X,X) + \sigma^2I)\) となる。また、\(\textrm{cov}(\mathbf{f}_*, \mathbf{y}) = \textrm{cov}(\mathbf{y},\mathbf{f}_*)^{\top} = K(X_*,X)\) であることも示せる。ここで \(K(X_*,X)\) は、テスト入力と訓練入力のすべての組に対してカーネルを評価して得られる \(m \times n\) 行列である。

(18.3.3)\[\begin{split}\begin{bmatrix} \mathbf{y} \\ \mathbf{f}_* \end{bmatrix} \sim \mathcal{N}\left(0, \mathbf{A} = \begin{bmatrix} K(X,X)+\sigma^2I & K(X,X_*) \\ K(X_*,X) & K(X_*,X_*) \end{bmatrix} \right)\end{split}\]

その後、標準的なガウスの恒等式を使って、同時分布から条件付き分布を求めることができる（例えば Bishop 第2章を参照）。

\(\mathbf{f}_* | \mathbf{y}, X, X_* \sim \mathcal{N}(m_*,S_*)\) であり、\(m_* = K(X_*,X)[K(X,X)+\sigma^2I]^{-1}\textbf{y}\)、\(S = K(X_*,X_*) - K(X_*,X)[K(X,X)+\sigma^2I]^{-1}K(X,X_*)\) である。

通常、予測共分散行列 \(S\) 全体を使う必要はなく、各予測の不確実性として \(S\) の対角成分だけを使う。そのため、しばしばテスト点の集合ではなく、単一のテスト点 \(x_*\) に対する予測分布を書くことが多いである。

カーネル行列には、上で述べた RBF カーネルの振幅 \(a\) や長さ尺度 \(\ell\) のように、推定したいパラメータ \(\theta\) もある。これらの目的のために、周辺尤度 \(p(\textbf{y} | \theta, X)\) を使う。これは、\(\mathbf{y},\mathbf{f}_*\) の同時分布を求めるために周辺分布を導出する過程ですでに得たものである。後で見るように、周辺尤度はモデル適合度とモデル複雑性の項に分解され、ハイパーパラメータ学習におけるオッカムの剃刀の概念を自動的に組み込みる。詳しくは MacKay Ch. 28 (MacKay, 2003) および Rasmussen and Williams Ch. 5 (Rasmussen and Williams, 2006) を参照しよ。

from d2l import torch as d2l
import numpy as np
from scipy.spatial import distance_matrix
from scipy import optimize
import matplotlib.pyplot as plt
import math
import torch
import gpytorch
import os

d2l.set_figsize()

18.3.2. GP 回帰における予測とカーネルハイパーパラメータ学習のための式

ここでは、ガウス過程回帰でハイパーパラメータを学習し予測を行う際に使う式をまとめる。繰り返しになるが、入力 \(X = \{x_1,\dots,x_n\}\) によって添字付けられた回帰目標 \(\textbf{y}\) のベクトルがあり、テスト入力 \(x_*\) で予測したいとする。分散 \(\sigma^2\) をもつ独立同分布な加法ゼロ平均ガウス雑音を仮定する。潜在的なノイズのない関数には、平均関数 \(m\) とカーネル関数 \(k\) をもつガウス過程事前分布 \(f(x) \sim \mathcal{GP}(m,k)\) を使う。カーネル自体には、学習したいパラメータ \(\theta\) がある。例えば RBF カーネル \(k(x_i,x_j) = a^2\exp\left(-\frac{1}{2\ell^2}||x-x'||^2\right)\) を使うなら、\(\theta = \{a^2, \ell^2\}\) を学習したいことになる。\(K(X,X)\) は、\(n\) 個の訓練入力のあり得るすべての組に対してカーネルを評価して得られる \(n \times n\) 行列を表す。\(K(x_*,X)\) は、\(k(x_*, x_i)\) を \(i=1,\dots,n\) について評価して得られる \(1 \times n\) ベクトルを表す。\(\mu\) は、すべての訓練点 \(x\) で平均関数 \(m(x)\) を評価して得られる平均ベクトルである。

通常、ガウス過程を扱うときは、2段階の手順に従う。 1. 周辺尤度をこれらのハイパーパラメータに関して最大化することで、カーネルのハイパーパラメータ \(\hat{\theta}\) を学習する。 2. 予測平均を点予測器として使い、予測標準偏差の 2 倍を用いて 95% の信用集合を作る。ここでは、学習済みハイパーパラメータ \(\hat{\theta}\) に条件付ける。

対数周辺尤度は単なる対数ガウス密度であり、次の形をしている:

(18.3.4)\[\log p(\textbf{y} | \theta, X) = -\frac{1}{2}\textbf{y}^{\top}[K_{\theta}(X,X) + \sigma^2I]^{-1}\textbf{y} - \frac{1}{2}\log|K_{\theta}(X,X)| + c\]

予測分布は次の形をとります:

(18.3.5)\[p(y_* | x_*, \textbf{y}, \theta) = \mathcal{N}(a_*,v_*)\]

(18.3.6)\[a_* = k_{\theta}(x_*,X)[K_{\theta}(X,X)+\sigma^2I]^{-1}(\textbf{y}-\mu) + \mu\]

(18.3.7)\[v_* = k_{\theta}(x_*,x_*) - K_{\theta}(x_*,X)[K_{\theta}(X,X)+\sigma^2I]^{-1}k_{\theta}(X,x_*)\]

18.3.3. 学習と予測の式の解釈

ガウス過程の予測分布について、いくつか重要な点がある。

• モデルクラスは柔軟であるが、GP 回帰では 閉形式 で 厳密な ベイズ推論を行うことができる。カーネルのハイパーパラメータを学習することを除けば、訓練 はない。予測に使う式を正確に書き下せる。この点でガウス過程はかなり例外的であり、その便利さ、多用途性、そして今なお高い人気に大きく貢献している。

• 予測平均 \(a_*\) は、訓練目標 \(\textbf{y}\) の線形結合であり、重みはカーネル \(k_{\theta}(x_*,X)[K_{\theta}(X,X)+\sigma^2I]^{-1}\) である。後で見るように、カーネル（とそのハイパーパラメータ）はモデルの汎化特性において極めて重要な役割を果たする。

• 予測平均は目標値 \(\textbf{y}\) に明示的に依存するが、予測分散は依存しない。代わりに、予測不確実性は、カーネル関数によって決まるように、テスト入力 \(x_*\) が目標位置 \(X\) から離れるにつれて増大する。ただし、不確実性は、データから学習されるカーネルハイパーパラメータ \(\theta\) を通じて、間接的に目標値 \(\textbf{y}\) に依存する。

• 周辺尤度は、モデル適合度とモデル複雑性（行列式の対数）に分解される。周辺尤度は、データと整合的でありつつ最も単純な適合を与えるハイパーパラメータを選びがちである。

• 主な計算ボトルネックは、\(n\) 個の訓練点に対する \(n \times n\) の対称正定値行列 \(K(X,X)\) について、線形方程式を解くことと対数行列式を計算することである。素朴に行うと、これらの操作はそれぞれ \(\mathcal{O}(n^3)\) の計算量と、カーネル（共分散）行列の各要素に対して \(\mathcal{O}(n^2)\) のメモリを要し、しばしばコレスキー分解から始める。歴史的には、これらのボトルネックのために GP はおよそ 10,000 点未満の問題に限られてきており、GP は「遅い」という評判を長年持っていたが、これは今ではほぼ 10 年近く不正確である。高度な話題では、GP を何百万点もの問題にスケールさせる方法を議論する。

• よく使われるカーネル関数では、\(K(X,X)\) はしばしば特異に近く、コレスキー分解や線形方程式を解くための他の操作で数値的問題を引き起こすことがある。幸い、回帰ではしばしば \(K_{\theta}(X,X)+\sigma^2I\) を扱うため、ノイズ分散 \(\sigma^2\) が \(K(X,X)\) の対角に加わり、条件数が大幅に改善される。ノイズ分散が小さい場合、あるいはノイズなし回帰を行う場合は、条件数を改善するために対角へ \(10^{-6}\) 程度の小さな “jitter” を加えるのが一般的である。

18.3.4. ゼロからの実例

回帰データを作成し、そのデータを GP で当てはめる。すべての手順をゼロから実装する。 次の式からデータをサンプルする。

(18.3.8)\[y(x) = \sin(x) + \frac{1}{2}\sin(4x) + \epsilon,\]

ただし \(\epsilon \sim \mathcal{N}(0,\sigma^2)\) である。求めたいノイズのない関数は \(f(x) = \sin(x) + \frac{1}{2}\sin(4x)\) である。まずはノイズの標準偏差を \(\sigma = 0.25\) とする。

def data_maker1(x, sig):
    return np.sin(x) + 0.5 * np.sin(4 * x) + np.random.randn(x.shape[0]) * sig

sig = 0.25
train_x, test_x = np.linspace(0, 5, 50), np.linspace(0, 5, 500)
train_y, test_y = data_maker1(train_x, sig=sig), data_maker1(test_x, sig=0.)

d2l.plt.scatter(train_x, train_y)
d2l.plt.plot(test_x, test_y)
d2l.plt.xlabel("x", fontsize=20)
d2l.plt.ylabel("Observations y", fontsize=20)
d2l.plt.show()

ここでは、ノイズのある観測が丸印で、求めたい青色のノイズのない関数が見える。

では、潜在的なノイズのない関数 \(f(x)\sim \mathcal{GP}(m,k)\) に対して GP 事前分布を指定しよう。平均関数 \(m(x) = 0\) と、RBF 共分散関数（カーネル）

(18.3.9)\[k(x_i,x_j) = a^2\exp\left(-\frac{1}{2\ell^2}||x-x'||^2\right).\]

を使う。

mean = np.zeros(test_x.shape[0])
cov = d2l.rbfkernel(test_x, test_x, ls=0.2)

長さ尺度は 0.2 から始めている。データを当てはめる前に、妥当な事前分布を指定できているかを考えることが重要である。この事前分布からのサンプル関数と、95% の信用集合（真の関数がこの領域内にある確率が 95% だと考える）を可視化してみよう。

prior_samples = np.random.multivariate_normal(mean=mean, cov=cov, size=5)
d2l.plt.plot(test_x, prior_samples.T, color='black', alpha=0.5)
d2l.plt.plot(test_x, mean, linewidth=2.)
d2l.plt.fill_between(test_x, mean - 2 * np.diag(cov), mean + 2 * np.diag(cov),
                 alpha=0.25)
d2l.plt.show()

これらのサンプルは妥当に見えるだろうか。関数の高レベルな性質は、モデル化したいデータの種類と整合しているだろうか。

では、任意のテスト点 \(x_*\) における事後予測分布の平均と分散を求める。

(18.3.10)\[\bar{f}_{*} = K(x, x_*)^T (K(x, x) + \sigma^2 I)^{-1}y\]

(18.3.11)\[V(f_{*}) = K(x_*, x_*) - K(x, x_*)^T (K(x, x) + \sigma^2 I)^{-1}K(x, x_*)\]

予測を行う前に、カーネルのハイパーパラメータ \(\theta\) とノイズ分散 \(\sigma^2\) を学習する必要がある。事前分布の関数が、当てはめるデータに比べて変動が速すぎるように見えたので、長さ尺度の初期値を 0.75 にする。また、ノイズの標準偏差 \(\sigma\) も 0.75 と仮定する。

これらのパラメータを学習するために、これらのパラメータに関して周辺尤度を最大化する。

(18.3.12)\[\log p(y | X) = \log \int p(y | f, X)p(f | X)df\]

(18.3.13)\[\log p(y | X) = -\frac{1}{2}y^T(K(x, x) + \sigma^2 I)^{-1}y - \frac{1}{2}\log |K(x, x) + \sigma^2 I| - \frac{n}{2}\log 2\pi\]

おそらく事前分布の関数は変動が速すぎた。長さ尺度を 0.4 と仮定してみよう。ノイズの標準偏差も 0.75 と仮定する。これらは単なるハイパーパラメータの初期値であり、周辺尤度からこれらのパラメータを学習する。

ell_est = 0.4
post_sig_est = 0.5

def neg_MLL(pars):
    K = d2l.rbfkernel(train_x, train_x, ls=pars[0])
    kernel_term = -0.5 * train_y @ \
        np.linalg.inv(K + pars[1] ** 2 * np.eye(train_x.shape[0])) @ train_y
    logdet = -0.5 * np.log(np.linalg.det(K + pars[1] ** 2 * \
                                         np.eye(train_x.shape[0])))
    const = -train_x.shape[0] / 2. * np.log(2 * np.pi)

    return -(kernel_term + logdet + const)


learned_hypers = optimize.minimize(neg_MLL, x0=np.array([ell_est,post_sig_est]),
                                   bounds=((0.01, 10.), (0.01, 10.)))
ell = learned_hypers.x[0]
post_sig_est = learned_hypers.x[1]

この例では、長さ尺度 0.299 とノイズの標準偏差 0.24 を学習した。学習されたノイズが真のノイズに非常に近いことに注意しよ。これは、この問題に対して GP が非常によく適合していることを示している。

一般に、カーネルの選択とハイパーパラメータの初期化には慎重な検討が不可欠である。周辺尤度の最適化は初期値に対して比較的頑健であるが、悪い初期化の影響を受けないわけではない。上のスクリプトをさまざまな初期値で実行し、どのような結果になるか試してみよ。

では、これらの学習済みハイパーパラメータで予測を行おう。

K_x_xstar = d2l.rbfkernel(train_x, test_x, ls=ell)
K_x_x = d2l.rbfkernel(train_x, train_x, ls=ell)
K_xstar_xstar = d2l.rbfkernel(test_x, test_x, ls=ell)

post_mean = K_x_xstar.T @ np.linalg.inv((K_x_x + \
                post_sig_est ** 2 * np.eye(train_x.shape[0]))) @ train_y
post_cov = K_xstar_xstar - K_x_xstar.T @ np.linalg.inv((K_x_x + \
                post_sig_est ** 2 * np.eye(train_x.shape[0]))) @ K_x_xstar

lw_bd = post_mean - 2 * np.sqrt(np.diag(post_cov))
up_bd = post_mean + 2 * np.sqrt(np.diag(post_cov))

d2l.plt.scatter(train_x, train_y)
d2l.plt.plot(test_x, test_y, linewidth=2.)
d2l.plt.plot(test_x, post_mean, linewidth=2.)
d2l.plt.fill_between(test_x, lw_bd, up_bd, alpha=0.25)
d2l.plt.legend(['Observed Data', 'True Function', 'Predictive Mean', '95% Set on True Func'])
d2l.plt.show()

オレンジ色の事後平均は、真のノイズのない関数とほぼ完全に一致していることがわかる。ここで示している 95% の信用集合は、データ点ではなく、潜在的な ノイズのない（真の）関数に対するものである。この信用集合は真の関数を完全に含んでおり、広すぎも狭すぎもしないように見える。データ点を含むことは期待していないし、含む必要もない。観測に対する信用集合が欲しいなら、次を計算する。

lw_bd_observed = post_mean - 2 * np.sqrt(np.diag(post_cov) + post_sig_est ** 2)
up_bd_observed = post_mean + 2 * np.sqrt(np.diag(post_cov) + post_sig_est ** 2)

不確実性には2つの源がある。認識論的不確実性 は 削減可能 な不確実性を表し、偶然的不確実性 または 不可約不確実性 もある。ここでの 認識論的不確実性 は、ノイズのない真の関数値に関する不確実性である。データから離れるにつれてこの不確実性は増大すべきである。なぜなら、データから離れるほど、データと整合する関数値の候補が増えるからである。より多くのデータを観測するにつれて、真の関数に対する信念はより確信を増し、認識論的不確実性は消えていく。この例における 偶然的不確実性 は観測ノイズである。データはこのノイズ付きで与えられており、減らすことはできない。

データにおける 認識論的不確実性 は、潜在的なノイズのない関数の分散 np.diag(post_cov) によって捉えられる。偶然的不確実性 はノイズ分散 post_sig_est**2 によって捉えられる。

残念ながら、不確実性の表し方について人々はしばしば不注意である。多くの論文では、エラーバーがまったく定義されていなかったり、認識論的不確実性と偶然的不確実性のどちらを可視化しているのか、あるいは両方なのかが不明瞭だったり、ノイズ分散とノイズ標準偏差、標準偏差と標準誤差、信頼区間と信用集合などを混同していたりする。不確実性が何を表しているのかを正確にしなければ、それは本質的に無意味である。

不確実性が何を表しているのかに注意を払うという観点から、ここではノイズのない関数の分散推定値の 平方根 に 2倍 を掛けていることが重要である。予測分布はガウス分布なので、この量によって 95% の信用集合を作ることができ、真の関数を 95% の確率で含むと考えられる区間に対する信念を表す。ノイズの 分散 はまったく別のスケールにあり、はるかに解釈しにくいである。

最後に、20 個の事後サンプルを見てみよう。これらのサンプルは、事後的にどのような関数がデータに適合しうると考えているかを示す。

post_samples = np.random.multivariate_normal(post_mean, post_cov, size=20)
d2l.plt.scatter(train_x, train_y)
d2l.plt.plot(test_x, test_y, linewidth=2.)
d2l.plt.plot(test_x, post_mean, linewidth=2.)
d2l.plt.plot(test_x, post_samples.T, color='gray', alpha=0.25)
d2l.plt.fill_between(test_x, lw_bd, up_bd, alpha=0.25)
plt.legend(['Observed Data', 'True Function', 'Predictive Mean', 'Posterior Samples'])
d2l.plt.show()

基本的な回帰アプリケーションでは、事後予測平均と標準偏差を、それぞれ点予測器と不確実性の指標として使うのが最も一般的である。モンテカルロ獲得関数を用いたベイズ最適化や、モデルベース RL のためのガウス過程のような、より高度な応用では、事後サンプルを取る必要があることがよくある。しかし、基本的な応用で厳密には必要でなくても、これらのサンプルはデータへの当てはまりについての直感を与えてくれ、可視化に含めると有用なことが多いである。

18.3.5. GPyTorch で簡単にする

見てきたように、基本的なガウス過程回帰は実際にはゼロからでもかなり簡単に実装できる。しかし、さまざまなカーネルを試したくなったり、近似推論を考えたり（これは分類でも必要である）、GP とニューラルネットワークを組み合わせたり、あるいは 10,000 点程度を超えるデータセットを扱ったりすると、ゼロからの実装は扱いにくく煩雑になる。SKI（KISS-GP とも呼ばれる）のようなスケーラブルな GP 推論の有力な手法の中には、高度な数値線形代数ルーチンを実装する数百行のコードを必要とするものもある。

このような場合、GPyTorch ライブラリは大いに役立ちる。GPyTorch については、ガウス過程の数値計算や高度な手法に関する今後のノートブックでさらに詳しく扱う。GPyTorch ライブラリには 多くの例 がある。パッケージの雰囲気をつかむために、単純な回帰の例 を見て、上の結果を GPyTorch で再現するようにどう適応できるかを示す。これは、上の基本的な回帰を再現するだけにしてはコードが多く見えるかもしれないし、ある意味ではその通りである。しかし、数千行の新しいコードを書く代わりに、下の数行を変えるだけで、さまざまなカーネル、スケーラブルな推論手法、近似推論をすぐに使えるようになる。

# First let's convert our data into tensors for use with PyTorch
train_x = torch.tensor(train_x)
train_y = torch.tensor(train_y)
test_y = torch.tensor(test_y)

# We are using exact GP inference with a zero mean and RBF kernel
class ExactGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ZeroMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.RBFKernel())

    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

このコードブロックでは、データを GPyTorch で使える形式に変換し、厳密推論を使うこと、そして使いたい平均関数（ゼロ）とカーネル関数（RBF）を指定している。例えば gpytorch.kernels.matern_kernel() や gpyotrch.kernels.spectral_mixture_kernel() を呼び出すだけで、他のカーネルも簡単に使える。ここまでで扱ってきたのは厳密推論だけであり、近似を行わずに予測分布を推定できる場合である。ガウス過程では、ガウス尤度がある場合にのみ厳密推論が可能である。より具体的には、観測がガウス過程で表されるノイズのない関数にガウス雑音が加わって生成されると仮定する場合である。今後のノートブックでは、これらの仮定が成り立たない分類などの他の設定を扱う。

# Initialize Gaussian likelihood
likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = ExactGPModel(train_x, train_y, likelihood)
training_iter = 50
# Find optimal model hyperparameters
model.train()
likelihood.train()
# Use the adam optimizer, includes GaussianLikelihood parameters
optimizer = torch.optim.Adam(model.parameters(), lr=0.1)
# Set our loss as the negative log GP marginal likelihood
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

ここでは、使いたい尤度（ガウス）、カーネルのハイパーパラメータを学習するために使う目的関数（ここでは周辺尤度）、そしてその目的関数を最適化するために使う手順（この場合は Adam）を明示的に指定している。Adam は「確率的」最適化手法であるが、この場合はフルバッチ Adam であることに注意しよ。周辺尤度はデータごとに因数分解できないため、データの「ミニバッチ」に対する最適化器を使っても収束が保証されない。L-BFGS のような他の最適化器も GPyTorch でサポートされている。標準的な深層学習とは異なり、周辺尤度をうまく最適化できることは良い汎化と強く結びついているため、計算コストが高すぎない限り、L-BFGS のような強力な最適化器を使いたくなることが多いである。

for i in range(training_iter):
    # Zero gradients from previous iteration
    optimizer.zero_grad()
    # Output from model
    output = model(train_x)
    # Calc loss and backprop gradients
    loss = -mll(output, train_y)
    loss.backward()
    if i % 10 == 0:
        print(f'Iter {i+1:d}/{training_iter:d} - Loss: {loss.item():.3f} '
              f'squared lengthscale: '
              f'{model.covar_module.base_kernel.lengthscale.item():.3f} '
              f'noise variance: {model.likelihood.noise.item():.3f}')
    optimizer.step()

Iter 1/50 - Loss: 0.982 squared lengthscale: 0.693 noise variance: 0.693
Iter 11/50 - Loss: 0.696 squared lengthscale: 0.519 noise variance: 0.313
Iter 21/50 - Loss: 0.441 squared lengthscale: 0.535 noise variance: 0.128
Iter 31/50 - Loss: 0.336 squared lengthscale: 0.521 noise variance: 0.057
Iter 41/50 - Loss: 0.353 squared lengthscale: 0.517 noise variance: 0.042

ここで実際に最適化手順を実行し、10 イテレーションごとに損失の値を出力している。

# Get into evaluation (predictive posterior) mode
test_x = torch.tensor(test_x)
model.eval()
likelihood.eval()
observed_pred = likelihood(model(test_x))

上のコードブロックにより、テスト入力に対して予測を行えるようになる。

with torch.no_grad():
    # Initialize plot
    f, ax = d2l.plt.subplots(1, 1, figsize=(4, 3))
    # Get upper and lower bounds for 95\% credible set (in this case, in
    # observation space)
    lower, upper = observed_pred.confidence_region()
    ax.scatter(train_x.numpy(), train_y.numpy())
    ax.plot(test_x.numpy(), test_y.numpy(), linewidth=2.)
    ax.plot(test_x.numpy(), observed_pred.mean.numpy(), linewidth=2.)
    ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.25)
    ax.set_ylim([-1.5, 1.5])
    ax.legend(['True Function', 'Predictive Mean', 'Observed Data',
               '95% Credible Set'])

最後に、当てはまりを描画する。

当てはまりは事実上同一であることがわかる。いくつか注意点がある。GPyTorch は 二乗された 長さ尺度と観測ノイズを扱っている。例えば、ゼロから書いたコードで学習されたノイズ標準偏差は約 0.283 でした。GPyTorch が見つけたノイズ分散は \(0.81 \approx 0.283^2\) である。GPyTorch の図では、潜在関数空間ではなく 観測空間 における信用集合も示しており、実際に観測データ点を覆っていることを示している。

18.3.6. まとめ

ガウス過程事前分布とデータを組み合わせて事後分布を作り、それを予測に使うことができる。また、周辺尤度を作ることもでき、これはガウス過程の変動速度などの性質を制御するカーネルハイパーパラメータの自動学習に役立つ。回帰における事後の構成とカーネルハイパーパラメータの学習の仕組みは単純で、コードはおよそ十数行で済む。このノートブックは、ガウス過程をすぐに「使い始めたい」読者にとって良い参考資料である。また、GPyTorch ライブラリも紹介した。基本的な回帰のための GPyTorch コードは比較的長いものの、他のカーネル関数や、今後のノートブックで扱うより高度な機能――スケーラブル推論や分類のための非ガウス尤度など――へは、数行変えるだけで簡単に拡張できる。

18.3.7. 演習

1. カーネルハイパーパラメータを 学習する ことの重要性と、ハイパーパラメータやカーネルがガウス過程の汎化特性に与える影響を強調してきた。ハイパーパラメータを学習する手順を飛ばし、代わりにさまざまな長さ尺度とノイズ分散を仮定して、予測への影響を確認してみせよ。長さ尺度が大きいとどうなるか。小さいとどうなるか。ノイズ分散が大きいとどうなるか。小さいとどうなるか。

2. 周辺尤度は凸最適化問題ではないが、長さ尺度やノイズ分散のようなハイパーパラメータは GP 回帰で信頼性高く推定できる、と述べた。これは一般に正しいである。実際、周辺尤度は、経験的自己相関関数（「コビアログラム」）を当てはめる空間統計学の従来手法よりも、長さ尺度ハイパーパラメータの学習に はるかに 優れている。少なくとも最近のスケーラブル推論の研究以前において、機械学習がガウス過程研究に与えた最大の貢献は、ハイパーパラメータ学習のための周辺尤度の導入だったと言えるだろう。

しかし、これらのパラメータの組み合わせが異なるだけで、多くのデータセットに対して解釈可能でもっともらしい説明が異なり、目的関数に局所最適が生じる。長さ尺度が大きい場合、真の基礎関数はゆっくり変化すると仮定していることになる。観測データが 実際に 大きく変動しているなら、大きな長さ尺度を正当化できるのは、大きなノイズ分散がある場合だけである。逆に、長さ尺度が小さい場合、当てはまりはデータの変動に非常に敏感になり、ノイズ（偶然的不確実性）で変動を説明する余地がほとんどなくなる。

これらの局所最適を見つけられるか試してみよ。大きな長さ尺度と大きなノイズ、そして小さな長さ尺度と小さなノイズで初期化してみよう。異なる解に収束するか？

3. ベイズ法の根本的な利点の1つは、認識論的不確実性 を自然に表現できることだと述べた。上の例では、認識論的不確実性の効果を完全には見ることができない。代わりに test_x = np.linspace(0, 10, 1000) で予測してみよ。予測がデータを超えて進むにつれて、95% の信用集合はどうなるか。その区間で真の関数を覆いますか。その領域で偶然的不確実性だけを可視化するとどうなるか。

4. 上の例を、訓練点数を 10,000、20,000、40,000 にして実行し、実行時間を測定してみよ。訓練時間はどのようにスケールするか。あるいは、テスト点数に対して実行時間はどうスケールするか。予測平均と予測分散で違いはあるか。理論的に訓練・テスト時間の計算量を求めることと、上のコードを異なる点数で実行することの両方で答えよ。

5. GPyTorch の例を、Matern カーネルなど異なる共分散関数で実行してみよ。結果はどう変わりますか。GPyTorch ライブラリにある spectral mixture カーネルはどうだろうか。周辺尤度で学習しやすいものとそうでないものはあるか。長距離予測と短距離予測で有用性に違いはあるか。

6. GPyTorch の例では観測ノイズを含めた予測分布を描いたが、ゼロからの例では認識論的不確実性だけを含めた。GPyTorch の例をやり直し、今度は認識論的不確実性だけを描画して、ゼロからの結果と比較せよ。予測分布は同じように見えるか？（同じはずである。）