GPML：ノート1（1章、2.1節）

以下の本を読みます。私の誤りは私に帰属します。お気付きの点がありましたらご指摘いただけますと幸いです。

www.gaussianprocess.org

キャラクターの原作とは無関係です。

f:id:cookie-box:20180305232608p:plain:w60

「訓練データ $\mathcal{D} = \bigl\{ (x_i, y_i), \cdots, (x_n, y_n)\bigr\}$ から任意の未知の点 $x$ の値をどう予測すればよいか」という問を考えましょう。本質的に有限個の点の情報から無限個の点の情報を得ることはできません。考える関数 $f(x)$ が訓練データの外でどうあるべきかについて僕たちは何も知らない。無限にある関数候補に優劣を付けられない。なのでそこは何か決め打つしかない。それには主に2つのアプローチがあるといっています。

考慮するべき関数のクラスを制限する。
すべての関数に事前確率を割り当てる（Ex. より滑らかなものが好ましいなど）。

実際に用いられる方法は 1. の関数のクラスを制限することでしょう。テキストでは例えば線形モデルといっていますが、特定の構造のニューラルネットや決定木のことも多いかもしれません。その予め決めたクラスの関数の中で訓練データに最もフィットするものを選ぶわけですが、この方法では表現力が不十分だと予測性能が悪くなり、といって表現力をもたせると訓練データにオーバーフィットするとあります。

なので 2. が実現できないかという流れなのでしょう。こちらで関数のクラスを制限することは避けたい、ただある程度滑らかであることだけ要請しておきたい、と。しかし、こちらには即座にセルフ突っ込みが入っていますね。「有限の時間でどうやって無限の関数候補を検討するんだ」と。

―そこで、「ガウス過程が助けにきてくれる」と、ガウス過程が初登場します。

f:id:cookie-box:20180305231302p:plain:w60

ごめん、「助けにきてくれる」っていわれても、意味がわからなすぎて「熱い展開！」とか「これで勝てる！」とかならない…。

そうですね、まだイントロダクションなので続く記述はいくぶん抽象的ですが、ガウス過程があると助かるという心は以下だと思います。この箇所にかいていないことも適当に想像でかきましたが。

訓練データ外のことはわからないので、あらゆる関数の上の確率分布（＝確率過程）を考えたい。
ここでその確率過程はガウス過程であることを仮定する。そうすれば、距離が近い点どうしは共分散が大きいような分散共分散行列を選ぶことで、滑らかな関数の事前確率を大きくすることができる。そして、未知の点上の値の分布は直接求まるので、無限の関数候補を個々に検討する必要はない。

何というか、2. のアプローチは未知の点を予測するのに「あるクラスの関数の中から尤もらしいものを特定する」ことによってそうするのではなく、「点どうしに相関関係を入れる」ことによってそうするというとしっくりくる気がします。ここではそれを「すべての関数に事前確率を割り当てる（滑らかな関数の事前確率が大きいような）」といっていますが。

3ページから4ページの前半は Figure 1.1 で1次元→1次元の回帰のデモンストレーションをしていますね。(a) の事前分布が、2点を観測すると (b) の事後分布になると。そして4ページの中ほどにはまさに「ガウス過程回帰するとは共分散を見出すことである」というような記述があります。

4ページの後半からは分類問題ならどうかという話ですね。座標 $x$ に検出された天体が星であるのか銀河であるのかを分類したいというシチュエーションのようです。座標 $x$ の天体が星である確率 $\pi(x) \in [0,1]$ を予測したいと。しかし、ガウス過程のある座標での実現値は $[0,1]$ には収まりませんから、ロジスティック関数で変換するのが常套手段であるようですね。2次元平面上でのこの2値分類の事後分布が Figure 1.2 の (d) です。

後はこの本の章立てが紹介されていますね。各章の内容は概して以下でしょうか。

ガウス過程の定義や、回帰問題の予測値の計算方法（ノイズがガウシアンなら解析的）。
2値/多値分類問題（非線形な活性化関数を用いるためにもはや解析的でない）。
様々な共分散関数とその性質、組み合わせ方。
共分散関数のパラメータ推定方法。
ガウス過程回帰はカーネルマシンの1手法に位置付けられるが、SVMなど他のカーネルマシンの紹介と、ガウス過程回帰との関係。
理論的な話（漸近理論や、学習曲線や、PACベイズ推測の枠組みについて）。
$n \times n$ 行列の逆行列への対処法。
その他の問題設定。

