GPyTorch 1.2 Tutorials : GPyTorch 回帰チュートリアル (翻訳/解説)
翻訳 : (株)クラスキャットセールスインフォメーション
作成日時 : 09/03/2020 (1.2)

* 本ページは、GPyTorch 1.2 ドキュメントの以下のページを翻訳した上で適宜、補足説明したものです：

GPyTorch Regression Tutorial

* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。

Tutorials : GPyTorch 回帰チュートリアル

イントロダクション

この notebook では、最も単純な例を使用して単純な関数上で RBF カーネル・ガウス過程を訓練し、GPyTorch の多くの設計の特徴を実演します。次の関数を 100 訓練サンプルでモデル化して :

$\begin{align} y &= \sin(2\pi x) + \epsilon \\ \epsilon &\sim \mathcal{N}(0, 0.04) \end{align}$

そして 51 テストサンプル上でテストしていきます。

Note: このノートブックはガウス過程の数学的背景を教えることを必ずしも意図していません、むしろ GPyTorch でどのように単純なものを訓練して予測を行なうかです。数学的扱いについては、Gaussian Processes for Machine Learning の 2 章が GP 回帰への非常に完全なイントロダクションを提供しています (このテキスト全体は高く推奨されます) :

http://www.gaussianprocess.org/gpml/chapters/RW2.pdf

import math
import torch
import gpytorch
from matplotlib import pyplot as plt
 
%matplotlib inline
%load_ext autoreload
%autoreload 2

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

訓練データのセットアップ

次のセルで、この例のために訓練データをセットアップします。[0, 1] 上の 100 の規則的な間隔のポイントを使用していきます、その上で関数を評価して訓練ラベルを得るためにガウスノイズを追加します。

# 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()) * math.sqrt(0.04)

モデルをセットアップする

次のセルは GPyTorch のユーザ定義ガウス過程モデルの最も重要な特徴を実演します。GPyTorch の GP モデルの構築は様々な方法で異なります。

まず、多くの既存の GP パッケージとは対照的に、ユーザのために full GP モデルは提供しません。むしろ、素早く一つを構築するために必要なツールを提供します。これは何故ならば、標準的な PyTorch でニューラルネットワークを構築することに類似して、どのようなコンポーネントを含む必要があるとしても柔軟性を持つことが重要だと私達は信じるからです。より複雑な例で見られるように、これはカスタムモデルの設計においてユーザに素晴らしい柔軟性を可能にします。

殆どの GP 回帰モデルについて、以下の GPyTorch オブジェクトを構築することが必要です :

GP モデル (gpytorch.models.ExactGP) – これは推論の殆どを扱います。
尤度 (gpytorch.likelihoods.GaussianLikelihood) – これは GP 回帰のために使用される最も一般的な尤度です。
平均 (= Mean) – これは GP の事前平均を定義します。(どの平均を使用するべきか知らない場合には、gpytorch.means.ConstantMean() で始めてください。)
カーネル – これは GP の 事前共分散 を定義します。(どのカーネルを使用するべきか知らない場合には、gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel()) で始めてください。)
MultivariateNormal 分布 (gpytorch.distributions.MultivariateNormal) – これは多変量正規分布を表わすために使用されるオブジェクトです。

GP モデル

GPyTorch のユーザ構築 (Exact, i.e. 非変分) GP モデルのコンポーネントは、大まかに言えば :

__init__ メソッド、これは訓練データと尤度を取り、そしてどのようなオブジェクトであれモデルの forward メソッドのために必要なものを構築します。これは通常は mean モジュールとカーネルモジュールのようなものを含みます。
forward メソッドは、ある $n \times d$ データ x を取りそして x で評価された事前平均と共分散を持つ MultivariateNormal を返します。
換言すれば、GP の事前平均と共分散行列を表わすベクトル $\mu(x)$ と $n \times n$ 行列 $K_{xx}$ を返します。

この仕様はモデルを定義するとき巨大な総量の柔軟性をそのままにしておきます。例えば、2 つのカーネルを加算を通じて組み立てるために、カーネルモジュールを直接加算することもできます :

1	`self.covar_module` `=` `ScaleKernel(RBFKernel()` `+` `WhiteNoiseKernel())`

