GPyTorch入門

はじめに

今回は、ガウス過程のためのPythonライブラリ「GPyTorch」を紹介したい。このページに簡単なチュートリアルがある。このチュートリアルを詳細な理論式で補いながら説明したい。ガウス過程自体の説明は本ブログ(ここここ)でもしているので、興味があれば参照して欲しい。

観測値の作成

観測値(訓練データ)を次式を用いてN(=100)個作成する。

(1)    \begin{eqnarray*} y_n&=&f(x_n)+\epsilon\nonumber \\ f(x_n)&=&\sin{\left(2\pi x_n\right)}\nonumber \\ \epsilon&\sim&\mathcal{N}(\epsilon|0,\sigma^2)\nonumber \\ \sigma&=&0.2 \end{eqnarray*}

\mathcal{N}(\epsilon|0,\sigma^2)は平均0、分散\sigma^2の正規分布を表す。これを実装したのが以下である。

ガウス過程の導入

いま、{\bf x}=\left(x_1,\cdots,x_N\right)\in\mathbb{R}^N {\bf y}=\left(y_1,\cdots,y_N\right)\in\mathbb{R}^Nとおき、ベクトル{\bf f}=\left(f(x_1),\cdots,f(x_N)\right)が次のガウス過程で作られると仮定する。

(2)    \begin{equation*} {\bf f}\sim p({\bf f}|{\bf x})=\mathcal{N}\left({\bf f}|\boldsymbol{\mu}\left({\bf x}\right),\boldsymbol{K}\left({\bf x}\right)\right) \end{equation*}

ここで、\boldsymbol{\mu}\left({\bf x}\right)\in\mathbb{R}^N\boldsymbol{K}\left({\bf x}\right)\in\mathbb{R}^{N\times N}である。\boldsymbol{\mu}\left({\bf x}\right)は平均ベクトル、\boldsymbol{K}\left({\bf x}\right)は共分散行列である。\boldsymbol{K}の成分K_{nm}=k(x_n,x_m)として、次のRBFカーネルを仮定する。

(3)    \begin{equation*} k(x_n,x_m)=\theta_1\exp{\left(-\frac{\left(x_n-x_m\right)^2}{2\theta_2^2}\right)} \end{equation*}

これらを実装したのが次のコードである。

ExactGPModelは式(2)の正規分布を表すクラスである。コンストラクタ内で、平均ベクトルmean_moduleと共分散行列covar_moduleを定義している。ここでは平均ベクトルとして定数値を仮定する(\boldsymbol{\mu}\left({\bf x}\right)=\boldsymbol{\mu}_0)。クラスRBFKernelは式(3)の指数関数部分を表す。\theta_2は、このクラスが持つ変数lengthscaleに相当する。式(3)のパラメータ\theta_1は関数ScaleKernel(6行目)により付与される。関数forwardが返す値は、クラスMultivariateNormalのインスタンスである。

尤度の導入

{\bf f}が式(2)から生成されるとき、{\bf y}を生成する分布(尤度)は次式から計算される。

(4)    \begin{eqnarray*} p({\bf y}|{\bf x}) &=& \int d{\bf f} p({\bf y}|{\bf f})p({\bf f}|{\bf x}) \nonumber \\ &=& \int d{\bf f} \mathcal{N}({\bf y}|{\bf f},\sigma^2)\;\mathcal{N}\left({\bf f}|\boldsymbol{\mu}\left({\bf x}\right),\boldsymbol{K}\left({\bf x}\right)\right)\nonumber \\ &=& \mathcal{N}\left({\bf y}|\boldsymbol{\mu}\left({\bf x}\right),\boldsymbol{K}\left({\bf x}\right)+\sigma^2\boldsymbol{I}\right) \end{eqnarray*}

\sigma^2=\theta_3として、\boldsymbol{\theta}=(\theta_1,\theta_2,\theta_3)とおき、パラメータ依存性を明示的に書くと、p({\bf y}|{\bf x})p({\bf y}|{\bf x},\boldsymbol{\theta})と書くことができる。p({\bf y}|{\bf x},\boldsymbol{\theta})は以下のGaussianLikelihoodで定義される。

ExactGPModelを使う場合、その引数likelihoodGaussianLikelihoodでなければならない。

最適化

