[GPML] 線形回帰モデルからガウス過程回帰モデルを導く

前回執筆した線形回帰モデルを変形することにより、ガウス過程回帰モデルを導き出すことができるということを、後のためにきちんとまとめておこうと思います。

線形回帰モデル

線形回帰モデルについて簡単に要約しておきます。

ある関数$f$に$\boldsymbol x$を入力した結果が、特徴ベクトル$\boldsymbol \phi(\boldsymbol x) \in \R^{N+1}$と重みベクトル$\boldsymbol w \in \mathbb R^{N+1}$を用いて

$$f(\bm x) = \bm{\phi}(\bm{x})^\mathsf{T} \bm{w}$$

と書けるとし、これを繰り返し観測して出力$\boldsymbol y$を得るという行為が

$$\bm y \sim \mathcal N(\bm \Phi^\mathsf{T} \bm w, \sigma^2\bm I)$$

で書けるとします。ここで、

$$\begin{align*} \bm \Phi &= \begin{bmatrix} \bm \phi(\bm x_1) & \bm \phi(\bm x_2) & \cdots & \bm \phi(\bm x_n) \end{bmatrix} \\ &= \begin{bmatrix} 1 & 1 & \cdots & 1 \\ \phi_1(\bm x_1) & \phi_1(\bm x_2) & \cdots & \phi_1(\bm x_n) \\ \phi_2(\bm x_1) & \phi_2(\bm x_2) & \cdots & \phi_2(\bm x_n) \\ \vdots&\vdots&\ddots&\vdots \\ \phi_N(\bm x_1) & \phi_N(\bm x_2) & \cdots & \phi_N(\bm x_n) \end{bmatrix} \end{align*}$$

$$\bm y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}$$

です。重み$\boldsymbol w$は直接観測することができないので確率変数と仮定し、事前知識なしの状態では

$$p(\bm w) = \mathcal N(\bm 0, \bm V_\textrm{pr})$$

であると仮定します。すると、$\boldsymbol x_1, \boldsymbol x_2, \cdots, \boldsymbol x_n$を観測して$y_1, y_2, \cdots, y_n$を得たという知識を利用した上での$\boldsymbol w$の分布は、

$$p(\bm{w}|\bm{\Phi}, \bm{y}) = \mathcal N( \bm{m}_\textrm{pos}, \bm{V}_\textrm{pos} )$$

$$\begin{align*} \bm{V}_\textrm{pos} &= \left( \frac{1}{\sigma^2} \bm{\Phi} \bm{\Phi}^\mathsf{T} + \bm{V}_\textrm{pr}^{-1} \right)^{-1} , \\ \bm{m}_\textrm{pos} &= \frac{1}{\sigma^2} \bm{V}_\textrm{pos} \bm{\Phi} \bm{y} \\ &= \frac{1}{\sigma^2} \left( \frac{1}{\sigma^2} \bm{\Phi} \bm{\Phi}^\mathsf{T} + \bm{V}_\textrm{pr}^{-1} \right)^{-1} \bm{\Phi} \bm{y} \end{align*}$$

と表すことができます。これを用いることで、未知の点$\boldsymbol x_*$を入力したときの$f_* := f(\boldsymbol x_*)$の予測分布は

$$p(f_*|\bm \phi(\bm x_*), \bm \Phi, \bm y) = \mathcal N( \bm \phi(\bm x_*)^\mathsf{T} \bm m_\textrm{pos}, \bm \phi(\bm x_*)^\mathsf{T} \bm V_\textrm{pos} \bm \phi(\bm x_*))$$

$$\begin{align*} \mathbb E[f_*] &= \bm \phi (\bm x_*)^\mathsf{T} \bm m_\textrm{pos}, \\ \mathbb{V}[f_*] &= \bm \phi (\bm x_*)^\mathsf{T} \bm V_\textrm{pos} \bm \phi (\bm x_*) \end{align*}$$

