ガウス過程（Gaussian Process）

6. ガウス過程（Gaussian Process）#

ガウス過程（Gaussian Process, GP）は、関数の分布を扱うベイズ的なモデルである。主に回帰、ベイズ最適化などで利用される。

特に以下では、関数の回帰を行う関数の分布を推定・生成する方法としての観点から解説する。

理論的な詳細等については下記のガウス過程の名著(通称: GPML)を参照されたい。著者のページから無償でPDFが入手可能である:
"Gaussian Processes for Machine Learning" by Carl Edward Rasmussen and Christopher K. I. Williams

6.1. 概要#

../_images/fb1984c68ee78a755e4fa62536b5e416ef4bd95225a60d3c12adf0b2feba04cd.png

今、図のようにデータが観測されているとする。このデータを説明する関数は無数に存在する。そんな中、単一の関数を選ぶのではなく、「もっともらしい」関数の分布を考えることができれば、より柔軟にデータを説明できると考えられる。

ガウス過程は、このような関数の分布を直接扱うための枠組みの１つである。背後にある考え方については、以下のとおりである。

任意の近傍二点に対して、その点での関数値は似ている(≒連続)であろう。
"似ている"度合いを非負の値を持つカーネル関数で定量化する。
既知のデータと未知の点での関数値の同時分布を考え、それが多変量ガウス分布に従うと仮定する。

とくに最後の点がガウス過程の名前の由来であると理解してもらってよい。

6.2. ガウス過程の定義#

観測データが \(\mathcal{D} = \{(\mathbf{x}_i, y_i)\}_{i=1}^{n}\) で与えられるとする。ここで、\(\mathbf{x}_i \in \mathbb{R}\) は説明変数、\(y_i \in \mathbb{R}\) は目的変数である。説明変数・関数値とも、多次元でも構わないが、ここでは単純化のために1次元としよう。

自分の興味のある\(m\)個の点 \(\mathbf{x}^* \) に対して、関数値 \(\mathbf{y}^*\) を予測したいとする。このときの\(\mathbf{x}^* = \{x_1^*, x_2^*, \ldots, x_m^*\}\) は、観測データ \(\mathcal{D}\) の説明変数 \(\mathbf{x}_i\) (データ点)とは異なる点であるとし、予測点とでも呼ぶことにする。

そのとき、\(\mathbf{x}, \mathbf{x}^*\) に対する関数値 \(\mathbf{y}, \mathbf{y}^*\) の同時分布が以下のように多変量ガウス分布に従うと仮定する。

\[\begin{split} \begin{pmatrix} \mathbf{y} \\ \mathbf{y}^* \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} \mu_\mathbf{x} \\ \mu_{\mathbf{x}^*} \end{pmatrix}, \begin{pmatrix} K(\mathbf{x}, \mathbf{x}) & K(\mathbf{x}, \mathbf{x}^*) \\ K(\mathbf{x}^*, \mathbf{x}) & K(\mathbf{x}^*, \mathbf{x}^*) \end{pmatrix} \right) \equiv \mathcal{N} \left( \boldsymbol{\mu}, \Sigma \right) \end{split}\]

ひとたび多次元正規分布の平均ベクトルと共分散行列が与えられれば、条件付き分布もまた多次元正規分布になるため、以下のように、予測点での関数値の分布を得ることができる。

\[ p(\mathbf{y}^*|\mathbf{y}) = \mathcal{N}\left( \mu_{\mathbf{x}^*} + K(\mathbf{x}^*, \mathbf{x}) K(\mathbf{x}, \mathbf{x})^{-1} (\mathbf{y} - \mu_{\mathbf{x}}), K(\mathbf{x}^*, \mathbf{x}^*) - K(\mathbf{x}^*, \mathbf{x}) K(\mathbf{x}, \mathbf{x})^{-1} K(\mathbf{x}, \mathbf{x}^*) \right) \]

多次元正規分布の条件付き確率分布の計算(3章)を用いた。この意味するところは、予測点での関数値 \(\mathbf{y}^*\) は、観測データ \(\mathcal{D}\) に基づいて平均と分散が決定される正規分布に従う、ということである。それぞれの予測点での関数値が分布に従うため、予測点での関数値の不確実性を考慮することができる。

6.3. カーネルとハイパーパラメータ#

カーネル関数 \(k(x, x')\) は、関数の滑らかさや周期性などの事前知識を反映する。代表的なカーネルは以下の通りである：

RBF（ガウス）カーネル:

\[ k(x, x') = \tau \exp\left(-\frac{(x-x')^2}{2l^2}\right) \]

ハイパーパラメータ：\(\tau\)（出力のスケール）、\(l\)（長さ尺度, 相関長:correlation lengthとも言う）
Matérnカーネル:

\[ k(x, x') = \tau \frac{2^{1-\nu}}{\Gamma(\nu)} \left( \sqrt{2\nu} \frac{|x-x'|}{l} \right)^\nu K_\nu \left( \sqrt{2\nu} \frac{|x-x'|}{l} \right) \]

ここで、\(K_\nu\) は第\(\nu\)種の変形ベッセル関数、\(\Gamma(\nu)\) はガンマ関数である。

ハイパーパラメータ：\(\tau\)（出力のスケール）、\(l\)（長さ尺度）、\(\nu\)（平滑度）

なかでも\(\nu=3/2\)や\(\nu=5/2\)がよく使われる。\(\nu=1/2\)のときは、指数カーネルとも呼ばれ、\(\nu \to \infty\)のときはRBFカーネルとなる。 RBFカーネルは無限回微分可能で非常に滑らかであるのに対し(それはデータが持つ連続性としては過剰な仮定と考えられ)、Matérnカーネルは\(\nu\)によって微分可能な回数が制御される。実用上においても、RBFカーネルで距離の近いデータ点を取りすぎると、数値的に不安定になることがある。

ハイパーパラメータは、対数尤度最大化などで推定されることが多い。上の、\(\tau, l\)をまとめて \(\theta = (\tau, l)\) と書くことにして、対数尤度の勾配は以下のように計算される。

\[ \frac{\partial}{\partial \theta_i} \log p(\mathbf{y}|\mathbf{x}, \theta) = \frac{1}{2} \mathbf{y}^\top K(\mathbf{x}, \mathbf{x})^{-1} \frac{\partial K(\mathbf{x}, \mathbf{x})}{\partial \theta_i} K(\mathbf{x}, \mathbf{x})^{-1} \mathbf{y} - \frac{1}{2} \mathrm{tr} \left(K(\mathbf{x}, \mathbf{x})^{-1} \frac{\partial K(\mathbf{x}, \mathbf{x})}{\partial \theta_i}\right) \]

また、カーネルが持つべき重要な性質としては、

対称性: \(k(x, x') = k(x', x)\)
正定値性: 任意の異なる点 \(\{x_1, x_2, \ldots, x_n\}\) に対して、対応するカーネル行列 \(K\) が正定値行列であること
非負性: \(k(x, x) \geq 0\) であること

が挙げられる。

6.4. ノイズの扱い#

ガウス過程では、観測データのノイズを自然に扱うことができる。

\[ y_i = f(x_i) + \epsilon_i, \quad \epsilon_i \sim \mathcal{N}(0, \sigma_n^2) \]

ここで、\(\epsilon_i\) は独立同分布に従うガウスノイズであるとする。このとき、観測データの関数値 \(\mathbf{y}\) の共分散行列は、以下のようにノイズの分散を加えたもの

\[ K(\mathbf{x}, \mathbf{x}) + \sigma_n^2 I \]

となる。ここで、\(I\) は単位行列である。このように、ノイズの分散を共分散行列に加えることで、ノイズを考慮した予測が可能となる。

一点注意すべきは、ノイズの分散 \(\sigma_n^2\) は、カーネルのハイパーパラメータとして推定されることが多いことである。既存のGPyなどのライブラリでは、データを与えた際に、ノイズの分散を含むハイパーパラメータを対数尤度最大化で推定する機能が備わっている。データは(多くの場合)標準化されていることが多いため、ノイズの分散は1よりもかなり小さい値になることが期待されるが、対数尤度最大化の結果、非常に大きな値になることがある。この場合、データが本来持つ不定性を過大評価してしまうため、予測が過度に平滑化されてしまうことがある。筆者の経験上、このあたりの問題を適切に理解した応用がなされていないものも多く散見するため、注意が必要である。

実装の背後にある理論や、データのドメイン知識が必要となる好例である。

6.5. 実装#

Pythonでは、GPyやscikit-learnなどのライブラリがガウス過程の実装を提供している。以下では、自前でMatérnカーネル(\(\nu=5/2\))を用いたガウス過程を実装してみる。

とりあえず、ハイパーパラメータを固定して、予測点での関数値の分布を計算してみよう。

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
import copy
from scipy import special
from matplotlib.animation import FuncAnimation
from IPython.display import HTML


## ガウス過程のカーネル(共分散行列)の設計
def Mat52(Thetas,r):
    tau,sigma = Thetas
    thetar = r * np.sqrt(5.0)/sigma
    return tau * (1.0 + thetar + (thetar**2) /3.0) * np.exp(-thetar)

def KernelMat(Thetas,xt,xp):
    lt = len(xt); lp=len(xp)
    Ktt = np.zeros((lt,lt)); Kpt = np.zeros((lp,lt));  Kpp = np.zeros((lp,lp))
    for j in range(lt):
        for i in range(j,lt):
            r = abs(xt[i]-xt[j])
            tmp  = Mat52(Thetas,r) 
            Ktt[i,j] = tmp; Ktt[j,i] = tmp 
        for i in range(lp):
            r= abs(xp[i]-xt[j])
            Kpt[i,j] = Mat52(Thetas,r)
    for j in range(lp):
        for i in range(j,lp):
            r= abs(xp[i]-xp[j])
            tmp  = Mat52(Thetas,r)
            Kpp[i,j] = tmp; Kpp[j,i] = tmp
    return Ktt,Kpt,Kpp

## 事後共分散行列の計算
def calcSj(cLinv,Kpt,Kpp,yt,mu_yt,mu_yp):
    tKtp= np.dot(cLinv,Kpt.T)
    return mu_yp + np.dot(Kpt,np.dot(cLinv.T,np.dot(cLinv,yt-mu_yt))), Kpp - np.dot(tKtp.T,tKtp)

## Cholesky分解
def Mchole(tmpA,ln) :
    cLL = np.linalg.cholesky(tmpA)
    logLii=0.0
    for i in range(ln):
        logLii += np.log(cLL[i,i])
    return np.linalg.inv(cLL), 2.0*logLii

## データ点＆予測点の設定
np.random.seed(0)
x = np.linspace(2, 8, 10)
y = f(x)
y = (np.array(y) - np.mean(y)) / np.std(y)  # 標準化
xp = np.linspace(0, 10, 200)

## ガウス過程の実装
for Thetas in [[1.0, 0.5, 1e-1], [1.0, 1.0, 1e-2]]:
    lt=len(x); lp=len(xp)
    Ktt,Kpt,Kpp = KernelMat(Thetas[:2], x, xp)
 
    mu_yt= np.mean(y)*np.ones(lt)
    mu_yp= np.mean(y)*np.ones(lp)

    cLinv,logdetK = Mchole(Ktt,lt) 
    mujoint,Sjoint = calcSj(cLinv, Kpt, Kpp, y, mu_yt, mu_yp)
    sigmaj=[ Sjoint[j][j] for j in range(lp)]
    ysamples = [np.random.multivariate_normal(mujoint,Sjoint) for i in range(10)]

    fig = plt.figure(figsize=(8, 3))
    plt.scatter(x, y, color="k", label="data")
    plt.plot(xp, mujoint, color="b", label="mean")
    plt.fill_between(xp, mujoint-1.96*np.sqrt(sigmaj), mujoint+1.96*np.sqrt(sigmaj), color="b", alpha=0.2, label="95% CI")
    for ys in ysamples:
        plt.plot(xp, ys, alpha=0.3, linestyle="--")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.title(f"GP regression with Mat52 kernel (tau={Thetas[0]:6.2f}, sigma={Thetas[1]})")
    plt.legend()
    plt.show()
    plt.close()

../_images/600c9ba899c4b229c02538502f1eff3fc848ae855f07300865000ef37226d421.png

../_images/57d6e0e48569dd620bb544f660835d04373a992d1fb7bd3a1ced86a9465ff934.png

青線が、事後分布の平均、青色のバンドが、95%信頼区間を示している。破線で示したのは、事後確率から得られたサンプルである。

\(d\)個のデータ点と相関した\(m\)個の予測点を同時にサンプルすれば、各点で独立した値ではなく、連続的な関数のサンプルを得ることができる。

また\(K(\mathbf{x}, \mathbf{x})\)の逆行列を計算する際に、数値的に不安定になる。上の実装では、コレスキー分解を用いた。

\[ K(\mathbf{x}, \mathbf{x}) = LL^T \]

ここで、\(L\)は下三角行列である。これを用いると、逆行列の計算は以下のように安定に行うことができる。

\[ K(\mathbf{x}, \mathbf{x})^{-1} = (L^T)^{-1} L^{-1} \]

念の為、同じハイパーパラメータで、GPyのガウス過程の実装と比較してみよう。

../_images/8544eea230ca032b08d8956c93ddb9dfe5c3937e9e59e0c2e287757b580bc6cc.png

両者がピッタリ重なっている事がわかる。ただし、GPyやscikit-learnなどのライブラリでは、データの分散や予測分散を別々にカスタマイズするのがやや面倒だったりもする。

6.6. ベイズ最適化#

ガウス過程は、ベイズ最適化の核としても広く使われる。関数の振る舞いをガウス過程でモデル化し、その不確実性を考慮しながら最適化を進めていく。

とくに、１回の観測(や数値実験など)にコストがかかり、勾配が利用できないような関数の最適化に有効である。また、探索と活用のトレードオフを考慮しながら、効率的に最適解に近づくことができることも多く機械学習のハイパーパラメータチューニングなどにも利用される。探索というのは「井の中の蛙状態が怖いので、未知の領域を探索してみよう」というイメージで、活用というのは「今までの経験から良さそうなところを掘り下げてみよう」というイメージである。

例えば今の関数の振る舞いを、パラメータ空間に対するコスト関数だと考え、その最小化や最大化を目指すとする。すでに得られているデータ点をもとに、次に評価する点を獲得関数で決定する。言い換えると、情報利得が高そうかつ、良さそうな点を次に評価するための、戦略を立てるための枠組みとも言える。

以下に、実装例を示す。なお、獲得関数には、期待改善量(Expected Improvement, EI)を用いている。

\[ EI(x) = (\mu(x) - y_\mathrm{best}) \Phi(z) + \sigma(x) \phi(z), \quad z = \frac{\mu(x) - y_\mathrm{best}}{\sigma(x)} \]

ここで\(\Phi, \phi\)はそれぞれ標準正規分布の累積分布関数と確率密度関数である。