「ガウス過程の機械学習」のノートと誤植訂正(3章)

$3$章からガウス過程に入り，かなり難しくなる印象です。数式も強そうな見た目になってきます。丁寧に読まないと混乱する場所が多そうです。誤植は付録のところくらいなようです。(素人目には)議論に少し怪しいところがあるような気もしますが，それも無視できるレベルです。ここでの記述に不正確なところなどあればご指摘ください(だんだん理解できているのかわからなくなってきた)。

#3.1 線形回帰モデルと次元の呪い
この章の初めの方では「次元」という言葉が二つの意味で用いられているので，ここで明確に区別しておこうと思います。

- 空間の次元$\cdots$$1$つのデータ$\boldsymbol{x}$の成分数のことです。数式を使うと，$\boldsymbol{x} \in \mathbb{R}^n$の$n$のことです。
- データ数の次元$\cdots$単にデータ数のことです。後で出てくる「無限次元のガウス分布」の無限次元はこちらのことを指します。与えられるデータ数は有限ですが，そこから無限のデータ点に対する出力を得ます(つまりこれは関数を得るということです)。

この節では線形回帰モデルの基底関数として動径基底関数
$$
\phi_h (x) = \exp \left( - \dfrac{ (\boldsymbol{x} - \boldsymbol{\mu}_h )^2} { \sigma^2} \right)
$$
を選ぶことを考えています。しかしこのとき空間の次元が大きいとき(つまり$\boldsymbol{x} \in R^n$として$n \gg 1$のとき)，必要なパラメータ数は刻み幅を$H$とすると$H^n$のように増加してしまいます。これを指して次元の呪いと言っていますが，厳密にはここでの次元は空間の次元ですね。本書でも次元に対して指数的にパラメータが増加すると述べられています。

#3.2　ガウス過程
この節では線形回帰モデルのパラメータ$\boldsymbol{w}$に確率分布を導入し，得られる出力$\boldsymbol{y}$の確率分布を導きます。
テキストでは厳密に線形回帰できるときを考えていますが、誤差を考慮に入れても同じように考えることができます。というか後の方を見るとそれについても扱っていますが，同じ形で導出した方がわかりやすいと思うので一応示しておきます。誤差は平均$0$で分散$\sigma$の正規分布に従うものとします。

```math
\begin{eqnarray}
\langle \hat{\boldsymbol{y}} \rangle  &=&  \langle \Phi \boldsymbol{w} + \boldsymbol{\epsilon} \rangle  = 0 \\
\Sigma &=& \langle \boldsymbol{y} \boldsymbol{y} ^T \rangle -\langle \boldsymbol{y}\rangle \langle \boldsymbol{y} ^T \rangle   \\
 &=& \Phi \langle \boldsymbol{w} \boldsymbol{w}^T \rangle \Phi ^T + \Phi \langle \boldsymbol{w} \boldsymbol{\epsilon}^T \rangle
+ \langle \boldsymbol{\epsilon} \boldsymbol{w}^T \rangle \Phi^T
+ \langle \boldsymbol{\epsilon} \boldsymbol{\epsilon}^T \rangle \\
 &=& \lambda \Phi \Phi^T + \sigma^2
\end{eqnarray}
```
縦ベクトル$\times$横ベクトルは行列になることに注意します。また$ \langle \boldsymbol{w} \boldsymbol{\epsilon}^T \rangle = 0$が成立するのは、これらの確率変数が独立ということを利用しています。
結局のところ、共分散行列の対角成分に余計な項が加わっただけという感じです。つまりあまり気にする必要はなかったということです。

パラメータ$\boldsymbol{w}$が消えてしまって困惑するかもしれません。通常の機械学習では$\boldsymbol{w}$をフィッティングすることが重要です。$\boldsymbol{w}$が消去されてしまったことによって，出力$\hat{\boldsymbol{y}}$を決定づけるのは共分散行列$K = \lambda^2 \Phi \Phi^T$のみになります。そしてこの$K$には観測値$\boldsymbol{y}$の情報は全く含まれておらず，得られる出力$\hat{\boldsymbol{y}}$にはまったく意味がないように思われます。
ここで重要になるのが前回扱った，ガウス分布に対する条件付確率です。つまり観測値$\boldsymbol{y}$を既知の条件とすることで，未知のデータ$\boldsymbol{x}^*$に対する予測値$\hat{\boldsymbol{y}^ *}$の確率分布が変化します。このような予測をたくさん(理想的には無限大のデータ数)することによって，関数形を得ようというのがガウス過程による回帰の大まかな手続きです。最後の予測には既知の観測値を用いているため，得られているデータに沿うような形になることが期待されます。そして$K$の選び方によって様々な特徴を持つということをこの後で見ます。関数形を得る手続きについては少し後回しです。