と書くことができます。以降、簡単のために$\boldsymbol\phi_*:=\boldsymbol\phi(\boldsymbol x_*)$を用いて

$$\begin{align*} \mathbb E[f_*] &= \bm \phi_*^\mathsf{T} \bm m_\textrm{pos}, \\ \mathbb{V}[f_*] &= \bm \phi_*^\mathsf{T} \bm V_\textrm{pos} \bm \phi_* \end{align*}$$

と表すことにします。

予測分布の変形

上記の線形回帰モデルの$\boldsymbol w$の事後分布の共分散行列と平均ベクトル

$$\begin{align*} \bm{V}_\textrm{pos} &= \left(\frac{1}{\sigma^2}\bm{\Phi}\bm{\Phi}^\mathsf{T} + \bm{V}_\textrm{pr}^{-1} \right)^{-1} , \\ \bm{m}_\textrm{pos} &= \frac{1}{\sigma^2} \bm{V}_\textrm{pos} \bm{\Phi} \bm{y} \\ &= \frac{1}{\sigma^2} \left(\frac{1}{\sigma^2}\bm{\Phi}\bm{\Phi}^\mathsf{T} + \bm{V}_\textrm{pr}^{-1} \right)^{-1} \bm{\Phi} \bm{y} \end{align*}$$

において、$\boldsymbol V_\textrm{pos}^{-1} = \boldsymbol A$と書き直すと、

$$\begin{align*} \bm{A} &= \frac{1}{\sigma^2} \bm{\Phi}\bm{\Phi}^\mathsf{T} + \bm{V}_\textrm{pr}^{-1} , \\ \bm{V}_\textrm{pos} &= \bm{A}^{-1} , \\ \bm{m}_\textrm{pos} &= \frac{1}{\sigma^2} \bm{A}^{-1} \bm{\Phi} \bm{y} \end{align*}$$

と書けます。これを、$f_*$の予測分布の期待値と分散

$$\begin{align*} \mathbb E[f_*] &= \bm \phi_*^\mathsf{T} \bm m_\textrm{pos}, \\ \mathbb{V}[f_*] &= \bm \phi_*^\mathsf{T} \bm V_\textrm{pos} \bm \phi_* \end{align*}$$

にそれぞれ代入すると、

$$\begin{align*} \mathbb E[f_*] &= \bm \phi_*^\mathsf{T} \bm m_\textrm{pos} = \frac{1}{\sigma^2} \bm \phi_*^\mathsf{T} \bm{A}^{-1} \bm{\Phi} \bm{y} \\ \mathbb V[f_*] &= \bm \phi_*^\mathsf{T} \bm{V}_\textrm{pos} \bm \phi_* = \bm \phi_*^\mathsf{T} \bm{A}^{-1} \bm \phi_* \end{align*}$$

と書けます。

ここから、$\boldsymbol{A}^{-1}$を計算する代わりに、別の方法で期待値と分散を計算できるように変形してみます。

期待値

$\boldsymbol{A} \boldsymbol{V}_\textrm{pr} \boldsymbol{\Phi}$を計算します。

$$\begin{align*} \bm{A} \bm{V}_\textrm{pr} \bm{\Phi} &= \left(\frac{1}{\sigma^2}\bm{\Phi}\bm{\Phi}^\mathsf{T}+\bm{V}_\textrm{pr}^{-1}\right)\bm{V}_\textrm{pr}\bm{\Phi} \\ &= \frac{1}{\sigma^2} \bm{\Phi} \bm{\Phi}^\mathsf{T} \bm{V}_\textrm{pr} \bm{\Phi} + \bm{\Phi} \\ &= \frac{1}{\sigma^2} \bm{\Phi} \left( \bm{\Phi}^\mathsf{T} \bm{V}_\textrm{pr} \bm{\Phi} + \sigma^2\bm{I} \right) \end{align*}$$

ここに、$\boldsymbol{A}^{-1}$を左から、$(\boldsymbol{\Phi}^\mathsf{T}\boldsymbol{V}_\textrm{pr}\boldsymbol{\Phi}+\sigma^2\boldsymbol{I})^{-1}$を右から掛けると、

