「Pythonではじめるベイズ機械学習入門」を読む　～第3章回帰モデル⑤(ガウス過程回帰モデル：ガウス尤度)～

はじめに

Pythonによるベイズモデルの実装をきちんと学ぼうと思い，森賀新・木田悠歩・須山敦志著「Pythonではじめるベイズ機械学習入門」を読むことにした。

本記事は，第3章「回帰モデル」のうちガウス過程回帰(ガウス尤度)に関する読書メモである。

本書の紹介ページ

www.kspub.co.jp

3.6 ガウス過程回帰モデル：ガウス尤度

本節では，ノンパラメトリック回帰モデルと呼ばれる非常に柔軟な回帰モデルであるガウス過程回帰(Gaussian Process Regression Model)について説明している。

3.6.1. モデル概要

ガウス過程(Gaussian Process)の定義は以下の通りである。

ある確率変数の集合を $F$ とする。任意の自然数 $N$ に対して， $f$ を関数として $F$ から選んだ $N$ 個の確率変数 $\{ f(x_1), f(x_2),\cdots, f(x_N) \}$ がガウス分布に従うとき， $F$ をガウス過程と呼び， $f$ はガウス過程に従うという。

$\begin{align} f \sim \mathrm{GP}(m(x), k(x, x')) \\ \end{align}$

つまりガウス過程は，関数からなるベクトル $f =(f(x_1), f(x_2),\cdots, f(x_N))$ を出力するシステムとみなすことができる。

カーネル関数

ガウス過程の直感的な性質として，「入力 $x$ が近ければ出力 $y$ も近くなる」というものがある。この近さを表現するために，カーネル関数を用いる。
カーネル関数 $k(x, x')$ は，入力 $x$ と $x'$ の距離を表現するものである。代表的なものとしてRBFカーネル(Radial Basis Function Kernel; 動径基底関数)がある。

$\begin{align} k(x, x'; \theta) = \exp \left( - \frac{ |x-x'|^2 }{ \theta } \right) \end{align}$

先述の通りガウス過程は，平均関数 $m(x)$ とカーネル関数 $k(x, x')$ によって，生成される関数 $f$ の性質が決まる。データの平均値を0に正規化しておくことで，平均関数はゼロベクトルと仮定できるので，ガウス過程の性質はカーネル関数によって特徴づけられることになる。

具体的な入力値の集合 $\mathrm{X}=\{x_1, \cdots, x_N \}$ を用意すると，これらの入力値を用いて，次のような共分散行列

$\begin{align} K = \begin{bmatrix} k(x_1, x_1) & \cdots & k(x_1, x_N) \\ \vdots & \ddots & \vdots \\ k(x_N, x_1) & \cdots & k(x_N, x_N) \\ \end{bmatrix} \tag{1} \end{align}$

を計算すると，ガウス過程は

$\begin{align} f &\sim \mathrm{GP}(m(x), k(x, x')) \\ \\ \Rightarrow f &\sim \mathcal{N}(0, K) \\ \end{align}$

のように $N$ 次元の多次元ガウス分布 $\mathcal{N}(0, K)$ できる。また $f$ は， $\mathcal{N}(0, K)$ から $N$ 次元ベクトルをサンプリングしたものとして表現できる。

なお，本書の参考文献( ベイズ深層学習 | 書籍情報 | 株式会社講談社サイエンティフィク )によると，入力ベクトルの集合 $\mathrm{X}=\{x_1, \cdots, x_N \}$ から2つの要素を取り出してカーネル関数の値を算出し行列にしたものを，特に $K_{\mathrm{X}\mathrm{X}}$ と表現している。
すなわち(1)式においては， $K = K_{\mathrm{X}\mathrm{X}}$ である。

カーネル関数の組合せ

カーネル関数には，RBFカーネル以外に，

Linearカーネル : 線形関係
Polynomialカーネル : 多項式
Periodicカーネル : 周期性

など，近さの関係性ごとに様々なカーネルが存在する。

また，2種類のカーネル関数 $k_1, k_2$ の和や積もまたカーネル関数になるので，新たなカーネル関数を設計することもできる。

データにもとづき，複数のカーネル関数の候補から適切なカーネル関数を選択するには，対数周辺尤度(log marginal likelihood)の比較を行なう。

入力 $x_i$ と，これに対応する出力 $y_i$ の組合せが $N$ 個存在し， $\mathrm{X} = \{x_1, \cdots, x_N \}$ および $\mathrm{Y} = (y_1, \cdots, y_N)^T$ と表すと，周辺対数尤度は，

$\begin{align} \ln (p(\mathrm{Y} | \mathrm{X}, \theta)) = \ln (\mathcal{N} (0, K_{\mathrm{X}\mathrm{X}})) \\ \end{align}$

となる。

ハイパーパラメータ $\theta$ も選択する必要があるが，データに対する過剰適合を防ぐために，ハイパーパラメータに対しても事前分布を設定し，ハイパーパラメータの事後分布を同時に推論することで，過剰適合を抑制できる。

推論

ガウス過程回帰では，新たなデータに対する予測分布の出力が可能である。
すなわち，学習データ $D = \{ \mathrm{X} , \mathrm{Y} \}$ を用いて学習を行なった後に，新規の入力値 $\mathrm{X_{\mathrm{new}}}$ が得られたときの予測値 $\mathrm{Y_{\mathrm{new}}}$ の予測分布を算出することができる。

問題設定を簡単にするために，入力値は1次元で $x^*$ で，これに対応する予測値は $y^*$ とする。
ガウス過程回帰の学習を行なうと， $\mathrm{Y} \sim \mathcal{N}(0, K_{\mathrm{X}\mathrm{X}} )$ のように表現できる。

予測分布を算出するには $y^*$ が従う分布を求めればよいのだが，これを求めるために $\mathrm{Y}' = (\mathrm{Y}^T, y^*)^T$ という $(N+1)$ 次元のベクトルを導入する。また， $\mathrm{X} = \{x_1, \cdots, x_N \}$ および $x^*$ から算出される $(N+1) \times (N+1)$ 行列を $K'$ とすると，ガウス分布の性質より $\mathrm{Y}'$ もまたガウス分布に従い，

$\begin{align} \mathrm{Y}' \sim \mathcal{N} (0, K') \\ \end{align}$

となる。すなわち，

$\begin{align} \begin{bmatrix} \mathrm{Y} \\ y^* \end{bmatrix} \sim \mathcal{N} \left( 0, \begin{bmatrix} K & k_* \\ k_*^T & k_{**} \end{bmatrix} \right) \end{align}$

となる。ただし，

$\begin{align} k_* &= (k(x^*, x_1), \cdots, k(x^*, x_N)^T) \\ k_{**} &= k(x^*, x^*) \\ \end{align}$

である。

ここで，ガウス分布の条件付き分布の公式を用いる。ガウス分布に従うベクトルが，

$\begin{align} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} \sim \mathcal{N} \left( \begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix} \begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{bmatrix} \right) \end{align}$

のように分割されているとすると，条件付分布 $p(y_2 | y_1)$ は，

$\begin{align} p(y_2 | y_1) \sim \mathcal{N} (\mu_2 + \Sigma_{21} \Sigma_{11}^{-1} (y_1 - \mu_1), \Sigma_{22} - \Sigma_{21} \Sigma_{11}^{-1} \Sigma_{12}) \\ \end{align}$

であるので，この式を用いるとガウス過程分布の予測分布は以下のようになる。

$\begin{align} p(y^* | x^* , D) &= \mathcal{N} (y^* | \mu^*, \Sigma^*) \\ \\ \mu^* &= k_*^T K_{\mathrm{X} \mathrm{X} }^{-1} \mathrm{Y} \\ \\ \Sigma^* &= k_{**} - k_*^T K_{\mathrm{X} \mathrm{X} }^{-1} k_* \\ \end{align}$

予測の入力が複数ある場合でも考え方は同じで，結果は本書のP130の通りとなる。

尤度関数

各関数の出力 $f_n = f(x_n)$ に対して独立なノイズが付与されることによって，データ $y_n$ が観測されていると仮定する。今回は観測分布としてガウス分布を仮定する。

$\begin{align} y_n = f_n + \varepsilon_n, \quad \varepsilon_n \sim \mathcal{N}(0, \sigma_y^2) \\ \\ \end{align}$

このとき，次のような尤度関数を定義していることになる。

$\begin{align} p(y_n | f_n) = \mathcal{N}(y_n | f_n, \sigma_y^2) = \frac{1}{\sqrt{2 \pi \sigma_y^2}} \exp \left( -\frac{ (y_n - f_n)^2 }{ 2 \sigma_y^2} \right) \\ \end{align}$

3.6.2 実装

サンプルコードを動かしながら，挙動を確認した。
github.com

データの準備

データは下図に示すような人工データを用いる。

GPyTorchの特徴

実装にはPyTorchおよびGPyTorchを用いる。

ガウス過程の最大の問題点は，予測分布の計算において $N \times N$ サイズの行列の逆行列計算が含まれることである。これにより計算量のオーダーが $\mathcal{O}(N^3)$ になる。
GPyTorchでは，blackbox matrix-matrix multiplication (BBMM)推論を用いることで，計算量を実質的に $\mathcal{O}(N^2)$ に減らしている。

モデルの定義

ガウス過程回帰のモデルは，以下の様に実装される。

サンプルコード

# ガウス過程回帰モデルの実装
class ExactGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)
        # 平均関数
        self.mean_module = gpytorch.means.ConstantMean()
        # カーネル関数
        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())

    # ガウス過程の生成過程
    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