あるいは forward メソッドでカーネルの出力を加算できます :

1	`covar_x` `=` `self.rbf_kernel_module(x)` `+` `self.white_noise_module(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)
 
# initialize likelihood and model
likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = ExactGPModel(train_x, train_y, likelihood)

モデル・モード

殆どの PyTorch モジュールのように、ExactGP は .train() と .eval() モードを持ちます。

.train() モードはモデルのハイパーパラメータを最適化するためです。
.eval() モードはモデル事後分布を通して予測を計算するためです。

モデルを訓練する

次のセルでは、ガウス過程のハイパーパラメータを訓練するために Type-II MLE を使用して処理します。

多くの他の GP 実装と比較してここでの最も明白な違いは、標準的な PyTorch でのように、コア訓練ループがユーザにより書かれることです。GPyTorch では、torch.optim からの標準的な PyTorch optimizer を利用し、そしてモデルの総ての訓練可能なパラメータはタイプ torch.nn.Parameter であるべきです。GP モデルは torch.nn.Module を直接拡張していますので、model.parameters() や model.named_parameters() 関数のようなメソッドへの呼び出しは (想像できるように) PyTorch 由来です。

殆どの場合、下のボイラープレートなコードは上手く動作します。それは標準的な PyTorch 訓練ループと同じ基本的なコンポーネントを持ちます :

総てのパラメータ勾配をゼロにします。
モデルを呼び出して損失を計算します。
勾配を満たすために損失上で backward を呼び出します。
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(model.parameters(), lr=0.1)  # Includes GaussianLikelihood parameters
 
# "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()

Iter 1/50 - Loss: 0.939   lengthscale: 0.693   noise: 0.693
Iter 2/50 - Loss: 0.908   lengthscale: 0.644   noise: 0.644
Iter 3/50 - Loss: 0.874   lengthscale: 0.598   noise: 0.598
Iter 4/50 - Loss: 0.837   lengthscale: 0.555   noise: 0.554
Iter 5/50 - Loss: 0.795   lengthscale: 0.514   noise: 0.513
Iter 6/50 - Loss: 0.749   lengthscale: 0.476   noise: 0.474
Iter 7/50 - Loss: 0.699   lengthscale: 0.440   noise: 0.437
Iter 8/50 - Loss: 0.649   lengthscale: 0.405   noise: 0.402
Iter 9/50 - Loss: 0.600   lengthscale: 0.372   noise: 0.369
Iter 10/50 - Loss: 0.556   lengthscale: 0.342   noise: 0.339
Iter 11/50 - Loss: 0.516   lengthscale: 0.315   noise: 0.310
Iter 12/50 - Loss: 0.480   lengthscale: 0.291   noise: 0.284
Iter 13/50 - Loss: 0.448   lengthscale: 0.270   noise: 0.259
Iter 14/50 - Loss: 0.413   lengthscale: 0.254   noise: 0.237
Iter 15/50 - Loss: 0.380   lengthscale: 0.241   noise: 0.216
Iter 16/50 - Loss: 0.355   lengthscale: 0.231   noise: 0.197
Iter 17/50 - Loss: 0.314   lengthscale: 0.223   noise: 0.179
Iter 18/50 - Loss: 0.292   lengthscale: 0.218   noise: 0.163
Iter 19/50 - Loss: 0.262   lengthscale: 0.214   noise: 0.148
Iter 20/50 - Loss: 0.236   lengthscale: 0.214   noise: 0.135
Iter 21/50 - Loss: 0.201   lengthscale: 0.216   noise: 0.122
Iter 22/50 - Loss: 0.176   lengthscale: 0.220   noise: 0.111
Iter 23/50 - Loss: 0.158   lengthscale: 0.224   noise: 0.102
Iter 24/50 - Loss: 0.125   lengthscale: 0.231   noise: 0.093
Iter 25/50 - Loss: 0.101   lengthscale: 0.239   noise: 0.085
Iter 26/50 - Loss: 0.078   lengthscale: 0.247   noise: 0.077
Iter 27/50 - Loss: 0.066   lengthscale: 0.256   noise: 0.071
Iter 28/50 - Loss: 0.052   lengthscale: 0.265   noise: 0.065
Iter 29/50 - Loss: 0.036   lengthscale: 0.276   noise: 0.060
Iter 30/50 - Loss: 0.036   lengthscale: 0.286   noise: 0.056
Iter 31/50 - Loss: 0.031   lengthscale: 0.297   noise: 0.052
Iter 32/50 - Loss: 0.028   lengthscale: 0.306   noise: 0.048
Iter 33/50 - Loss: 0.030   lengthscale: 0.315   noise: 0.045
Iter 34/50 - Loss: 0.035   lengthscale: 0.322   noise: 0.043
Iter 35/50 - Loss: 0.039   lengthscale: 0.326   noise: 0.041
Iter 36/50 - Loss: 0.043   lengthscale: 0.329   noise: 0.039
Iter 37/50 - Loss: 0.047   lengthscale: 0.327   noise: 0.038
Iter 38/50 - Loss: 0.052   lengthscale: 0.323   noise: 0.037
Iter 39/50 - Loss: 0.048   lengthscale: 0.317   noise: 0.036
Iter 40/50 - Loss: 0.051   lengthscale: 0.309   noise: 0.036
Iter 41/50 - Loss: 0.051   lengthscale: 0.302   noise: 0.036
Iter 42/50 - Loss: 0.047   lengthscale: 0.295   noise: 0.036
Iter 43/50 - Loss: 0.048   lengthscale: 0.288   noise: 0.036
Iter 44/50 - Loss: 0.047   lengthscale: 0.281   noise: 0.037
Iter 45/50 - Loss: 0.047   lengthscale: 0.276   noise: 0.037
Iter 46/50 - Loss: 0.040   lengthscale: 0.273   noise: 0.038
Iter 47/50 - Loss: 0.037   lengthscale: 0.271   noise: 0.039
Iter 48/50 - Loss: 0.040   lengthscale: 0.270   noise: 0.040
Iter 49/50 - Loss: 0.033   lengthscale: 0.269   noise: 0.042
Iter 50/50 - Loss: 0.032   lengthscale: 0.269   noise: 0.043

モデルで予測する

次のセルでは、モデルで予測を行ないます。これを行なうには、単純にモデルと尤度を eval モードにして、テストデータ上で両者のモジュールを呼び出します。

ちょうどユーザ定義の GP モデルが forward から事前平均と共分散を含む MultivariateNormal を返すように、eval モードの訓練された GP モデルは事後平均と共分散を含む MultivariateNormal を返します。こうして、予測平均と分散を得て、それから与えられたテストポイントにおける GP からのサンプリング関数は次のような呼び出しで達成されるでしょう :

f_preds = model(test_x)
y_preds = likelihood(model(test_x))
 
f_mean = f_preds.mean
f_var = f_preds.variance
f_covar = f_preds.covariance_matrix
f_samples = f_preds.sample(sample_shape=torch.Size(1000,))

gpytorch.settings.fast_pred_var コンテキストは必要ではありませんが、ここでは LOVE を使用してより高速な予測分布を得て、クールな特徴の一つを使用するプレビューを与えています :

# 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))

モデル fit をプロットする

次のセルでは、ガウス過程モデルの平均と信頼 (= confidence) 領域をプロットします。confidence_region メソッドは平均の上と下の 2 標準偏差を返すヘルパー・メソッドです。

with torch.no_grad():
    # Initialize plot
    f, ax = plt.subplots(1, 1, figsize=(4, 3))
 
    # Get upper and lower confidence bounds
    lower, upper = observed_pred.confidence_region()
    # Plot training data as black stars
    ax.plot(train_x.numpy(), train_y.numpy(), 'k*')
    # Plot predictive means as blue line
    ax.plot(test_x.numpy(), observed_pred.mean.numpy(), 'b')
    # Shade between the lower and upper confidence bounds
    ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)
    ax.set_ylim([-3, 3])
    ax.legend(['Observed Data', 'Mean', 'Confidence'])

以上

2020年7月
月	火	水	木	金	土	日
		1	2	3	4	5
6	7	8	9	10	11	12
13	14	15	16	17	18	19
20	21	22	23	24	25	26
27	28	29	30	31