読者です 読者をやめる 読者になる 読者になる

naoya_t@hatenablog

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

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

グラフ描きたいだけ PRML Octave

「ガウス過程」- PRML §6.4

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

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

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

だから仕方がない。

本題に入ろう。

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

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

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

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

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

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

足掛かりになるのは shuyo さんの3年前の記事。
PRML6章「ガウス過程による回帰」を 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できるという事だよなー、と気づいて目の前が少し開けた。ビジョンが1次元→2次元に広がった感じ。分からないことその1が解決。

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

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

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

式(6.63)に与えるパラメータを(手動で)いろいろ試して、それっぽく描かれる値を見つけた。やたー。
ぬか喜びだった。この影部分って標準偏差の2倍と明記してあるわけだし当然\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 ガウス過程による回帰」は攻略できたが、まだ1つ攻略すべき高い山が・・・

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

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

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

PRML 読書会 #8 「カーネル法」 - Mi manca qualche giovedì?

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

そもそも「関数 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軸を1つの次元と考えるよりは、無限次元、ないし linspace(-1,1,100) なら独立した100の次元と考えるべきなのではないか!?*2 ビジョンが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じゃないぞと。

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

ここで下巻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

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

関連度自動決定 (ARD)

ガウス過程というものが何なのかが見えてきた勢いで、図6.9「ガウス過程における ARD事前分布からのサンプルを示したもの」を再現できるのでは?と思って試してみた。
式(6.71) って x が二次元になっただけだし計算難しくないし。

Octaveに3次元プロットをさせるのが初めてだったのでそちらで嵌った。*3

\eta_1=\eta_2=1
f:id:n4_t:20121113213118p:plain
\eta_1=1,\eta_2=0.01
f:id:n4_t:20121113213132p:plain

あとがき+余談

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

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

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

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

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