ガウス過程回帰入門

こんにちは、ぐぐりら(<a href="https://twitter.com/guglilac">@guglilac</a>)です。

最近、機械学習プロフェッショナルシリーズのガウス過程と機械学習を購入してみたのですが、とてもわかりやすくて良い本でした。

Twitterでもいい評判が流れていたのもありずっと気になっていたのですが、時間がなくてすぐには買えず、ようやく買って勉強しています。

最後の方は難しそうで読めていないのですが、とりあえずガウス過程の大枠が学べたので基本的なところを自分の言葉で説明して実装することで定着させようと思います。

勉強し始めたばかりなので間違っているところあればコメントください。

### ガウス過程
ガウス過程とは、なにか。というところでまず自分は昔つまづいていました。
以前ネットに転がっている記事で勉強しようとした時は「ガウス過程とは無限次元のガウス分布である」とか「ガウス過程は関数上の確率分布だ」とかが最初に書いてあって、よくわかりませんでした。

勉強してみて、ガウス過程とは何かと聞かれれば確かに「ガウス過程とは無限次元のガウス分布である」は理解できるのでそう書かれているのはよく分かるのですが、いきなり書かれてもわからないよ、という気持ちにもなりました。

あと、ガウス過程の過程という単語は時系列の分野から生じたものですが、必ずしもインデックスは時刻を表すものではないので時刻は一旦忘れて考えてみるといいと思います。

今回は導出等の細かい話はせずに、ガウス過程ではどういうことを実際にはやっているのかを書いていきたいと思います。

$x_1,\cdots, x_N$という入力データ点があり、これらに対応するモデルの出力を$f(x_1),\cdots,f(x_N)$とします。
 
 機械学習ではこの関数$f$をいい感じにデータから見つけてこよう、というのが目標になります。
 
 ガウス過程では$f(x_1),\cdots,f(x_N)$がN次元ガウス分布$N(\mu,\Sigma)$に従うとします。
 
 イメージとしては、入力$x_1,\cdots, x_N$が与えられると、それによって決まる$N(\mu,\Sigma)$から出力$f(x_1),\cdots,f(x_N)$がぽんっとサンプルできるかんじです。
 
 
 普通に線形回帰などをするときは、入力点$x_1,\cdots, x_N$に対して出力$f(x_1),\cdots,f(x_N)$は確定します。関数の形はいつも一つです。
 
 一方ガウス過程では、入力が与えられても関数の形は一つに定まりません。
何回もサイコロをふるたびに違う出力$f(x_1),\cdots,f(x_N)$が$N(\mu,\Sigma)$から得られます。

ガウス過程での関数の形と言っていますが、「$f(x_1),\cdots,f(x_N)$はN点しかないのに関数ってなに？」という疑問を持つかもしれません。というか以前自分も思っていました。

実際には、ガウス過程は無限に入力$x_1,\cdots$をとってきた時に出力$f(x_1),\cdots$もガウス分布に従う、という定義になっています。とても細かく入力をとってきた時に$f(x_1),\cdots$がガウス分布に従って得られるので、これは関数と言えるわけです。

実際には無限を考えることはなく、有限の入力に対して有限個の出力がガウス分布に従うのでそれだけ考えていても大丈夫です。

どの入力をとってきてもそれに対応する出力が多次元ガウス分布に従う、というところがポイントです。

なので、ガウス過程は関数が確率変数になっている確率分布と解釈できます。サイコロをふると関数が一つ得られるみたいなイメージです。

例えば、平均0,分散共分散行列が単位行列の多次元ガウス分布から関数を(すなわち$x_i$を細かくとった$f(x_1,\cdots,x_N)$を)3つランダムにサンプルしてくるとこんな感じです。

この状態では完全にランダムに関数がサンプルされているように見えますが、学習データを用いることでサンプルしてくる分布がデータに近づきます。

ガウス過程では、具体的には出力はガウス分布$N(0,K)$に従います。
ただし出力は平均をとることで0にできるのでここでは平均は0と考えます。

$K$はグラム行列などとよばれ、

$$
K= \Phi\Phi^T
$$

と定義します。ただし、$\Phi$は計画行列とします。

$K_{ij}=\phi(x_i)^T\phi(x_j)$
$\phi$は入力から特徴量空間への写像です。

$K_{ij}$が大きい、つまり$\phi(x_i)^T\phi(x_j)$が大きい、つまり$x_i$と$x_j$が特徴量空間で似ている点のとき、$f(x_i),f(x_j)$の共分散も大きくなるので出力が近くなります。

似ている点は似ている出力になってほしいので、これは納得できると思います。

実際は、SVMなどのように特徴量写像$\phi$を具体的に計算する必要はなくて、特徴量空間での内積のみが必要なので、$x_i,x_j$間の特徴量空間での内積をカーネル関数$k(x_i,x_j)$として定義してこれを計算します。カーネルトリックというやつです。

今回の実装ではガウスカーネルを使っています。
ガウスカーネルと使うと得られる関数は無限回微分可能なとてもなめらかな関数になるという性質があるらしいです。

問題によってカーネルを変えることでデータ間の「似ている」という状態を指定することができます。これは数値データでなくても文字列や木構造といった構造にもカーネルを与えればいいのでいろんなデータ構造を扱うことができて良さげです。