$$\bm{V}_\textrm{pr}\bm{\Phi}(\bm{\Phi}^\mathsf{T}\bm{V}_\textrm{pr}\bm{\Phi}+\sigma^2\bm{I})^{-1} = \frac{1}{\sigma^2}\bm{A}^{-1}\bm{\Phi}$$

更に、$\boldsymbol{x}_*^\mathsf{T}$を左から、$\boldsymbol{y}$を右から掛けると、

$$\bm{\phi}_*^\mathsf{T} \bm{V}_\textrm{pr} \bm{\Phi}( \bm{\Phi}^\mathsf{T} \bm{V}_\textrm{pr} \bm{\Phi} + \sigma^2\bm{I} )^{-1} \bm{y} = \frac{1}{\sigma^2} \bm{\phi}_*^\mathsf{T} \bm{A}^{-1} \bm{\Phi} \bm{y}$$

このことから、

$$\mathbb E[f_*] = \frac{1}{\sigma^2} \bm \phi_*^\mathsf{T} \bm{A}^{-1} \bm{\Phi} \bm{y} = \bm{\phi}_*^\mathsf{T} \bm{V}_\textrm{pr} \bm{\Phi}( \bm{\Phi}^\mathsf{T} \bm{V}_\textrm{pr} \bm{\Phi} + \sigma^2 \bm{I} )^{-1} \bm{y}$$

と書けます。

共分散行列

Sherman & Morisonの公式

$$( \bm{Z} + \bm{U} \bm{W} \bm{V}^\mathsf{T} )^{-1} = \bm{Z}^{-1} - \bm{Z}^{-1} \bm{U} ( \bm{W}^{-1} + \bm{V}^\mathsf{T} \bm{Z}^{-1} \bm{U} )^{-1} \bm{V}^\mathsf{T} \bm{Z}^{-1}$$

に、$\boldsymbol{Z} = \boldsymbol{V}_\textrm{pr}^{-1}$、$\boldsymbol{U} =\boldsymbol{V} =\boldsymbol{\Phi}$、および$\boldsymbol{W} =\boldsymbol{I}/\sigma^2$を代入すると、

$$\left( \bm{V}_\textrm{pr}^{-1} + \frac{1}{\sigma^2} \bm{\Phi} \bm{\Phi}^\mathsf{T} \right)^{-1} = \bm{V}_\textrm{pr}-\bm{V}_\textrm{pr} \bm{\Phi} ( \sigma^2 \bm{I} + \bm{\Phi}^\mathsf{T} \bm{V}_\textrm{pr} \bm{\Phi} )^{-1} \bm{\Phi}^\mathsf{T} \bm{V}_\textrm{pr}$$

左辺は$\boldsymbol{A}^{-1}$に一致するので、

$$\bm{A}^{-1} = \bm{V}_\textrm{pr}-\bm{V}_\textrm{pr} \bm{\Phi} ( \bm{\Phi}^\mathsf{T} \bm{V}_\textrm{pr} \bm{\Phi} + \sigma^2 \bm{I})^{-1} \bm{\Phi}^\mathsf{T} \bm{V}_\textrm{pr}$$

$\boldsymbol{\phi}_*^\mathsf{T}$を左から、$\boldsymbol{\phi}_*$を右から掛けると、

$$\bm \phi_*^\mathsf{T} \bm{A}^{-1} \bm \phi_* = \bm \phi_*^\mathsf{T} \bm{V}_\textrm{pr} \bm \phi_* - \bm \phi_*^\mathsf{T} \bm{V}_\textrm{pr} \bm{\Phi}( \bm{\Phi}^\mathsf{T} \bm{V}_\textrm{pr} \bm{\Phi}+\sigma^2 \bm{I})^{-1} \bm{\Phi}^\mathsf{T} \bm{V}_\textrm{pr} \bm \phi_*$$

