甲斐性なしのブログ

うどんくらいしか食べる気がしない

ガウス過程回帰で学習データにない任意の入力点に対する関数出力値をMCMCサンプリングする方法

ガウス過程 ベイズ推論 機械学習

はじめに

表題についてちょっと調べた＆やってみたのでメモ。 やりたいことは表題の通りなんだけど、少し数式で表現して整理。

  • 今、学習データとして \left\{{\bf X}, {\bf y} \right\} = \left\{({\bf x}_1, y_1), \cdots, ({\bf x}_N, y_N)\right\}のN個得られているとする。
  • 入力{\bf x}と関数fが与えられたとき、出力yはp(y | f({\bf x}))に従うとする（ハイパーパラメータを持つ場合もあるが今回は固定値として扱うので表記から省略）。
  • そして、この関数fはガウス過程p(f)={\mathcal gp}(0, k({\bf x}, {\bf x}'))に従う確率変数とする（ k({\bf x}, {\bf x}')は共分散関数）。
    • p(y | f({\bf x}))がガウス分布でないなら、関数の事後分布p(f|{\bf X}, {\bf y}) \propto \prod_{i = 1}^{N}p( y_i | f({\bf x}_i))p(f)は解析的に得られないのでそういう場合は近似するかMCMC等でサンプリングを行う。
  • やりたいのは、M個の任意の未知入力\left\{{\bar {\bf x}}_1, \cdots, {\bar{\bf x}}_M \right\}が与えられて、事後分布からサンプリングされた関数fの出力\left\{f({\bar {\bf x}}_1), \cdots, f({\bar{\bf x}}_M) \right\}を得たい。
    • 予測分布からのサンプリングも得たいという話もあるけど今回は対象外。
      • とりあえず、事後分布からの出力値f({\bar {\bf x}})が得られれば、これからp(y |f({\bar {\bf x}}))が計算できるので、それっぽいyは得られそう。

ようするに、学習入力点N個に対する関数fの出力を並べたベクトルを{\bf f}_N = \left[f({\bf x})_1 \cdots, f({\bf x}_N)\ \right]^Tとし、未知入力点M個に対する関数fの出力を並べたベクトルを{\bar {\bf f}}_M = \left[f({\bar {\bf x}})_1 \cdots, f({\bar {\bf x}}_M)\ \right]^Tとした場合、

\displaystyle p(f|{\bf X}, {\bf y}) = p({\bf f}_N | {\bf y}) \propto p({\bf y}| {\bf f}_N) p({\bf f}_N)

からのMCMCサンプリングはできるのは分かるけど、 p({\bar {\bf f}}_M | {\bf y}) からのサンプリングはどうやったらいいんだ？って話。

方法1 ~MCMCサンプリングの時、同時分布 p({\bf f}_N, {\bar {\bf f}}_M )からサンプリング~

非常にお世話になっている書籍

ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)

ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)

のサポートページに、コーシー分布を観測モデルとしたガウス過程回帰のpythonコードgpr-cauchy.py, elliptical.pyがある。これは楕円スライスサンプリングというMCMCサンプリングの一手法を用いて関数出力値の推定を行っているが、これを見るとサンプルを得るところで

p({\bf f}_N, {\bar {\bf f}}_M ) = {\mathcal N} \left( {\bf f}_N, {\bar {\bf f}}_M| 0, \begin{pmatrix} {\bf K}_{NN} & {\bf K}_{NM} \\ {\bf K}_{NM}^T & {\bf K}_{MM} \\ \end{pmatrix} \right)

（ただし、{\bf K}_{NN}はK_{i,j} = , k({\bf x}_i, {\bf x}_j) \; (i,j \in \left\{1, \cdots, N \right\})という要素を持つN \times Nの行列、{\bf K}_{NM}はK_{i,j} = , k({\bf x}_i, {\bar {\bf x}}_j) \; (i \in \left\{1, \cdots, N \right\}, j \in \left\{1, \cdots, M \right\})という要素を持つN \times Mの行列、{\bf K}_{MM}はK_{i,j} = , k({\bar{\bf x}}_i, {\bar{\bf x}}_j) \; (i,j \in \left\{1, \cdots, M \right\})という要素を持つM \times Mの行列）

からのサンプルを使って更新している。そして、そのサンプルを採択するか否かの尤度計算はp({\bf y}| {\bf f}_N) p({\bf f}_N)で計算、つまり学習入力点だけで行っている*1。

なるほど、確かにこうすればp({\bf f}_N, {\bar {\bf f}}_M | {\bf y}) からサンプリングしたことになり、そのサンプルから{\bf f}_Nを捨てれば周辺化することになるからp({\bar {\bf f}}_M | {\bf y}) からサンプリングしたことになる。

方法2 ~ p({\bf f}_N | {\bf y}) からのサンプルを用いて p({\bar {\bf f}}_M | {\bf y}) からサンプリング~

方法1だと、学習時に評価したい未知入力\left\{{\bar {\bf x}}_1, \cdots, {\bar{\bf x}}_M \right\}がわかっている必要があるが、場合によってはそれが困難だったりする。そういう場合はどうするか？

p({\bf f}_N | {\bf y})からのサンプルがS個 {\bf f}_N^{(1)}, \cdots, {\bf f}_N^{(S)}得られているとすると、