### 3.2.1　ガウス過程の意味
与えられた$2$つのデータ点があったときに、これらの特徴を表す特徴量空間への写像$\phi$を考えて、特徴量空間での$2$つの"近さ"というものを測ります。それが近いとき、相関が強くなり似たような出力をするということです。
ここを読むと、確かに自然言語処理などと相性がいいのかもしれないと感じました。自然言語はそのままだと"距離"などを考えるのは難しく，特徴量空間に飛ばしてしまって考えるということですかね。

### 3.2.2 カーネルトリック
特徴ベクトル$\boldsymbol{\phi}$を用いて，カーネル関数は以下のように定義できます：

$$
k(\boldsymbol{x} _n ,\boldsymbol{x} _n' ) = \boldsymbol{\phi} (\boldsymbol{x} _n)^T \boldsymbol{\phi} (\boldsymbol{x} _n')
$$
そして特徴ベクトルのことは無視して，カーネル関数だけを適切に選ぶことで計算の手間を減らすというのがカーネルトリックです。ここから得られる共分散行列(カーネル行列)には対称行列であること、逆行列の存在と正定値行列であるということが要請されます。
カーネル関数とカーネル行列の関係はどのようになっているんでしょうか。$k(\boldsymbol{x},\boldsymbol{x}' ) = k(\boldsymbol{x}',\boldsymbol{x} )$ならば対称行列であることは直ちにわかります。正定値行列であることととかは安直に考えれば$\int \boldsymbol{x}^T k(\boldsymbol{x},\boldsymbol{x}' ) \boldsymbol{x}' d\boldsymbol{x} d\boldsymbol{x}'>0$ですかね。しかしカーネル関数を指定した段階でカーネル行列が良い性質を持っているかは判定できたほうが良いはずなので，おそらくそのような判定法は存在するのだろうと期待します。

### 3.2.3 ガウス過程の定義
ガウス過程のイメージをつかむうえでは図3.5がとてもわかりやすいですね。つまりカーネル関数$k(\cdot,\cdot)$を指定すると，いくつかの入力に対する出力たちが，そのカーネル$k$によって決まる共分散行列を持つガウス分布に従うということです。
カーネル法とガウス過程に書いてあることは、カーネル法について無知なのでよくわかりませんでした。勉強します。

### 3.2.4 ガウス過程からのサンプル
動径基底関数をカーネルに選んだときを考えます。
"近い"入力データを持ってくるとカーネル関数が大きな値を持ち，よってより大きな相関を持つため，出力も似たような値になります。逆に"遠い"入力データに対しては出力はあまり相関を持ちません。よって観測値などを使わずにサンプルを取り出すと滑らかにつながった出力$\boldsymbol{y}$が得られるはずです。しかし，繰り返しますが，このままでは観測値の情報を使っていないため，てんでばらばらの謎の関数です。そこで観測値のデータを与えた条件付確率を用いてサンプルを取り出すことで，回帰が可能になるということをしばらく後で見ます。

# 3.3 ガウス過程とカーネル
### 3.3.1 RBFカーネルと基底関数
動径基底関数を基底関数とした線形モデルを考え，細かい刻み幅で並べると動径基底関数カーネルになるということが述べられています。しかし(3.31)式から(3.32)式の変形が微妙なように感じます。$dh$に対応するものがないというのはどうも不自然だと思います。
この場合
$$
\phi_h (x) = \dfrac{\tau}{\sqrt H} \exp \left( - \dfrac{ (x - h/H)^2} { r^2} \right)
$$
のような修正をすればうまく変形ができるはずです。
なお，$\boldsymbol{x}$の空間次元が$n$次元のときも全く同様に証明することができます。$n$次元空間でのガウス積分が出てきますが，結局定数として$\theta_1$に押し込められるので問題ありません。

カーネルと関数形に載っている内容は難しいですが面白いですね。カーネル関数$k$を持つガウス過程からのサンプル(関数)の空間が，その期待値の空間を含んでいるのは当たり前のように感じますが，その差が無限小になるというのはよくわかりません。期待値をとってしまうため，後者の空間はかなり小さくなるように感じるのですがどうでしょうか。カーネル法について勉強しようと思いました。

### 3.3.2さまざまなカーネル
いろいろなカーネルの特徴が述べられています。その中でもマターンカーネルの話に感心しました。確率過程で選び出したサンプルの解析性についてまで調べることができるというのはすごいことだと思います。少し特殊関数への興味が沸いたような気がします。

#3.4 ガウス過程回帰モデル
ようやく回帰の話になりました。ここも少しややこしいです。
問題を明確にします。既知の入力データ$X = (\boldsymbol{x} _1 ,\cdots,\boldsymbol{x} _N)$と出力データ$\boldsymbol{y}$があります。また適当にカーネル関数$k(\boldsymbol{x},\boldsymbol{x} ')$を選びます。このとき新しい入力データ$\boldsymbol{x} ^* $に対する出力$y ^*$の確率分布を求めることが目標です。


まず$\boldsymbol{y}$を知らないものとします。いま$y$と$\boldsymbol{x}$の間に$y = f(\boldsymbol{x})$の関係があり，この関数$f$が平均$0$でカーネル関数$k$のガウス過程によって生成されているものとします。このとき$K _{n n'} = k(\boldsymbol{x} _n,\boldsymbol{x} _{n'})$で与えられるカーネル行列$K$を用いて，$\boldsymbol{y} \sim  \mathcal{N} (0,K)$となります。

### 3.4.1 ガウス過程回帰の予測分布
ここで$\boldsymbol{y}' = (y_1,\cdots,y_N,y^*)^T$の分布を考えると，これもまたガウス分布に従います。つまり$\boldsymbol{y}' \sim \mathcal{N} (0,K')$となります。ただし$K'$は$K$に$N+1$行と$N+1$列に$\boldsymbol{x}^ *$による成分を追加したものです。ここで$\boldsymbol{y}$の情報を利用します。つまり$\boldsymbol{y}$が既知としたときの$y^ 
 *$の条件付確率分布は多変量ガウス分布についての知識から

$$
p(y ^* | \boldsymbol{x} ^ * \,X,\boldsymbol{y}) = \mathcal{N} (\boldsymbol{k} _*^T K^{-1}\boldsymbol{y} , k _{ * * } - \boldsymbol{k} _{ * } ^T K^{-1} \boldsymbol{k} _ * )
$$
となります。ここでようやくこの分布に従ってサンプルを生成することになります。この手前の分布では既知のデータ$\boldsymbol{y}$の情報が活かされていないため，全く謎の曲線を引くことになります。

予測値が複数ある場合も同様です。しかし(3.78)には注意です。目立たないように書いてありますが，$k _ { * * }$には対角成分しかありません。つまり$m \neq n$のとき，$y_ m ^ * $と$y_n^ * $の相関がないということです。
$$
\boldsymbol{k} _ { * * } (n,m) = k(\boldsymbol{x}_n^*,\boldsymbol{x}_m ^ *) \delta(n,m)
$$
つまり$1$つ$1$つ新しい予測をしても，それは他の予測に影響を及ぼさないため問題ないということを意味しています。これを式で表したのが(3.80)ですが，出力$ y _  m ^ * $と$ y _ n ^ * $の間に相関がないということを考えれば自明です。

### 3.4.2 ガウス過程回帰の計算
図3.17はここまでに示した手順をそのまま実行するだけのアルゴリズムです。図3.18はforループを使用せず，大きな行列を一気に計算してしまうことで計算時間を短縮しようしています。大まかにいうと，以下のようになっています：

```math
\begin{eqnarray}
 X &=& (\boldsymbol{x}_1,\cdots,\boldsymbol{x} _N) \\
 \boldsymbol{z} &=& (\boldsymbol{x}_1^2,\cdots,\boldsymbol{x}_N^2) \\
 (Z_1) _ { i j } &=& \boldsymbol{x} _ i ^2 \\
 (Z_2) _ { i j } &=& \boldsymbol{x} _ j ^2 \\
 (X^T X) _ { i j } &=& \boldsymbol{x}^T_i \boldsymbol{x}_j \\
 J _ { i j } &=& (Z_1) _ { i j }+(Z_2) _ { i j } -2 (X^T X) _ { i j } = (\boldsymbol{x} _ i - \boldsymbol{x} _ j)^2 \\
 K _ {i j } &=& k(\boldsymbol{x} _ i ,\boldsymbol{x} _ j ) \\
\end{eqnarray}
```
数式で書くとむしろわかりずらくなったような気もしますが，要するにfor分を使わないために組み込みの関数を利用して処理しているということです。

### 3.4.3 ガウス過程回帰の要素表現
カーネル法についての記述はとりあえず飛ばします。
ガウス過程はNNの中間ノード無限の極限として得られるそうです。コラムに詳しい記述があります。

# 3.5 ガウス過程回帰のハイパーパラメータ推定
しれっと(3.87)にデルタ関数がありますが，対角成分にだけ効いてきた誤差$\boldsymbol{\epsilon}$の分散に対応します。しかしそのときの注釈では$n \neq n'$だが$\boldsymbol{x} _ n = \boldsymbol{x} _ {n '}$のときには効いてきてはいけないと述べています。デルタ関数で実装してしまうとこれを区別できなくなってしまうためまずいような気もします。まぁ大した影響はなさそうですが一応。
ちなみにここについても，
$$
p(\boldsymbol{y} | X,\theta) = \dfrac{p(\boldsymbol{y} ,\theta | X)}{p(\theta|X)} = \dfrac{p(\boldsymbol{y} | X)}{p(\theta|X)} p(\theta | X,\boldsymbol{y})
$$
という関係があります。実際にしたいことは$p(\theta | X,\boldsymbol{y})$を最大化する$\theta$を求めるということ(のはず)ですが，$p(\theta|X)$が未知のためできません。$p(\theta | X)$が$\theta$について一様分布と信じればテキストと同じ条件が得られますし，$\theta$について正規分布していると信じればリッジ回帰のときと同じようになります。ここではテキスト通り一様分布しているものと考えましょう。
エビデンスを最大化するパラメータ$\theta$を計算するためには勾配法を利用します。そのためには行列$K_\theta$を含んだ式を微分する必要がありますが，付録にある公式を利用するとそれほど手間なく計算できます。数値微分せずとも解析的に計算できるということがうれしいですね。しかもパラメータの数は通常の機械学習に比べて圧倒的に少ないため，おそらくこのハイパーパラメータ推定はそれほど大変ではないのではないかと思います。
細かいですが付録の誤植です。公式A.3のところですが，余因子行列の説明で$C^{-1} = |C| \tilde{C}$となっているのは$C^{-1} = |C|^{-1} \tilde{C}$の間違いです。公式A.3の証明はすっかり忘れてましたが，こんなに簡単にできるんですね。

# 3.6　ガウス過程回帰の一般化
まず大前提として，考えている入力$\boldsymbol{x}$と出力$y$は$y = f(\boldsymbol{x}) + \epsilon$という関係で結ばれています。ここで$f$は適当な滑らかな関数で，$\epsilon$は誤差を表します。今までは誤差が正規分布に従っていたため，$p(\boldsymbol{y} | \boldsymbol{f})$も正規分布に従っていました。ベイズの公式から
$$
p( \boldsymbol{f} | \boldsymbol{y} ) = \dfrac{p(\boldsymbol{f})}{p(\boldsymbol{y})} p(\boldsymbol{y} | \boldsymbol{f})
$$
が成立します。$f$がガウス過程であるため$p(\boldsymbol{f})$はガウス分布ということを考えると，$p( \boldsymbol{f} | \boldsymbol{y} ) $もガウス分布に従っていました。ここで$\epsilon$の分布として別のものを考えるとどのようになるかを考えます。

### 3.6.1 ロバストなガウス過程回帰
$p(\boldsymbol{f})$は相変わらずガウス分布ですが，誤差$\epsilon$がコーシー分布に従うものとします。このとき誤差の項をカーネルの中に含めることができず，$\boldsymbol{f}$の推定に条件付分布を用いた解析的手法のようなものは使うことができなくなります。ただコーシー分布を用いるご利益としては外れ値に対しロバストになるということです。これはコーシー分布自体がロバストであることを考えれば自然です。

### 3.6.2 ガウス過程識別モデル
ガウス過程回帰はロジスティック関数などを利用すると識別に利用することもできます。正規分布を積分したプロビット関数の近似としてロジスティック関数をとらえるということは初めて知りました。ものすごい近似精度でほとんど区別できません。

### 3.6.3 ポアソン回帰モデル
ある期間の間に機械が故障する回数などの従う分布のモデルとしてポアソン分布があります：
$$
p(y | \lambda) = \dfrac{\lambda^y e^{-\lambda}}{y!}
$$
この分布に従うと$y$の期待値は$\lambda$になります。ここで$\lambda$がある$x$に依存する場合を考えます。このままだとイメージしづらいので，$y$をある期間での機械の故障回数，$x$を時間(たとえば年度など)としましょう。このときいくつかのデータを利用して平均故障回数の推移$\lambda(x)$を確率分布も含めて推測するということを目標とします。
手続きは今までとまったく同様です。つまり様々な年度$x$の値に対応する平均値を$\boldsymbol{\lambda}$としたとき，
$$
p( \boldsymbol{\lambda} | \boldsymbol{y} ) = \dfrac{p(\boldsymbol{\lambda})}{p(\boldsymbol{y})} p(\boldsymbol{y} | \boldsymbol{\lambda})
$$
によって$\boldsymbol{\lambda}$の分布がわかります。$p(\boldsymbol{\lambda})$として今まで通りガウス過程(ガウス分布)，$p(\boldsymbol{y} | \boldsymbol{\lambda})$としてポアソン分布を用いれば$\lambda(x)$の概形が分布付きで得られるという感じです。

#コラム ニューラルネットワークとガウス過程
$1$層のNNの隠れ層のサイズを無限大にするとガウス過程と等価になることの説明があります。さらっと読んだだけですが，それほど難しい議論はしてなさそうです。ここを読むとガウス過程を使うことの利点がかなり感じられました。また単純に理論的な観点からもここは面白いです。