はじめに
今回はガウス過程の説明を最低限に行いつつ、GPyTorchを使ってガウス過程をさらっと見てみようという記事です。 関連記事に以下のようなものがあるので参考にしてみてください。
ガウス過程（GP）
GPは関数を生成する確率分布（確率過程）であり、平均関数 $m$ と共分散関数 $K$ をパラメータとして関数 $f$ を生成する。
$$ f \sim {\mathcal GP} (m, K) $$
GP回帰では上記で生成された関数を用いて、入力データ $x$ と出力データ $y$ を下記でモデル化する。
$$ \begin{align} y &= f(x) + \epsilon \\ \epsilon &\sim {\mathcal N}(0, \sigma) \end{align} $$
ここに、$\epsilon$ は平均 $0$ で標準偏差 $\sigma$ に従うノイズ項である。 一般の機械学習でも、モデルに対して何らかのノイズ項を設けて出力にバラつきがあることを表現する。 加えてGP回帰では、モデルとして使われる関数を通常の線形回帰で使われる $f(x) = wx + b$ などのように仮定せずに、関数自体も確率的に揺らぐ何かであるという緩い仮定しか設けない。パラメータを用いて特定の形状の関数を作り、パラメータを調整することでデータにフィッティングさせる方法を「パラメトリックな手法」と呼び、GPのような手法は「ノンパラメトリック」な手法と呼ばれる。
しかし、実際には全くパラメータを置かずにモデルを表現するということはできない。学習がパラメータの調整によって実施されていないというだけで、モデルには手で設計するパラメータ（ハイパーパラメータ）が存在するし、一段上の立場に立って手で設計したハイパーパラメータ自体もデータから決めようという試みもある（経験ベイズ法や第二種最尤推定、変分ベイズ法などを用いることができる）。下記でGPの中身を概観する。
GPyTorchを使ったモデリング
さっそくコードを見ながらガウス過程を見ていく。
class GPModel(gpytorch.models.ExactGP): def __init__(self, train_x, train_y, likelihood, lengthscale_prior=None, outputscale_prior=None): super(GPModel, self).__init__(train_x, train_y, likelihood, ) self.mean_module = gpytorch.means.ConstantMean() self.covar_module = gpytorch.kernels.ScaleKernel( gpytorch.kernels.RBFKernel(lengthscale_prior=lengthscale_prior), outputscale_prior=outputscale_prior ) def forward(self, x): mean_x = self.mean_module(x) covar_x = self.covar_module(x) return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)
上記のコードでガウス過程のモデルクラスを簡単に作成することができる。
コード概要
GPyTorchでは決まりきったGPでの計算処理を隠蔽して、学習で重要な項目となるモデルの構造とハイパーパラメータ、確率分布を柔軟に設計できるようになっている。
__init__メソッドの引数であるtrain_xとtrain_yは学習に用いる入力と出力のペアである。後述するが、GP（含めカーネル法を用いる手法）ではモデル自体が学習データを保持して、学習データそれ自体を使って新しいデータに対する予測を記述する（通常の線形回帰などでは学習データで学習パラメータを調整し、その後は学習データは不要で、調整したパラメータを使って予測を行う）。
また、likelihoodはGPから生成される関数 $f$ によって定まる $x,y$ の関係性に関しての確率分布を定める。回帰ではガウス分布を仮定した場合の尤度を用いる。lengthscale_priorとoutputscale_prior はGPから生成される関数 $f$ 自体がどのように生成されるかを決めるハイパーパラメータの振る舞いを定める（ハイパーパラメータの事前分布）。
superはGPyTorchで実装されている ExactGP クラスを継承し、学習データと尤度関数を引数として初期化している。この中でデータと尤度を使った計算は隠蔽されているため、残りのモデリングに集中できるという仕組みになっている。
その後の
self.mean_module=gpytorch.means.ConstantMean()
self.covar_module = gpytorch.kernels.ScaleKernel(
gpytorch.kernels.RBFKernel(lengthscale_prior=lengthscale_prior),
outputscale_prior=outputscale_prior
)
によって $f$ が生成されるGPの平均関数と共分散関数を定めている。具体的には
$$ f \sim {\mathcal GP} (m, K) $$
に関して $m = \rm Const$ としているのが一行目である。二行目が重要であり、共分散関数を定めているこの部分にGPの重要な要素が含まれている。共分散関数は、2つのデータを引数に取るカーネル関数 $k(x_i, x_j)$ と呼ばれるものによって規定される。ここではRadial Basis Function（RBF）カーネルを利用している。RBFカーネルは
$$ k(x _ i, x _ j) = \sigma ^ 2 \exp \left\{ -\frac{(x _ i - x _ j) ^ 2}{2l ^ 2} \right\} $$
と書かれる。ここで $l$ は lengthscale と呼ばれるハイパーパラメータである。GPyTorchでは、ハイパーパラメータにも事前分布を与えることができる（もちろん与えないこともできる）ので、コードでは事前分布lengthscale_priorを与えるようなコードの書き方をしている。同様に $\sigma$ は outputscale と呼ばれるハイパーパラメータである。
カーネル関数は各データのインデックスを2つ $i, j$ と指定し、それに対応するデータを受け取った時に、上記で計算されるスカラー値を返す関数である。直感的にはカーネル関数で計算されるスカラー値は、引数として渡した2つのデータの類似度のようなものに対応する。共分散関数は実体としては行列の形式をしており、$i$ 行 $j$ 列に、$k(x_i,x_j)$ の値が格納されている。すなわち、全てのデータから2つのデータを取り出した時に、カーネル関数によって計算される値を、インデックスをちゃんとつけて保存しているということある。仮に $x_i$ と $x_j$ のカーネル関数による値が大きい (=類似度が高い）という場合には、共分散関数で対応する $(i, j)$ の要素が大きな値を持ち、結果として $f(x_i), f(x_j)$ は相関が大きいというようなことを表現できる。
今、outputscaleをRBFカーネルのハイパーパラメータのごとく説明したが、実際のところはRBFカーネルとして必須なのは lengthscale と呼ばれる $l$ の方だけかもしれない（例えば、SVMのようにカーネル法を利用する他の手法での説明を見てほしい。SVMでは特徴空間での値の配置関係だけが問題になるからだろうか）。なぜここで、追加で $\sigma$ なるハイパーパラメータを考えたのかを説明する。仮に lengthscale しかない場合には、カーネルの値は $1$ より大きな値を取らない。すなわち、分散共分散の各成分が$1$ より大きな値を取らないということである。分散が本当にそうであるかというと、そんなわけもないので、大きさの具合を調整するパラメータを追加したのだ。しかも都合の良いことに
$$ \begin{align} y &= f(x) + \epsilon \\ \epsilon &\sim {\mathcal N}(0, \sigma) \end{align} $$
において、式で現れる $\epsilon$ をカーネルに包含してしまえたりもする。というのも $\epsilon$ は要するにただ分散に影響を与えるだけであるので、そうであるならば、共分散関数の対角成分にだけ、ノイズ分の分散を加えておけば良い。そういうわけで、RBFカーネルが gpytorch.kernels.ScaleKernel なるクラスにわたされているのだ。
学習コード
GPでは手元のデータと定められたカーネル関数から、新たなデータに対する予測を計算できる。また、新たな入力データとその出力のセットが得られたならば、単にGPのモデルにデータを追加してやれば予測モデル自体をアップデートできる（ループ計算無しで！！）。なぜなら、手元の各入力データの類似度と各データに対応する出力との関係性から、新たなデータに対する予測が解析的に定まるからである。
するとGPにとって重要なのは、そもそもそのカーネル関数の設定は合っていますか？ということになる。先程までハイパーパラメータと呼んでいたものを、何らかの方法でデータから決めてやりたいということである。これをモデルの対数周辺尤度を最大化するようなハイパーパラメータの選定（第二種最尤推定）によって決める。通常、対数周辺尤度がどんな形式であるのか、計算を実行するのは積分を含むために困難である。ところがGPの場合は内部の分布の性質が良いために数式として得ることが可能である。また、GPyTorchは作ったモデルと尤度関数を与えてやれば、対数周辺尤度を作ってくれるモジュールgpytorch.mlls.ExactMarginalLogLikelihood(likelihood, gpr) が準備されている。
勘の良い人ならば、対数周辺尤度と言っても、一体なにを周辺化したのか？と疑問に思うかもしれない。確率モデル $p(y \mid x)$ があるときに、 カーネル関数のハイパーパラメータ $\theta$ を明示すれば $p(y \mid x, \theta)$ である。ところでガウス過程がノンパラメトリックであるというのを差し置いて、本来、学習パラメータ $w$ があり、 $p(y \mid x, w)$ であったとする。そしてハイパーパラメータも明示するならば $p(y \mid x, w, \theta)$ である。この時、$w$ についての周辺化 $$ p(y \mid x,\theta) =\int _ w p(y \mid x, w, \theta) p(w) {\rm d}w$$ というものを考えることができる。この $w$ を綺麗に解析的に消去してしまったことでパラメータが数式上から消えてなくなっているのである。GPはそうした後にカーネル法云々によって定式化されていく。なので、 $\log p(y \mid x, \theta)$ はまぎれもなく対数周辺尤度なのである。
そうなれば、あとは勾配法の手続きによって対数周辺尤度のハイパーパラメータによる最適化が実施可能で、PyTorchのコードで簡単に実現できる。以下にTrainer クラスを実装しておく。
class Trainer(object): def __init__(self, gpr, likelihood, optimizer, mll): self.gpr = gpr self.likelihood = likelihood self.optimizer = optimizer self.mll = mll def update_hyperparameter(self, epochs): self.gpr.train() self.likelihood.train() for epoch in range(epochs): self.optimizer.zero_grad() output = self.gpr(self.gpr.train_inputs[0]) loss = - self.mll(output, self.gpr.train_targets) loss.backward() self.optimizer.step() if epoch % (epochs//10) == 0: print('Epoch %d/%d - Loss: %.3f ' % ( epoch + 1, epochs, loss.item(), ))
データとモデルの準備
GPの学習に必要になる、モデル、ハイパーパラメータ（の事前分布）、尤度関数、対数周辺尤度を準備する。
またモデルの構築時点で学習データが必要なため、適当なToyDataを準備しておく。あとはPyTorchでの学習のためにoptimizerも準備しておく。
lengthscale_priorにはl_priorとして平均 $1$ で標準偏差 $10$ の正規分布を用いる。 outputscaleに対しても同様である。標準偏差を大きくとっておけば、事前分布が学習に大きく寄与することはなくなってくる。特にデータに対する知識を持ち合わせていない場合は、ハイパーパラメータに偏見を持たせないような無情報事前分布を用いる場合が多い。
## GPでは入力は多次元前提なので (num_data, dim) という shape ## 一方で出力は一次元前提なので (num_data) という形式にする train_inputs = torch.linspace(0, 1, 10).reshape(10, 1) train_targets = torch.sin(2*np.pi*train_inputs).reshape(10) + 0.3*torch.randn(10) likelihood = gpytorch.likelihoods.GaussianLikelihood() l_prior = gpytorch.priors.NormalPrior(loc=torch.tensor(1.), scale=torch.tensor(10.)) s_prior = gpytorch.priors.NormalPrior(loc=torch.tensor(1.), scale=torch.tensor(10.)) gpr = GPModel(train_inputs, train_targets, likelihood, lengthscale_prior=l_prior, outputscale_prior=s_prior) mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, gpr) optimizer = torch.optim.RMSprop(params=gpr.parameters(), lr=1e-2)
学習と結果
下記のコードで学習を実行できる。
trainer = Trainer(gpr, likelihood, optimizer, mll)
trainer.update_hyperparameter(2000)
学習実施後は下記のようにテストデータを流してみる。テストデータは1点でも良いのだが、どうせならどんな曲線を得たのか、学習データより広い領域のあらゆる点を見てみよう。
test_inputs = torch.linspace(-0.2, 1.2, 100) gpr.eval() likelihood.eval() with torch.no_grad(): predicts = likelihood(gpr(test_inputs)) predicts_mean = predicts.mean predicts_std = predicts.stddev
ここで、gpr の __call__ の結果を likelihood に渡してやることを忘れないようにしたい。少し独特だが、最終的な確率モデルとしての振る舞いを司っているのが likelihood なのだ。こいつの戻り値にはいろいろなプロパティが含まれているが mean と stddev を今回は見てみる。
plt.plot(test_inputs.numpy(), predicts_mean.numpy())
plt.fill_between(test_inputs.numpy(),
predicts_mean.numpy() - 0.9*predicts_std.numpy(),
predicts_mean.numpy() + 0.9*predicts_std.numpy(), alpha=0.4)
plt.plot(train_inputs.numpy(), train_targets.numpy(), "ro")
上記では学習データが赤い丸、青い線が回帰曲線（予測分布の平均曲線）、薄い青の領域が予測された各点における標準偏差の大きさの $0.9$ 掛けを上下塗りつぶした。
重要なのは仮にパラメトリックに上記のデータを回帰したい場合は、三次関数なり三角関数なりを仮定しなければならない。しかしGPではデータを得た時に適応的に回帰モデルを柔軟に変えていくことができる（モデルは共分散関数によって規定され、共分散関数の正体はデータの類似度を埋め込んだものであった。すなわちデータが追加されたらGPが共分散関数を通して変化する）。
また、確率モデルであるため予測が分布で得られ、したがってバラつきを評価することができる。これは言い換えれば、ある点でバラつきが大きいというのは予測の自信のなさを表している。分からないことは分からないと正直に言うのである。だから、データが得られていない領域（0より小さく、1より大きい領域）は薄く塗りつぶされた領域が広がっている。確率モデルでない場合は、単に回帰曲線だけが得られ、それがどの程度モデルにとって自信のある結果なのかはわからない。
ハイパーパラメータの振る舞い
lengthscale
ここで強い事前分布を仮定することで、ハイパーパラメータの値を強く制限してみる。事前分布に引きづられた学習結果になることの確認と、ハイパーパラメータがどういう意味を持ち、どんな振る舞いをするのかを確認する。まずはlengthscaleが非常に小さくなるような事前分布を仮定する。こういう場合には正規分布で作る事前分布の平均で値を指定し、標準偏差を異様に小さくすれば良い。
## GPでは入力は多次元前提なので (num_data, dim) という shape ## 一方で出力は一次元前提なので (num_data) という形式にする train_inputs = torch.linspace(0, 1, 10).reshape(10, 1) train_targets = torch.sin(2*np.pi*train_inputs).reshape(10) + 0.1*torch.randn(10) likelihood = gpytorch.likelihoods.GaussianLikelihood() l_prior = gpytorch.priors.NormalPrior(loc=torch.tensor(0.01), scale=torch.tensor(0.01)) s_prior = gpytorch.priors.NormalPrior(loc=torch.tensor(1.), scale=torch.tensor(10.)) gpr = GPModel(train_inputs, train_targets, likelihood, lengthscale_prior=l_prior, outputscale_prior=s_prior) mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, gpr) optimizer = torch.optim.RMSprop(params=gpr.parameters(), lr=1e-2) trainer = Trainer(gpr, likelihood, optimizer, mll) trainer.update_hyperparameter(2000)
こうして学習したものを、先ほどと同様に図で見てみる。
lengthscaleが小さいと、GPの平均関数は $0$ の定数にされているため、基本は $0$ なのだが、データがあるところだけ値が飛び出ている。これは、周囲のデータとの類似度をどれだけ尊重するかを調整しているパラメータであると言える。今回、周りのデータを参考にする度合いが極端に小さいため、各点各点の丸暗記に走ったのだ。逆にlengthscaleが大きければ近くのデータを参考に予測を実施する。あまりに大きければ、全く関係ない場所のデータにも引きづられてしまうこととなる。
outputscale
次はoutputscale を見てみる。
## GPでは入力は多次元前提なので (num_data, dim) という shape ## 一方で出力は一次元前提なので (num_data) という形式にする train_inputs = torch.linspace(0, 1, 10).reshape(10, 1) train_targets = torch.sin(2*np.pi*train_inputs).reshape(10) + 0.1*torch.randn(10) likelihood = gpytorch.likelihoods.GaussianLikelihood() l_prior = gpytorch.priors.NormalPrior(loc=torch.tensor(1.), scale=torch.tensor(10.)) s_prior = gpytorch.priors.NormalPrior(loc=torch.tensor(0.01), scale=torch.tensor(0.01)) gpr = GPModel(train_inputs, train_targets, likelihood, lengthscale_prior=l_prior, outputscale_prior=s_prior) mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, gpr) optimizer = torch.optim.RMSprop(params=gpr.parameters(), lr=1e-2) trainer = Trainer(gpr, likelihood, optimizer, mll) trainer.update_hyperparameter(2000)
outputscaleを小さくすると回帰曲線がほとんど平坦になっている。名前から想像できる通り、曲線が変動できる度合いを指定しているもので、これが小さいために、学習の中で値を大きく取ることができなかったのだ。逆にこの値を大きくしすぎれば、仮にどこかに異常に大きなデータが混入してしまった場合に、無理やりその値を追いかけるようになる。
最後に
今回はGPyTorchを使ってGPの基本を見た。
実は今回、勾配法によって対数周辺尤度を $\log p(y\mid x, \theta)$ を最大化することでハイパーパラメータを獲得したのだが、通常、この関数は複数の峰を持つとされている。すなわち学習によって得られるハイパーパラメータの推定値は、初期値に依存するということである。ハイパーパラメータはほぼ決め打ちのつもりで、強い事前分布を仮定したならば、そもそも最適解を決め込んでるも同然であるから関係はなさそうだが、事前分布を緩めた場合には、本来の対数周辺尤度の形の影響が強く出る。
もっとモデルとデータを信用し、慎重にハイパーパラメータについて議論する必要がある場合は勾配法ではなくMCMCなどの利用を検討しなければならない。
また、共分散関数はデータに応じて膨らんでいき、内部の計算ではデータの3乗のオーダーで計算量が増えていく。これは近年のビッグデータを扱う場合には無視できない欠点だ。従って、あえて学習データを制限する（例えば関数があまり複雑でない領域はデータの密度を少なくしたりするなど）工夫が必要である。学習データを少なくするのはもったいないので、全てのデータを使って、性質を知っていると嬉しいような領域を推定する試みもできる。そういう場合には変分ベイズ法などより高度な手法を混ぜ込む必要が出てくる。
またGP自体は多入力1出力のシステムである。近年は多入力多出力が当たり前になってきているので、これに対応するためにマルチタスクGPなどの利用が必要かもしれない（これは、GPがそれぞれの次元を個別に対応する、すなわち出力の数だけ別のGPを用意する、という素朴な利用方法を発展させ、各GPがうまい具合にカーネルを通じて情報を共有する仕組みである）。
