PyTorch : Pyro examples : ガウス過程 (翻訳)

翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 10/28/2018 (v0.2.1)

* 本ページは、Pyro のドキュメント Examples : Gaussian Processes を翻訳した上で適宜、補足説明したものです：

• http://pyro.ai/examples/gp.html

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

 

ガウス過程

イントロダクション

ガウス過程 は教師あり、教師なし、そして強化学習問題においてさえも使用されてきました、そして洗練された数学的理論で記述されます (主題の概要については [1, 4] 参照)。それらはまた概念的に非常に魅力的です、何故ならばそれらは関数に渡る事前分布を定義する直感的な方法を与えるからです。そして最後に、ガウス過程はベイジアン設定で定式化されますので、それらは不確かさのパワフルな考えを装備するようになります。

幸い、Pyro は pyro.contrib.gp モジュールでガウス過程のための何某かのサポートを提供します。このチュートリアルの目標はこのモジュールのコンテキストでガウス過程 (GP) への簡単なイントロダクションを与えることです。殆どは Pyro の GP インターフェイスをどのように使用するかについて焦点を与えてそして一般的な GP についてのより詳細については読者にリファレンスを参照してもらいます。

私達が興味あるモデルは次で定義されます :
\[
f \sim \mathcal{GP}\left(0, \mathbf{K}_f(x, x’)\right)
\]

そして
\[
y = f(x) + \epsilon,\quad \epsilon \sim \mathcal{N}\left(0, \beta^{-1}\mathbf{I}\right)
\]

ここで \(x, x’ \in\mathbf{X}\) は入力空間の点で \(y\in\mathbf{Y}\) は出力空間の点です。\(f\) はカーネル \(\mathbf{K}_f\) で指定された GP 事前分布からのドローで \(\mathbf{X}\) から \(\mathbf{Y}\) への関数を表します。最後に、\(\epsilon\) はガウス観測ノイズを表します。

私達の GP のカーネルとして 動径基底関数カーネル (RBF カーネル) を使用します :
\[
k(x,x’) = \sigma^2 \exp\left(-\frac{\|x-x’\|^2}{2l^2}\right)
\]

ここで \(\sigma^2\) と \(l\) はカーネルを指定するパラメータです; 特に、\(\sigma^2\) は分散または二乗振幅 (= amplitude squared) そして \(l\) は長さスケールです。下でこれらのパラメータのための何某かの直感を得るでしょう。

 

インポート

最初に、必要なモジュールをインポートします。

from __future__ import print_function
import os
import matplotlib.pyplot as plt
import torch

import pyro
import pyro.contrib.gp as gp
import pyro.distributions as dist
from pyro.infer import SVI, Trace_ELBO
from pyro.optim import Adam

smoke_test = ('CI' in os.environ)  # ignore; used to check code integrity in the Pyro repo
pyro.enable_validation(True)       # can help with debugging
pyro.set_rng_seed(0)

チュートリアルを通して GP を可視化することを望みます。そのためプロットのためのヘルパー関数を定義します :

# note that this helper function does three different things:
# (i) plots the observed data;
# (ii) plots the predictions from the learned GP after conditioning on data;
# (iii) plots samples from the GP prior (with no conditioning on observed data)

def plot(plot_observed_data=False, plot_predictions=False, n_prior_samples=0,
         model=None, kernel=None, n_test=500):

    plt.figure(figsize=(12, 6))
    if plot_observed_data:
        plt.plot(X.numpy(), y.numpy(), 'kx')
    if plot_predictions:
        Xtest = torch.linspace(-0.5, 5.5, n_test)  # test inputs
        # compute predictive mean and variance
        with torch.no_grad():
            if type(model) == gp.models.VariationalSparseGP:
                mean, cov = model(Xtest, full_cov=True)
            else:
                mean, cov = model(Xtest, full_cov=True, noiseless=False)
        sd = cov.diag().sqrt()  # standard deviation at each input point x
        plt.plot(Xtest.numpy(), mean.numpy(), 'r', lw=2)  # plot the mean
        plt.fill_between(Xtest.numpy(),  # plot the two-sigma uncertainty about the mean
                         (mean - 2.0 * sd).numpy(),
                         (mean + 2.0 * sd).numpy(),
                         color='C0', alpha=0.3)
    if n_prior_samples > 0:  # plot samples from the GP prior
        Xtest = torch.linspace(-0.5, 5.5, n_test)  # test inputs
        noise = (model.noise if type(model) != gp.models.VariationalSparseGP
                 else model.likelihood.variance)
        cov = kernel.forward(Xtest) + noise.expand(n_test).diag()
        samples = dist.MultivariateNormal(torch.zeros(n_test), covariance_matrix=cov)\
                      .sample(sample_shape=(n_prior_samples,))
        plt.plot(Xtest.numpy(), samples.numpy().T, lw=2, alpha=0.4)

    plt.xlim(-0.5, 5.5)

 

データ

