この記事では、主にガウス過程回帰モデルについて説明します。ガウス過程の導出については、詳しくは参考文献[1]を参照してください。

ガウス過程とは

入力$\mathbf{x} \rightarrow y$を予測する回帰関数（regressor)の確率モデル

データ $D=\{ (\mathbf{x_{n}}, y_{n}) \}_{n=1}^{N}$

が与えられたとき、新しい$\mathbf{x_{n+1}}$に対する$y_{n+1}$を予測する。

$y$の期待値だけではなく分散も予測できるので、ベイズ最適化による$\mathbf{x}$の探索に用いられる。

図1. ガウス過程回帰の例。ばつ印はデータ点、実線は_y_の期待値、青帯は_y_の$\pm 2 \sigma$の範囲を示す。

ガウス過程の導出

入力__x__の特徴ベクトルを$\phi (\mathbf{x}) = (\phi_{0}(\mathbf{x}),\phi_{1}(\mathbf{x}),...,\phi_{H}(\mathbf{x}))^{T}$としたとき、線形回帰モデルは

\begin{align}
\hat{y} &= \mathbf{w}^{T}\mathbf{\phi}(\mathbf{x}) \\
&= w_{0}\phi_{0}(\mathbf{x}) + w_{1}\phi_{1}(\mathbf{x}) \cdots + w_{H}\phi_{H}(\mathbf{x})\\
\end{align}

と表せる。_N_個の訓練データに対して、まとめて書くと

\begin{pmatrix}
\hat{y_{1}} \\
\hat{y_{2}} \\
\vdots \\
\hat{y_{N}} \\
\end{pmatrix}
=
\begin{pmatrix}
\phi_{0}(x_{1}) & \phi_{1}(x_{1}) & \cdots & \phi_{H}(x_{1}) \\
\phi_{0}(x_{2}) & \phi_{1}(x_{2}) & \cdots & \phi_{H}(x_{2}) \\
\vdots          &                 &        & \vdots \\
\phi_{0}(x_{N}) & \phi_{1}(x_{N}) & \cdots & \phi_{H}(x_{N}) \\
\end{pmatrix}
\begin{pmatrix}
w_{0} \\
w_{1} \\
\vdots \\
\vdots \\
w_{H} \\
\end{pmatrix}

となる。計画行列$\mathbf{\Phi}$を使って書くと

$$ \mathbf{\hat{y}} = \mathbf{\Phi}\mathbf{w} $$

となる。__計画行列$\mathbf{\Phi}$は定数行列であることに注意する。以降、簡単のため$\mathbf{\hat{y}} = \mathbf{y}$とする。ここで、リッジ回帰の場合と同様に、重み$\mathbf{w}$が、平均__0、分散$\lambda^{2}\mathbf{I}$のガウス分布

$$ \mathbf{w} \sim \mathcal{N}(\mathbf{0},\lambda^{2}\mathbf{I}) $$

から生成されたとものとする。ガウス分布の定数倍はガウス分布となるため$\mathbf{y}=\mathbf{\Phi w}$もガウス分布に従う。すると__y__の期待値は$E[\mathbf{y}]=E[\mathbf{\Phi w}]=\mathbf{\Phi}E[\mathbf{w}]=\mathbf{0}$となる。また$\mathbf{y}$の共分散行列は、

\begin{align}
\Sigma = E[\mathbf{yy}^{T}]-E[\mathbf{y}]E[\mathbf{y}]^{T}&=E[(\mathbf{\Phi w})(\mathbf{\Phi w})^{T}] = \mathbf{\Phi}E[\mathbf{ww}^{T}]\mathbf{\Phi}^T \\
&= \lambda^{2} \mathbf{\Phi \Phi}^{T}
\end{align}

となる。結果として、__y__の分布は多変量ガウス分布

$$ \mathbf{y} \sim \mathcal{N}(\mathbf{0}, \lambda^{2} \mathbf{\Phi \Phi}^{T})$$

となる。線形回帰モデルにあった重み$\mathbf{w}$が__期待値を取ることによって消えていることに注意する。__

ガウス過程の定義

どんな_N_個の入力の集合$(\mathbf{x_{1}},\mathbf{x_{2}},\cdots,\mathbf{x_{N}})$についても、対応する出力$\mathbf{y}=(y_{1},y_{2},\cdots,y_{N})$の同時分布$p(\mathbf{y})$が多変量ガウス分布に従うとき、x__と__y__の関係は__ガウス過程（Gaussian process）に従う

と言う。

カーネルトリック

__y__の共分散行列を

$$ \mathbf{K} = \lambda^{2} \mathbf{\Phi \Phi}^{T} $$

とおく。この$(n,n')$要素は、

$$ K_{nn'} = \lambda^{2}\mathbf{\phi}(\mathbf{x_{n}})\mathbf{\phi}(\mathbf{x_{n'}})^{T}$$

