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戦略

JAX

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

バイナリー・クロスエントロピー

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

NeRFを用いた点群生成の試み

Google Vision APIでOCR

アーカイブ

カテゴリー