データは次からサンプリングされた $20$ 点から成ります :
\[
y = 0.5\sin(3x) + \epsilon, \quad \epsilon \sim \mathcal{N}(0, 0.2)
\]

$x$ は区間 \([0, 5]\) から一様にサンプリングされています。

N = 20
X = dist.Uniform(0.0, 5.0).sample(sample_shape=(N,))
y = 0.5 * torch.sin(3*X) + dist.Normal(0.0, 0.2).sample(sample_shape=(N,))

plot(plot_observed_data=True)  # let's plot the observed data

 

モデルを定義する

最初に 2 つのハイパーパラメータ variance と lengthscale の値を指定して、RBF カーネルを定義します。それから GPRegression オブジェクトをコンストラクトします。ここでもう一つのハイパーパラメータ noise を供給します、これは上の \(\epsilon\) に相当します :

kernel = gp.kernels.RBF(input_dim=1, variance=torch.tensor(5.),
                        lengthscale=torch.tensor(10.))
gpr = gp.models.GPRegression(X, y, kernel, noise=torch.tensor(1.))

この GP 関数事前分布からのサンプルがどのようなものか見てみましょう。これはデータ上で条件付ける前であることに注意してください。これらの関数が取る形状 — それらの滑らかさ、垂直スケール, etc. — は GP カーネルにより制御されます。

plot(model=gpr, kernel=kernel, n_prior_samples=2)

例えば、分散とノイズをより小さくすればより小さい垂直振幅を持つ関数サンプルを見ます :

kernel2 = gp.kernels.RBF(input_dim=1, variance=torch.tensor(0.1),
                         lengthscale=torch.tensor(10.))
gpr2 = gp.models.GPRegression(X, y, kernel2, noise=torch.tensor(0.1))
plot(model=gpr2, kernel=kernel2, n_prior_samples=2)

 

推論

上でカーネル・ハイパーパラメータを手動で設定しました。ハイパーパラメータをデータから学習することを望む場合、推論を行なう必要があります。最も単純な (共役な) ケースでは対数周辺尤度上で勾配上昇を行ないます。pyro.contrib.gp では単に Pyro SVI インターフェイスを使用します (この場合どのような潜在確率変数もサンプリングさえしませんが)。

optim = Adam({"lr": 0.005})
svi = SVI(gpr.model, gpr.guide, optim, loss=Trace_ELBO())
losses = []
num_steps = 2500 if not smoke_test else 2
for i in range(num_steps):
    losses.append(svi.step())

# let's plot the loss curve after 2500 steps of training
plt.plot(losses);

何か合理的なものを学習したかを見てみましょう。

plot(model=gpr, plot_observed_data=True, plot_predictions=True)

ここで太い赤いカーブは平均予測で青いバンドは平均のまわりの 2-sigma 不確かさを表しています。合理的なカーネル・ハイパーパラメータを学習したようです、何故ならば平均と不確かさの両者がデータへの合理的な適合を与えるからです。(e.g. 巨大過ぎる学習率を選択したり悪い初期ハイパーパラメータを選択すれば学習は簡単に悪い方向に進んだ可能性があることに注意してください。)

分散と長さスケールが正であるときに限りカーネルは well-defined であることに注意してください。ハイパーパラメータが適切な領域に制約されることを確実にするために Pyro は裏では PyTorch constraints (ドキュメント 参照) を使用してます。ハイパーパラメータの制約された値を見るために get_param(…) メソッドを使用します。

gpr.kernel.get_param("variance").item()

0.20291300117969513

gpr.kernel.get_param("lengthscale").item()

0.5022150278091431

gpr.get_param("noise").item()

0.04273035749793053

データを生成する正弦曲線の周期は \(T = 2\pi/3 \approx 2.09\) ですので、1/4 周期におおよそ等しい長さスケールの学習は意味があります。

 

MAP を使用してモデルを fit する

ハイパーパラメータのための事前分布を定義する必要があります。

# note that our priors have support on the positive reals
gpr.kernel.set_prior("lengthscale", dist.LogNormal(0.0, 1.0))
gpr.kernel.set_prior("variance", dist.LogNormal(0.0, 1.0))
# we reset the param store so that the previous inference doesn't interfere with this one
pyro.clear_param_store()
optim = Adam({"lr": 0.005})
svi = SVI(gpr.model, gpr.guide, optim, loss=Trace_ELBO())
losses = []
num_steps = 2500 if not smoke_test else 2
for i in range(num_steps):
    losses.append(svi.step())
plt.plot(losses);

plot(model=gpr, plot_observed_data=True, plot_predictions=True)

学習したハイパーパラメータを調査しましょう :

for param_name in pyro.get_param_store().get_all_param_names():
    print('{} = {}'.format(param_name, pyro.param(param_name).item()))

RBF$$$variance_MAP = 0.24471372365951538
GPR$$$noise = 0.04222385957837105
RBF$$$lengthscale_MAP = 0.5217667818069458

