はじめに

本記事はMLPシリーズ本「ガウス過程と機械学習」を読んで私なりに理解したガウス過程回帰とその御利益を説明します。
（本や文献には”関数の確率分布”とか”無限次元のガウス分布”などと説明がされていますが、自分にはあまりしっくりこず、今回実際に実装してみることで理解が進みましたので、その感覚で本記事を書こうと思いました。ですので、専門家からするとガウス過程の直接的な説明でなかったり、誤解を与えるような表現もあるかもしれませんが、その点ご容赦ください。）

参考図書
ガウス過程と機械学習

本記事で作成したコード

#想定している問題
既知のデータ$\boldsymbol{D} \{ (\boldsymbol{x_1},y_1),(\boldsymbol{x_2},y_2),\cdots ,(\boldsymbol{x_N},y_N)\}$を持っているとします。このとき、$\boldsymbol{x}$は入力、$y$は出力です。
このとき、新しい入力$\boldsymbol{x^*}$に対する$y^*$を予想するという問題が今回の想定です。

例えば、スイカの肥料の量、水の量、ビニールハウスの室温を入力$\boldsymbol{x}$として、半年後のスイカの重量を$y$とします。過去の経験から、ある条件$\boldsymbol{x^*}$のときに半年後のスイカの重量はどれくらいになるか予想するといった感じです。
さらに、半年後のスイカの重量を最大化するためにはどのような条件$\boldsymbol{x}$が良いか確率的に求める問題がベイズ最適化です。（本記事に関連させて別途投稿したいと思います。）

ガウス過程回帰の数式表現

