5. ガウス過程回帰
5.1. 動機付け
式 (1.13) において、データ数とパラメータ数が等しい、すなわち \(\boldsymbol{\Phi}\) が正方行列となるケースを考える
(5.5)\[\begin{split}\boldsymbol{y}
\equiv
\begin{pmatrix}
y_1 \\
y_2 \\
\vdots \\
y_N
\end{pmatrix}
=
\begin{pmatrix}
\phi_0(\boldsymbol{x}_1) & \phi_1(\boldsymbol{x}_1) & \cdots & \phi_{N-1}(\boldsymbol{x}_1) \\
\phi_0(\boldsymbol{x}_2) & \phi_1(\boldsymbol{x}_2) & \cdots & \phi_{N-1}(\boldsymbol{x}_2) \\
\vdots & \cdots & \cdots & \vdots \\
\phi_0(\boldsymbol{x}_N) & \phi_1(\boldsymbol{x}_N) & \cdots & \phi_{N-1}(\boldsymbol{x}_N) \\
\end{pmatrix}
\begin{pmatrix}
w_0 \\
w_1 \\
\vdots \\
w_{N-1}
\end{pmatrix}
\equiv
\boldsymbol{\Phi}\boldsymbol{w} .\end{split}\]
さらに \(\boldsymbol{\Phi}\) が正則であるとして、上式を \(\boldsymbol{w}\) について解く
(5.6)\[\boldsymbol{w}
= \boldsymbol{\Phi}^{-1} \boldsymbol{y} .\]
パラメータ \(\boldsymbol{w}\) に対する事前分布として
(5.7)\[p(\boldsymbol{w})
\propto
\exp{\left(-\frac{\beta}{2}\boldsymbol{w}^\top\boldsymbol{w}\right)}\]
を考え、ここに式 (5.6) を代入すると次式が得られる
(5.8)\[\exp{\left(-\frac{\beta}{2}\boldsymbol{w}^\top\boldsymbol{w}\right)}
=
\exp{\left(-\frac{\beta}{2}\boldsymbol{y}^\top\boldsymbol{\Phi}^{-\top}\boldsymbol{\Phi}^{-1}\boldsymbol{y}\right)}
=
\exp{\left(-\frac{\beta}{2}\boldsymbol{y}^\top
\left(\boldsymbol{\Phi}\boldsymbol{\Phi}^\top\right)^{-1}\boldsymbol{y}\right)} .\]
上式はパラメータ \(\boldsymbol{w}\) が消去され、代わりに \(\boldsymbol{y}\) に対する事前分布を与えた式となっている。式 (5.5) から
(5.9)\[\begin{split}\boldsymbol{\Phi}\boldsymbol{\Phi}^\top
=
\begin{pmatrix}
\phi_0(\boldsymbol{x}_1) & \phi_1(\boldsymbol{x}_1) & \cdots & \phi_{N-1}(\boldsymbol{x}_1) \\
\phi_0(\boldsymbol{x}_2) & \phi_1(\boldsymbol{x}_2) & \cdots & \phi_{N-1}(\boldsymbol{x}_2) \\
\vdots & \cdots & \cdots & \vdots \\
\phi_0(\boldsymbol{x}_N) & \phi_1(\boldsymbol{x}_N) & \cdots & \phi_{N-1}(\boldsymbol{x}_N) \\
\end{pmatrix}
\begin{pmatrix}
\phi_0(\boldsymbol{x}_1) & \phi_0(\boldsymbol{x}_2) & \cdots & \phi_0(\boldsymbol{x}_N) \\
\phi_1(\boldsymbol{x}_1) & \phi_1(\boldsymbol{x}_2) & \cdots & \phi_1(\boldsymbol{x}_N) \\
\vdots & \cdots & \cdots & \vdots \\
\phi_{N-1}(\boldsymbol{x}_1) & \phi_{N-1}(\boldsymbol{x}_2) & \cdots & \phi_{N-1}(\boldsymbol{x}_N) \\
\end{pmatrix}\end{split}\]
であり、
(5.10)\[k(\boldsymbol{x}_i, \boldsymbol{x}_j)
\equiv
\sum_{n=0}^{N-1} \phi_n(\boldsymbol{x}_i) \phi_n(\boldsymbol{x}_j)\]
という関数を定義して、これを要素とする行列を
(5.11)\[\begin{split}\boldsymbol{K}
\equiv
\begin{pmatrix}
k(\boldsymbol{x}_1, \boldsymbol{x}_1) & k(\boldsymbol{x}_1, \boldsymbol{x}_2) & \cdots & k(\boldsymbol{x}_1, \boldsymbol{x}_N) \\
k(\boldsymbol{x}_2, \boldsymbol{x}_1) & k(\boldsymbol{x}_2, \boldsymbol{x}_2) & \cdots & k(\boldsymbol{x}_2, \boldsymbol{x}_N) \\
\vdots & \cdots & \cdots & \vdots \\
k(\boldsymbol{x}_N, \boldsymbol{x}_1) & k(\boldsymbol{x}_N, \boldsymbol{x}_2) & \cdots & k(\boldsymbol{x}_N, \boldsymbol{x}_N) \\
\end{pmatrix}\end{split}\]
と置くと、\(\boldsymbol{y}\) の事前分布は次のように書ける
(5.12)\[p(\boldsymbol{y})
=
\frac{1}{\sqrt{(2\pi)^N\left| K \right|}}
\exp{\left(
-\frac{1}{2}
\boldsymbol{y}^\top \boldsymbol{K}^{-1} \boldsymbol{y}
\right)} .\]
以上より、パラメータ\(\boldsymbol{w}\) を経由することなく、 \(\boldsymbol{y}\) についての確率分布を直接的に扱う定式化の可能性が示唆される。そうすることで、基底関数やパラメータの数の増大で計算が困難になる、いわゆる「次元の呪い」からの解放される事になる。
5.2. ガウス過程との関係
いま、多数の（あるいは連続的な）データから、任意に \(N\) 個のデータを取り出すことを考える
(5.13)\[\{y_1, y_2, \dots, y_N\} .\]
この\(N\) 個のデータの同時分布が、式 (5.12) のような多変量ガウス分布であるとき、これは「ガウス過程」と呼ばれる確率過程の定義そのものになっている。もちろん今はデータのラベルに時間の意味はないが、この対応から、式 (5.12) の確率分布を利用した回帰は「ガウス過程回帰」と呼ばれる。
5.3. 出力にノイズがある場合
出力データにノイズがあり、次のような形でデータが生成される場合を考える
(5.14)\[\boldsymbol{t} = \boldsymbol{y} + \boldsymbol{\varepsilon} .\]
\(\boldsymbol{y}\) が与えられた時の \(\boldsymbol{t}\) の分布を
(5.15)\[p(\boldsymbol{t} | \boldsymbol{y})
=
\frac{1}{\sqrt{(2\pi)^N(\sigma_\varepsilon^2)^N}}
\exp{\left(
-\frac{1}{2\sigma_\varepsilon^2}
(\boldsymbol{t}-\boldsymbol{y})^\top (\boldsymbol{t}-\boldsymbol{y})
\right)}\]
でモデル化するとき、\(\boldsymbol{t}\) の周辺分布は
(5.16)\[p(\boldsymbol{t})
=
\int
p(\boldsymbol{t} | \boldsymbol{y})
p(\boldsymbol{y})
d\boldsymbol{y}\]
と計算される。これを (5.12) および (5.15) を用いて計算すると
(5.17)\[p(\boldsymbol{t})
=
\frac{1}{\sqrt{(2\pi)^N\left| \boldsymbol{K} + \beta^{-1}\boldsymbol{I} \right|}}
\exp{\left(
-\frac{1}{2}
\boldsymbol{t}^\top
\left( \boldsymbol{K} + \beta^{-1}\boldsymbol{I} \right)^{-1}
\boldsymbol{t}
\right)}\]
が得られる。これは、共分散行列を
(5.18)\[\boldsymbol{K} \longrightarrow \boldsymbol{K} + \beta^{-1}\boldsymbol{I}\]
と置き換えれば (5.12) と本質的には同じものであることがわかる。
5.4. 予測分布
いま、\(N_1\) 個の入出力データのセットがあるとして、新たな \(N\) 個の入力と対応する出力を
(5.19)\[\begin{split}\boldsymbol{X}
&\equiv
\{\boldsymbol{x}_{1}, \boldsymbol{x}_{2}, \dots, \boldsymbol{x}_{N}\} ,\\\\
\boldsymbol{t}
&\equiv
\{t_{1}, t_{2}, \dots, t_{N}\}\end{split}\]
と表し、入出力がそろっている \(N_1\) 個の（既知の）データを
(5.20)\[\begin{split}\boldsymbol{X}_{N_1}
&\equiv
\{\boldsymbol{x}_{N+1}, \boldsymbol{x}_{N+2}, \dots, \boldsymbol{x}_{N+N_1}\} ,\\\\
\boldsymbol{t}_{N_1}
&\equiv
\{t_{N+1}, t_{N+2}, \dots, t_{N+N_1}\} \\\end{split}\]
と表すことにする。関数 \(k(\boldsymbol{x}, \boldsymbol{x}')\) の形も与えられているとして、 共分散行列をブロック行列の形で以下のように定義する
(5.21)\[\begin{split}\begin{pmatrix}
\boldsymbol{K}_{NN} & \boldsymbol{K}_{NN_1} \\
\boldsymbol{K}_{N_1N} & \boldsymbol{K}_{N_1N_1} \\
\end{pmatrix}
\equiv
\begin{pmatrix}
k(\boldsymbol{x}_1, \boldsymbol{x}_1) & \cdots & k(\boldsymbol{x}_1, \boldsymbol{x}_{N}) & \cdots & k(\boldsymbol{x}_1, \boldsymbol{x}_{N+N_1}) \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
k(\boldsymbol{x}_N, \boldsymbol{x}_1) & \cdots & k(\boldsymbol{x}_N, \boldsymbol{x}_{N}) & \cdots & k(\boldsymbol{x}_N, \boldsymbol{x}_{N+N_1}) \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
k(\boldsymbol{x}_{N+N_1}, \boldsymbol{x}_1) & \cdots & k(\boldsymbol{x}_{N+N_1}, \boldsymbol{x}_N) & \cdots & k(\boldsymbol{x}_{N+N_1}, \boldsymbol{x}_{N+N_1}) \\
\end{pmatrix} .\end{split}\]
ここで \(\boldsymbol{t}\) と \(\boldsymbol{t_{N_1}}\) の同時分布を考える。指数部の中身だけ書き出すと
(5.22)\[\begin{split}-\frac{1}{2}
\left(
\boldsymbol{t}^\top \;\;
\boldsymbol{t}_{N_1}^\top
\right)
\begin{pmatrix}
\boldsymbol{K}_{N N} & \boldsymbol{K}_{N N_1} \\
\boldsymbol{K}_{N_1 N} & \boldsymbol{K}_{N_1 N_1} \\
\end{pmatrix}^{-1}
\begin{pmatrix}
\boldsymbol{t} \\
\boldsymbol{t}_{N_1} \\
\end{pmatrix}
\equiv
-\frac{1}{2}
\left(
\boldsymbol{t}^\top \;\;
\boldsymbol{t}_{N_1}^\top
\right)
\begin{pmatrix}
\boldsymbol{L}_{N N} & \boldsymbol{L}_{N N_1} \\
\boldsymbol{L}_{N_1 N} & \boldsymbol{L}_{N_1 N_1} \\
\end{pmatrix}
\begin{pmatrix}
\boldsymbol{t} \\
\boldsymbol{t}_{N_1} \\
\end{pmatrix} .\end{split}\]
ここで、共分散行列の逆行列（精度行列）の各ブロックは、ブロック行列の逆行列の公式を用いれば次のようになる
\[\begin{split}\boldsymbol{L}_{NN}
&=
\left(
\boldsymbol{K}_{NN}
-\boldsymbol{K}_{NN_1} \boldsymbol{K}^{-1}_{N_1N_1} \boldsymbol{K}_{N_1N}
\right)^{-1} ,\\\\
\boldsymbol{L}_{NN_1}
&=
-\left(
\boldsymbol{K}_{NN}
-\boldsymbol{K}_{NN_1} \boldsymbol{K}^{-1}_{N_1N_1} \boldsymbol{K}_{N_1N}
\right)^{-1}
\boldsymbol{K}_{NN_1} \boldsymbol{K}^{-1}_{N_1N_1} ,\\\\
\boldsymbol{L}_{N_1N}
&=
-\boldsymbol{K}^{-1}_{N_1N_1} \boldsymbol{K}_{N_1N}
\left(
\boldsymbol{K}_{NN}
-\boldsymbol{K}_{NN_1} \boldsymbol{K}^{-1}_{N_1N_1} \boldsymbol{K}_{N_1N}
\right)^{-1} ,\\\\
\boldsymbol{L}_{N_1N_1}
&=
\boldsymbol{K}^{-1}_{N_1N_1}
-
\boldsymbol{K}^{-1}_{N_1N_1}
\boldsymbol{K}_{N_1N}
\left(
\boldsymbol{K}_{NN}
-\boldsymbol{K}_{NN_1} \boldsymbol{K}^{-1}_{N_1N_1} \boldsymbol{K}_{N_1N}
\right)^{-1}
\boldsymbol{K}_{NN_1}
\boldsymbol{K}^{-1}_{N_1N_1}\end{split}\]
となる。これらを用いて (5.22) は次のように書ける
\[\begin{split}-\frac{1}{2}
\left(
\boldsymbol{t}^\top \;\;
\boldsymbol{t}_{N_1}^\top
\right)
\begin{pmatrix}
\boldsymbol{L}_{N N} & \boldsymbol{L}_{N N_1} \\
\boldsymbol{L}_{N_1 N} & \boldsymbol{L}_{N_1 N_1} \\
\end{pmatrix}
\begin{pmatrix}
\boldsymbol{t} \\
\boldsymbol{t}_{N_1} \\
\end{pmatrix}
&=
-\frac{1}{2}
\left(
\boldsymbol{t}^\top \boldsymbol{L}_{NN} \boldsymbol{t}
+ \boldsymbol{t}^\top \boldsymbol{L}_{NN_1} \boldsymbol{t}_{N_1}
+ \boldsymbol{t}_{N_1}^\top \boldsymbol{L}_{N_1N} \boldsymbol{t}
+ \boldsymbol{t}_{N_1}^\top \boldsymbol{L}_{N_1N_1} \boldsymbol{t}_{N_1}
\right) \\
&=
-\frac{1}{2}
\left(
\boldsymbol{t}^\top + \boldsymbol{t}_{N_1}^\top \boldsymbol{L}_{N_1N} \boldsymbol{L}^{-1}_{NN}
\right)
\boldsymbol{L}_{NN}
\left(
\boldsymbol{t} + \boldsymbol{L}^{-1}_{NN} \boldsymbol{L}_{NN_1} \boldsymbol{t}_{N_1}
\right)
-\frac{1}{2}
\boldsymbol{t}_{N_1}^\top
\left(
\boldsymbol{L}_{N_1N_1} - \boldsymbol{L}_{N_1N} \boldsymbol{L}^{-1}_{NN} \boldsymbol{L}_{NN_1}
\right)
\boldsymbol{t}_{N_1} .\end{split}\]
\(\boldsymbol{K}\) に戻して式を整理する
\[\begin{split}\boldsymbol{L}_{N_1N} \boldsymbol{L}^{-1}_{NN}
&=
-\boldsymbol{K}^{-1}_{N_1N_1} \boldsymbol{K}_{N_1N} ,\\
\\
\boldsymbol{L}^{-1}_{NN} \boldsymbol{L}_{NN_1}
&=
\left( \boldsymbol{L}_{N_1N} \boldsymbol{L}^{-1}_{NN} \right)^\top
=
-\boldsymbol{K}_{NN_1} \boldsymbol{K}^{-1}_{N_1N_1} ,\\
\\
\boldsymbol{L}_{N_1N_1} - \boldsymbol{L}_{N_1N} \boldsymbol{L}^{-1}_{NN} \boldsymbol{L}_{NN_1}
&=
\boldsymbol{K}^{-1}_{N_1N_1} .\\\end{split}\]
以上をまとめると
(5.23)\[\begin{split}-\frac{1}{2}
\left(
\boldsymbol{t}^\top \;\;
\boldsymbol{t}_{N_1}^\top
\right)
\begin{pmatrix}
\boldsymbol{L}_{N N} & \boldsymbol{L}_{N N_1} \\
\boldsymbol{L}_{N_1 N} & \boldsymbol{L}_{N_1 N_1} \\
\end{pmatrix}
\begin{pmatrix}
\boldsymbol{t}\\
\boldsymbol{t}_{N_1}\\
\end{pmatrix}
&=
-\frac{1}{2}
\left(
\boldsymbol{t}^\top - \boldsymbol{t}_{N_1}^\top \boldsymbol{K}^{-1}_{N_1N_1} \boldsymbol{K}_{N_1N}
\right)
\boldsymbol{L}_{NN}
\left(
\boldsymbol{t} -\boldsymbol{K}_{NN_1} \boldsymbol{K}^{-1}_{N_1N_1} \boldsymbol{t}_{N_1}
\right)
-\frac{1}{2}
\boldsymbol{t}_{N_1}^\top
\boldsymbol{K}^{-1}_{N_1N_1}
\boldsymbol{t}_{N_1}\end{split}\]
となる。確率の同時分布と条件付き分布の関係式
(5.24)\[p(\boldsymbol{t}, \boldsymbol{t}_{N_1})
=
p(\boldsymbol{t} | \boldsymbol{t}_{N_1}) p(\boldsymbol{t}_{N_1})\]
より、
(5.25)\[p(\boldsymbol{t} | \boldsymbol{t}_{N_1})
\propto
\exp{
\left(
-\frac{1}{2}
\left(
\boldsymbol{t}^\top - \boldsymbol{t}_{N_1}^\top \boldsymbol{K}^{-1}_{N_1N_1} \boldsymbol{K}_{N_1N}
\right)
\boldsymbol{L}_{NN}
\left(
\boldsymbol{t} -\boldsymbol{K}_{NN_1} \boldsymbol{K}^{-1}_{N_1N_1} \boldsymbol{t}_{N_1}
\right)
\right)
}\]
が得られ、これが求めようとしていた \(\boldsymbol{t}_{N_1}\) が既知である状況での、\(\boldsymbol{t}\) の条件付き確率分布となる。
5.5. 予測関数
式 (5.25) より、平均
(5.26)\[\boldsymbol{m}
\equiv
\boldsymbol{K}_{NN_1} \boldsymbol{K}^{-1}_{N_1N_1} \boldsymbol{t}_{N_1}\]
が、(5.19) の入力に対する出力の予測となり、(5.25) は、その誤差を与える事になる。式 (5.26) を書き下すと次のようになる
(5.27)\[\begin{split}\begin{pmatrix}
m(\boldsymbol{x}_{1}) \\
m(\boldsymbol{x}_{2}) \\
\vdots \\
m(\boldsymbol{x}_{N}) \\
\end{pmatrix}
=
\begin{pmatrix}
\sum_{i=N+1}^{N+N_1} \sum_{j=N+1}^{N+N_1}
k(\boldsymbol{x}_{1}, \boldsymbol{x}_i) k^{-1}(\boldsymbol{x}_i, \boldsymbol{x}_j) t(\boldsymbol{x}_j) \\
\sum_{i=N+1}^{N+N_1} \sum_{j=N+1}^{N+N_1}
k(\boldsymbol{x}_{2}, \boldsymbol{x}_i) k^{-1}(\boldsymbol{x}_i, \boldsymbol{x}_j) t(\boldsymbol{x}_j) \\
\vdots \\
\sum_{i=N+1}^{N+N_1} \sum_{j=N+1}^{N+N_1}
k(\boldsymbol{x}_{N}, \boldsymbol{x}_i) k^{-1}(\boldsymbol{x}_i, \boldsymbol{x}_j) t(\boldsymbol{x}_j)
\end{pmatrix}\end{split}\]
となる。これは「平均」として得られる関数の予測が
(5.28)\[m(\boldsymbol{x}) = \sum_{i=1}^{N_1} c_i k(\boldsymbol{x}, \boldsymbol{x}_i)\]
という形で表せることを示している。ここで
(5.29)\[c_i
\equiv
\sum_{j=1}^{N_1} k^{-1}(\boldsymbol{x}_i, \boldsymbol{x}_j) t(\boldsymbol{x}_j)\]
とおき、また既知データの添え字を \(1, 2, \dots, N_1\) と振り直した。式 (5.28) は「再生核ヒルベルト空間」というものと密接な関係がある。 詳細はヒルベルト空間の項に譲り、ここでは主要な結果を列挙するだけにとどめる。
再生核（定義）
関数 \(f(x)\) と、ヒルベルト空間の要素としての \(f\) をあまり区別せず、おおらかに定義を述べると、
\[\int k(x,x') f(x') dx' = f(x)\]
となるような関数 \(k(x,x')\) を「再生核」と呼ぶ。
ムーア・アロンシャインの定理
正定値対称な関数 \(k(x,x')\) があれば、これに付随するヒルベルト空間（再生核ヒルベルト空間）が一意に定まる
Representer定理
次のような正則化項付きの損失関数を考える
\[L
=
E\left[ (x_1, y_1, f(x_1)), (x_2, y_2, f(x_2)), \dots, (x_N, y_N, f(x_N)) \right]
+
\Phi(\|f\|^2) ,\]
ここで正則化項 \(\Phi\) は単調増加関数であるとする。このとき、上式を最小化する \(f(x)\) は、次の形に書ける
\[f(x) = \sum_{i=1}^N c_i k(x, x_i) .\]