じゃあまずはともかく Chapter 2 の「回帰」だな…序文に「ガウス過程モデルの解釈は色々ある」ってあるな。このテキストで扱うのは「関数空間」派と「重み空間」派ってことかな？とりあえず後者を先にみていくのか。えっと、訓練データが $\mathcal{D} = \bigl\{ (x_i, y_i), \cdots, (x_n, y_n)\bigr\}$ であるのは変わらなくて、入力は $D$ 次元、出力は1次元の実数で、入力ベクトルを束ねた $D \times n$ 行列 $X$ を計画行列とよぶと（これの転置を計画行列とする流儀の方が多いけど意図的に）。そしていまの状況は入力が与えられた下での出力の分布に興味がある（入力の分布自体には興味がない）のか。

それでまず、 $y=f(x)+\varepsilon, \; f(x) = x^\top w, \; \varepsilon \sim N(0, \sigma_n^2)$ というモデルを考えるのかな。これは入力の線形和に分散 $\sigma_n^2$ の独立なガウスノイズがのっているというモデルか。そうすると訓練データに対するモデルの尤度 $p(y|X,w)$ は式 (2.3) になって、 $w$ の事前分布を $w \sim N(0,\Sigma_p)$ として事後分布を求めるのかな？それがなんで式 (2.5) になるんだっけ？

式から $X$ による条件付けの部分を取った方がわかりやすいかもしれません。

$\displaystyle p(w|y)= \frac{p(y|w)p(w)}{p(y)}$

これなら $p(y|w)p(w)=p(w|y)p(y)$ を変形しただけですね。しかしいま考えているモデルは $w$ のみから $y$ を出すモデルではありませんから、 $X$ による条件付けを入れましょう。

$\displaystyle p(w|y, X)= \frac{p(y|w, X)p(w|X)}{p(y|X)}$

ここで $w$ の事前分布は計画行列 $X$ に応じて決めているのではないので $p(w|X)=p(w)$ とするだけです。それでは実際に事後分布を求めておきましょうか。

$\begin{split} \displaystyle p(w|X, y) &\propto \exp \left[-\frac{(y - X^\top w)^\top (y - X^\top w)}{2\sigma_n^2} - \frac{w^\top \Sigma_p^{-1} w}{2} \right] \\ &\propto \exp \left[-\frac{- y^\top X^\top w - w^\top X y + w^\top X X^\top w}{2\sigma_n^2} - \frac{w^\top \Sigma_p^{-1} w}{2} \right] \end{split}$

よって、これは以下の形に平方完成できますから、係数比較で $\bar{w}$ を求めましょう。

$\begin{split} \displaystyle & \exp \left[-\frac{1}{2}(w - \bar{w})^\top \left( \frac{X X^\top}{\sigma_n^2} + \Sigma_p^{-1} \right) (w - \bar{w}) \right]\\ &\propto \exp \left[-\frac{1}{2}w^\top \left( \frac{X X^\top}{\sigma_n^2} + \Sigma_p^{-1} \right) w +\frac{1}{2} \bar{w}^\top \left( \frac{X X^\top}{\sigma_n^2} + \Sigma_p^{-1} \right) w +\frac{1}{2} w^\top \left( \frac{X X^\top}{\sigma_n^2} + \Sigma_p^{-1} \right) \bar{w}\right] \end{split}$

見比べると $\bar{w} = \sigma_n^{-2}(\sigma_n^{-2} X X^\top +\Sigma_p^{-1} )^{-1}Xy$ であることがわかりますね。つまり、事後分布は平均が $\bar{w}$ で分散共分散行列が $(\sigma_n^{-2} X X^\top +\Sigma_p^{-1} )^{-1}$ のガウス分布です。ときにハヤト、線形モデル $f(x) = x^\top w$ を最小2乗フィッティングしたときの $\hat{w}$ ってどうなりましたっけ。

それは $\hat{w} = (X X^\top)^{-1}Xy$ だろ（下スライド）。リッジ正則化するなら $\hat{w} = (X X^\top + \lambda I)^{-1}Xy$ だっけ。あれ？これって $\bar{w}$ とすごい似てる？

似ているが混同するなと 10 ページにあります。

ええ…じゃあ訊くなよ…。

$\hat{w}$ はペナルティ項付きの尤度を最大にする $w$ であり、リッジ回帰の最終的な答えともいうべきものでしょう。他方、 $\bar{w}$ はベイズ推測の事後分布の平均であり最大点ですが、それだけなんです。事後分布の特徴的な点ではあるがこの点での予測値が特別な役割をもつのではないといいたいようです。このケースではたまたまモデルも事後分布も対称なので予測値の分布の平均値が $\bar{w}$ による予測値と一致するんですが。…続く 11 ページにも「未知の点を予測するときは、すべてのありうるパラメータでの予測値を事後分布による重みで平均します」とありますから、ベイズ推測とはあくまでそういうもので、事後分布のある点での予測値を使用するものではないということなのかと。しかし、渡辺ベイズ本の1章にも事後確率最大化推測とか平均プラグイン推測とかありましたから「特別な役割をもたせる人もいるのでは」となりそうなんですが、あくまで「ベイズ推測」というとそういうことではないということなのかもしれません。