(数式は参考図書と同じ表現です。太字の場合はベクトルもしくは行列です。）　　
　　

既知データの集合$\boldsymbol{D} \{ (\boldsymbol{x_1},y_1),(\boldsymbol{x_2},y_2),\cdots ,(\boldsymbol{x_N},y_N)\}$があるとします。
このとき、ガウス過程では、入力$\boldsymbol{x}$に対する出力$\boldsymbol{y}$が多変量ガウス分布に従うとします。
数式で書くと、

\boldsymbol{y} \sim \mathcal{N}(\boldsymbol{0},\boldsymbol{K})

確率密度関数は、

\begin{bmatrix}
y_1 \\
\vdots\\
 y_N\\
\end{bmatrix}
=　\frac{1}{(\sqrt{2\pi})^N \sqrt{|\boldsymbol{K}|}}
exp \biggl(
-\frac{1}{2} \boldsymbol{x}^T \boldsymbol{K}^{-1} \boldsymbol{x}
\biggr)

です。

つまり、出力$\boldsymbol{y}$は平均$\boldsymbol{0}$の周辺に分散$\boldsymbol{K}$で分布する多変量ガウスモデルに従うと事前に想定し、実際に得られた出力$\boldsymbol{y}$はこの確率モデルの実現値だと解釈します。

このもとで、新しいデータ $\boldsymbol{x}^*$が与えられたとき、$y^*$ の予測分布は、以下の通りに計算できます。

p(y^*|\boldsymbol{x}^*,\boldsymbol{D})=\mathcal{N}(\boldsymbol{k}_*^T \boldsymbol{K}^{-1} \boldsymbol{y},\quad k_{**}-\boldsymbol{k}_*^T \boldsymbol{K}^{-1} \boldsymbol{k}_*)\\

つまり、$y^*$ は、平均 $\boldsymbol{k}_*^T \boldsymbol{K}^{-1} \boldsymbol{y}$、分散　$k_{**}-\boldsymbol{k}_*^T \boldsymbol{K}^{-1} \boldsymbol{k}_*$　のガウス分布ということです。
　　
　　
ここに出てきた変数の意味は、
　　
$\qquad\boldsymbol{K}$ ：カーネル行列 (N x N次元)

\boldsymbol{K} = \begin{bmatrix}
k(\boldsymbol{x_1},\boldsymbol{x_1}) & \cdots & k(\boldsymbol{x_1},\boldsymbol{x_N})\\
\vdots & & \vdots \\
k(\boldsymbol{x_N},\boldsymbol{x_1}) & \cdots & k(\boldsymbol{x_N},\boldsymbol{x_N})
\end{bmatrix}

$\qquad k$: カーネル関数(スカラ)

いろいろなカーネルがありますが、ガウスカーネルの場合は以下の通りです。\\
k(\boldsymbol{x},\boldsymbol{x'})=\theta_1exp \biggl(
-\frac{|\boldsymbol{x}-\boldsymbol{x'}|^2}{\theta_2}
\biggr)\\
ここで、\theta_1と\theta_2はパラメータです。

$\qquad \boldsymbol{k}^*$: (N次元ベクトル)

\boldsymbol{k^*} = \begin{bmatrix}
k(\boldsymbol{x_1},\boldsymbol{x^*})\\
\vdots\\
k(\boldsymbol{x_N},\boldsymbol{x^*})
\end{bmatrix}

$\qquad k^{**}$: (スカラ)

k^{**}=k(\boldsymbol{x^*},\boldsymbol{x^*})

数値実験

適当な１次元の関数を作り上記の式を計算してみて、ガウス過程回帰がどのような挙動になるか見てみます。
（１次元ですので、この場合データ$\boldsymbol{x}$はスカラ$x$です。）

適当な関数を用意

y=e^{x}sin(2\pi x)

データ準備

既知データ $D$ として３つの点を用意しました。

\begin{align}
(x_1, y_1) &= (0.79,-2.14)\\
(x_2, y_2) &= (-0.34,-0.61)\\
(x_3, y_3) &= (0.64,-1.48)\\
\end{align}

##カーネル行列を計算する

\begin{align}
\boldsymbol{K} &= \begin{bmatrix}
k(x_1,x_1) & k(x_1,x_2) & k(x_1,x_3)\\
k(x_2,x_1) & k(x_2,x_2) & k(x_2,x_3)\\
k(x_3,x_1) & k(x_3,x_2) & k(x_3,x_3)\\
\end{bmatrix}\\
\\
&=
\begin{bmatrix}
1.00000 & 0.00000 & 0.66333\\
0.00000 & 1.00000 & 0.00000\\
0.66333 & 0.00000 & 1.00000]\\
\end{bmatrix}\\
\end{align}

　　　
上述したように、カーネル行列$K$は$\boldsymbol{y}$の確率モデルである多変量ガウス分布の共分散行列で用いられます。したがって、カーネル行列の要素$(i,j)$の数値が大きいということは、$y_iとy_j$の相関が大きいということです。
さらに、カーネル行列の要素$(i,j)$は、上述したとおり、（ガウスカーネルを用いる場合）

k(x_i, x_j)=\theta_1exp \bigl(
-\frac{|x_i-x_j|^2}{\theta_2}\bigr)

です。この関数は、$x_i$と$x_j$が近い時に大きな値になります。

以上をまとめると、
　　

x_iとx_jが近い\\
↓\\
k(x_i, x_j)が大きい\\
↓\\
多次元ガウス分布の共分散行列の(i,j)要素が大きい\\
↓\\
y_iとy_jが近い

　　
　　
先ほど実際に計算したカーネル行列を見ると,　$k(x_1,x_2)=0$です。これは$x_1$と$x_2$が離れているためです。一方で、$k(x_1,x_3)=0.66$です。これは$x_1$と$x_3$が近いためです。
繰り返しになりますが、$\boldsymbol{y}$は$\boldsymbol{K}$を共分散行列とした多変量ガウス分布に従うため、以上の結果から、$y_1$と$y_2$は離れた値になるだろう、$y_1$と$y_3$は近しい値になるだろうと予想していることになります。

回帰する場合、つまり、既知データ$(x_1,y_1),(x_2,y_2),\cdots$を持っていて、新しいデータ$x^*$に対する$y^*$を予想する場合、ガウス過程回帰は$x^*$と近い$x_i$に対する$y_i$周辺を予想します。

回帰をしてみる

$x^*=0.4$に対してyを予想してみます。
再び、既知のデータ$D$のもとで、新しいデータ$x^*$に対する$y^*$の分布は、以下の通り表される。

p(y^*|\boldsymbol{x}^*,\boldsymbol{D})=\mathcal{N}(\boldsymbol{k}_*^T \boldsymbol{K}^{-1} \boldsymbol{y},\quad k_{**}-\boldsymbol{k}_*^T \boldsymbol{K}^{-1} \boldsymbol{k}_*)\\

したがって、$y^*$の期待値 ($\boldsymbol{k}_*^T \boldsymbol{K}^{-1} \boldsymbol{y}$)は、

\begin{align}
\boldsymbol{k}_*^T \boldsymbol{K}^{-1} \boldsymbol{y}
&=
\begin{bmatrix}
0.051 & 0.000 & 0.308\\
\end{bmatrix}
\begin{bmatrix}
1.00000 & 0.00000 & 0.66333\\
0.00000 & 1.00000 & 0.00000\\
0.66333 & 0.00000 & 1.00000]\\
\end{bmatrix}^{-1}
\begin{bmatrix}
-2.14\\
-0.61\\
-1.48\\
\end{bmatrix}\\
&=-0.141
\end{align}

です。

さらに、分散　$k_{**}-\boldsymbol{k}_*^T \boldsymbol{K}^{-1} \boldsymbol{k}_*$　は、

\begin{align}
k_{**}-\boldsymbol{k}_*^T \boldsymbol{K}^{-1} \boldsymbol{k}_*
&=1.0-
\begin{bmatrix}
0.051 & 0.000 & 0.308\\
\end{bmatrix}
\begin{bmatrix}
1.00000 & 0.00000 & 0.66333\\
0.00000 & 1.00000 & 0.00000\\
0.66333 & 0.00000 & 1.00000]\\
\end{bmatrix}^{-1}
\begin{bmatrix}
0.051\\
0.000\\
0.308\\
\end{bmatrix}\\
&=0.863
\end{align}

です。

以上より、$y^*$は、$-0.141\pm \sqrt{0.863}=-1.069 \sim 0.787$あたりと予想されます。

さらに、同じことを繰り返し、全ての$x$に対する$y$を予想すると以下のようになります。
データがある部分は不確かさ０で予想しています。一方で、データが少ない領域は予想の不確かさが大きいです。
　　

ガウス過程回帰では、$\boldsymbol{x}$ に対する $y$の分布形状を何も仮定していません。つまり、$y$がどのような形状でもきちんとその分布を表現し予測することができます。このように$y$の分布を仮定しない方法をノンパラメトリックといい、ガウス過程回帰の優れている点です。

次の投稿では、$y$を最適化するための最適解$\boldsymbol{x}$ を求める方法としてベイズ最適化について書こうと思います。