はじめに
この記事は『Multi-Class Gaussian Process Classification Made Conjugate: Efficient Inference via Data Augmentation』を理解するために勉強しながら書いたものである。本記事は和訳ではないので必要な部分のみをピックアップしながら数式の行間を埋めていく形になるが、原論文は背景もしっかりまとまっていて数理的にもエレガントなのでぜひ原論文も読んでほしい。
本記事の数式は原論文に倣うが、一部変更点としてベクトルを表すために太字は使わないことに注意する（こだわりがあるわけではなく単に面倒くさい）。数式番号は論文のものをそのまま用いる。
また、この論文では data augmentation という概念が出てくるが、これはモデルに潜在変数を追加する操作のことを指しており、深層学習の文脈で画像を切り抜き・拡大縮小・回転などしてデータセットの画像をカサ増しさせる data augmentation とは異なると思われることに注意する。
ざっくり背景の補足
ガウス過程は回帰に用いる場合は予測分布が closed form で書けるが、分類に用いる場合は予測分布を closed form で書けないので
- 変分推論（Variational Inference）
- 信念伝播法（EP: Expectation Propagation method）
- ラプラス近似（Laplace approximation）
などで解く必要があるとされていた。その後、マルコフ連鎖モンテカルロ法（MCMC: Markov Chain Monte Carlo methods）を使った解法なども提案されたが、ガウス過程回帰の簡潔さに比べると、個人的にはあまり勉強する気にもなれないくらい複雑なモデルと解法になるという時代が 20 年ほど続いていた。
『Multi-Class Gaussian Process Classification Made Conjugate: Efficient Inference via Data Augmentation』はそのガウス過程分類をついに共役（conjugate）にしたぞ、という論文である。正確には平均場近似による変分推論の文脈で事前分布と事後分布を条件付き共役（conditionally conjugate）にした、という内容になっていて、簡単に言えばパラメータのアップデートは closed form の式での交互更新で行えるようになったということである。
従来のアルゴリズムに比べても高速で、精度も過去のガウス過程分類の SOTA に匹敵するらしいので、ようやくガウス過程分類に光が差してきた感じがする。
1. Introduction
省略。
与えられたN点のデータセットの入力をX = (x _ 1, \ldots, x _ N)、ラベルをy = (y _ 1, \ldots, y _ N)とおく。ただしy _ i \in \{ 1, \ldots, C \}でCは分類したいクラスの個数である。
マルチクラスガウス過程分類ではガウス過程に従う潜在変数
f = (f ^ 1, \ldots, f ^ C) \quad \text{where} \,\, f ^ c \sim \operatorname{GP}(0, k ^ c)
を導入する。ただしk ^ cはそのガウス過程に対応するカーネル関数である。このときラベルのカテゴリカル尤度は
p(y _ i = k | x _ i, f _ i) = g ^ k \bigl( f(x _ i) \bigr)
とモデル化される。ただしg ^ k (f)は実数ベクトルを確率ベクトルに変換する関数である。f(x _ i) \in \mathbb{R} ^ Cなので、おそらくこれをg ^ k \bigl( f(x _ i) \bigr) \in \mathbb{P}、ただし
\mathbb{P} := \{ p \in [0, 1] ^ C \mid {\textstyle \sum _ {c = 1} ^ C p ^ c = 1 } \}
に変換したいということだと思う。
もっとも一般的なカテゴリカル尤度の構成方法は
p(y _ i = k | f _ i) = \frac{\exp (f _ i ^ k)}{{\textstyle \sum _ {c=1} ^ C \exp (f _ i ^ c)}}
である。ただしf _ i ^ c = f ^ c (x _ i)であり、わかりやすさのためにx _ iは条件から省いてある。
しかしラベルのカテゴリカル尤度についてこのモデル化を採用するとモデルは共役ではなくなる。ガウス過程分類の歴史はこのカテゴリカル尤度をどう工夫するかという歴史であり、原論文には過去の論文で提案された他のいくつかのモデル化方法が載っている。
3. Conjugate multi-class Gaussian process classification
提案手法ではカテゴリカル分布のモデル化をロジスティックソフトマックス（logistic-softmax）
p(y _ i = k | f _ i) = \frac{\sigma (f _ i ^ k)}{{\textstyle \sum _ {c=1} ^ C \sigma (f _ i ^ c)}} \tag{5}
によって行う。ただし\sigmaはロジスティック関数で\sigma(z)=(1 + \exp(-z)) ^ {-1}である。あとは
- 補助変数\lambda _ iを用いて分母にあるガウス過程潜在変数f _ i ^ kをばらばらにする
- ポアソン確率変数を導入することでモデル尤度を簡単にする
- モデルの共役表現を得るためにシグモイド関数の Pólya–Gamma 表現を用いる
の三段階の data augmentation により条件付き共役モデルを得る。
Augmentation 1: Gamma augmentation
積分等式
\frac{1}{z} = \int _ 0 ^ \infty \exp(- \lambda z) d \lambda
を(5)式に用いると
\begin{aligned}
p(y _ i = k | f _ i) &= \frac{\sigma (f _ i ^ k)}{{\textstyle \sum _ {c=1} ^ C \sigma (f _ i ^ c)}} \\
&= \sigma (f _ i ^ k) \int _ 0 ^ \infty \exp \Bigl(- \lambda _ i \sum _ {c=1} ^ C \sigma (f _ i ^ c) \Bigr) d \lambda _ i
\end{aligned}
となる。これはギブスサンプリング界隈ではよく知られたテクニックらしいが、変分推論ではあまり使わないらしい。ここで\lambda _ iを新たな潜在変数とみなすと、拡張された尤度は
p(y _ i = k | f _ i, \lambda _ i) = \sigma (f _ i ^ k) \prod _ {c=1} ^ C \exp \Bigl(- \lambda _ i \sigma (f _ i ^ c) \Bigr) \tag{6}
と求まり、同時に不適切事前分布p(\lambda _ i) \propto \mathbb{1} _ {[0,\infty)}(\lambda _ i)をモデルに導入することになる。
行間
ここは若干の行間があるので埋めていく。まず不適切事前分布（improper prior）とは領域の全体\Omegaで積分したときに値が1にならない確率分布、つまり
\int _ \Omega p(x) dx \neq 1
となる密度関数p(x)をいう。まぁ確率分布の定数倍とか、どうせ計算中は気にしないで最後につじつま合わせするのでどうでもいいんじゃない、という場合には使っても問題ないのではなかろうか。平均場近似とかまさにそんな感じなので。
p(\lambda _ i) \propto \mathbb{1} _ {[0,\infty)}(\lambda _ i)という表記はあまり見かけないが、文脈から察するに区間[0, \infty)で定義された定数関数（つまり一様分布）っぽい。あとの議論でわかることだが、ここではおそらく
という定数関数を想定していると思われる（この定数関数を区間[0, \infty)で積分すると積分した値は+\inftyに発散するけどとりあえず目をつぶろう）。
あとは確率分布の周辺化と条件付き確率分布の定義より
\begin{aligned}
p(y _ i = k | f _ i) &= \int p(y _ i = k, \lambda _ i | f _ i) d \lambda _ i \\
&= \int p(y _ i = k | f _ i, \lambda _ i)p(\lambda _ i) d \lambda _ i \\
&= \int p(y _ i = k | f _ i, \lambda _ i) d \lambda _ i \quad (\,\because \, p(\lambda _ i) = 1) \\
\end{aligned}
が成り立つ。ここでp(y _ i = k | f _ i, \lambda _ i)が(6)式で表されるとすると、
\begin{aligned}
\int p(y _ i = k | f _ i, \lambda _ i) d \lambda _ i &= \int _ 0 ^ \infty \sigma (f _ i ^ k) \prod _ {c=1} ^ C \exp \Bigl(- \lambda _ i \sigma (f _ i ^ c) \Bigr) d \lambda _ i \\
&= \int _ 0 ^ \infty \sigma (f _ i ^ k) \exp \Bigl(- \lambda _ i \sum _ {c=1} ^ C \sigma (f _ i ^ c) \Bigr) d \lambda _ i \\
&= \sigma (f _ i ^ k) \int _ 0 ^ \infty \exp \Bigl(- \lambda _ i \sum _ {c=1} ^ C \sigma (f _ i ^ c) \Bigr) d \lambda _ i
\end{aligned}
となって、これは(5)式におけるp(y _ i = k | f _ i)の定義に一致するから、p(y _ i = k | f _ i, \lambda _ i)には(6)式を割り当ててもよいことがわかる。
Augmentation 2: Poisson augmentation
ポアソン分布\operatorname{Po}(\cdot | \lambda)のモーメント母関数は
\exp(\lambda(z - 1)) = \sum _ {n = 0} ^ \infty z ^ n \operatorname{Po}(z | \lambda)
である。\sigma(f) = 1 - \sigma(-f)という事実を踏まえてz = \sigma(-f)を代入すると、
\begin{aligned}
\exp (- \lambda _ i \sigma (f _ i ^ c)) &= \exp (\lambda _ i(\sigma(- f _ i ^ c) - 1)) \\
&= \sum _ {n _ i ^ c = 0} ^ \infty (\sigma (- f _ i ^ c)) ^ {n _ i ^ c} \operatorname{Po}(n _ i ^ c | \lambda _ i)
\end{aligned}
となるので、拡張された尤度は
p(y _ i = k | f _ i, \lambda _ i, n _ i) = \sigma (f _ i ^ k) \prod _ {c=1} ^ C (\sigma(- f _ i ^ c)) ^ {n _ i ^ c} \tag{7}
とし、n _ i = (n _ i ^ 1, \ldots, n _ i ^ C)の事前分布はポアソン分布p(n _ i ^ c | \lambda _ i) = \operatorname{Po}(n _ i ^ c | \lambda _ i)とする。
行間
この拡張の行間の埋め方はさっきと同じである。周辺化と条件付き確率の定義から
\begin{aligned}
p(y _ i = k | f _ i, \lambda _ i) &= \int p(y _ i = k, n _ i | f _ i, \lambda _ i) d n _ i \\
&= \int p(y _ i = k | f _ i, \lambda _ i, n _ i) p(n _ i | \lambda _ i) d n _ i
\end{aligned}
であり、(7)式とn _ iに関する事前分布から
\begin{aligned}
\int p(y _ i = k | f _ i, \lambda _ i, n _ i) p(n _ i | \lambda _ i) d n _ i &= \sum _ {n _ i ^ 1 = 0} ^ \infty\cdots \sum _ {n _ i ^ C = 0} ^ \infty p(y _ i = k | f _ i, \lambda _ i, n _ i) p(n _ i | \lambda _ i) \\
&= \sum _ {n _ i ^ 1 = 0} ^ \infty\cdots \sum _ {n _ i ^ C = 0} ^ \infty \Bigl( \sigma (f _ i ^ k) \prod _ {c=1} ^ C (\sigma(- f _ i ^ c)) ^ {n _ i ^ c} \prod _ {c = 1} ^ C \operatorname{Po}(n _ i ^ c | \lambda _ i) \Bigr) \\
&= \sum _ {n _ i ^ 1 = 0} ^ \infty\cdots \sum _ {n _ i ^ C = 0} ^ \infty \Bigl( \sigma (f _ i ^ k) \prod _ {c=1} ^ C (\sigma(- f _ i ^ c)) ^ {n _ i ^ c} \operatorname{Po}(n _ i ^ c | \lambda _ i) \Bigr) \\
&= \sigma (f _ i ^ k) \prod _ {c=1} ^ C \Bigl( \sum _ {n _ i ^ c = 0} ^ \infty (\sigma(- f _ i ^ c)) ^ {n _ i ^ c} \operatorname{Po}(n _ i ^ c | \lambda _ i) \Bigr) \\
&= \sigma (f _ i ^ k) \prod _ {c=1} ^ C \exp (- \lambda _ i \sigma (f _ i ^ c))
\end{aligned}
となって、これは(6)式におけるp(y _ i = k | f _ i, \lambda _ i)の定義に一致するから、p(y _ i = k | f _ i, \lambda _ i, n _ i)は(7)式で定義してもよいことがわかる。
Augmentation 3: Pólya-Gamma augmentation
最後に、シグモイド関数の Pólya-Gamma 表現によってスケール混合ガウシアン表現に変換する。
\sigma(z) ^ n = \int _ 0 ^ \infty 2 ^ {-n} \exp \Bigl( \frac{nz}{2} - \frac{z ^ 2}{2} \omega \Bigr) \operatorname{PG}(\omega | n, 0) d \omega \tag{8}
上式は原論文だと誤植でd \omegaが抜けているが、上式の根拠は参考文献『Bayesian inference for logistic models using Pólya-Gamma latent variables』の(7)式、
\frac{(e ^ \psi) ^ a}{(1 + e ^ \psi) ^ b} = 2 ^ {-b} e ^ {\kappa \psi} \int _ 0 ^ \infty e ^ {- \omega \psi ^ 2 / 2} \operatorname{PG}(\omega | b, 0) d \omega
であろうから抜けているのはd \omegaで間違いないであろう。上式はa = 0, b=-n, \psi =zで置き換えると\kappa = \frac{n}{2}から(8)式に一致する。
Pólya-Gamma 分布\operatorname{PG}(\omega | n, b)はそのモーメントを解析的に求めることができ、効率的なサンプリングも可能なので augmentation に適している。このデータ拡張を(7)式に適用すると
p(y _ i = k | f _ i, \lambda _ i, n _ i, \omega _ i) = \prod _ {c = 1} ^ C 2 ^ {-(y _ i ^ {\prime c} + n _ i ^ c)} \exp \left( \frac{(y _ i ^ {\prime c} - n _ i ^ c)f _ i ^ c}{2} - \frac{(f _ i ^ c) ^ 2}{2} \omega _ i ^ c \right) \tag{9}
とできる。ただし\omega _ i = (\omega _ i ^ 1, \ldots, \omega _ i ^ C)は Pólya-Gamma 分布
p(\omega _ i | n _ i, y _ i) = \prod _ {c = 1} ^ C \operatorname{PG}(\omega _ i ^ c | y _ i ^ {\prime c} + n _ i ^ c, 0)
にしたがうものとする。ここでy ^ \primeはラベルについてのN \times C次元の one-hot エンコーディングである。
行間
(9)式を導く行間の埋め方はここまでの 2 回と同じ手順に加えて Appendix A.1 に書かれた操作が必要になる。Appendix A.1 を見ると実際はさらに\tilde{\omega} _ iという変数をモデルに導入していることがわかる。
以上で、(9)式がf _ iに関するガウシアン形式になったのでガウス過程潜在モデルが共役になった。モデル全体としても条件付き共役になったということを確認していく。
The final model
以上のデータ拡張をもってモデルの条件付き分布は closed form で与えられる。
ガウス過程のf ^ cは
p(f ^ c | y, \omega ^ c, n ^ c) = \mathcal{N} \Bigl( f ^ c | \frac{1}{2} A ^ c (y ^ {\prime c} - n ^ c), A ^ c \Bigr)
で与えられる。ただし共分散行列A ^ cは
A ^ c = \bigl( \operatorname{diag}(\omega ^ c) + K _ c \bigr) ^ {-1}
であり、K _ cはガウス過程f ^ cのカーネル行列である。さらに\lambda _ iの条件付き分布は
p(\lambda _ i | n _ i) = \operatorname{Ga} \Bigl( \lambda _ i | 1 + \sum _ {c = 1} ^ C n _ i ^ c, C \Bigr)
と表される。ただし\operatorname{Ga}(\cdot | a, b)は形状パラメータaと比率パラメータbを持つガンマ分布である。\lambda _ iの事前分布は不適切事前分布になっているが、計算上問題にならないのでよしとする。また、n, \omegaに関しても条件付き分布は
\begin{aligned}
& p(n _ i ^ c | f _ i ^ c, \lambda _ i) = \operatorname{Po}(n _ i ^ c | \lambda _ i \sigma (f _ i ^ c)) \\
& p(\omega _ i | n _ i ^ c, f _ i ^ c, y _ i) = \operatorname{PG}(\omega _ i | y _ i ^ {\prime c} + n _ i ^ c, |f _ i ^ c|)
\end{aligned}
と求められる。
行間
ベイズモデルを作るときの基本は同時分布を設計することである。今回の場合であればそれ以外の分布は同時分布からベイズの定理や適切な周辺化によって解析的に closed form で求まるようにしてくれているはずだから、同時分布
p(y, f, \lambda, n, \omega, \tilde{\omega})
が分かればあとの分布は自然に求められるはずである（Appendix A.1 の適切な式変形で\tilde{\omega}は\omegaに吸収されるが一応書いておく）。
同時分布を書いてくれていれば誤解なく計算を進めることができたのだが、残念ながら書いていないようなので頑張るしかない。さて、こういったモデルでは各データに対して i.i.d. を仮定することが多いのでまず
p(y, f, \lambda, n, \omega, \tilde{\omega}) = \prod _ {i = 1} ^ N p(y _ i, f _ i, \lambda _ i, n _ i, \omega _ i, \tilde{\omega} _ i)
だろう。さらに、ここまでの議論から
\begin{aligned}
&p(y _ i, f _ i, \tilde{\omega} _ i, \omega _ i, n _ i, \lambda _ i) \\
= &p(y _ i | f _ i, \tilde{\omega} _ i, \omega _ i, n _ i, \lambda _ i) p(f _ i) p(\tilde{\omega} _ i) p(\omega _ i | n _ i)p(n _ i | \lambda _ i) p(\lambda _ i)
\end{aligned}
とする（式を見やすくするために変数は並べ替えた）。ここで、論文には書いていないがおそらく
p(y _ i | f _ i, \tilde{\omega} _ i, \omega _ i, n _ i, \lambda _ i) = \prod _ {k = 1} ^ C p(y _ i = k | f _ i, \tilde{\omega} _ i, \omega _ i, n _ i, \lambda _ i) ^ {y _ i ^ {\prime k}}
である。以降は論文にしたがって
\begin{aligned}
&p(y _ i = k | f _ i, \tilde{\omega} _ i, \omega _ i, n _ i, \lambda _ i) = 1 \\
&p(f _ i) = \prod _ {c = 1} ^ C \operatorname{GP}(f _ i ^ c | 0, k ^ c) \\
&p(\tilde{\omega} _ i) = \operatorname{PG}(\tilde{\omega} _ i | 1, 0) \\
&p(\omega _ i | n _ i) = \prod _ {c=1} ^ C \operatorname{PG}( \omega _ i ^ c | n _ i ^ c, 0) \\
&p(n _ i | \lambda _ i) = \prod _ {c = 1} ^ C \operatorname{Po}(n _ i ^ c | \lambda _ i ^ c) \\
&p(\lambda _ i) = \text{const.}
\end{aligned}
である。あとはベイズの定理と周辺化で頑張る。たとえば
\begin{aligned}
p(f | y, \tilde{\omega}, \omega, n) &= \int d \lambda \cdot p(f, \lambda | y, \tilde{\omega}, \omega, n) \\
&= \int d \lambda \cdot \frac{p(y, f, \tilde{\omega}, \omega, n, \lambda)}{p(y, \tilde{\omega}, \omega, n)} \\
&= \int d \lambda \cdot p(y, f, \tilde{\omega}, \omega, n, \lambda) \cdot \text{const.} \\
& \propto \int d \lambda \cdot p(y, f, \tilde{\omega}, \omega, n, \lambda) \\
&= \int d \lambda \cdot p(y | f, \tilde{\omega}, \omega, n, \lambda)p(f) p(\tilde{\omega}) p(\omega | n)p(n | \lambda) p(\lambda) \\
&= p(f) p(\tilde{\omega}) p(\omega | n) \int p(y | f, \tilde{\omega}, \omega, n, \lambda) p(n | \lambda) p(\lambda) d \lambda \\
&= p(f) p(\tilde{\omega}) p(\omega | n) \int p(y, n, \lambda | f, \tilde{\omega}, \omega) d \lambda \\
&= p(f) p(\tilde{\omega}) p(\omega | n) p(y, n | f, \tilde{\omega}, \omega)
\end{aligned}
となり、論文中でp(f ^ c | y, \omega ^ c, n ^ c)を求めている理由がわかる（f ^ cは\lambdaに依存しないので周辺化して積分消去している）。実際に積分計算して求めるのは、やってもいいがどうせ大変だろうし今回は無料記事なのでやりたくない。論文を信じることにする。あとは同様にして
p(\tilde{\omega}, \omega | y, f, n) = \int d \lambda \cdot p(\tilde{\omega}, \omega, \lambda | y, f, n)
p(n | f, \lambda) = \int dy \int d\tilde{\omega} \int d\omega \cdot p(n, y, \tilde{\omega}, \omega | f, \lambda)
p(\lambda | n) = \int dy \int df \int d\tilde{\omega} \int d\omega \cdot p(\lambda, y, f, \tilde{\omega}, \omega | n)
からp(\tilde{\omega}, \omega _ i | n _ i ^ c, f _ i ^ c, y _ i)、p(n _ i ^ c | f _ i ^ c, \lambda _ i)、p(\lambda _ i | n _ i)を求めればよいことがわかる。
4. Inference
『Variational Learning of Inducing Variables in Sparse Gaussian Processes』のスパースガウス過程回帰を使って計算量を削減している。
疲れたのでいったん休ませてほしい。続きは気力があったらあとで書きます。
まとめ
修士のときに私も平均場近似で欠損値がある行列の特異値分解モデルを共役にしたのだが、そのときだってA4紙何枚使ったか覚えてないくらいモデルをいじって試行錯誤した（一度の計算で10〜20枚は使うので確実に100枚以上は使った）。今回の論文で使われているテクニックは高度なものばかりだし、ガウス過程分類を共役にするためにどれだけの時間と労力を注ぎ込んだのだろうと著者たちの苦労が偲ばれる。ガウス過程分類が登場してから 20 年近く誰も思いつかなかったモデルを発見したわけだから、このエレガントな論文の裏にはエゲツない執念があったに違いない。著者たちに敬意と賞賛を送りたい。
Discussion