元バイオ系

元バイオウェット系がデータサイエンスやらを勉強していくブログ。 基本自分用のまとめ。

ガウス過程を理解したい1

時間が空きましたが、ガウス過程を理解するために多変量ガウス分布について学習してきました(以下を参照)。

hotoke-x.hatenablog.com

hotoke-x.hatenablog.com

ガウス過程と機械学習の第5章、補助変数法でよくわからなくなってしまったので原著を読みました。

原著論文の紹介をいきなりしても良いのですが、論文のモチベーションがわかるようにまずはガウス過程の導入から行っていきます。

登場する用語の整理

ガウス過程に入る前に、用語の整理をしておきます。

  1. 動径基底関数

    $$ \begin{align} \exp \left(- \frac{\left(x - \mu \right)^{2}}{\sigma^{2}} \right) \end{align} $$

    のような形の関数のこと。(オイラーの公式で三角関数に分解できることを想像すると、名前に納得感があると思った。)

  2. 動径基底関数(radial basis function, RBF)回帰
    線形回帰モデルの基底関数に、動径基底関数を用いたモデル。

    $$ \begin{align} y = \sum_{h=-H}^{H} w_{h} \exp \left(- \frac{\left(x - \mu_h \right)^{2}}{\sigma^{2}} \right) \label{reg:rbf} \end{align} $$

    素直にこの式に従うとデータが低次元じゃないと計算できません(次元の呪い)。たとえば、1次元なら任意の曲線を表現するために \displaystyle xの値を細かくとれば良いですが、2次元で任意の局面を表現するためには細かいメッシュを切る必要があります。

ガウス過程とは

任意次元の入力に対する出力の同時分布が多変量ガウス分布に従うもののこと。

ガウス過程のポイント

  1.  wに関数を導入したことで、線形回帰モデルの任意の点(つまり無限点)の基底関数を表現可能になった
  2. 式\eqref{reg:rbf}において、 \displaystyle wについて全範囲 \displaystyle [-\infty, \infty]を解析的に積分することで、観測された次元に限って他を積分消去(周辺化)してしまう(扱うのは有限次元で済む)
  3. 周辺化によってデータ点がある次元だけで事後分布を表現できる

無限点とは何のことか、周辺化が何をしているのかのイメージをつかみたければガウス過程と機械学習 p66の図3.5を参照。ちゃんと理解したい人は再生核ヒルベルト空間について調べましょう。

ちゃんと式で考えていきます。

ちゃんと式で考えていきます。入力 \displaystyle \boldsymbol xに対して

$$ \begin{align} \left(\begin{array}{c}{\widehat{y}_{1}} \\ {\widehat{y}_{2}} \\ {\vdots} \\ {\widehat{y}_{N}}\end{array}\right)&=\left(\begin{array}{cccc}{\phi_{0}\left(\boldsymbol{x}_{1}\right)} & {\phi_{1}\left(\boldsymbol{x}_{1}\right)} & {\cdots} & {\phi_{H}\left(\boldsymbol{x}_{1}\right)} \\ {\phi_{0}\left(\boldsymbol{x}_{2}\right)} & {\phi_{1}\left(\boldsymbol{x}_{2}\right)} & {\cdots} & {\phi_{H}\left(\boldsymbol{x}_{2}\right)} \\ {\vdots} & {\vdots} & {\ddots} & {\vdots} \\ {\phi_{0}\left(\boldsymbol{x}_{N}\right)} & {\phi_{1}\left(\boldsymbol{x}_{N}\right)} & {\cdots} & {\phi_{H}\left(\boldsymbol{x}_{N}\right)}\end{array}\right)\left(\begin{array}{c}{w_{0}} \\ {w_{1}} \\ {\vdots} \\ {w_{H}}\end{array}\right) \\ \boldsymbol{\widehat y} &= \boldsymbol{\Phi w} \end{align} $$

とすると、

