ガウス過程回帰を考え方から導出してみた


今回は図解がまだ準備できてないが、図つきの記事は他の記事で色々書いてあるので、今回は数式ゴリゴリで導出する記事にしてみた。図解系は解説は後日少しは書き足すつもり。

線形代数と直和とかの概念が理解できていて、PCAの導出、ベイズ推定の記事を読める人が見れば、ある程度導出の筋道が分かりやすくなっているのではと思う。

別記事に、やはり単回帰と重回帰系、カーネル回帰系の記事は書いた方がいいかもしれない。今回は考え方の基本からガウス過程回帰を個人的に導出したかったので先にガウス過程回帰を記載したが、次に書こうと思う。


## 確率過程とガウス過程

確率過程とガウス過程について、ざっくり要約すると以下になる。

- ベクトルをランダムサンプリングするときに、添字$i$同士に相関をもたせてサンプリングした場合は、iを横軸とした際の各成分はランダムにカーブっぽい形状を描く。こういうサンプルを**確率過程**とよぶ。
- ベクトルの添字$i$の代わりに連続値$x$をもちいた無限次元ベクトルV(x)は関数と等価である。このことから、確率過程は関数をランダムサンプリングする過程だ、ともいえる。
- 有名な確率過程にマルコフ過程がある。これは現在の状態が直前の状態にのみ依存して決まる、という過程である。これもiに依存したサンプリングになっている。
$$p(y_i|i)=p(y_i|y_{i-1})$$
- 今回回帰に使う**ガウス過程**とは、yの分布形状が正規分布となる過程だといえる。
$$p(y|x)=\mathcal N(y|m(x), v(x,x))$$
正規分布のパラメータは平均$m$, 分散$v$であり、これらはxに依存して決まる、と定義すればガウス過程の式は立てられる。(ここで、$m$と$v$はx自身の平均や分散をとる訳ではなく、解こうとする問題に応じて決まることに注意する).
- **ガウス過程回帰**とは、上記ガウス過程の式を回帰問題にあてはめて、$m,v$の関係式を導くことでy(x)を予測しよう、という試みである。
- ガウス過程によって求まるyは正規分布として求まるので、誤差範囲も込みで予測できる。
- 回帰した分布($(x'_i,y'_i)$が既知のときの**事後分布**)は、全く既知のデータがないときの分布(**事前分布**)を仮定しないと求まらない。事前分布は
$$m(x)=0, v(x_i,x_j)=\exp(-(x_i-x_j)^2)$$
を使うことがよくある。

さて、これを踏まえてガウス過程回帰の着想から導出までを完走しよう。

## ガウス過程回帰の着想と確率過程の利用

まずは思想を簡単にするために、
x-yの二次元平面上の回帰問題として考える。
いま、
$$(x_1,y_1),...,(x_n,y_n)$$
のn個のデータが得られた前提を考える。

このとき、新しいデータ点$x_{n+1}$に対して$y_{n+1}$はどんな値をとるのが確からしいのだろうか？
これを推定する方法はいくつかあるが、ガウス過程回帰では(これは厳密には正しい表現ではないが)、$y_1,y_2,y_3,...$のすべてについて、バラつきが各々$x_i$に依存した正規分布に従うと考える。

いま、$$i=\in [1,N]$$となる一つの$i$について着目しよう。$y_i$は$x_i$に依存することを前提とするので、$y_i$の分布はxに関する適当な関数$m$,$k$を用いて

$$p(y_i)=\mathcal N(y_i|m(x_i), k(x_i))$$

と書ける。y_iはここからサンプリングされることで得られる、と考える。

いま、$y=x$を真の傾向とする母集団について回帰することを考えよう。$$$x_j=1+\frac{j}{100}, j\in[1,100]$$の100点についてサンプリングしてみる。
すると、たとえ有利な条件(m(x)=x, k(x)=1)を何らかの手法で求めたとしても、とても回帰できたとは言えない結果がでてしまう。


これは、2点$x_i,x_j$のサンプル間の相関を考慮してないことが原因である。

yはxの値に依存した関数の形になるような答えを求めようとしている。また、xの値が似ているならばyの値も似るような関数で回帰したいのである。それを考慮すると$x_i,x_j$は2点のxの距離に応じて同じ値がでるようなサンプリングが必要である。

相関を議論するは**共分散**を考慮すればよい。


### 共分散おさらい

ここまでの話を一旦忘れて共分散について議論しよう。
値が一つのみで構成されるサンプル$x_i, i\in [1,N]$についての平均と分散(標本分散)は以下で表される.

