機械学習に詳しくなりたいブログ

機械学習や数学について勉強した内容を中心に書きます。100%趣味です。記事は数学的に厳密でなかったり誤りを含んでいるかもしれません。ご指摘頂ければ幸いです。

ガウス過程

概要

今回はガウス過程の導出と、線形回帰における出力\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}です。ここでカーネル法による正則化最小二乗法(1)の式(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}が生成されるのですから、なめらかな曲線が得られることがわかります。

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

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

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

ガウス過程による回帰の記事1つ目:

www.iwanttobeacat.com

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

www.iwanttobeacat.com


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

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()