ガウス過程

概要

今回はガウス過程の導出と、線形回帰における出力 $\mathbf{y}$ の事前分布はカーネル関数を使って表すことができるということの確認です。計算自体は難しくないのですが、ガウス過程という名称から確率過程のように時間軸を意識してしまうと戸惑ってしまいます。

線形回帰における出力 $\mathbf{y}$ の事前分布

まずは線形回帰で何度もやってきた通り、回帰式のモデルを

$y = \mathbf{w}^{T}\mathbf{x} \tag{1}$

と設定します。線形回帰をMAP推定で解くで計算したのと同じく、係数 $\mathbf{w}$ の事前分布を

$p(\mathbf{w}) = N(\mathbf{w}|\mathbf{0},\alpha^{-1}\mathbf{I}) \tag{2}$

とします。平均が0で、各パラメータには相関がなく、全て等しい分散 $\alpha^{-1}$ であるということです。線形回帰をMAP推定で解くの式(5)で $\mathbf m_0=0$ 、 $\mathbf S_0 = \alpha^{-1}\mathbf I$ とすれば同じことです。

ここで最小二乗法の解の導出の式(4)と同じように計画行列 $\mathbf{X}$ を用いて

$\mathbf{y} = \mathbf{X}\mathbf{w} \tag{3}$

とします。 $\mathbf{y}$ は訓練データ $x_1, \cdots , x_N$ に対する値 $y_1, \cdots , y_N$ をまとめたベクトルになっています。冒頭で書いた通り、今回はこの $\mathbf{y}$ の分布を求めることが目標です。

さて、正規分布に従う独立な確率変数の和は、やはり正規分布に従うという性質があります。正規分布に従うと仮定した係数 $\mathbf{w}$ と入力 $\mathbf{x}$ の線形和が $y$ ですから、この性質より $y$ は正規分布に従うことがわかります。*1

今回の目的は $\mathbf{y}$ の分布を求めることでしたから、正規分布に従うということは平均と共分散がわかれば良いということになります。まずは平均を求めますと、

$\begin{eqnarray*} E[\mathbf{y}] &=& E[ \mathbf{X}\mathbf{w}] \tag{4} \\ &=& \mathbf{X} E[\mathbf{w}] \tag{5} \\ &=& \mathbf{0} \tag{6} \end{eqnarray*}$

です。 $\mathbf{w}$ は式(2)で平均 $\mathbf{0}$ としていますから式(5)は $\mathbf{0}$ と計算できます。

次に $\mathbf{y}$ の共分散を求めます。多変量正規分布の式(2)の通り、

$\mathrm{Cov}(y_a,y_b) = E[y_a y_b] - E[y_a]E[y_b] \tag{7}$

です。ここで式(6)より $E [ \mathbf{y} ] = \mathbf{0}$ でしたから、

$\mathrm{Cov}(y_a,y_b) = E[y_a y_b] \tag{8}$

です。共分散行列は多変量正規分布の式(4)ですので、

$\begin{eqnarray*} \mathrm{Cov} [ \mathbf{y} ] &=& E[ \mathbf{y}\mathbf{y}^{T} ] \tag{9} \\ &=& E[ \mathbf{X}\mathbf{w} \mathbf{w}^{T}\mathbf{X}^{T} ] \tag{10} \\ &=&\mathbf{X} E[ \mathbf{w} \mathbf{w}^{T} ] \mathbf{X}^{T} \tag{11} \\ \displaystyle &=& \mathbf{X} \frac{1}{\alpha} \mathbf{I} \mathbf{X}^{T} \tag{12} \\ \displaystyle &=& \frac{1}{\alpha} \mathbf{X} \mathbf{X}^{T} \tag{13} \\ \end{eqnarray*}$

と計算できます。式(9)から式(10)は転置行列の定理を使っています。また、 $E[ \mathbf{w} \mathbf{w}^{T}$ ]は前提条件より $\alpha^{-1}\mathbf{I}$ です。ここでカーネル法による正則化最小二乗法（１）の式(10)より、

$\mathbf{K}= \mathbf{X}\mathbf{X}^{T} \tag{14}$

でした。以上をまとめると、求めたい分布は