$$\bar x=\frac{1}{N}\sum_i x_i$$
$$V_x=\frac{1}{N}\sum_i(x_i-\bar x)^2$$

値が複数ある場合(1サンプルをベクトルXとしよう)、複数の値の間で相関がある場合も考えて、共分散という考え方が出てきたのだった。

$x_{ji}$をi番目のサンプル$X$に対するj成分だとすると、平均と分散、共分散は以下で計算できる。

$$\bar x_j=\frac{1}{N}\sum_i x_{ji}$$
$$V_{jk}=\frac{1}{N}\sum_i (x_{ji}-\bar x_j)(x_{ki}-\bar x_k)$$

この$V$のことを**分散共分散行列**というのだった。これは対角成分が分散であり、非対角成分は$j,k$の間の相関を表していた。
(より詳しくは主成分分析の記事を参照)

### 確率過程: Xの値に相関するYの分布

さて、冒頭の話を思い出してみよう。
今回ちょっと捻っているのは、サンプル点一つ一つが独立の特徴と考えていることだ。
「一連のサンプルY」というベクトルYはy=f(x)という関数のもと、「一連のサンプルX」というベクトルXをランダムサンプルした結果得られたもの、と考える。

したがって、$(X_i,X_j)$が同じもしくは近い値であるほど$(Y_i,Y_j)$が相関しているべきだし、逆に$(X_i,X_j)$が遠いならば$(Y_i,Y_j)$は独立していてよい、ということになる。

つまり、Yの分散共分散行列の各成分$v(x_i,x_j)$は$x_i$と$x_j$の距離に応じて決まる非負値であり、$x_i,x_j$の距離が大きいほど0に漸近するような関数だと良い。これを**カーネル関数**、$k(x_i,x_j)$として、

$$p(Y)=\mathcal N(Y|m(X), K(X,X))$$

ただし、
$$K(X,X)_{ij}=k(x_i,x_j)$$

として定義できる。(Kを**グラム行列**という)
mとkは適当に決めたものを**事前分布**として定義する。
あとは、データ$X^\prime$に対して$Y^\prime$が得られた条件で新しい$X$に対する事後分布を得れば対応する$Y$が得られる、という算段である。

ベクトルYの構成サンプル次元が十分に大きかった場合、任意のXに対するYが、**Xの値に依存する形で**定義できるようになるだろう。そういう意味で

$$p(Y)=\mathcal(Y|m(X),K(X,X))$$

はXに関する関数をサンプリングしていると言い換えることができる。このような過程を**確率過程**と呼び、今回のように分布関数が正規分布のときには特に**ガウス過程**と呼ぶ。


### 問題のまとめと、ガウス過程回帰のゴール

いま、現時点で得られたデータ点$(x^\prime_i,y^\prime_i)$からyをxの値に依存する前提として、新しい$x_j$に対する$y_j$を求める回帰を行いたいのだった。

複数のxに依存した複数のyの値をサンプリングするのはガウス過程の式

$$p(Y|X)=\mathcal N(Y|m(X), K(X,X))$$

で定義できる。事前分布として$m$は適当に決めることとして、$k$は以下の条件を満たすような式の中から適当に決めればいいだろう。

$$k(x_i,x_j)\geq 0$$
$$\lim_{|x_i-x_j|\rightarrow\infty} k(x_i,x_j)=0$$
$$k(x_i,x_j)=k(x_j,x_i)$$

なお、事前分布の例として
$$m(X)=0$$
$$k(x_i,x_j)=\exp(-(x_i-x_j)^2)$$
を選ぶことはよく見る。
今回、汎用性を考えて事前分布を特定の形に仮定せず、$m(x),k(x_i,x_j)$のまま完走しよう。

いま、既知のサンプル群$(X',Y')$を用いて、Xが与えられているがYが未知のサンプルの$(X,Y)$を推定することを考える。
これは、条件付き分布$p(Y|X,X',Y')$を求めることと同じである。
これを計算するのは同時分布を用いて