MAP 値は事前分布のために MLE 値とは異なることに注意してください。パラメータ名が変わっていることも見て取れます。これはカーネルのハイパーパラメータ上に事前分布を置いたからで、それらを確率変数に変換しました (それらはもはや torch.nn.Parameters ではありません)。内部では、GP モジュールはこの違いを反映するために末尾の (= trailing) _MAP を各ハイパーパラメータ名に追加しています。

 

スパース GP

巨大なデータセットについては対数周辺尤度の計算は伴う高価な行列演算によりコスト高です ([1] のセクション 2.2 参照)。いわゆる「スパース」変分法は GP をより巨大なデータセットのために実行可能とするために開発されてきました。これは研究の大きな領域で詳細の総てには進みません。代わりにこれらのメソッドを活用するために pyro.contrib.gp で SparseGPRegression をどのように利用できるかを素早く示します。

最初に、より多くのデータを生成します。

N = 1000
X = dist.Uniform(0.0, 5.0).sample(sample_shape=(N,))
y = 0.5 * torch.sin(3*X) + dist.Normal(0.0, 0.2).sample(sample_shape=(N,))
plot(plot_observed_data=True)

スパース GP の使用は上で使用した基本的な GP の使用に非常に類似しています。単に特別なパラメータ \(X_u\) (誘導点) を追加する必要があるだけです。

# initialize the inducing inputs
Xu = torch.arange(20.) / 4.0

# initialize the kernel and model
kernel = gp.kernels.RBF(input_dim=1)
# we increase the jitter for better numerical stability
sgpr = gp.models.SparseGPRegression(X, y, kernel, Xu=Xu, jitter=1.0e-5)

# the way we setup inference is similar to above
pyro.clear_param_store()
optim = Adam({"lr": 0.005})
svi = SVI(sgpr.model, sgpr.guide, optim, loss=Trace_ELBO())
losses = []
num_steps = 2500 if not smoke_test else 2
for i in range(num_steps):
    losses.append(svi.step())
plt.plot(losses);

# let's look at the inducing points we've learned
print("inducing points:\n{}".format(pyro.param("SGPR$$$Xu").data.numpy()))
# and plot the predictions from the sparse GP
plot(model=sgpr, plot_observed_data=True, plot_predictions=True)

inducing points:
[0.0295011  0.24607515 0.44147182 0.7126421  1.0303956  1.3245093
 1.5803136  1.877679   2.1162455  2.4676545  2.4189417  2.8084621
 3.070498   3.3058653  3.577411   3.8880596  4.1536465  4.4360414
 4.709121   4.9547815 ]

モデルがデータへの合理的な適合を学習することを見て取れます。現在 Pyro で実装されている 3 つの異なるスパース近似があります :

• “DTC” (Deterministic Training Conditional)
• “FITC” (Fully Independent Training Conditional)
• “VFE” (変分自由エネルギー, Variational Free Energy)

デフォルトでは、SparseGPRegression は推論法として “VFE” を使用します。SparseGPRegression への異なる approx フラグを渡すことにより他のメソッドを使用できます。

 

より多くのスパース GP

上の GPRegression と SparseGPRegression はガウス尤度に制限されます。GP で他の尤度を使用することもできます — 例えば、分類問題のために Bernoulli を使用することもできます — しかし推論問題はより難しくなります。このセクションでは、VariationalSparseGP モジュールをどのように使用するかを示します、これは非ガウス尤度を扱うことができます。そこで上で行なったことと比較することができ、私達は依然としてガウス尤度を使用していきます。ポイントは内部で成された推論は他の尤度もサポートできるということです。

# initialize the inducing inputs
Xu = torch.arange(10.) / 2.0

# initialize the kernel, likelihood, and model
kernel = gp.kernels.RBF(input_dim=1)
likelihood = gp.likelihoods.Gaussian()
# turn on "whiten" flag for more stable optimization
vsgp = gp.models.VariationalSparseGP(X, y, kernel, Xu=Xu, likelihood=likelihood, whiten=True)

pyro.clear_param_store()
# instead of defining our own training loop, we will
# use the built-in support provided by the GP module
num_steps = 1500 if not smoke_test else 2
losses = vsgp.optimize(optimizer=Adam({"lr": 0.01}), num_steps=num_steps)
plt.plot(losses);

plot(model=vsgp, plot_observed_data=True, plot_predictions=True)

これだけのことです。pyro.contrib.gp モジュールのより詳細については ドキュメント を見てください。そして分類のための GP を使用するサンプルコードについては こちら を見てください。

 

Reference

• [1] Deep Gaussian processes and variational propagation of uncertainty, Andreas Damianou
• [2] A unifying framework for sparse Gaussian process approximation using power expectation propagation, Thang D. Bui, Josiah Yan, and Richard E. Turner
• [3] Scalable variational Gaussian process classification, James Hensman, Alexander G. de G. Matthews, and Zoubin Ghahramani
• [4] Gaussian Processes for Machine Learning, Carl E. Rasmussen, and Christopher K. I. Williams
• [5] A Unifying View of Sparse Approximate Gaussian Process Regression, Joaquin Quinonero-Candela, and Carl E. Rasmussen

以上