$p(\mathbf{y}) = N(\mathbf{y}|\mathbf{0},\alpha^{-1}\mathbf{K}) \tag{15}$

となることがわかります。カーネル関数さえ定まればいきなり $y$ の事前分布が求まってしまうっていうことなんですが、なんか驚きです。係数 $\mathbf w$ の事前分布を設定し、計算の過程で $\mathbf w$ の前提条件から $\mathbf w$ を式からなくしてしまい、 $y$ の事前分布になるという。

さて、ここでようやくタイトルのガウス過程の話です。ガウス過程とは、任意の個数の入力を選んだとき、それに対する出力が常に多次元正規分布に従うことをいいます。ちょっと説明が正確でないかもしれません、、、Wikiを参考にしました：ガウス過程 - Wikipedia

$p(\mathbf{y})$ は式(15)の通り多次元正規分布ですから、これを満たしています。つまりガウス過程となっています。確率過程のように時間軸が関係しそうな名前ですが、とにかく多数個入力の出力が多次元正規分布に従っていればガウス過程と呼べます。

式(15)を使って、 $\mathbf y$ を出力してみた実験結果です。いくつかサンプリングしてみました。

上のグラフは、5次の多項式モデルを考え、 $\mathbf{K}$ を計算しています。下のグラフはガウスカーネルで $\mathbf{K}$ を計算しました。ガウスという名前がたくさん出てきてややこしいですが、ガウスカーネルを使っていない上の図もガウス過程です。

これらの結果はなめらかな曲線になっていますが、これは $\mathbf{K}$ が実際どのような値になっているかを見れば明らかです。例えば-1～1の範囲で均等に5点 $x$ を取り、この $x$ に対してガウスカーネルで $\mathbf{K}$ を計算すれば、以下のようになります。（ $\sigma=0.3$ としました。）

$\mathbf{K} = \left( \begin{array}{ccccc} 1& 0.249& 0.00387& 3.73E-06& 2.23E-10 \\ 0.249& 1& 0.249& 0.00387& 3.73E-06 \\ 0.00387& 0.249& 1& 0.249& 0.00387 \\ 3.73E-06& 0.00387& 0.249& 1& 0.249 \\ 2.23E-10& 3.73E-06& 0.00387& 0.249& 1 \end{array} \right) \tag{16}$

$x$ の値が近いほうが共分散が高い、つまり同じような値を取る確率が高いことがわかります。また、逆に $x$ の値が遠ければ共分散は0に近く無相関に近くなります。この共分散行列によって $\mathbf{y}$ が生成されるのですから、なめらかな曲線が得られることがわかります。

入力は多次元になっても同じようにサンプリングすることができます。こちらもガウスカーネルを用いました。

ガウス過程を回帰や分類で使う

このガウス過程を回帰や分類問題に使っていきます。

ガウス過程による回帰の記事１つ目：

www.iwanttobeacat.com

回帰の最終的な実験結果はこちらです。

www.iwanttobeacat.com

ガウス過程による分類の記事１つ目

www.iwanttobeacat.com

分類の最終的な実験結果はこちらです。

www.iwanttobeacat.com

今回のコードです。

# ガウス過程からのサンプリング
import matplotlib.pyplot as plt
import numpy as np


# カーネル
def kernel(x1, x2):
    # 5次の多項式モデルを考えたときのカーネル関数（内積）
    ret = 0
    for i in range(5):
        ret += x1**i * x2**i
    return ret
    # ガウスカーネル
    # sigma = 0.3
    # return np.exp(-(x1 - x2) ** 2 / (2 * sigma ** 2))



# グラム行列
def gram_matrix(x):
    gram = np.empty([N, N])
    for i in range(N):
        for j in range(N):
            gram[i, j] = kernel(x[i], x[j])
    return gram


# x
N = 100
x = np.linspace(-1, 1, N)

# α
ALPHA = 1.0

# グラム行列K
K = gram_matrix(x) / ALPHA

for i in range(5):
    y = np.random.multivariate_normal(np.zeros(N), K)
    plt.plot(x, y)


# 結果の表示
plt.xlim(-1, 1)
plt.show()