このことから、

$$\mathbb V[f_*] = \bm{\phi}_*^\mathsf{T} \bm{A}^{-1} \bm{\phi}_* = \bm{\phi}_*^\mathsf{T} \bm{V}_\textrm{pr} \bm{\phi}_* - \bm{\phi}_*^\mathsf{T} \bm{V}_\textrm{pr} \bm{\Phi}( \bm{\Phi}^\mathsf{T} \bm{V}_\textrm{pr} \bm{\Phi} + \sigma^2 \bm{I} )^{-1} \bm{\Phi}^\mathsf{T} \bm{V}_\textrm{pr} \bm{\phi}_*$$

まとめると

以上のことから、

$$\begin{align*} \mathbb E[f_*] &= \boldsymbol{\phi}_*^\mathsf{T} \boldsymbol{V}_\textrm{pr} \boldsymbol{\Phi} ( \boldsymbol{\Phi}^\mathsf{T} \boldsymbol{V}_\textrm{pr} \boldsymbol{\Phi} + \sigma^2 \boldsymbol{I} )^{-1} \boldsymbol{y} \\ \mathbb V[f_*] &= \boldsymbol{\phi}_*^\mathsf{T} \boldsymbol{V}_\textrm{pr} \boldsymbol{\phi}_* - \boldsymbol{\phi}_*^\mathsf{T} \boldsymbol{V}_\textrm{pr} \boldsymbol{\Phi} ( \boldsymbol{\Phi}^\mathsf{T} \boldsymbol{V}_\textrm{pr} \boldsymbol{\Phi} + \sigma^2 \boldsymbol{I} )^{-1} \boldsymbol{\Phi}^\mathsf{T} \boldsymbol{V}_\textrm{pr} \boldsymbol{\phi}_* \end{align*}$$

と変形できました。

さらなる変形

ここで改めて、予測分布と計画行列が、

$$p(f_*|\bm{\phi}_*, \bm X, \bm y) = \mathcal N( \bm{\phi}_*^\mathsf{T} \bm m_\textrm{pos} , \bm{\phi}_*^\mathsf{T} \bm V_\textrm{pos} \bm{\phi}_*)$$

$$\mathbb V[f_*] = \bm \phi_*^\mathsf{T} \bm V_\textrm{pos} \bm \phi_*$$

$$\bm \Phi = \begin{bmatrix} \bm{\phi}_1 & \bm{\phi}_2 & \cdots & \bm{\phi}_n \end{bmatrix}$$

であることに注意して、導いた式を見てみましょう:

$$\begin{align*} \mathbb E[f_*] &= \boldsymbol{\phi}_*^\mathsf{T} \boldsymbol{V}_\textrm{pr} \boldsymbol{\Phi} ( \boldsymbol{\Phi}^\mathsf{T} \boldsymbol{V}_\textrm{pr} \boldsymbol{\Phi} + \sigma^2 \boldsymbol{I} )^{-1} \boldsymbol{y} \\ \mathbb V[f_*] &= \boldsymbol{\phi}_*^\mathsf{T} \boldsymbol{V}_\textrm{pr} \boldsymbol{\phi}_* - \boldsymbol{\phi}_*^\mathsf{T} \boldsymbol{V}_\textrm{pr} \boldsymbol{\Phi} ( \boldsymbol{\Phi}^\mathsf{T} \boldsymbol{V}_\textrm{pr} \boldsymbol{\Phi} + \sigma^2 \boldsymbol{I} )^{-1} \boldsymbol{\Phi}^\mathsf{T} \boldsymbol{V}_\textrm{pr} \boldsymbol{\phi}_* \end{align*}$$

すると、右辺のほとんどが、

$$\bm{\phi}^\mathsf{T} \bm V_\textrm{pr} \bm{\phi}'$$

で構成されていることに気づくと思います。実際、

