カーネル法とガウス過程（PRML下巻第6章）

「ガウス過程」- PRML §6.4

３年前、本レーンでこの辺りを読んだ頃には線形代数力が弱すぎて、というか行列式と式展開から何もイメージがつかめなくて、皆さんの~~空中戦~~発表を聞きながら式を書き写してるだけ、みたいな感じだったのを覚えている。あたかも理解できてるような事を喋れていれば「中国語の部屋」的アナロジーが使えるのだけれどそれにも到底及ばないような。

PRML読書会も３周目にもなると、突如PRMLが日本語で書いてあるように読めるようになる瞬間がある（実話）。※１周目、遅くとも２周目でそれが来ていれば復々習レーンは無かった。

とかちょっと前置きをしてみたのだけれど、本質的にこれは「グラフ描きたいだけのエントリ」である。

PRML読んでて一番楽しいのは図表の再現プログラムを書くこと
— naoya tさん (@naoya_t) 11月 12, 2012

だから仕方がない。

本題に入ろう。

そもそも「××過程」とは何なんだろう

ガウス過程、ディリクレ過程、のような「××過程」というのがそもそも意味不明だった。「過程（process）」って何？みたいなところから。

より一般的には、確率過程 (stochastic process) y(x) とは、任意の有限な値集合 y(x_1), ..., y(x_n) に対して、矛盾のない同時分布を与えるものである。

日本語は分かるのだけれど脳内で像を結ばない。
図6.4とか見て、ぼんやりMCMC*1みたいなのを想像していた。書かれている式の意味するところも理解できずに、そんな想像していたまま３年が過ぎたある秋。

復々習レーンのMCMC（二代目司会進行の意）であるところの@sleepy_yoshi大先生が、そろそろHackathonやろう的な事を言い出した。ニューラルネットワークの回（第5章）に会場でNNの図表再現コードを書いてて楽しかったので、図表再現Hackathonでもしたいなと思い３年の間にずいぶん風化した我がPRMLをパラパラと捲ってみる。自分的謎度の高かった６章の図表の数々。これが再現できたら少し自信がつけられるかもしれない。

shuyoさんのコードの写経から

足掛かりになるのは shuyo さんの3年前の記事。
PRML６章「ガウス過程による回帰」を R で試す - Mi manca qualche giovedì?

とりあえずRをOctaveに直訳してみる。

3引数の、というか第3引数に関数を取る outer() 関数に該当するものがないので別途準備
diag(10) の挙動がRとOctaveで違う。Octaveでは diag(c(10)) = [10] と解釈されてしまうので eye(10) で代用。

f:id:n4_t:20121113125626p:plain
出来た。分かった気になってきた。というか数式とコードの対応が今ならはっきりと見える。

x_1〜x_N は与えられたデータ点、 $x_{N+1}$ は新しいデータ点なのだけれど、新しいデータ点 $x_{N+1}$ をx上のぽつんとした固定の1点で考えていて、なんでそれでグラフが描けるんだろう意味不明だなあと思っていた。でも新しいデータ点として任意の値をとりうる、というかとっていいのだから linspace(-1,1,100) とか取ったら t_N+1 を曲線としてplotできるという事だよなー、と気づいて目の前が少し開けた。ビジョンが１次元→２次元に広がった感じ。分からないことその１が解決。

標準偏差影付きグラフ（図6.8）を描きたい

さて。図6.8（下巻p.21）みたいに標準偏差エリアの影付きのグラフが描きたい。
$m(x_N+1)$ に加え、 $\sigma^2(x_N+1)$ も求めたい。これは式(6.67) 相当の計算をすればよい。ベクトル化して綺麗に計算する書き方が分からなかったので forループ。

Octaveで影付き、というか塗りつぶしをするにはどうしたら良いか。行きにμ+2σ側、帰りにμ-2σを通るような多角形と考えればfillで行けそう

式(6.63)に与えるパラメータを（手動で）いろいろ試して、それっぽく描かれる値を見つけた。やたー。
ぬか喜びだった。この影部分って標準偏差の２倍と明記してあるわけだし当然 $\mu\pm2\sigma$ で引くべきなのに $\mu\pm2\sigma^2$ で計算していた。

図表で使われているパラメータの探索そのものが機械学習の練習問題...

描画コードはsqrt(σ2)を使うだけだから修正は少ない。問題はパラメータ探索。
手動でパラメータを試すのめんどくさいよね。
制御点をいくつか選んで、目標（＝テキストに印刷してあるグラフの実測値w）からどのくらいずれているかのスコアリングを使ったsimulated annealing（焼きなまし法）でそれっぽいのが出来た。 $(\theta_0=0.33,\theta_1=25,\theta_2=1,\theta_3=0.03)$
f:id:n4_t:20121113201726p:plain

SAのコードも長いけど貼っておきますか

これでとりあえず「§6.4.2 ガウス過程による回帰」は攻略できたが、まだ１つ攻略すべき高い山が・・・

ガウス過程からのサンプリング

ガウス過程による事前分布からの「サンプル」ってさらっと言うけど何なの？

図6.4、図6.5 の見方がわからん。「ガウス過程による事前分布からのサンプル」ってキャプションに書いてあるけど、線形回帰ならともかく、ガウス過程からそんなんとれないでしょ？　うーん。それぞれのグラフ４色あるけど、同じ色は同じような傾向があるように見えない？　もしかして、それぞれの色ごとに何かサンプルが決まっていて、そのサンプルに対してこのカーネル使った場合のガウス過程回帰を描いているんじゃあない？　それだ！！！
PRML 読書会 #8 「カーネル法」 - Mi manca qualche giovedì?