### ガウス過程で回帰
ガウス過程を使って回帰モデルを作ります。

細かいことを書くと大変なのでここは省略していまいますが、

学習データ$x_1,\cdots, x_N$に対応する出力$f(x_1),\cdots,f(x_N)$がガウス分布に従っているとき、
新しい入力データ集合$X$に対応する出力$y$を予測したいです。

学習データ$x_1,\cdots, x_N$と新しい入力データ$X$を合体させて、出力も合体させて$f(x_1),\cdots,f(x_N),y$を考えると、これら全体がまたガウス分布に従います。

ガウス分布の一部の変数が条件づけられた場合、残りの変数もまたガウス分布に従うので、

学習データが与えられたもとでの新しいデータ点での条件つきガウス分布$p(y|X,D)$は頑張って計算すると次のようになります。

$$
p(y|X,D)= N(k^T_xK^{-1}y,k_{xx}-k_{x}^TK^{-1}k_{x})
$$

ただし、$k_x$は新しいデータ$X$と学習データのグラム行列(かならずしも正方行列にはならないのでそう呼ぶのかわからないけど。。雰囲気で察してください)、$k_{xx}$は新しいデータ$X$のグラム行列です。

各点における平均や分散が得られるので、この平均を回帰結果とすればいいし、分散もあるのでその回帰結果がどれくらい信頼できるかも求めることができます。

学習データに近い点の回帰にはある程度自信をもっているものの、学習データがあまりない領域のデータに対して回帰すると分散は大きくなってしまう、なんてことが表現できます。

この「予測結果の確信度がわかる」という性質はベイズ最適化などでも使われています。

また、計算量についてですが、$K$の逆行列を計算する際に$O(N^3)$($N$はデータ数)かかるのがガウス過程の一番の短所です。

計算量を落とす工夫もたくさんあるらしいです。ガウス過程と機械学習の参考書にもいろいろ載っています。

### 実装
pythonで実装しました。
ガウス過程回帰モデルです。

```python
class GaussianProcessRegressor:
    def __init__(self, kernel, c):
        self.kernel = kernel
        self.c = c

def fit(self, X, y):
        self.X_train = X
        self.y_mean = np.mean(y)
        y_train = y - self.y_mean
        self.y_train = np.reshape(y_train, (len(y_train), 1))
        self.K = self.get_gram_matrix(X, X, flag=True)

def predict(self, X):
        mu, _ = self.calc_mu_sigma(X)
        return mu

def calc_mu_sigma(self, X):
        k_star = self.get_gram_matrix(self.X_train, X, flag=False)
        k_starstar = self.kernel.run(X, X) + self.c
        K_inv = np.linalg.inv(self.K)
        mu = k_star.T @ K_inv @ self.y_train
        sigma = k_starstar - k_star.T@K_inv@k_star
        return mu + self.y_mean, sigma

def get_gram_matrix(self, A, B, flag):
        """AとBのgram行列を返す.
        return's shape (len(A),len(B))
        flag: クロネッカーつけるかどうか
        """
        return np.array([[self.kernel.run(A[i], B[j]) + self.c * self.kronecker(i, j, flag) for j in range(len(B))] for i in range(len(A))])

def kronecker(self, i, j, flag):
        return 1 if i == j and flag else 0
```

今回使ったガウスカーネルです。
```python
class GaussianKernel:
    def __init__(self,h):
        self.h = h
        
    def run(self,x,y):
        return np.exp(-np.linalg.norm(x-y)**2/self.h)
```

### 実験と結果

今回は$y= x\cos x$からデータを生成しています。ノイズは乗せられる実装にしてはいますが、簡単のため実験ではノイズなしで行っています。

こんな感じ。

このデータを学習データとしてガウス過程回帰モデルをtrainすると、結果はこうなりました。

赤い曲線が予測結果(推定したガウス分布の平均)で、緑色の領域は推定したガウス分布の共分散行列を使って計算した2σ区間を表していて、予測の確信度が表現されています。

この学習した結果のガウス分布から関数を3つサンプルしてくるとこうなります。

先ほど完全にランダムにサンプルしてきた時と比較して、真の関数の形に近い関数がサンプルされていることがわかります。

次のgifは一つずつデータを増やしていったときの予測結果の変化を表しています。

データが少ないと分散(緑色の幅)も大きくなっていて、データが増えてくると観測された点の近くは信頼区間が小さくなっていきます。

赤色の線がガウス過程の予測結果です。(ガウス分布の平均)

### おわりに

今回はガウス過程回帰を実装してみました。
細かい導出とかは省略して、自分が勉強した時に(する前に)詰まっていたところをなるべく日本語にしてみたのですが、うまく説明するにはもう少し勉強が必要ですね。。。

実装も、本当はカーネルのハイパーパラメータを勾配法やMCMCで最適化することができるのですが、今回は時間がなく断念。。

詳しく勉強したい方はぜひ参考書買って勉強するといいと思います。わかりやすかったです。

ありがとうございました。

このブログを検索

ぐぐりらにっき

ガウス過程回帰入門

コメント

コメントを投稿

ガウス過程回帰入門

この記事をシェアする

コメント

コメントを投稿