$$\begin{align*} \bm{\phi}_*^\mathsf{T} \bm V_{\textrm{pr}} \bm \Phi &= \bm{\phi}_*^\mathsf{T} \bm V_{\textrm{pr}} \begin{bmatrix} \bm{\phi}_1 & \bm{\phi}_2 & \cdots & \bm{\phi}_n \end{bmatrix} \\ &= \begin{bmatrix} \bm{\phi}_*^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_1 & \bm{\phi}_*^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_2 & \cdots & \bm{\phi}_*^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_n \end{bmatrix} , \\ \\ \bm \Phi^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_* &= \begin{bmatrix} \bm{\phi}_1 \\ \bm{\phi}_2 \\ \vdots \\ \bm{\phi}_n \end{bmatrix} \bm V_{\textrm{pr}} \bm{\phi}_* \\ &= \begin{bmatrix} \bm{\phi}_*^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_1 \\ \bm{\phi}_*^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_2 \\ \vdots \\ \bm{\phi}_*^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_n \end{bmatrix} , \\ \\ \bm \Phi^\mathsf{T} \bm V_{\textrm{pr}} \bm \Phi &= \begin{bmatrix} \bm{\phi}_1 \\ \bm{\phi}_2 \\ \vdots \\ \bm{\phi}_n \end{bmatrix} \bm V_{\textrm{pr}} \begin{bmatrix} \bm{\phi}_1 & \bm{\phi}_2 & \cdots & \bm{\phi}_n \end{bmatrix} \\ &= \begin{bmatrix} \bm{\phi}_1^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_1 & \bm{\phi}_1^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_2 & \cdots & \bm{\phi}_1^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_n \\ \bm{\phi}_2^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_1 & \bm{\phi}_2^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_2 & \cdots & \bm{\phi}_2^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_n \\ \vdots & \vdots & \ddots & \vdots \\ \bm{\phi}_n^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_1 & \bm{\phi}_n^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_2 & \cdots & \bm{\phi}_n^\mathsf{T} \bm V_{\textrm{pr}} \bm{\phi}_n \end{bmatrix} \end{align*}$$

となっています。そこで