$$ \begin{align} \mathbb{E}\left[\boldsymbol{y}\right] &= \mathbb{E}\left[\boldsymbol{\Phi w}\right] = \boldsymbol{\Phi} \mathbb{E}\left[\boldsymbol{w}\right] = \boldsymbol{0} \\ \Sigma &=\mathbb{E}\left[\boldsymbol{y y}^{\mathrm T}\right]-\mathbb{E}[\boldsymbol{y}] \mathbb{E}[\boldsymbol{y}]^{\mathrm T} \\ &= \mathbb{E}\left[(\boldsymbol{\Phi} \boldsymbol{w})(\boldsymbol{\Phi} \boldsymbol{w})^{\mathrm T}\right] \\ &= \boldsymbol{\Phi} \mathbb{E}\left[\boldsymbol{w} \boldsymbol{w}^{\mathrm T}\right] \boldsymbol{\Phi}^{\mathrm T} = \lambda^{2} \boldsymbol{\Phi} \Phi^{\mathrm T} \end{align} $$

より、

$$ \begin{align} y \sim \mathcal{N} \left(\boldsymbol{0}, \lambda^{2}\boldsymbol{\Phi \Phi^{\mathrm T}} \right) \end{align} $$

となります。ここでは、 \boldsymbol{w} \sim \mathcal{N}(\boldsymbol{0}, \lambda^{2} \boldsymbol{I})としました。 \displaystyle \boldsymbol{w}について知らなくてよい点で強力です。

分散共分散行列(今回のような \displaystyle A^{\mathrm T}Aの形の行列をグラム行列ともいう) \displaystyle \boldsymbol{K} = \lambda^{2}\Phi\Phi^{\mathrm T} \displaystyle n, n'成分を見てみると、

$$ \begin{align} \boldsymbol{K}_{n, n'} = \lambda^{2} \phi(\boldsymbol{x_{n}})^{\mathrm T}\phi(\boldsymbol{x_{n'}}) \end{align} $$

となっており、基底関数の内積の定数倍となっています。 \displaystyle \boldsymbol{K}は分散共分散行列なので、正則な対称行列であることに注意しましょう。

ちなみに、データ点が増えすぎると計算量とメモリが爆発するので一般には計算方法を工夫する必要があります。その手法として補助変数法があります(後述)。

カーネルトリック

上記の内積を、 \displaystyle \boldsymbol{x_{n}} \displaystyle \boldsymbol{x_{n}}のカーネル関数(kernel function)と呼び、以下のように定義する。

$$ \begin{align} k(\boldsymbol{x}_{n}, \boldsymbol{x}_{n '}) = \phi(\boldsymbol{x_{n}})^{\mathrm T}\phi(\boldsymbol{x_{n'}}) \end{align} $$

 \boldsymbol{w}は積分消去されましたが、このままだと特徴ベクトル \phi(\boldsymbol x)が無限次元の場合表現できません。 そこで、 \phi(\cdot)を直接表現することを避け、内積の結果の式だけを利用します(カーネルトリック)。

ガウス過程回帰

カーネル関数

以下のような観測モデルを考えます。