と与えられる。ここで明示的に特徴ベクトル$\mathbf{\phi}(\mathbf{x})$を与える必要はなく、$K_{nn'}$の値がわかれば良い。この値を与える関数をカーネル関数と呼び、

$$ k(\mathbf{x_{n}},\mathbf{x_{n'}}) = \lambda^{2}\mathbf{\phi}(\mathbf{x_{n}})\mathbf{\phi}(\mathbf{x_{n'}})^{T} $$

と書く。カーネル関数は入力__x__の類似度を測る関数と見ることができる。

ガウスカーネル

カーネル関数として、ガウスカーネル、または動径基底関数がよく用いられる。

k(\mathbf{x}, \mathbf{x'}) = \theta_{1} \exp \biggl( -\frac{|\mathbf{x} - \mathbf{x'}|^{2}}{\theta_{2}} \biggr)

補足

ガウスカーネルは、入力__x__の各座標軸に沿って無限個の動径基底関数をグリッド上に配置し、無限個の要素からなる特徴ベクトル$\phi(\mathbf{x})$を求め、その内積を取ったものに等しい。

観測ノイズ

_y_にノイズがのっている場合、

$$ y_{n} = f(\mathbf{x_{n}}) + \epsilon_{n} $$

と表せる。観測ノイズ$\epsilon$が平均0、分散$\sigma^{2}$のガウス分布に従う場合、__y__の分布は

$$ p(\mathbf{y}) = \mathcal{N}(\mathbf{0}, \mathbf{K} + \sigma^{2}\mathbf{I}) $$

となる。各観測におけるノイズは相関がないので、共分散行列の対角成分が__0__となることから、このようになる。カーネル関数は、

