はじめに
この記事は以下の知識を前提とします。
余談ですが時間が経つのは速いもので、上記の記事を書いてから一年以上も経つのですね…。
変分ベイズ法あらまし
事後分布とその近似分布
まず、変分ベイズ法で近似推論しようとしている事後分布を、確率変数の実現値として得られているデータ $D$ と確率変数として扱うパラメータ $\theta$ を用い、確率変数として扱わないハイパーパラメータをまとめて $\alpha$ として下記のように記します。
$$ p _ \alpha (\theta \mid D) = \frac{p _ \alpha (D, \theta )} {p _ \alpha(D)} = \frac{p _ \alpha (D \mid \theta)p _ \alpha ( \theta )} {p _ \alpha(D)} $$
変分ベイズ法では、この事後分布を解析的に求めることができない場合、代わりに
$$ q _ \eta (\theta) \simeq p _ \alpha (\theta \mid D) $$
という状況を期待し、 $\eta$ を確率変数として扱わない変分パラメータとした $\theta$ の近似事後分布を仮定するのでした。そして $\eta$ を例えば下記のELBOを最大化することで決定し、事後分布を（近似的に）求めることができたとする流れです。
$$ {\rm ELBO(\eta)} = \mathbb{E} _ {q _ {\eta}({\theta})} \left [ \log p _ {\alpha}(D\mid \theta) + \log p _ \alpha (\theta) - \log q _ {\eta}(\theta) \right] $$
この式では目的関数 $\rm ELBO$ の変数は $\eta$ です。すなわち近似事後分布 $q _ \eta$ をどのような形に設定すれば、$[]$ 内の期待値を最大化できるのだろうか？という問題になっています。
第二種最尤推定
ところで、今、$\alpha$ という近似前の元のモデルのハイパーパラメータは決め打ちという状況です。しかし、せっかくなので、このハイパーパラメータも同時に最適化してしまおうということを考え
$$ {\rm ELBO(\eta, \alpha)} = \mathbb{E} _ {q _ {\eta}({\theta})} \left [ \log p _ {\alpha}(D\mid \theta) + \log p _ \alpha (\theta) - \log q _ {\eta}(\theta) \right] $$
としてしまうことも多いです。このことは、下記の対数周辺尤度 $\log p _ \alpha (D)$ 最大化によるハイパーパラメータの決定方法（第二種最尤推定）に基づいています。
$$ \begin{align} \alpha _ {ML} &= {\rm argmax} _ \alpha \log p _ \alpha (D)\\ &= {\rm argmax} _ \alpha \log \int _ \theta p _ \alpha (D \mid \theta) p _ \alpha (\theta) \mathrm d \theta \end{align} $$
どういうことかと言うと、本来は上記の式を使って $\alpha$ を最適化していきたいのですが、積分計算がかなり難しいという事実があります。そこで
$$ \log p _ {\alpha}({D}) \ge {\rm ELBO} $$
という関係を使い、 $\log p _ \theta$ を直接最大化する代わりに$ \rm ELBO$ を最大化して事を済ましてしまおうということになります。$\rm ELBO$ には変分近似分布 $q _ \eta$ とモデル自身である $p _ \alpha(D \mid \theta) p _ \alpha (\theta)$ の両方が含まれています。
したがって、第二種最尤推定という試みと、事後分布の近似という試みを同時に一本の式の最大化問題として扱ってしまおうということができるのです（しかし元の問題に戻れば、対数周辺尤度を$\rm ELBO$よりも良く評価できる関数があるなら、そちらを目的関数にすべきであるし、事後分布を変分近似したいという観点でも、もっと良い関数があるならそうすべきである）。
ローカルパラメータとグローバルパラメータ
グローバルパラメータ
ここで、下記の端折られていた事後分布
$$ p _ \alpha (\theta \mid D) = \frac{p _ \alpha (D, \theta )} {p _ \alpha(D)} = \frac{p _ \alpha (D \mid \theta)p _ \alpha ( \theta )} {p _ \alpha(D)} $$
について少し見ていきましょう。事後分布 $\theta$ は確率変数パラメータを観測データ $D$ の条件付き分布で表したものでした。これは真ん中の式を見てわかる通り、同時分布 $p _ \alpha (D, \theta)$ というものを考え、観測データの分布 $p _ \alpha (D)$ で除したものになっています（そしてこれは、単に条件付き確率の定義である）。分母は既に観測されているデータ （＝確率変数の実現値）であるので $p _ \alpha(D)$ が計算可能かどうかは別として、分布の形状が決まった時点で定数になっています。
さて、上記の式ではパラメータ $\theta$ と観測データ $D$ の同時分布を
$$ p _ \alpha (D, \theta) = p _ \alpha (D \mid \theta)p _ \alpha (\theta) $$
と表していますが、これも同時分布と条件付き分布の間で常に成り立つ式変形をしているだけに過ぎません。いま $\theta$ に条件付けられた状況で、$D = \{ x _ 1, ..., x _ N\}$ の中の $x _ i $ が（条件付き）独立に同分布から生成されているとすると、
$$ p _ \alpha(D, \theta) = p _ \alpha(\theta) \prod _ i ^ N p _ \alpha(x _ i \mid \theta) $$
と書き表すことができます。このようにしてモデルを設計したのです。式を見るとわかる通り、$\theta$ は $\prod$ の外にいることがわかるはずです。これを観測データ $x _ i$ にとって $\theta$ が グローバルパラメータ であると表現することにしましょう。
一旦脱線：同時分布のモデリング
一度、話を脱線します。当然、後でローカルパラメータの話をするのですが、ここで同時分布をモデリングするという感覚を理解して置かなければ「グローバルパラメータ」だの「ローカルパラメータ」だのを考える意味がわからなくなってしまいます。
先程までは、データ $D = \{x _ i\} _ {i = 1, ..., N}$ という観測値に対して、確率変数パラメータ $\theta$ によってモデリングすることを考え、いわゆる尤度関数 $p _ \alpha(D \mid \theta)$ と事前分布 $p _ \alpha(\theta)$ を扱ってきました。
ここで肝になっているのは、結局のところ
$$ p _ \alpha(x, \theta) = p _ \alpha(\theta) p _ \alpha(x \mid \theta) $$
という同時分布を考え、$x _ i$ が独立同分布に生成されるとすることで
$$ p _ \alpha(D, \theta) = \prod _ i p _ \alpha(x _ i ,\theta) = p _ \alpha (\theta) \prod _ i ^ N p _ \alpha(x _ i \mid \theta) $$
なる同時分布のモデリングを実施したということです。
同時分布のモデリングができてしまえば、単に同時分布と条件付き分布の関係式から事後分布は形式的に導けます。
さて、ここではもう少し手の込んだ同時分布を考えてみましょう。パラメータ $z$ を導入し、$p _ \alpha(x, z, \theta)$。言うならばデータ $x$ の生成仮定にはまた別のパラメータが一枚噛んでいることをモデルに組み込むということです。その噛ませ方を下記のように表現してみます。
$$ p _ \alpha(x, z, \theta) = p _ \alpha(\theta)p _ \alpha(z \mid \theta)p _ \alpha(x \mid z) $$
これは、パラメータ $\theta$ に条件付けられて $x$ が決まるというこれまでのモデルとは異なり、一旦 $z$ を挟んでから $x$ が最終的に観測されているとモデリングしていることになります。ここで、$z$ なるものは特に潜在変数と呼ばれることが多いようです。（$\theta$ にしてもそうであるが）$z$ というのは観測データ $x$ を生成しているであろう物理的実体と紐付いている必要は必ずしもありません。統計モデリングを実施している時点で、多くの場合は物理的な実体であるとか、世の中の真実とは異なってはいるが「便利なモデル」を作ろうとしているということに注意が必要です。
ローカルパラメータ
話を戻しましょう。このように同時分布の設計を自由に実施し（もちろんいい加減に…という意味ではない）、上記で見たように
$$ p _ \alpha(x, z, \theta) = p _ \alpha(\theta)p _ \alpha(z \mid \theta)p _ \alpha(x \mid z) $$
というモデルを考えることになったとしましょう。ここで観測データが $D = \{ x _ 1, ..., x _ N \}$ と得られており、この観測を同時分布に当てはめるならどうすればいいでしょうか。例えば、下記のようにできます。
$$ p _ \alpha(D, z, \theta) = p _ \alpha(\theta)p _ \alpha(z \mid \theta)\prod _ i p _ \alpha(x _ i \mid z) $$
これは「グローバルパラメータ $\theta$ に条件付けられたグローバルパラメータ $z$ に対して条件付き独立で $x$ が生成する」というようなモデルになっています。グラフィカルモデルでは下記のように表現します。
白丸は観測されていないことを表し、灰色が観測されている事を表します。黄色の四角で右下に $N$ とあるのは、N個の変数が独立に存在することを意味します。実を言うとこれはつまらないモデリングです。説明したかったのはこれではありません。
代わりに「グローバルパラメータ $\theta$ に対して条件付き独立で生成される $z _ i$ に条件付けられた $x$ が生成されている」というようなモデルを見ましょう。グラフィカルモデルは
となっており、そのような同時分布モデルの式は下記になります。
$$ p _ \alpha(D, z _ 1, ..., z _ N, \theta) = p _ \alpha(\theta)\prod _ i p _ \alpha(z _ i \mid \theta)p _ \alpha(x _ i \mid z _ i) $$
違いは、どの段階で条件付き独立になっているかです（グラフィカルモデルを見るとわかりやすい）。前者では $x$ が生成される段階になって初めて条件付き独立を導入しています。後者は $\theta$ によって $z$ が生成される段階で条件付き独立が導入され、その後の $z$ と $x$ の関係はデータ点毎に関係性を持っている状況となっています。そして、同時分布は $N$ 個の $x _ i$ を考えているのと同様で $N$ 個の $z _ i$ を考慮しなければならなくなりました。
今ここに、$\theta$ はグローバルパラメータで $z$ はローカルパラメータとなっているモデルを導入したことになります。
償却推論
ローカルパラメータの事後分布
ここまでで、ローカルパラメータを含む同時分布のモデルを把握できました。さて、ベイズ推論では観測データからパラメータの事後分布求めようと試みますが、今事後分布を書き下すことにすれば
$$ p _ \alpha (z _ 1, ..., z _ N, \theta \mid D) = \frac{ p _ \alpha(\theta)\prod _ i p _ \alpha(z _ i \mid \theta)p _ \alpha(x _ i \mid z _ i) }{p _ \alpha(D)} $$
となっています。当然、この事後分布も解析的に求めるのは無理そうという前提で、変分推論を用いて行くことにしましょう。 準備するのは事後分布を近似できそうな変分近似分布 $q _ \eta (z _ 1, ..., z _ N, \theta)$ です。我々は、同時分布を設計した段階で、 $z _ i$ は互いに $\theta$ の条件付き独立であることを組み込んでいるのですから、変分近似分布でも条件付き独立を入れてしまって構わなく、
$$ \begin{align} q _ \eta (z _ 1, ..., z _ N, \theta) &= q _ \eta (z _ 1 \mid \theta)...q _ \eta (z _ N \mid \theta)q _ \eta(\theta)\\ & = q _ \eta(\theta)\prod _ i q _ \eta(z _ i \mid \theta ) \end{align} $$
という形式で変分近似分布を準備することにしてみましょう。さて、こうすれば、$\eta$ を使って上手なパラメトリックな分布を準備し、変分推論を実行することができるでしょうか。素朴には、まとめて表記していた変分パラメータ $\eta$ を個別に
$$ q _ {\eta _ \theta} (\theta)\prod _ i q _ {\eta _ {z _ i}}(z _ i \mid \theta ) $$
と設定してしまうことができます（すなわち、$z _ i$ のインデックス毎に変分パラメータを個別に準備するということである）。
実を言うと、この推論はかなり難しい問題になってしまいます。なぜなら $z _ i$ が $N$ 個もあるせいで変分パラメータもその数だけ準備してしまうことになり、かつ $\eta _ {z _ i}$ を知る手がかりは対応するインデックスを持つ $x _ i$ だけということになってしまうからです。これはあたかも事後分布の推論において $z _ i$ の振る舞いを記述するには $x _ i$ しか頼りにできないということです。
変分パラメータを予測するモデルの導入
$z _ i$ を推論する上で $x _ i$ しか使えないというのは妥当なのでしょうか。もともとは $z$ に条件付けられて $x$ が生成されるというモデルだったのです。そして今事後分布を推論するためにデータ $\{x _ i \} _ {i=1,...,N} $ を使って行こうという状況です。$x, z$ に関係性があるのだとすれば、$x, z$ の関係性を表現できるような仕組みがあれば良いはずです。
例えば $\eta _ {z _ i}$ という変分パラメータと $x _ i$ との関係性を $\eta _ {z _ i} = f(x _ i) $ などと結べたらどうでしょうか。そして $f$ がうまく作れた暁には $\eta = f(x)$ と対応する$x$ に対する変分パラメータ $\eta$ 自体を予測できるとしたらどうでしょうか。
今、改めて変分近似分布を表記するために $f (x)$ という関数を導入することにして
$$ q _ \eta (z _ 1, ..., z _ N, \theta) = q _ \eta (\theta) \prod _ i q _ {f(x _ i)}(z _ i \mid \theta ) $$
としてしまいます。すると求めるべきはグローバルパラメータ $\theta$ の近似分布を記述する変分パラメータ $\eta$ と、$z$ の方に関しては変分パラメータ $\eta _ {z _ i}$ を $N$ 個考える代わりに、$x$ と変分パラメータ $\eta _ z$ 間の関係性を記述しているであろう関数 $\eta = f (x)$ になります。当然のことながら、関数 $f$ として何かを準備せねばならず、これは $\eta _ {z _ i}$ をたくさん考えなければならないという問題を関数 $f(x)$ に押し付けたに過ぎません。
しかし、この押しつけが推論を効率的なものにしてくれます。
なぜなら、$\eta _ {z _ i}$ を求めるためには対応するインデックスである $x _ i $ しか使えないような状況になっていた中で、その状況を回避し、$\eta = f(x)$ という $x$ の値に依存して変分パラメータ $\eta$ は変わることを反映できているからです。
これは要するにインデックス $i$ に基づいて個々にパラメータ $\eta _ i$ を作るのではなく、$i$ 毎に異なる $\eta _ i$ を$x$ のデータをすべて使って予測してしまおうということです。
このように、ローカルパラメータの近似事後分布を推論するために、変分パラメータ自体データから回帰してしまおうという方法を償却推論といいます。
結果として、償却推論を用いた変分近似はデータから変分パラメータを出力する関数 $f$ を最適化するという問題に落ち着いたのです。
