はじめに
ガウス過程回帰は機械学習領域ではかなりメジャーなので、皆様御存知かと思いますが、
ここで本当に初学者がガウス過程回帰にいきつくような解説と、多次元まで対応したPythonのガウス過程回帰クラスを提示したいと思います。実装については第2回で書きます。
まずは線形回帰モデルを考える
まずは泣く子も黙る線形回帰モデルを考えることにします。
$w$はパラメータベクトル、$\phi(x)$は非線形の基底関数です。
$$
y(\boldsymbol{x}, \boldsymbol{w})=\boldsymbol{w}^T\phi(\boldsymbol{x}),
$$
これは以前の投稿ベイズ線形クラスのPython実装で紹介したものです。
少し先走って書くと、$w$の事前分布を決めれば、関数$y(w,x)$に対する事前分布が決まり、訓練集合が与えられると$w$の事後分布が決まり、結果として関数の事後分布が求まり、そこから新しいデータ$x$に対する予測分布$p(t|x)$が導出されます。
wの事前分布
wの事前分布としてガウス分布を考えます。
$$
p(\boldsymbol{w})=\mathcal{N}(\boldsymbol{w}|\boldsymbol{0},\alpha^{-1}\boldsymbol{I})
$$
ある要素$\boldsymbol{x}_n$に対する関数の値$y_n$を持つ集合を$\boldsymbol{y}$で表すと、下記のように表されます。
$$
\boldsymbol{y}=\Phi\boldsymbol{w}
$$
ここで$\Phi$は要素$\Phi_{nk}=\phi_k(\boldsymbol{x}_n)$を持つ計画行列です。
yの分布
$\boldsymbol{y}$はガウス分布$\boldsymbol{w}$の線形結合になるので、ガウス分布に従うことがすぐに理解できます。
\begin{array}
\mathbb{E}[\boldsymbol{y}]=\Phi\mathbb{E}[\boldsymbol{w}]=0 \\
\mathbb{V}[\boldsymbol{y}]=\mathbb{E}[\boldsymbol{y}\boldsymbol{y}^T]-\mathbb{E}[\boldsymbol{y}]^2=\Phi\mathbb{E}[\boldsymbol{w}\boldsymbol{w}^T]\Phi^T=\frac{1}{\alpha}\Phi\Phi^T
\end{array}
この$\frac{1}{\alpha}\Phi\Phi^T$は、$k$はカーネル関数とした、下記のような要素を持つ$n \times m$のグラム行列です。
$$
K_{nm}=k(\boldsymbol{x}_n, \boldsymbol{x}_m)=\frac{1}{\alpha}\phi(\boldsymbol{x}_n)^T\phi(\boldsymbol{x}_m)
$$
とてもシンプルにカーネル関数のみを用いて$y(x)$の同時分布を記述することが出来ました。
結局、ガウス過程って何？
任意の点集合$\boldsymbol{x}_1,...,\boldsymbol{x}_N$に対する$y(\boldsymbol{x})$の同時分布が、ガウス分布に従うもの。
つまり$w$の事前分布をガウス分布と定義した上の線形回帰の例も、結果としてガウス過程にほかならないわけです。
ガウス分布に従うということは、$y$の同時分布がたかだか2次までの統計量(平均と共分散など)で完全に記述されることがポイントです。逆に言えば、その仮定なしにガウス過程は使えないということです。
線形回帰の例でも見ましたが、$y(x)$の平均がゼロであるとすれば、ガウス過程はもはや下記のようなカーネル関数だけを考えれば良くなります。これは任意の2つの入力$\boldsymbol{x}$に対する$y(x)$の共分散です。
$$
\mathbb{E}[y(\boldsymbol{x}_n)y(\boldsymbol{x}_m)]=k(\boldsymbol{x}_n, \boldsymbol{x}_m)
$$
線形回帰の例では、基底関数として$\phi$を使っていたので、カーネル関数は下記のようになりました。
$$
k(\boldsymbol{x}_n, \boldsymbol{x}_m)=\frac{1}{\alpha}\phi(\boldsymbol{x}_n)^T\phi(\boldsymbol{x}_m)
$$
カーネル関数は基底関数から間接的に求めても良いですし、陽に定義することも可能です。
有名なものはガウスカーネルや指数カーネルなどがあります。
ガウス過程回帰を考える
ガウス過程の定義がわかったところで、一般的なガウス過程回帰をやってみましょう。
問題の準備をしていきます。
ガウス過程回帰の準備
観測される目的変数にはノイズが含まれることを定義しましょう。
$$
t_n=y_n+\epsilon_n
$$
ここで$\epsilon_n$はn番目の観測値に加えられるノイズで、それぞれの観測値に対して独立に決定されるとします。さらにこの$\epsilon$もガウス分布に従うとします。
$$
p(t_n|y_n)=\mathcal{N}(t_n|y_n,\beta^{-1})
$$
$\beta$は精度パラメータで、ノイズの精度を表します。これを同時分布として表すと、下記のようになります。ここで$\boldsymbol{I}_N$は$N\times N$の単位行列です。
$$
p(\boldsymbol{t}|\boldsymbol{y})=\mathcal{N}(\boldsymbol{t}|\boldsymbol{y},\beta^{-1}\boldsymbol{I}_N)
$$
ガウス過程は周辺分布$p(y)$がグラム行列で表されるガウス分布に従うという定義なので、周辺分布は下記のように表せます。
$$
p(\boldsymbol{y})=\mathcal{N}(\boldsymbol{y}|0, K)
$$
これらの式を用いて$\boldsymbol{t}$の周辺分布を計算します。
ガウス分布の周辺分布の計算は比較的簡単です。
$$
p(\boldsymbol{t})=\int p(\boldsymbol{t}|\boldsymbol{y})p(\boldsymbol{y})d\boldsymbol{y}
=\mathcal{N}(\boldsymbol{t}|0,C) \tag{1}
$$
$$
C(\boldsymbol{x}_n, \boldsymbol{x}_m)=k(\boldsymbol{x}_n, \boldsymbol{x}_m)+\beta^{-1}\delta _{nm} \tag{2}
$$
$y$と$\epsilon$を足し合わせた分布の共分散は、独立であるという前提で、それぞれ$k(\boldsymbol{x}_n, \boldsymbol{x}_m)$と、$\beta^{-1}\delta _{nm}$を足したものとしています。以上によって、データ点集合の上での同時分布をモデル化することができました。
予測分布の計算
ここから訓練データではなく、新しいデータへの予測はどうなるかを考えていきます。
notationを整理すると、下記のようになります。
- 訓練データ：
- 入力：$\boldsymbol{x}_1,...,\boldsymbol{x}_N$
- 目的変数：$\boldsymbol{t}_N=(t_1,...,t_N)^T$
- 新しいデータ
- 入力：$\boldsymbol{x}_{N+1}$
- 目的変数：$t_{N+1}$
このときに、予測分布$p(t_{N+1}|\boldsymbol{t}_ N)$を求めたい。
そのために、まずは同時分布$p(\boldsymbol{t}_{N+1})$を計算します。これは(1)式から直ちにわかります。
$$
p(\boldsymbol{t}_{N+1})=\mathcal{N}(\boldsymbol{t} _{N+1}|0, \boldsymbol{C} _{N+1})
$$
ここまで読んでいただければ御存知の通り、$\boldsymbol{C} _ {N+1}$は$N+1\times N+1$の共分散行列であり、その要素は（2）式のとおり、カーネル関数と精度パラメータで表されます。
ここから予測分布$p(t_{N+1}|\boldsymbol{t}_ N)$を計算するために、共分散行列を分割します。
\boldsymbol{C} _ {N+1}=
\begin{pmatrix}
\boldsymbol{C}_N & \boldsymbol{k} \\
\boldsymbol{k}^T & c
\end{pmatrix}
共分散行列を分割したのは、多変量ガウス分布の同時分布がガウス分布に従うとき、その変数を分割して、一方の変数集合が与えられたときに、もう一方もガウス分布に従うことはよく知られています。 この周辺分布の計算を応用するためです。この辺はまたどこかの記事でまとめておこうかと。こちらの記事にまとめました。
（具体的には、一方の変数集合が訓練データ、もう一方がN+1番目のデータとする）
結果として、$p(t_{N+1}|\boldsymbol{t}_ N)$は次の様な平均と分散で表されるガウス分布になります。
$$
\mu(\boldsymbol{x}_{N+1})=\boldsymbol{k}^T\boldsymbol{C}_N^{-1}\boldsymbol{t}
$$
$$
\sigma^2(\boldsymbol{x}_{N+1})=c-\boldsymbol{k}^T\boldsymbol{C}_N^{-1}\boldsymbol{k}
$$
ようやくよく見るガウス過程回帰の平均と分散が出てきました。
補足：ハイパーパラメータの学習
ここまで見てくると、観測点とカーネル関数が与えられれば、パラメータ無しに予測ができるようにみえます。逆に言えば、ガウス過程による予測は、共分散関数（カーネル関数）の選択に依存している部分が多いです。そのため、共分散関数を固定として持つのではなく、パラメトリックな関数として定義し、そのパラメータを学習することを考えます。
いま、尤度関数の最大化（最尤推定）を行うために、ガウス分布の尤度関数$p(\boldsymbol{t}|\boldsymbol{\theta})$を考えましょう。一般的な多次元ガウス分布の尤度関数を計算してみると、下記のような式が導かれます。
$$
\ln p(\boldsymbol{t}|\boldsymbol{\theta})=-\frac{1}{2}\ln|\boldsymbol{C}_N|-\frac{1}{2}\boldsymbol{t}^T\boldsymbol{C}_N^{-1}\boldsymbol{t}-\frac{N}{2}\ln(2\pi)
$$
少し計算すると、この勾配は、下記のようになります。
$$
\frac{\partial}{\partial\theta_i}\ln p(\boldsymbol{t}|\boldsymbol{\theta})=-\frac{1}{2}\mathrm{Tr}\left(\boldsymbol{C}_N^{-1}\frac{\partial\boldsymbol{C}_N}{\partial\theta_i}\right)+\frac{1}{2}\boldsymbol{t}^T\boldsymbol{C}_N^{-1}\frac{\partial\boldsymbol{C}_N}{\partial\theta_i}\boldsymbol{C}_N^{-1}\boldsymbol{t}
$$
具体的にカーネル関数を定義します。PRMLに従って、ガウスカーネルと線形項と定数項を加えたような下記の形で計算してみましょう。
$$
k(\boldsymbol{x}_n,\boldsymbol{x}_m)=\theta_0\exp\left(-\frac{\theta_1}{2}||\boldsymbol{x}_n-\boldsymbol{x}_m||^2\right)+\theta_2+\theta_3\boldsymbol{x}_n^T\boldsymbol{x}_m
$$
各パラメータの勾配は下記のようになります。
これらの勾配がゼロになり、尤度が最大化するようなするようなパラメータが最尤推定されたパラメータとなります。
$$
\frac{\partial\boldsymbol{C}_N}{\partial\beta}=-\frac{1}{\beta^2}
$$
$$
\frac{\partial\boldsymbol{C}_N}{\partial\theta_0}=\exp\left(-\frac{\theta_1}{2}||\boldsymbol{x}_n-\boldsymbol{x}_m||^2\right)
$$
$$
\frac{\partial\boldsymbol{C}_N}{\partial\theta_1}=-\frac{\theta_0}{2}||\boldsymbol{x}_n-\boldsymbol{x}_m||^2\exp\left(-\frac{\theta_1}{2}||\boldsymbol{x}_n-\boldsymbol{x}_m||^2\right)
$$
$$
\frac{\partial\boldsymbol{C}_N}{\partial\theta_2}=1
$$
$$
\frac{\partial\boldsymbol{C}_N}{\partial\theta_3}=\boldsymbol{x}_n^T\boldsymbol{x}_m
$$
まとめ
線形回帰の$w$というパラメータにガウス分布を導入する流れから、カーネルを一般化し、ガウス過程回帰を理解してきました。パラメトリックモデルを経由せず、関数に対して直接的に分布を定義していくように見えます。
要はデータが手に入って、その共分散行列さえ計算できれば、多彩なカーネルで回帰できちゃうぞ、ということです。
次回は実装編です。