$$k(\bm{x}, \bm{x}') = \bm{\phi}^\mathsf{T} \bm V_\textrm{pr} \bm{\phi}' = \bm{\phi}(\bm x)^\mathsf{T} \bm V_\textrm{pr} \bm{\phi}(\bm x')$$

と定義し、さらに

$$\begin{align*} \boldsymbol{\phi}_*^\mathsf{T} \boldsymbol V_{\textrm{pr}} \boldsymbol{\phi}_* &= k_{**} &&(1 \times 1) \\ \boldsymbol{\phi}_*^\mathsf{T} \boldsymbol V_{\textrm{pr}} \boldsymbol{\Phi} &= \boldsymbol k_*^\mathsf{T} &&(1 \times n) \\ \boldsymbol{\Phi}^\mathsf{T} \boldsymbol V_{\textrm{pr}} \boldsymbol{\phi}_* &= \boldsymbol k_* &&(n \times 1) \\ \boldsymbol{\Phi}^\mathsf{T} \boldsymbol V_{\textrm{pr}} \boldsymbol{\Phi} &= \boldsymbol K &&(n \times n) \end{align*}$$

と書き換えてみます。すると、

$$\begin{align*} \mathbb E[f_*] &= \bm{k}_*^\mathsf{T} (\bm{K} + \sigma^2 \bm{I})^{-1} \bm{y} \\ \mathbb V[f_*] &= k_{**} - \bm{k}_*^\mathsf{T} (\bm{K} + \sigma^2 \bm{I})^{-1} \bm{k}_* \end{align*}$$

と書けます。

ガウス過程へ

こうして導き出された

$$\begin{align*} \mathbb E[f_*] &= \bm{k}_*^\mathsf{T} (\bm{K} + \sigma^2 \bm{I})^{-1} \bm{y} \\ \mathbb V[f_*] &= k_{**} - \bm{k}_*^\mathsf{T} (\bm{K} + \sigma^2 \bm{I})^{-1} \bm{k}_* \end{align*}$$

を、多変量ガウス分布

$$p(\bm{x}_1, \bm{x}_2) = \mathcal N\!\left( \begin{bmatrix} \bm{x}_1 \\ \bm{x}_2 \end{bmatrix} \middle| \begin{bmatrix} \bm{0} \\ \bm{0} \end{bmatrix} , \begin{bmatrix} \bm{V}_1 & \bm{V}_{12} \\ \bm{V}_{12}^\mathsf{T} & \bm{V}_2 \end{bmatrix} \right)$$

の条件付き確率の公式

$$p(\bm{x}_2|\bm{x}_1) = \mathcal N(\bm x_2| \bm{V}_{12}^\mathsf{T}\bm{V}_1^{-1}\bm{x}_1, \bm{V}_2 - \bm{V}_{12}^\mathsf{T}\bm{V}_1^{-1} \bm{V}_{12})$$

$$\begin{align*} \mathbb E[\bm{x}_2|\bm{x}_1] &= \bm{V}_{12}^\mathsf{T}\bm{V}_1^{-1} \bm{x}_1 \\ \mathrm{Cov}[\bm{x}_2|\bm{x}_1] &= \bm{V}_2 - \bm{V}_{12}^\mathsf{T}\bm{V}_1^{-1} \bm{V}_{12} \end{align*}$$

と比較してみましょう。かなり似ていることに気付けると思います。この公式については次の記事で導出しているので参考にしてください:

https://zenn.dev/mory22k/articles/ec8ff69576fbbe

公式へのはめ込み

この公式に当てはめられるように、

$$p(\bm{y}, f_*|\bm \phi_*, \bm \Phi) = \mathcal N\!\left( \begin{bmatrix}\bm{y} \\f_* \end{bmatrix} \middle| \begin{bmatrix}\bm{0} \\\bm{0} \end{bmatrix},\begin{bmatrix}\bm{K}+\sigma^2\bm I &\bm{k}_*\\\bm{k}_*^\mathsf{T} &k_{**}\end{bmatrix} \right)$$

$$p(\boldsymbol w) = \mathcal N(\boldsymbol 0, \boldsymbol V_\textrm{pr})$$

$$k(\bm{x}, \bm{x}') = \bm \phi^\mathsf{T} \bm V_\textrm{pr} \bm \phi'$$

と定めてみます。これをもとに、改めて公式を適用してみると、

$$p(f_*|\bm \phi_*, \bm \Phi, \bm{y}) = \mathcal N( f_*| \bm{k}_*^\mathsf{T} (\bm{K} + \sigma^2 \bm{I})^{-1} \bm{y} , k_{**} - \bm{k}_*^\mathsf{T} (\bm{K} + \sigma^2 \bm{I})^{-1} \bm{k}_* )$$

$$\begin{align*} \mathbb E[f_*|\boldsymbol y] &= \boldsymbol{k}_*^\mathsf{T} (\boldsymbol{K} + \sigma^2 \boldsymbol{I})^{-1} \boldsymbol{y} \\ \mathbb V[f_*|\boldsymbol y] &= k_{**} - \boldsymbol{k}_*^\mathsf{T} (\boldsymbol{K} + \sigma^2 \boldsymbol{I})^{-1} \boldsymbol{k}_* \end{align*}$$

となります。線形回帰モデルを作成した際の予測分布の計算結果と一致する結果が得られましたね。

ガウス過程を仮定する

問題は、仮定した

$$p(\bm{y}, f_*|\bm \phi_*, \bm \Phi) = \mathcal N\!\left( \begin{bmatrix}\bm{y} \\f_* \end{bmatrix} \middle| \begin{bmatrix}\bm{0} \\\bm{0} \end{bmatrix},\begin{bmatrix}\bm{K}+\sigma^2\bm I &\bm{k}_*\\\bm{k}_*^\mathsf{T} &k_{**}\end{bmatrix} \right)$$

がもつ意味です。簡単のために、$f$がノイズなし、すなわち$\sigma = 0$であるような場合を考えてみましょう。

$$p(\bm{f}, f_*|\bm \phi_*, \bm \Phi) = \mathcal N\!\left( \begin{bmatrix} \bm{f} \\ f_* \end{bmatrix} \middle| \begin{bmatrix} \bm{0} \\ \bm{0} \end{bmatrix} , \begin{bmatrix} \bm{K} & \bm{k}_* \\ \bm{k}_*^\mathsf{T} & k_{**} \end{bmatrix} \right)$$

このうちの、特に$\boldsymbol f$だけに注目します:

$$p(\bm{f}|\bm X) = \mathcal N\!\left( \bm{f} \middle| \bm{0} , \bm{K} \right)$$

これを素直に読み解くと、どうやら

1. 関数$f$に$\boldsymbol x_1, \boldsymbol x_2, \cdots, \boldsymbol x_n$を入力して得られる値$f(\boldsymbol x_1), f(\boldsymbol x_2), \cdots, f(\boldsymbol x_n)$は、多変量ガウス分布に従う
2. 1. で述べたことは、任意の正の整数$n$に対して成り立つ
3. データを得る前の情報として、平均ベクトルに$\boldsymbol 0$、共分散行列に$[\boldsymbol K]_{pq} = k(\boldsymbol x_p, \boldsymbol x_q)$を与えている

ということを仮定しているようです。任意の個数の出力がすべて多変量ガウス分布に従うという仮定はなかなか面白いもので、後述するように、この仮定によって、関数$k(\boldsymbol x_p, \boldsymbol x_q)$の設計次第でかなりの種類の$f$ (定義域は連続的に限る) の回帰分析が行えます。

ガウス過程

さて、このような仮定が成り立つとき、$f$はガウス過程に従うといい、

$$f \sim \mathcal{GP}(0, k(\bm x, \bm x'))$$

などと書きます。$k(\boldsymbol x, \boldsymbol x')$は共分散関数と呼ばれ、

$$\bm K = \begin{bmatrix} k(\bm x_1, \bm x_1) & k(\bm x_1, \bm x_2) & \cdots & k(\bm x_1, \bm x_n) \\ k(\bm x_2, \bm x_1) & k(\bm x_2, \bm x_2) & \cdots & k(\bm x_2, \bm x_n) \\ \vdots & \vdots & \ddots & \vdots \\ k(\bm x_n, \bm x_1) & k(\bm x_n, \bm x_2) & \cdots & k(\bm x_n, \bm x_n) \end{bmatrix}$$

をグラム行列と呼びます。ガウス過程を回帰モデルとして捉える文脈では、それぞれカーネル関数、カーネル行列という呼び方もされます。また、もしも平均ベクトルとして$\boldsymbol 0$ではなく

$$\bm m = \begin{bmatrix} m(\bm x_1) \\ m(\bm x_2) \\ \vdots \\ m(\bm x_n) \end{bmatrix}$$

であるようなベクトル$\boldsymbol m$を指定した場合、

$$f \sim \mathcal{GP}(m(\bm x), k(\bm x, \bm x'))$$

と書くことになります。$m(\boldsymbol x)$は平均関数などと呼ばれます。

ガウス過程回帰モデルと予測分布

今までの話を別の角度から見てみましょう。まず、$f$の任意の関数値を並べたベクトル

$$\bm f = \begin{bmatrix} f(\bm x_1) \\ f(\bm x_2) \\ \vdots \\ f(\bm x_n) \end{bmatrix}$$

は、多変量ガウス分布に従うと仮定します。

$$p(\bm{f}|X) = \mathcal N\!\left( \bm{f} \middle| \bm{0} , \bm{K} \right)$$

そして、

1. $\boldsymbol f$のうちの最初の$k$個を観測済みの値$\boldsymbol y$、残りの$(n-k)$個を未観測の値$\boldsymbol f_*$に分ける
2. $\boldsymbol y$に対応する入力値の集合を$X = \{\bm x_1, \bm x_2, \cdots, \bm x_k\}$、カーネル行列を$\boldsymbol K_f$とする
3. $\boldsymbol f_*$に対応する入力値を$X_* = \{\bm x_{k+1}, \bm x_{k+2}, \cdots, \bm x_{n}\}$、カーネル行列を$\boldsymbol K_{**}$とする
4. 両者の間の共分散行列のうち右上のほうを$\boldsymbol K_{f*}$と書くことにする

つまり、

$$p(\bm{y}, \bm{f}_*|X_*, X) = \mathcal N\!\left( \begin{bmatrix} \bm{y} \\ \bm{f}_* \end{bmatrix} \middle| \begin{bmatrix} \bm{0} \\ \bm{0} \end{bmatrix} , \begin{bmatrix} \bm{K}_f & \bm{K}_{f*} \\ \bm{K}_{f*}^\mathsf{T} & \bm K_{**} \end{bmatrix} \right)$$

と表すことにします。すると、多変量ガウス分布の公式

$$p(\bm{x}_1, \bm{x}_2) = \mathcal N\!\left( \begin{bmatrix} \bm{x}_1 \\ \bm{x}_2 \end{bmatrix} \middle| \begin{bmatrix} \bm{0} \\ \bm{0} \end{bmatrix} , \begin{bmatrix} \bm{V}_1 & \bm{V}_{12} \\ \bm{V}_{12}^\mathsf{T} & \bm{V}_2 \end{bmatrix} \right)$$

$$p(\bm{x}_2|\bm{x}_1) = \mathcal N( \bm x_2 | \bm{V}_{12}^\mathsf{T} \bm{V}_1^{-1} \bm{x}_1 , \bm{V}_2 - \bm{V}_{12}^\mathsf{T} \bm{V}_1^{-1} \bm{V}_{12} )$$

$$\begin{align*} \mathbb E[\bm{x}_2|\bm{x}_1] &= \bm{V}_{12}^\mathsf{T}\bm{V}_1^{-1} \bm{x}_1 \\ \mathrm{Cov}[\bm{x}_2|\bm{x}_1] &= \bm{V}_2 - \bm{V}_{12}^\mathsf{T}\bm{V}_1^{-1} \bm{V}_{12} \end{align*}$$

により、$\boldsymbol y$が定まったもとでの$\boldsymbol f_*$の分布は

$$p(\bm{f}_*|\bm{y}, X_*, X) = \mathcal N( \bm{K}_{f*}^\mathsf{T} \bm{K}_f^{-1} \bm{y} , \bm{K}_{**} - \bm{K}_{f*}^\mathsf{T} \bm{K}_f^{-1} \bm{K}_{f*} )$$

と書くことができます。ここでさらに、$\bm y = [y_1, y_2, \cdots, y_n]^\mathsf{T}$と$X = \{\bm x_1, \bm x_2, \cdots, \bm x_n\}$を

$$\mathcal D_n = \{ (\bm x_1, y_1), (\bm x_2, y_2), \cdots, (\bm x_n, y_n) \}$$

とまとめて書くことにすると、

$$p(\bm{f}_*|\mathcal D_n, X_*) = \mathcal N( \bm{K}_{f*}^\mathsf{T} \bm{K}_f^{-1} \bm{y} , \bm{K}_{**} - \bm{K}_{f*}^\mathsf{T} \bm{K}_f^{-1} \bm{K}_{f*} )$$

となります。これが、ガウス過程回帰モデルを調べたときによく出てくる予測分布の式の正体です。

参考文献

1. Carl Edward Rasmussen and Christopher K. I. Williams: "Gaussian Process for Machine Learning," 2. Regression, https://gaussianprocess.org/gpml/chapters/, The MIT Press, 2006.
2. 大羽成征, 持橋大地: 「ガウス過程と機械学習」, 機械学習プロフェッショナルシリーズ, 講談社