Juliaのガウス過程ライブラリ GaussianProcesses.jl を触ってみる。 チュートリアルはまだ未整備のようなので、サンプルノートブックとガウス過程と機械学習」などを参考にしながらやっていく。 Juliaのはv1.1.0, GaussianProcesses.jlはv0.9.0を使っている。
インストールはいつものように
using Pkg
Pkg.add("GaussianProcesses")
カーネル関数
すでに多くの種類のカーネル関数が実装されており、 ガウスカーネル、線形カーネル、指数カーネル、周期カーネル、Matérnカーネルなど、 「ガウス過程と機械学習」に出てきたカーネル関数はほぼカバーされているようだ。
全て見ていくとキリがないのでいくつか基本的な部分だけ確認することにする。
ガウスカーネルのカーネル関数は「ガウス過程と機械学習」では $$ k(\mathbf{x}, \mathbf{x}^\prime) = \theta_1 \exp \left( - \frac{|\mathbf{x} - \mathbf{x}^\prime|^2}{\theta_2}\right) $$
となっているが、GaussianProcesses.jlでこれに対応するのは [GaussianProcesses.SEIso]
(http://stor-i.github.io/GaussianProcesses.jl/latest/kernels.html#GaussianProcesses.SEIso) である。
関数のパラメトライズは対数スケールとなっているため十分注意する必要があり、カーネル関数 SEIso(ll::T, lσ::T) は、
\(\sigma = \log (\) lσ \()\), \(l = \log (\) ll \()\) とした時に
$$ k(\mathbf{x}, \mathbf{x}^\prime) = \sigma^2 \exp \left( - \frac{(\mathbf{x} - \mathbf{x}^\prime)^\top(\mathbf{x} - \mathbf{x}^\prime))}{2l^2}\right) $$
だから、\( k(\mathbf{x}, \mathbf{x}^\prime) = \exp \left( - |\mathbf{x} - \mathbf{x}^\prime|^2\right) \) は
rbf = SEIso(log(1/sqrt(2)), 0.0)
に対応している。「ガウス過程と機械学習」の pp.69 のように、\( (x_1, x_2, x_3, x_4) = (1, 2, 3, 4) \) とした時の共分散行列を求める。
xs = [1, 2, 3, 4]
cov(rbf, xs, xs)
このようにすると xs の全体が1要素として認識されてしまい
1.0
が帰ってきてダメなので、reshape で 1×4 Arrayに変換して次のようにする。
xs = reshape([1.0, 2, 3, 4], 1, :)
cov(rbf, xs, xs)
4×4 Array{Float64,2}:
1.0 0.367879 0.0183156 0.00012341
0.367879 1.0 0.367879 0.0183156
0.0183156 0.367879 1.0 0.367879
0.00012341 0.0183156 0.367879 1.0
となり本と同様の結果が得られた。ちなみに
xs = reshape([1, 2, 3, 4], 1, :)
cov(rbf, xs, xs)
は
InexactError: Int64(0.3678794411714422)
となり動かないので、xs を Real のベクターにする。
ガウス過程からのサンプル
次はこのカーネル関数で定義したガウス過程からサンプリングを行う。GaussianProcesses.jl でガウス過程を呼び出すメソッド GPE は、
学習データからフィッティングを行うものと行わないものの両方が用意されている。
すでにカーネル関数を定義したので、使うのはフィッティングしない方。
ちなみに GPE は観測ノイズを考慮しないノイズフリーのガウス過程のフィッティングをするメソッドで、観測ノイズも考慮したフィッティングを行いたければ GP を使う。
この辺りは 公式ドキュメント を見た方が良いかもしれない。
先ほどのRBFカーネルで定まるガウス過程を定義し、サンプルを5つ取得する。
コード上ではカーネルのみを指定して平均は与えていないが、GPE の mean 引数のデフォルトが MeanZero() なので平均0のガウス過程を考えていることになる。
mean 引数を変更すれば様々な平均関数のガウス過程を考えられる(はず、未確認)。
サンプルは rand で行うことができ、
xs = -25:0.1:25
gp_rbf = GPE(kernel = rbf)
samples = rand(gp_rbf, xs, 5)
plot(xs, samples, label = "", lw = 2)
非常にお手軽である。
カーネルの和、積
カーネルの和と積も定義されている。線形カーネルと周期カーネルで和、積を使ってみる。
kernel_sum = Lin(2.0) + Periodic(0.0, 0.0, 2.0)
gp_sum = GPE(kernel = kernel_sum)
samples = rand(gp_sum, xs, 5)
plot(xs, samples, label = "", lw = 2)
kernel_prod = Lin(2.0) * Periodic(0.0, 0.0, 2.0)
gp_prod = GPE(kernel = kernel_prod)
samples = rand(gp_prod, xs, 5)
plot(xs, samples, label = "", lw = 2)
これまた簡単。
「やってみた」系の記事になってしまったが、今日はここまで。
まだ少し触っただけではあるが、手軽にガウス過程が使えて便利。
次は何か具体例を用いて GP, GPE のフィッティングを行いたい。
上の内容のjupyter notebook:
https://github.com/matsueushi/notebook_blog/blob/master/gaussianprocess-julia.ipynb
参考
FAIRBROTHER, Jamie, et al. GaussianProcesses. jl: A Nonparametric Bayes package for the Julia Language. arXiv preprint arXiv:1812.09064, 2018.
GitHubレポジトリ: STOR-i/GaussianProcesses.jl, サンプルノートブック
ドキュメント: http://stor-i.github.io/GaussianProcesses.jl/latest/
持橋大地, 大羽成征. ガウス過程と機械学習, 講談社 機械学習プロフェッショナルシリーズ, 2019