ふーん…？まあ何にせよ、未知の点 $x_\ast$ に対する予測値 $f_\ast$ の分布は (2.9) 式で…何でこうなるんだっけ？

渡辺ベイズ本 1.2.3 節の「計算できる例」を参照すれば計算できるでしょう。計算できるわけですから。

あー…事前分布から事後分布への変換がハイパーパラメータの更新で表せるのをつかうんだっけか。それで、今回の場合は予測値の分布はガウス分布で、その平均は $x_\ast^\top \bar{w}$ で分散は $x_\ast^\top (\sigma_n^{-2} X X^\top +\Sigma_p^{-1} )^{-1} x_\ast$ になるのか。ん？「予測値の不確かさは入力の大きさが大きくなるほど大きくなり、これは線形モデルに期待する性質と合致している」みたいにあるけど、どういう意味？そんなこと線形モデルに期待してたっけ？

線形モデル $f(x) = x^\top w$ は $w$ が不確かなわけですから、 $f(x)$ の不確かさはノルムが大きい $x$ ではそれに比例して大きいだろうということなんでしょうか。今回考えた「入力の線形和に分散決め打ちのガウスノイズをのせるモデルを考え、重みベクトルをベイズ更新する」という方針ではそれが満たされていますね、というくらいでは。

まあいいや。Figure 2.1 にこの推測の図示があって、図 (a) の等高線が事前分布で、図 (b) が訓練データで（3点だけなんだな）、図 (c) の等高線は尤度で、図 (d) の等高線は事後分布で、図 (b) に戻って実線と点線が予測分布の信頼区間って感じか。図のキャプションをみると、図 (a)(c)(d) の等高線は1シグマ、2シグマで、図 (b) の破線は2シグマか。というか図 (c) と図 (d) めっちゃ似てるな。

事前分布がかかっているかどうかの違いしかありませんから、それは似ているでしょう。図 (c)(d) で特筆すべきなのは、図 (c) で切片の不確かさはまだ大きいのに対し傾きの不確かさはずっと小さいこと、図 (d) では傾きの広がりは図 (c) とほぼ変わらないのに対し切片は平均も分散も図 (c) からやや変わっていることですね（とキャプションにあるんですが）。この例では切片より傾きがずっとよく特定されます。図 (b) の直線について、切片を 0.5 増やしたときの尤度の減り方と傾きを 0.5 増やしたときの尤度の減り方を想像すればこれはわかるでしょう。

ちなみにプロデューサーさんが Figure 2.1 を再現しました（gist）。

f:id:cookie-box:20210130162410p:plain:w590

暇か。

ただもちろん線形モデルでは表現力が乏しいので、簡単に表現力をもたせる手法として入力 $x$ を基底関数 $\phi(x)$ で適当な高次元空間に送ればよい、と続く 2.1.2 節にありますね。 $\phi(x) = (1, x, x^2, x^3)^\top$ などとしてしまえばいいわけです。後の方の5章ではこの基底関数をどのように選ぶべきかという話題が出てくるようですが、さしあたり適当な基底関数が既に手に入っているものとします。いまや計画行列は $N \times n$ 行列 $\Phi = \Phi(X)$ となり、モデルは $f(x) = \phi(x)^\top w$ となりました。予測分布はどうなりますか？

どうなるって、 $X$ が $\Phi$ になるだけだよな。平均が $\phi(x_\ast)^\top \sigma_n^{-2}(\sigma_n^{-2} \Phi \Phi^\top +\Sigma_p^{-1} )^{-1} \Phi y$ で分散共分散行列が $\phi(x_\ast)^\top (\sigma_n^{-2} \Phi \Phi^\top +\Sigma_p^{-1} )^{-1} \phi(x_\ast)$ のガウス分布なんじゃ。

ええしかし、 $N \times N$ 行列 $\sigma_n^{-2} \Phi \Phi^\top +\Sigma_p^{-1}$ の逆行列を求めなければならないのがネックです。より表現力の高い基底関数を利用しようとするほど $N$ は大きいはずですし。

そうかもだけど仕方ないだろ。

思い出してください。多変量正規分布の積を平方完成した多変量正規分布の分散共分散行列を計算するときに工夫したことがありませんでしたか？

雑記： Schur 補行列と Sherman-Morrison-Woodbury の公式とカルマンフィルタとガウス過程回帰の整理 - クッキーの日記

そんなことやった覚えな…うわあやってる。

以下を今回の場合にあてはめてみましょう。

