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$ の正規分布を表す。これを実装したのが以下である。

# Training data is 100 points in [0,1] inclusive regularly spaced
train_x = torch.linspace(0, 1, 100)
# True function is sin(2*pi*x) with Gaussian noise
train_y = torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * 0.2

# Training data is 100 points in [0,1] inclusive regularly spaced

train_x = torch.linspace(0, 1, 100)

# True function is sin(2*pi*x) with Gaussian noise

train_y = torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * 0.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*}$

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

# We will use the simplest form of GP model, exact inference
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)

# We will use the simplest form of GP model, exact inference

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)

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で定義される。

# initialize likelihood and model
likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = ExactGPModel(train_x, train_y, likelihood)

# initialize likelihood and model

likelihood = gpytorch.likelihoods.GaussianLikelihood()

model = ExactGPModel(train_x, train_y, likelihood)

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

最適化

尤度の対数をとったもの

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

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

# this is for running the notebook in our testing framework
import os
smoke_test = ('CI' in os.environ)
training_iter = 2 if smoke_test else 50

# Find optimal model hyperparameters
model.train()
likelihood.train()

# Use the adam optimizer
optimizer = torch.optim.Adam([
    {'params': model.parameters()},  # Includes GaussianLikelihood parameters
], lr=0.1)

# "Loss" for GPs - the marginal log likelihood
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

for i in range(training_iter):
    # Zero gradients from previous iteration
    optimizer.zero_grad()
    # Output from model
    output = model(train_x)
    # Calc loss and backprop gradients
    loss = -mll(output, train_y)
    loss.backward()
    print('Iter %d/%d - Loss: %.3f   lengthscale: %.3f   noise: %.3f' % (
        i + 1, training_iter, loss.item(),
        model.covar_module.base_kernel.lengthscale.item(),
        model.likelihood.noise.item()
    ))
    optimizer.step()

# this is for running the notebook in our testing framework

import os

smoke_test = ('CI' in os.environ)

training_iter = 2 if smoke_test else 50

# Find optimal model hyperparameters

model.train()

likelihood.train()

# Use the adam optimizer

optimizer = torch.optim.Adam([

{'params': model.parameters()}, # Includes GaussianLikelihood parameters

], lr=0.1)

# "Loss" for GPs - the marginal log likelihood

mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

for i in range(training_iter):

# Zero gradients from previous iteration

optimizer.zero_grad()

# Output from model

output = model(train_x)

# Calc loss and backprop gradients

loss = -mll(output, train_y)

loss.backward()

print('Iter %d/%d - Loss: %.3f lengthscale: %.3f noise: %.3f' % (

i + 1, training_iter, loss.item(),

model.covar_module.base_kernel.lengthscale.item(),

model.likelihood.noise.item()

))

optimizer.step()

7,8行目でmodelとlikelihoodをそれぞれ訓練モードに設定している。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}_{**}$ は対角行列である。これらを計算するのが次のコードである。

# Get into evaluation (predictive posterior) mode
model.eval()
likelihood.eval()

# Test points are regularly spaced along [0,1]
# Make predictions by feeding model through likelihood
with torch.no_grad(), gpytorch.settings.fast_pred_var():
    test_x = torch.linspace(0, 1, 51)
    observed_pred = likelihood(model(test_x))

# Get into evaluation (predictive posterior) mode

model.eval()

likelihood.eval()

# Test points are regularly spaced along [0,1]

# Make predictions by feeding model through likelihood

with torch.no_grad(), gpytorch.settings.fast_pred_var():

test_x = torch.linspace(0, 1, 51)

observed_pred = likelihood(model(test_x))

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

plt.rcParams["font.size"] = 18
plt.figure(figsize=(15,10))
with torch.no_grad():
    # Get upper and lower confidence bounds
    lower, upper = observed_pred.confidence_region()
    # Plot training data as black stars
    plt.plot(train_x.numpy(), train_y.numpy(), 'k*')
    # Plot predictive means as blue line
    plt.plot(test_x.numpy(), observed_pred.mean.numpy(), 'b')
    # Shade between the lower and upper confidence bounds
    plt.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)
    plt.ylim([-3, 3])
    plt.legend(['Observed Data', 'Mean', 'Confidence'])
    plt.show()

plt.rcParams["font.size"] = 18

plt.figure(figsize=(15,10))

with torch.no_grad():

# Get upper and lower confidence bounds

lower, upper = observed_pred.confidence_region()

# Plot training data as black stars

plt.plot(train_x.numpy(), train_y.numpy(), 'k*')

# Plot predictive means as blue line

plt.plot(test_x.numpy(), observed_pred.mean.numpy(), 'b')

# Shade between the lower and upper confidence bounds

plt.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)

plt.ylim([-3, 3])

plt.legend(['Observed Data', 'Mean', 'Confidence'])

plt.show()

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

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

まとめ

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

参考文献

ガウス過程と機械学習：とても良い本である。ベイズ推定による機械学習もおすすめ。

Kumada Seiya

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

「御社が第1志望」と伝えるべきなのか

サラリーマンはオワコンなのか！？

GPyTorch入門

はじめに

観測値の作成

ガウス過程の導入

尤度の導入

最適化

予測

まとめ

参考文献

Kumada Seiya

最近の記事

顔認証（DeepFace vs GPT4o）

時系列解析再訪

πの新しい公式

中心極限定理

DBSCANを用いたクラスタリング

レイアウト解析

GANのminmax戦略

PyTorch C++の導入〜その3〜

テンソルネットワークの入り口（その2）

ChatGPTがもたらす哲学の再構築

Variational Autoencoderのグレイ画像への適用

PyTorch C++の導入〜その2〜

Google Vision APIでOCR

アーカイブ

カテゴリー