３年前。全く訳がわからなかった。正直、一昨日まで全く訳がわかっていなかった。

そもそも「関数 y(x) の上での確率分布」って何だよ・・・

式(6.49)からまた読み直す。（何度目だこのサイクル）
任意の x_n, x_m の組、すなわちデータ集合がないとこのカーネル、というか共分散行列って定義できないんじゃないか？と思っていた。昨日までは。

ゴールは（この節の理解とかはもう諦めてて、ぶっちゃけ）図6.4, 図6.5を再現すること。

式(6.63) と、θ0〜θ3 の具体的な数値は与えられている。x の関数としてグラフを描くのは分かるのだが k(x_n, x_m) に何が来るんだろう。

あ。カーネル関数の計算に目標値 t_i は要らないんだった。xだけ適当に、linspace(-1,1,100)で与えたらどうなるんだろう。それって何を意味するんだろう。

そもそも共分散ってことは多次元だよね・・・
ハッ！x軸を１つの次元と考えるよりは、無限次元、ないし linspace(-1,1,100) なら独立した100の次元と考えるべきなのではないか！？*2 ビジョンが２次元→無限次元に広がった！

とすると、「μ と Σ さえ決まればよい」多次元ガウス分布としての y(x) が見えてきた。
μはとりあえず0で、Σは例の式(6.63)で求まる。
とすると、やるべき事はxの関数のプロットじゃない。多変量（100次元！）ガウス分布からのサンプリングだ。これはサンプリング法を扱ったPRML第11章に書いてある。

平均μ、共分散Σをもつ多変量ガウス分布に従うベクトル値の変数を生成するには、 $\sigma=LL^T$ の形を取るコレスキー分解 (Cholesky decomposition) を用いればよい (Press et al.,1992)。このとき、もし z がベクトル値の確率変数であり、その各要素が独立で、平均0、分散1のガウス分布に従うとすれば、y = μ + Lz は平均μ、共分散Σのガウス分布に従う。

コレスキー分解。何も考えないでライブラリ関数使うもんね楽勝じゃん！

・・・甘かった。chol(Σ) がエラーを返す。Σがpositive definiteじゃないぞと。

正定値？半正定値？なにそれ美味しいの？と思ってた３年前とは違うもんね

ここで下巻p.5の

関数 k(x, x') が有効なカーネルであるための必要十分条件は、任意の {x_n} に対して、要素が k(x_n, x_m) で与えられるグラム行列 K が半正定値であることである。なお、行列が半正定値であることと、行列のすべての要素が非負であることとは異なることに注意する。

という伏線を回収。というかそりゃそうだよねと。

でもカーネル関数もパラメータも与えられてしまっているのにどうしたら良いのか。p.19の最後〜p.21に

カーネル関数についての唯一の制約は、(6.62)で与えられる共分散行列が正定値でなければならないことである。λ_i を K の固有値とすると、C の対応する固有値は λ_i + β^-1 になる。したがって、カーネル行列 k(x_n, x_m) が、任意の x_n と x_m に対して半正定値、つまり、λ_i ≧0 であることを示せば十分である。これは、Kの固有値の中で、零であるような λ_i があった場合でも、β>0 であることから、Cの固有値はすべて正となるからである。この制約は、カーネル関数について以前議論したものと同じであり、6.2節で用いたテクニックを再び用いて、適切なカーネルを設計することができる。

ここで@naoya_tは姑息な手段に出る。

「λ_i が負なら、β^-1≧-λ_i にしちゃえばいいじゃない」

ランダムなサンプリングなのでPRMLの図表と全く同じにはならないけれど、雰囲気はかなり近い。

図6.4

左図: thetas = thetas_6_5(1,:); k = @my_kernel

f:id:n4_t:20121113211435p:plain

右図: k = @ornstein_uhlenbeck;

f:id:n4_t:20121113211447p:plain

図6.5

(1.00, 4.00, 0.00, 0.00)
f:id:n4_t:20121113211435p:plain
(9.00, 4.00, 0.00, 0.00)
f:id:n4_t:20121113211633p:plain
(1.00, 64.00, 0.00, 0.00)
f:id:n4_t:20121113212105p:plain

(1.00, 0.25, 0.00, 0.00)
f:id:n4_t:20121113212216p:plain
(1.00, 4.00, 10.00, 0.00)
f:id:n4_t:20121113212222p:plain
(1.00, 4.00, 0.00, 5.00)
f:id:n4_t:20121113212229p:plain

分からないことその２〜その∞が解決。風邪っぽいけれど今日はとても気分がよい。

あとがき+余談

図表の再現が出来て、おまけにガウス過程が何なのかが（実装に落とせるレベルにまで）見えてきたのが嬉しくて無駄に長いエントリになってしまった。

本記事の「ガウス過程からのサンプリング」は、深夜のカラオケルームでふと思いついてその場で実装してみた物。余談というかPRというか、来る12/15にPRML復々習レーン#7をやるのだけれど、その（前か）後にPRMLカラオケクラスタひみつ集会Vol.1というのをやります。PRML読みで、特に昨今のアニメ主題歌に造詣の深い方はぜひいらして下さい。

*1:PRML§11.2 マルコフ過程モンテカルロ

*2:この発想の転換が出来るようになったのはCourseraで量子力学を勉強したお陰だ。Observableとか

*3:3次元プロットはOctaveのやつよりmatplotlibの方がいいな

naoya_t@hatenablog

いわゆるチラシノウラであります