尤度の対数をとったもの

(5)    \begin{equation*} \ln{p({\bf y}|{\bf x},\boldsymbol{\theta})} \end{equation*}

を最大にするような\boldsymbol{\theta}を勾配降下法により求める。これを実現するコードが以下である。

7,8行目でmodellikelihoodをそれぞれ訓練モードに設定している。PyTorchでお馴染みの手順である。11行目で最適化器Adamを定義する。16行目にあるExactMarginalLogLikelihoodは式(5)に相当する。24行目のloss

(6)    \begin{equation*} -\ln{p({\bf y}|{\bf x},\boldsymbol{\theta})} \end{equation*}

の値である。これを最小にするような計算が行われている。ループの中身は、PyTorchを用いたニューラルネットワークの最適化手順と同じである。PyTorchのクラスAdamが使われていることに注意する。上記コードを実行すると以下の出力を得る。

Iter 1/50 – Loss: 0.947 lengthscale: 0.693 noise: 0.693
Iter 2/50 – Loss: 0.916 lengthscale: 0.644 noise: 0.644
Iter 3/50 – Loss: 0.882 lengthscale: 0.598 noise: 0.598
Iter 4/50 – Loss: 0.844 lengthscale: 0.555 noise: 0.554
Iter 5/50 – Loss: 0.801 lengthscale: 0.514 noise: 0.513
Iter 6/50 – Loss: 0.752 lengthscale: 0.476 noise: 0.474
Iter 7/50 – Loss: 0.699 lengthscale: 0.439 noise: 0.437

Iter 47/50 – Loss: -0.064 lengthscale: 0.287 noise: 0.029
Iter 48/50 – Loss: -0.066 lengthscale: 0.284 noise: 0.030
Iter 49/50 – Loss: -0.069 lengthscale: 0.282 noise: 0.030
Iter 50/50 – Loss: -0.071 lengthscale: 0.281 noise: 0.031

損失(Loss)とlengthscale(\theta_2)、noise(\theta_3)が最適値へ向かう様子が出力される。

予測

予測したいM(=51)個の位置を{\bf x}^*=(x_1^*,\cdots,x_M^*)とおく。対応するyの値{\bf y}^*=(y_1^*,\cdots,y_M^*)は次式で与えられる。

(7)    \begin{equation*} {\bf y}^*\sim \mathcal{N}({\bf y}^*|\boldsymbol{K}_*^T\boldsymbol{K}^{-1}{\bf y},\boldsymbol{K}_{**}-\boldsymbol{K}_*^T\boldsymbol{K}^{-1}\boldsymbol{K}_*) \end{equation*}

ここで、\boldsymbol{K}+\sigma^2\boldsymbol{I}を改めて\boldsymbol{K}と置いてある。行列\boldsymbol{K}_*\boldsymbol{K}_{**}の成分は次式で定義される。

(8)    \begin{eqnarray*} k_*(x_n,x_m)&=&k(x_n,x_m^*)\;\;(n=1,\cdots,N, m=1,\cdots,M) \nonumber \\ k_{**}(x_m,x_m)&=&k(x_m^*,x_m^*)\;\;(m=1,\cdots,M)\nonumber \end{eqnarray*}

行列\boldsymbol{K}_{**}は対角行列である。これらを計算するのが次のコードである。

2行目、3行目で評価モードに切り替えている。9行目で式(7)を計算している(MultivariateNormalのインスタンスが作られる)。次のコードでグラフを描く。

これを実行すると次の画像を得る。

黒の星印は観測値(訓練データ)を、青の曲線は式(7)の平均ベクトル\boldsymbol{K}_*^T\boldsymbol{K}^{-1}{\bf y}を、淡い青の領域は標準偏差の下限と上限を表す。

まとめ

今回は、ガウス過程のためのPythonライブラリ「GPyTorch」のチュートリアルを詳細な理論式で補いながら説明した。深層学習で使われるPyTorchと統一的なコーディングができることが大きなメリットである。

参考文献

Kumada Seiya

Kumada Seiya

仕事であろうとなかろうと勉強し続ける、その結果”中身”を知ったエンジニアになれる

最近の記事

  • 関連記事
  • おすすめ記事
  • 特集記事

アーカイブ

カテゴリー

PAGE TOP