# 尤度にガウス分布を設定
likelihood = gpytorch.likelihoods.GaussianLikelihood()
# モデルのインスタンス化
model = ExactGPModel(x_data, y_data, likelihood)

まずモデルの構造として，mean_moduleとcovar_moduleを設定している。
サンプルコードではRBFカーネルを用いているが，これ以外のカーネル関数を用いる場合はcovar_moduleの定義を変更すればよい。

カーネル関数に含まれるハイパーパラメータは，勾配法によって求める。損失関数は負の周辺対数尤度である。学習ステップごとの損失関数の推移は下図のようになる。

学習の後，予測分布を算出する。予測の入力は，-1から1までを50等分した値であり，予測結果は下図のようになる。

学習データが少ない領域の分散は大きくなり，予測の不確実性が高くなっていることが分かる。

まとめと感想

ノンパラメトリック回帰の1つである，ガウス過程回帰モデルについて学んだ。

本書は理論と実装の両方を紹介している良い本なのだが，紙面の関係でどうしても理論部分の説明を他書に譲る必要が出てくる。ただ，参考文献もきちんと書いてくれているので，行間はガウス過程と機械学習 | 書籍情報 | 株式会社講談社サイエンティフィクやベイズ深層学習 | 書籍情報 | 株式会社講談社サイエンティフィクで補うことができた。

ガウス過程は，定義がとっつきにくかったが，いくらか具体例を紹介していただけると理解できた。特に予測分布は，多変量正規分布の条件付き分布のよい復習になった。

機械学習に初めて触れた2000年代前半は，カーネル法やSVMが最盛期であり，当時からカーネル行列の逆行列計算やその近似計算については議論になっていた。アルゴリズムや実装の面で工夫がされてきているので，しっかり実践していきたい。

本記事を最後まで読んでくださり，どうもありがとうございました。