$$ \begin{align} \left\{\begin{array}{l} y_n &= f\left(\boldsymbol{x_n}\right) + \epsilon_n \\ \epsilon &\sim \mathcal{N}(0, \sigma^2) \end{array}\right. \end{align} $$

すると、 \displaystyle X=(\boldsymbol{x}_1, \boldsymbol{x}_2, \ldots, \boldsymbol{x}_3)が与えられたとき、 \displaystyle \boldsymbol{y}=(y_1, y_2, \ldots, y_n)の確率分布は

$$ \begin{align} p(\mathbf{y} | \mathbf{X}) &=\int p(\mathbf{y}, \mathbf{f} | \mathbf{X}) d \mathbf{f} \\ &=\int p(\mathbf{y} | \mathbf{X}) p(\mathbf{f} | \mathbf{X}) d \mathbf{f} \\ &=\int p(\mathbf{y} | \mathbf{f}) p(\mathbf{f} | \mathbf{X}) d \mathbf{f} \\ &=\int \mathcal{N}\left(\mathbf{y} | \mathbf{f}, \sigma^{2} \mathbf{I}\right) \mathcal{N}(\mathbf{f} | \boldsymbol{\mu}, \mathbf{K}) d \mathbf{f} \end{align} $$

なお、 \displaystyle p(\mathbf{y} | \mathbf{X}) = p(\mathbf{y} | \mathbf{f}) を利用しました。ガウス分布の再生性より

$$ \begin{align} p(\mathbf{y} | \mathbf{X}) &= \int \mathcal{N}\left(\boldsymbol{\mu},\boldsymbol{K} + \sigma^2 \boldsymbol{I} \right) \end{align} $$

となります。ガウス分布の再生性は積率母関数を使えば容易に証明可能です。以上より、 \boldsymbol{y}はガウス過程に従い、そのカーネル関数は

$$ \begin{align} k^{\prime}\left(\mathbf{x}_{n}, \mathbf{x}_{n^{\prime}}\right)=k\left(\mathbf{x}_{n}, \mathbf{x}_{n^{\prime}}\right)+\sigma^{2} \delta\left(n, n^{\prime}\right) \end{align} $$

となることがわかります。 \deltaはクロネッカーのデルタです。

予測(欠損値補完など)

非観測点 \boldsymbol{x}^*における出力 y^*を予測するためには、 \boldsymbol{x}^*を含めたカーネル行列

$$ \begin{align} \left(\begin{array}{l}{\mathbf{y}} \\ {y^{*}}\end{array}\right) \sim \mathcal{N}\left(\mathbf{0}, \left(\begin{array}{cc} {\boldsymbol{K}} & {\boldsymbol{k}_{*}} \\ {\boldsymbol{k}_{*}^{T}} & {k_{* *}} \end{array}\right)\right) \end{align} $$

を用いて y' = (y_1, \ldots, y_n, y^*)^{\mathrm{T}}をサンプリングすればよいです。ここで、

$$ \begin{align} \left\{\begin{array}{l}{\boldsymbol{k}_{*}=\left(k\left(\boldsymbol{x}^{*}, \boldsymbol{x}_{1}\right), k\left(\boldsymbol{x}^{*}, \boldsymbol{x}_{2}\right), \ldots, k\left(\boldsymbol{x}^{*}, \boldsymbol{x}_{N}\right)\right)^{T}} \\ {k_{* *}=k\left(\boldsymbol{x}^{*}, \boldsymbol{x}^{*}\right)}\end{array}\right. \end{align} $$

です。条件付き多変量ガウス分布の公式より、

$$ \begin{align} \left\{\begin{array}{l}{p\left(y^{*} | \boldsymbol{x}^{*}, \mathcal{D}\right)=\mathcal{N}\left(\boldsymbol{k}_{*}^{T} \boldsymbol{K}^{-1} \boldsymbol{y}, k_{* *}-\boldsymbol{k}_{*}^{T} \boldsymbol{K}^{-1} \boldsymbol{k}_{*}\right)} \\ {\mathcal{D}=\left\{(\boldsymbol{x}_1, y_1), (\boldsymbol{x}_2, y_2), \ldots, (\boldsymbol{x}_N, y_N) \right\}} \end{array}\right. \end{align} $$

が得られます。なお、予測値(欠損値など)が複数点ある場合でも同様です。詳細は書籍p85を参照のしてください。

基本はここまでです。しかし、このままではデータの量が増えるのに伴って計算量が爆発してしまいます。これを避けるために計算を効率化する補助変数法と呼ばれる計算方法があります。

補助変数法については次回まとめます。

参考

  1. 持橋 大地、大羽 成征、(2019)「ガウス過程と機械学習」講談社
  2. Quinonero-C, ela, J. & Rasmussen, C. E. A Unifying View of Sparse Approximate Gaussian Process Regression. 6, 1939–1959 (2005).