$$p(Y|X,X',Y')=\frac{p(Y,Y'|X,X')}{p(Y'|X')}$$
とかける。
ただし、$p(Y|X,X')=p(Y|X), p(Y'|X,X')=p(Y'|X')$すなわち、真のYは過去に選ばれたX'に依存しないし過去のY'は未来のXに依存しないものとする。

右辺分子の同時分布は$(X',Y')$と$(X,Y)$のベクトルをつなげて(直和して)算出すればよい。
$$\hat X=X'\oplus X$$
$$\hat Y=Y'\oplus Y$$
とおくと、

$$p(Y,Y'|X,X')=p(\hat Y|\hat X)$$
$$p(\hat Y|\hat X)=\mathcal N(\hat Y|m(\hat X), K(\hat X,\hat X))$$

となる。ゆえに、

$$p(Y|Y',\hat X)=\frac{p(\hat Y|\hat X)}{p(Y'|X')}$$

これが、
$$p(Y|Y',\hat X)=\mathcal N(Y|m_\ast(X), K_\ast(X,X))$$
なる形に変換できれば、(Y',X')をデータとして得たときの(X,Y)分布、つまりばらつき込みの回帰予測が求まったことになる。
このような$m_\ast,K_\ast$を求めることが、ここから先の変形の目的である。

## ガウス過程回帰の導出

式を組み立てる材料は揃った。あとは上で求めた条件付き分布

$$p(Y|Y',\hat X)=\frac{p(\hat Y|\hat X)}{p(Y'|X')}$$

の右辺を一旦expを使った式で直接計算して係数比較できる形に変形すれば導出完了となる。分子分母をexpに置き換えて両辺にlogをとると、

$$log(p(Y|Y',\hat X))=-(\hat Y-\hat m)^T\hat K^{-1}(\hat Y-\hat m)+(Y'-m')^TK'^{-1}(Y'-m')+Const.$$

となる。ただし、
$$\hat m=m(\hat X), m'=m(X')$$
$$\hat K=K(\hat X,\hat X), K'=K(X',X')$$
である。

目的の式を得るには、この式の右辺から$\hat X,\hat Y$を$X,Y,X',Y'$に戻して、$(Y-m)$について平方完成すればよい。
そのため、確率過程の定義から$\hat m,\hat K$を戻そう。

$$m(\hat X)=m(X')\oplus m(X)$$
$$K(\hat X,\hat X)=K(X'\oplus X, X'\oplus X)$$

$K(X,X)=k(x_i,x_j)$であり、
かつ$\hat X$の成分は番号の若い方からX'、Xの順につなげて定義するのだった。
つまり、これを$K(X',X'),K(X',X),K(X,X'),K(X,X)$
を使って置き換えると、以下のブロック行列
$$K(\hat X,\hat X)=\begin{pmatrix}
K(X',X') & K(X',X) \\
K(X,X') & K(X,X)
\end{pmatrix}$$
となる。

同様に、一旦
$$K(\hat X,\hat X)^{-1}=\hat G$$
とおくと、

$$\hat G=\begin{pmatrix}
G_a&G_b\\
G_c&G_d
\end{pmatrix}$$

とでき、$G_a,G_b,G_c,G_d$の行列形状は$K(X',X'),K(X',X),K(X,X'),K(X,X)$
の形状と対応させた定義ができる(逆行列をとることでは行列形状は変わらないので、こうできることは自明である)。
ただし、$G_\ast$の**各々の値は各逆行列に対応している訳ではない点に注意**する(そもそも$G_b,G_c$に至っては長方形行列の場合すらある)。この辺の対応関係はあとで議論しよう。
なお、$\hat K$は対称行列なので逆行列である$\hat G$も対称行列であり、ゆえに$G_b=G_c^T$であることや、$G_a=G_a^T,G_d=G_d^T$であることは、以降の変形で利用する。

ここまでの変形を踏まえてlog(p)の式の変形を試みる。(以降、前提データ(X',Y',X)をまとめてDとおく)

$$\log p(Y|D)=-(Y'\oplus Y-m'\oplus m)^T\hat G(Y'\oplus Y-m' \oplus m)+(Y'-m')^T[K(X',X')]^{-1}(Y'-m')+Const.$$

$$\begin{eqnarray}
=&-(Y'-m')^TG_a(Y'-m')\\
&-(Y'-m')^TG_b(Y-m)\\
&-(Y-m)^TG_b^T(Y'-m')\\
&-(Y-m)^TG_d(Y-m)\\
&+(Y'-m')^T[K(X',X')]^{-1}(Y'-m')\\
&+Const.
\end{eqnarray}
$$

ここで、
$$-(Y-m_\ast)^TG_\ast(Y-m_\ast)+...$$

で、省略した先には$Y$は含まれないようにする。展開すると、

$$-(Y-m)^TG_\ast(Y-m)-(Y-m)^TG_\ast(m-m\ast)-(m-m_\ast)^TG_\ast(Y-m)+...$$

となる(省略した部分はYを含まない)。係数比較を行うと、

$$G_\ast=G_d$$
$$G_\ast(m-m_\ast)=G_b^T(Y'-m')$$
$$(m-m_\ast)^TG_\ast=(Y'-m')^TG_b$$
となる。$G_\ast, m_\ast$に関する式が3本になるので、一見解けない可能性もあるが、3番目の式は両辺転置とると

$$G_d^T(m-m_\ast)=G_b^T(Y'-m')$$

となり、$G_d^T=G_d$により二番目の式と一致する。つまり解は存在するということになる。

二番目と三番目の式は$m_\ast$について解くと、
$$m_\ast=m-G_d^{-1}G_b^T(Y'-m')$$
となる。

まとめると、

$$K_\ast=G_d^{-1}$$
$$m_\ast=m-G_d^{-1}G_b^T(Y'-m')$$

で定義される

$$p(Y|D)=\mathcal N(Y|m_\ast,K_\ast)$$

がガウス過程回帰の結果となる。

### ブロック行列の変形

最後に、SをKの式で表現できればゴールである。

これは、2x2のブロック行列の逆行列が(逆行列とる部分が正則ならば)以下のようになることを利用する。

$$\begin{pmatrix}
A & B \\
C & D
\end{pmatrix}^{-1}=
\begin{pmatrix}
S_D^{-1} & -A^{-1}BS_A^{-1} \\
-S_A^{-1}CA^{-1} S_A^{-1}
\end{pmatrix}
$$
ただし、
$$S_A=D-CA^{-1}B,\quad S_D=A-BD^{-1}C$$

をSchur補行列と呼ぶ。

### 最後の変形

さきほどの
$$\hat G=\hat K^{-1}$$
であることを思い出せば、
$$G_a^{-1}=K_a-K_bK_d^{-1}K_b^T$$
$$G_d^{-1}=K_d-K_b^TK_a^{-1}K_b$$
$$G_b^T=-G_dK_b^TK_a^{-1}$$

となる。ただし、

$$\hat K=\begin{pmatrix}
K(X',X') & K(X,X') \\
K(X',X) & K(X,X)
\end{pmatrix}$$
$$=\begin{pmatrix}
K_a & K_b\\
K_b^T & K_d
\end{pmatrix}
$$

これをガウス過程回帰の各係数の式に適用すると、

$$K_\ast=G_d^{-1}=K_d-K_b^TK_a^{-1}K_b$$
$$m_\ast=m+K_b^TK_a^{-1}(Y'-m')$$

となり、

$$K_\ast=K(X,X)-K(X,X')[K(X',X')]^{-1}K(X',X)$$
$$m_\ast=m(X)+K(X,X')[K(X',X')]^{-1}(Y'-m')$$

と導ける。上記が新しいXに対応する平均と分散なので、これで導出は完走である。


## Yの誤差とノイズの設定

さて、ガウス過程の確率分布

$$p(Y)=\mathcal N(Y|m(X), K(X,X))$$
について考察を始めよう。まず、この式のままだと過学習が懸念される。
なぜなら、$$x_i=x_j,i\neq j$$な$x_j$に対して、$$k(x_i,x_i)=k(x_i,x_j)$$となるからである。つまり、$(x_i,x_j)$の共分散が分散と一致し、$y_i,y_j$は完全に相関する。つまり同じ値となることを要請するからである。

現実問題として、$y_i$は真値からの測定誤差があるだろうから、同じ$x_i$を繰り返し測定したところで、同じ$y_i$はすこしずつ違った値が出て当然である。**確率過程**は一本の関数$f(x)$をランダムサンプリングするので、関数が決まってからの繰り返し誤差を許容せず$x_i=x_j$なら$f(x_i)=f(x_j)$になってしまう。

繰り返し誤差を許容するには、

$$y=f(x)+\epsilon$$
とノイズを混入させればいい。ノイズが
$$p(\epsilon)=\mathcal N(\epsilon y|0, \lambda^2)$$
に従う(つまり反復誤差が\labmda)であるとき、yの分布は

$$p(Y)=\mathcal N(Y|m(X),K(X,X)+\lambda^2I)$$

となる。単位行列で定義するのは、繰り返し誤差は試行回数に無相関であるため共分散に影響を与えないためである。

一見導出式が変わってしまうように見えるが、導出の過程で
$$k(x_i,x_j)=k(x_i;i,x_j;j)$$
と置き直してもグラム行列$K$の対称性は失わないため、更新式に現れるKについて、既知のデータ成分について
$$k(x^\prime_ii,x^\prime_jj)=k(x^\prime_i,x^\prime_j)+\lambda^2\delta_{ij}$$
と置き換えてやればよい。

なお、上記例は既知のデータのみ誤差を与えたが、求める方の分布も誤差を与えてやってもよい。更新式はその場合、事後分布の分散にシンプルに誤差が加わる結果になる：

$$K_\ast=K(X,X)+\lambda^2 I - \cdots$$

これは、未知の$X$を与えたときに**測定誤差込みで**$Y$の分布を予測した結果になる。通常、回帰で予測したいのは$Y$の**真値の存在範囲**なので、そういう意味では予測対象には誤差を考慮しないほうが使いやすいといえるだろう。

なお、測定誤差成分はべつに単位行列である必要もなく、対角行列であれば係数の違いは問わない。もし$i$ごとに異なる係数を割り当てた場合は、係数が大きい成分ほど$y^\prime_i$に対する精度が低いといえる。
$$p(Y)=\mathcal N(y|m(X),K(X,X)+\Lambda)$$

$$\Lambda_{ij}=\lambda_i^2\delta_{ij}$$

とするとき、$\lambda_i$は各々の$y_i$の測定誤差に対応する。

## 事前分布の意味と平均関数の設定

ところで、事前分布のパラメータ関数

$$
m(X)=0$$
$$
K(X,X)=\exp(-\frac{(x_i-x_j)^2}{\beta^2})
$$

は何を意味しているのか？

- yは平均0分散1である
- xが近いほどyも同じ値をとる(相関係数が$1/e$となるときの距離を半径と呼ぼう)
- この半径は$\beta$である

となる。いまデータ$x_i,y_i$が得られると事後分布の平均と分散は更新されるが、$x_i$が遠いところの$x$付近は相関せず、相変わらず事前分布のまま平均0分散1が期待される。
要は事前分布平均0での回帰は「得られたデータが遠い点は0と見なしましょう」と言っているわけだ。
逆に、データの傾向的にあるべき姿がある程度予測ついている場合、事前分布はその予測する関数を使った方が「データがない点」への接続の確からしさはあがる、ともいえる。

たとえば、根本的には$y=x$に従うことが予想される中で、一部得られたデータで関数を微調節したい、という場合は事前分布として、
$$
m(X)=X
$$
とした方が都合のいい場合もある。だいたいの
場合は$y=ax+b$くらいしか形状予測ができないので、aとbがわかっている例というのは多くないが。
なお一応、aやbが未知の状態を事前分布として求めることも可能なはずである。これについては別記事で議論しよう。

## 事前分布の分散について

さて、事前分布の分散はもう少し拡張したら、

$$K(X,X)+\Lambda_{ij}=\kappa_{i}\kappa_{j}\exp(-\frac{(x_i-x_j)^2}{\beta_i\beta_j})+\lambda_{ij}$$

としてもカーネルの条件を満たすので、拡張として成り立つはずだ。では、この$\kappa$や$\beta$は何を意味するか？分散を増やすという効果なので正則化に似ていそうだが、共分散も変化する点に注意する。

これらに解釈を与えるとしたら、下記だと考える(もしかするとどこかの教材に書いてあるかもしれないが)

- $\beta_i$はデータ点$(x_i,y_i)$のX半径。そのデータ点が$x_i$の周辺どこまで影響するかを示した半径となる、$\beta_i,\beta_j$をそれぞれ掛けると、分散への影響は各点で$\beta_i^2$という別々の値がかかるので、わかりやすい。
- $\kappa_i$はデータ点$(x_i,y_i)$のY半径の扱いになる。ここで$\Lambda_ij$もy方向の誤差に相当するので紛らわしいが、後者は**偶然誤差**であり、似たサンプルをいくつか打てば真値の誤差は狭まっていく。一方、$\kappa_i$は、その場のサンプリングを何回増やしても誤差が小さくならないタイプの誤差である。得られたデータごとに異なる値は設定できるので、データごとに測定に使った方法自体にどのくらいのY方向誤差がのりそうか、を加味する際に利用する。
- $\beta_i$はX方向の誤差によって決まる値ともいえるが、どちらかというとデータ点$X_i$をどのくらいのx範囲の予測に利用するか、を定義していると考えた方がしっくりくる。

色々考察は可能だが長くなりそうなので、ここから先は気が向いたら別記事にしようとおもう。