以下のサイズの行列 $A, \, B, \, C, \, D$ があるとき、
$A \in \mathbb{R}^{n \times n}, \quad B \in \mathbb{R}^{n \times m}$
$C \in \mathbb{R}^{m \times n}, \quad D \in \mathbb{R}^{m \times m}$ 以下の恒等式で逆行列がとられている行列に逆行列が存在するならば、以下の恒等式が成り立ちます。 $(A + BDC)^{-1} = A^{-1}-A^{-1}B(D^{-1}+CA^{-1}B)^{-1}CA^{-1}$ これを Sherman-Morrison-Woodbury の公式といいます。

	カルマンフィルタのフィルタ操作	線形基底関数モデルのベイズ推定
更新する分布	現在の観測に基づいて状態 $x_t$ の分布を前ステップでの一期先予測から更新したい	訓練データに基づいてパラメータ $w$ の分布を事前分布から事後分布に更新したい
$(A + BDC)^{-1}$	フィルタ分布の分散共分散行列 $({V_{t\|t-1}}^{-1} + {H_t}^{\top} {R_t}^{-1} H_t)^{-1}$ $= V_{t\|t-1}$ $\quad - V_{t\|t-1} {H_t}^{\top} (H_t V_{t\|t-1} {H_t}^{\top} + R_t)^{-1} H_t V_{t\|t-1}$	事後分布の分散共分散行列 $(\Sigma_p^{-1} + \sigma_n^{-2} \Phi \Phi^\top)^{-1}$ $= \Sigma_p$ $\quad - \Sigma_p \Phi (\Phi^\top \Sigma_p \Phi + \sigma_n^{2} I)^{-1} \Phi^\top \Sigma_p$
$A$	前回の一期先予測の分散共分散行列の逆行列 ${V_{t\|t-1}}^{-1}$	事前分布の分散共分散行列の逆行列 $\Sigma_p^{-1}$
$B = C^\top$	システムの転置 ${H_t}^{\top}$	計画行列 $\Phi$
$D$	観測ノイズの分散共分散行列の逆行列 ${R_t}^{-1}$	モデルの分散共分散行列の逆行列 $\sigma_n^{-2} I$

なので $n \times n$ 行列の逆行列を求めるので済むようになるわけです。 $n$ もデータ数なので大きいイメージがありますが、それより $N$ がずっと大きいケースを想定しています。こうしてみると、計画行列 $\Phi$ はパラメータ $w$ の側からみれば自分を $y$ に変換してくれるシステムなんですね。

というわけで改めて予測分布の分散共分散行列をかき直すとこうです。教科書に倣って $K = \Phi^\top \Sigma_p \Phi$ としました。

$\begin{split} \phi(x_\ast)^\top (\sigma_n^{-2} \Phi \Phi^\top +\Sigma_p^{-1} )^{-1} \phi(x_\ast) &= \phi(x_\ast)^\top \Sigma_p \phi(x_\ast) -\phi(x_\ast)^\top \Sigma_p \Phi (\Phi^\top \Sigma_p \Phi + \sigma_n^{2} I)^{-1} \Phi^\top \Sigma_p \phi(x_\ast) \\ &= \phi(x_\ast)^\top \Sigma_p \phi(x_\ast) -\phi(x_\ast)^\top \Sigma_p \Phi (K + \sigma_n^{2} I)^{-1} \Phi^\top \Sigma_p \phi(x_\ast) \end{split}$ それで、予測分布の平均については Sherman-Morrison-Woodbury の公式は要りませんでした。逆行列の前後に付いているのを逆行列の中に押し込めてみましょう。初手で $\Phi \Phi^\top$ を破壊できるんですよね。そもそも $(A + BDC)^{-1}$ の $BDC$ の部分が邪魔で崩したくて Woodbury の公式を使おうとなるんですが、勝手に崩れました。 $\begin{split} \phi(x_\ast)^\top \sigma_n^{-2}(\sigma_n^{-2} \Phi \Phi^\top +\Sigma_p^{-1} )^{-1} \Phi y &= \phi(x_\ast)^\top (\Phi^\top +\sigma_n^{2} \Phi^{-1} \Sigma_p^{-1} )^{-1} y \\ &= \phi(x_\ast)^\top ( (\Phi^\top \Sigma_p +\sigma_n^{2} \Phi^{-1}) \Sigma_p^{-1} )^{-1} y \\ &= \phi(x_\ast)^\top ( (K +\sigma_n^{2}I ) \Phi^{-1} \Sigma_p^{-1})^{-1} y \\ &= \phi(x_\ast)^\top \Sigma_p \Phi ( K +\sigma_n^{2}I )^{-1} y \end{split}$ なので、元より逆行列の中身が $n \times n$ 行列で済んでいます。

つづいたらつづく