\displaystyle p({\bar {\bf f}}_M | {\bf y}) = \int p({\bar {\bf f}}_M | {\bf f}_N) p({\bf f}_N | {\bf y})d {\bf f} _N \sim \frac{1}{S} \sum_{s = 1}^{S}p({\bar {\bf f}}_M | {\bf f}_N^{(s)})

と近似できる。ここでガウス過程の性質より、

p({\bar {\bf f}}_M | {\bf f}_N^{(s)}) = {\mathcal N} ({\bar {\bf f}}_M | {\bf K}_{NM}^T {\bf K}_{NN}^{-1} {\bf f}_N^{(s)} , {\bf K}_{MM} - {\bf K}_{NM}^T {\bf K}_{NN}^{-1} {\bf K}_{NM} )

なので、p({\bar {\bf f}}_M | {\bf y})は全て同じ混合比の混合正規分布で近似できる。従ってs \in \left\{1 \cdots S \right\}の乱数を発生させた後、それに紐づく多次元正規分布 {\mathcal N} ({\bar {\bf f}}_M | {\bf K}_{NM}^T {\bf K}_{NN}^{-1} {\bf f}_N^{(s)} , {\bf K}_{MM} - {\bf K}_{NM}^T {\bf K}_{NN}^{-1} {\bf K}_{NM} )からサンプリングすればよい。

実装

方法2によるコーシー分布ガウス過程回帰を上記サポートページにあるgpr-cauchy.pyを改変する形でやってみた。

import sys
import putil
import numpy as np
from pylab import *
from numpy import exp,log
from elliptical import elliptical
from numpy.linalg import cholesky as chol

def gpr_cauchy (f, param):
    x,y,gamma,Kinv = param
    M = len(x)
    return gpr_cauchy_lik (y[0:M], f[0:M], gamma, Kinv)

def gpr_cauchy_lik (y, f, gamma, Kinv):
    return - np.sum (log (gamma + (y - f)**2 / gamma)) \
           - np.dot (f, np.dot(Kinv, f)) / 2

def kgauss (tau,sigma):
    return lambda x,y: exp(tau) * exp (-(x - y)**2 / exp(sigma))

def kernel_matrix (xx, kernel):
    N = len(xx)
    eta = 1e-6
    return np.array (
        [kernel (xi, xj) for xi in xx for xj in xx]
    ).reshape(N,N) + eta * np.eye(N)

def gpr_mcmc(x, y, iters, xmin, xmax, gamma):
    
    N = len(x)
    K_NN = kernel_matrix(x, kgauss(1,1))
    K_NNinv = inv(K_NN)

    #こうすると多次元正規分布からサンプリングすることになる
    P = chol (K_NN)
    f = P @ randn(N)
    
    #MCMCサンプルのうち最初の200個は使わない
    remove_num = 200
    S = iters - remove_num
    f_posterior = np.zeros((f.shape[0], S))
    
    for iter in range(iters):
        #楕円スライスサンプリング
        f,lik = elliptical (f, P, gpr_cauchy, (x,y,gamma,K_NNinv))
        print ('\r[iter %2d]' % (iter + 1),)
        if iter >= remove_num:
            f_posterior[:, iter - remove_num] = f
    
    #未知入力点
    x_new = np.linspace(xmin, xmax, 100)
    M = len(x_new)
    xx = np.hstack((x, x_new))

    K = kernel_matrix(xx, kgauss(1, 1))
    K_NM = K[:N, N:N + M]
    K_MM = K[N:N + M, N:N + M]

    g = np.zeros (len(x_new))
    #50個サンプリング
    sampling_num = 50
    for t in range(sampling_num):
        s = np.random.randint(S)
        f_n = f_posterior[:, s]
        mu = K_NM.T@ K_NNinv @ f_n
        sig = K_MM - K_NM.T @ K_NNinv @ K_NM
        P = chol(sig)
        f_m = P @ np.random.randn(M) + mu
        g = g + f_m
        plot (x_new, f_m)
    plot (x, y, 'bx', markersize=14)
    plot (x_new, g/sampling_num, 'k', linewidth=3)

def usage ():
    print ('usage: gpr-cauchy.py data.xyf iters [output]')
    sys.exit (0)

def main ():
    xmin = -5
    xmax =  5
    ymin = -7.5
    ymax = 12.5
    gamma = 0.2
    
    if len(sys.argv) < 3:
        usage ()
    else:
        [x, y, f] = np.loadtxt (sys.argv[1]).T
        iters = int (sys.argv[2])

    gpr_mcmc (x, y, iters, xmin, xmax, gamma)
    axis ([xmin, xmax, ymin, ymax])
    show ()

if __name__ == "__main__":
    main ()

楕円スライスサンプリングを行う関数ellipticalについては、 サポートページのelliptical.pyから変えてないので省略。また多次元正規分布からはコレスキー分解を用いてサンプリングしている（この辺が参考になる）。これを、

python .\gpr-cauchy_new.py  .\gpr-cauchy.dat 1000

として、実行すると以下のような結果が得られた。

f:id:YamagenSakam:20200123235146p:plain

うん、いい感じ。

今後やりたいこと

  • ハイパーパラメータも事前分布与えてMCMCサンプリング
  • 予測分布からのサンプリング
  • ガウス過程分類
  • 変分近似

*1:事前分布の方はp({\bf f}_N, {\bar {\bf f}}_M)の同時分布を評価するべきかと思ったがどうなんだろ？実際両方やってみたけど、結果はあまり変わらなかった