$$ k'(\mathbf{x_{n}}, \mathbf{x_{n'}}) = k(\mathbf{x_{n}}, \mathbf{x_{n'}}) + \sigma_{2} \delta(n, n') $$

となる。ここで$\delta$は$n=n'$の時は1、それ以外の時は0を返す関数である。

ガウス過程回帰モデル

データ$D=\{ (\mathbf{x_{n}}, y_{n}) \}_{n=1}^{N}$が与えられたときデータに含まれない$\mathbf{x^{*}}$に対する値$y^{*}$を予測するには、

$\mathbf{y'}=(y_{1},y_{2},\cdots,y_{N},y^{*})$として

$(x_{1},x_{2},\cdots,x_{N})$に$x^{*}$を加えた$(n+1)\times(n+1)$の共分散行列$\mathbf{K'}$とすれば、これら全体がまたガウス分布に従うことから、

$$ \mathbf{y'} \sim \mathcal{N}(\mathbf{0},\mathbf{K'}) $$

となる。つまり

\begin{pmatrix}
\mathbf{y} \\
y^{*}
\end{pmatrix}
\sim
\mathcal{N}
\begin{pmatrix}
\mathbf{0},
\begin{pmatrix}
\mathbf{K} & \mathbf{k_{*}} \\
\mathbf{k}_{*}^{T} & k^{**}
\end{pmatrix}
\end{pmatrix}

と表せる。ここで、

\left\{
\begin{align}
\mathbf{k_{*}} &= (k(\mathbf{x}^{*}, \mathbf{x}_{1}),k(\mathbf{x}^{*}, \mathbf{x}_{2}),\cdots,k(\mathbf{x}^{*}, \mathbf{x}_{N}))^{T} \\
k_{**} &=  k(\mathbf{x}^{*},\mathbf{x}^{*})
\end{align}
\right.

は、新しい入力$\mathbf{x}^{*}$と学習データの入力の類似度（カーネル関数の値）を並べたベクトル、

および$\mathbf{x}^{*}$自身との類似度である。

ガウス分布の公式より、

\begin{pmatrix}
\mathbf{y}_{1} \\
\mathbf{y}_{2}
\end{pmatrix}
=\mathcal{N}
\begin{pmatrix}
\begin{pmatrix}
\mu_{1} \\
\mu_{2}
\end{pmatrix}
,
\begin{pmatrix}
\Sigma_{11} & \Sigma_{12} \\
\Sigma_{21} & \Sigma_{22}
\end{pmatrix}
\end{pmatrix}

において、ベクトルの一部$\mathbf{y_{1}}$が与えられたとき、$\mathbf{y_{2}}$はガウス分布

p(\mathbf{y_{2}}|\mathbf{y_{1}}) = \mathcal{N}(\mu_2 + \Sigma_{21} \Sigma^{-1}_{11} (y_{1} - \mu_{1}), \Sigma_{22} - \Sigma_{21}  \Sigma^{-1}_{11}  \Sigma_{12})

に従う。ここでは$\mathbf{\mu_{1}}=\mathbf{0}, \mathbf{\mu_{2}}=\mathbf{0}$なので、

p(\mathbf{y_{2}}|\mathbf{y_{1}}) = \mathcal{N}(\Sigma_{21} \Sigma^{-1}_{11}y_{1}, \Sigma_{22} - \Sigma_{21}  \Sigma^{-1}_{11}  \Sigma_{12})

と簡単に表すことができる。よって、

p(y^{*}|\mathbf{x}^{*},D)=\mathcal{N}(\mathbf{k}^{T}_{*}\mathbf{K}^{-1}\mathbf{y},k_{**} - \mathbf{k}^{T}_{*}\mathbf{K}^{-1}\mathbf{k}_{*})

がガウス過程による予測分布として与えられる。上の式から$y^{*}$の期待値、分散を予測できる。

アルゴリズム

学習データのカーネル行列__K__を計算
__K__の逆行列を計算
$\mathbf{yy} = \mathbf{K}^{-1}\mathbf{y}$を計算
$\mathbf{x}^{*}$に対する学習データの入力の類似度$\mathbf{k_{*}}$を計算
$\mathbf{x}^{*}$の自身に対する類似度$s = k(\mathbf{x}^{*}, \mathbf{x}^{*})$を計算
$y^{*}$の平均$\mu = \mathbf{k_{*}} \mathbf{yy} $を計算
$y^{*}$の分散$\sigma^{2} = s - \mathbf{k_{*}} \mathbf{K}^{-1} \mathbf{k}^{T}_{*} $を計算
学習データに含まれない複数の入力に対する出力を予測する場合、4から7を繰り返す

サンプルコードによる数値実験

$\mathbf{x}$が一次元の場合のサンプルコード。$\mathbf{x}$が高次元の場合は、ガウスカーネルを計算する関数を改良すればいけるはず・・・。Jupyter notebook上で実装しています。

データはMLPシリーズ『ガウス過程回帰と機械学習』サポートページからダウンロードしました。

ライブラリのインポート

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

データのロード

data = np.loadtxt("data/simple2.dat")

x_train = data[:, 1]
y_train = data[:, 0]

plt.plot(x_train, y_train, "x")

plt.ylim(-1.3, 4.8)
plt.xlim(-3, 4)
plt.xlabel("x")
plt.ylabel("y")

plt.tight_layout()

plt.savefig("data.png", dpi=300)

図2. サンプルデータ

ガウス過程回帰の関数

パラメータ

ガウスカーネル$\theta_{1} = 1, \theta_{2} = 0.4$
ノイズの分散$\theta_{3} = 0.1$


def gaussian_kernel(x, x_train, t1=1, t2=0.4):
    xx = (x_train - x) ** 2
    return t1 * np.exp(-xx / t2)
    
def gpr(x_test, x_train, y_train, kernel, t3=0.1):
    
    mean = y_train.mean()
    
    K = []
    for x in x_train:
        K.append(kernel(x, x_train))
    K = np.array(K) + np.eye(len(x_train)) * t3
    A = np.linalg.inv(K)
    yy = np.dot(A, y_train - mean)

    mu = []
    var = [] 
    for x in x_test:
        k = kernel(x, x_train)
        s = gaussian_kernel(x, [x]) + t3
        mu.append(np.dot(k, yy))
        var.append(s[0] - np.dot(k, np.dot(A, k.T)))
    
    mu = np.array(mu) + mean
    var = np.array(var)
    
    return mu, var

実行結果

x_test = np.linspace(-3, 4, 71, endpoint=True)
mu, var = gpr(x_test, x_train, y_train, gaussian_kernel)

plt.plot(x_test, mu)
plt.plot(x_train, y_train, "x", c="C0")
plt.fill_between(x_test, mu - var * 2, mu + var * 2, facecolor='C0', alpha=0.5)
plt.title(r"$\theta_{1}=1, \theta_{2}=0.4, \theta_{3}=0.1$")
plt.ylim(-1.3, 4.8)
plt.xlim(-3, 4)
plt.xlabel("x")
plt.ylabel("y")
plt.tight_layout()

plt.savefig("GP.png", dpi=300)

図3. ガウス過程回帰の実行結果。実線は_y_の期待値、青帯は_y_の$\pm 2 \sigma$の範囲を示す。

考察

雑音ありの設定なので$(\theta_3 = 0.1)$、実線が学習データを通っていない。
学習データの入力付近では分散が小さくなっている。
ベイズ最適化では、期待値に標準偏差の定数倍を足した値が最大となる点を探索するため、図3の例では2から少し離れたところを探索するはずである。
滑らかな実線が得られており、少数のデータに対しても過学習せず、柔軟（？）にフィットできている。

まとめ

ガウス過程回帰モデルについて説明した。

ガウス過程回帰は、学習データの入力と新しい入力$x^{*}$との類似度に基づいて、ガウス分布の公式から、$y^{*}$の期待値と分散を予測するものである。

参考文献