</link>

図 1 の左上に示した Thomas Bayes (1701-1761) は，その名の通り，死後 3 年経った1764 年に発表した論文で初めてベイズの定理を発見した。 ただし，ベイズはその定理で，一様事前分布を用いている [1]。 図 1 の右下に示した Pierre-Simon Laplace (1749-1827) は，Bayes の研究 を知らなかったようだ。 だが 25 歳の時に書いた手記でより一般的な形で 同定理を発見し，その広い応用可能性を示した [2]。 これらの問題について S.M. Stiegler はこう書いている:

この手記が与えた影響は計り知れない。 Bayes 自身の論文は 1780 年まで無視され，20 世紀まで科学的議論において重要な役割を果たさなかったので，「Bayes 的」な考え方はここから数学界に初めて広まったのである。 また，現在も採用されている事後分布の漸近的解析の数学的手法を導入したのも，この Laplace の論文であった。 また，最適推定の最古の例として，事後期待損失を最小化する推定量の導出とその特徴づけが行われたのも，この論文である。 2 世紀以上の時を経て，我々数学者，統計学者は，この科学の傑作に自分たちのルーツを認めるだけでなく，そこから学ぶことができるのである [3]。

最尤法 (ML) は，現代の統計的信号処理で最もよく使われる手法の 1 つである。 EM アルゴリズムは，ML 推定のための反復纂法であり，多くの利点を持ち，統計的信号処理の問題を解くための標準的な方法論となっている。 しかし EM アルゴリズムにはある種の要件があり，複雑な問題への適用を著しく制限している。 最近 EM アルゴリズムの制約条件を緩和した変分ベイズ推論と呼ばれる新しい手法が登場し，急速に普及しつつある。 さらに EM アルゴリズムがこの方法論の特殊な場合と見なすことができることを示すことができる。 本稿では，まず，信号処理のコミュニティを対象としたベイズ変分推論のチュートリアルの紹介を行う。 線形回帰と Gauss 混合モデリングを例として，EM アルゴリズムと比較して，Bayes 変分推論が提供する追加的な機能を示す。

1. はじめに

ML の方法論は，現代の統計的信号処理の基本的な定番の 1 つである。 EM アルゴリズムは，ML 推定値を得るために多くの利点を提供する反復アルゴリズムである。 1977 年に Dempster らによって正式に発表されて以来[4]，EM アルゴリズムは ML 推定の標準的な方法論となった。 IEEE のコミュニティでも，EM は着実に人気を集めており，応用例も増えてきている。 IEEE のジャーナルに EM アルゴリズムが初めて掲載されたのは 1988 年であり，光子制限画像の断層像再構成の問題を扱っていた[5, 6]。 以来 EM アルゴリズムは，画像や映像の復元やセグメンテーション，画像のモデリング，通信や音声認識における搬送波周波数同期やチャネル推定など，幅広い用途で用いられる統計的信号処理のツールとして人気を博してきた。

EM アルゴリズムの背後にある概念は非常に直感的で自然なものである。 EM 的なアルゴリズムは[4] よりも前から統計学の文献に存在していたが，そのようなアルゴリズムは実際には特殊な文脈における EM アルゴリズムであった。 このような最初の文献は 1886 年に遡り，Newcomb が 2 つの一変量規範の混合物のパラメータの推定を考察している[7]。 しかし，そのようなアイデアが統合され，EM アルゴリズムの一般的な定式化が確立されたのは [4] であった。 [4] 以前の EM アルゴリズムの歴史に関する良いサーベイが [8] にある。

本論文は EM アルゴリズムに関するチュートリアルではない。 このようなチュートリアルは 1996 年に IEEE Signal Processing Magazine に掲載された[9]。 本論文の目的は EM アルゴリズムの欠点を改善する統計的推測のための新しい方法論を提示することである。 この方法は変分近似と呼ばれ[10]，EMアルゴリズムが適用できない複雑なベイズモデルを解くために用いることができる。 変分近似に基づく Bayes 推定は，それが最初に導入された 1990 年代半ば以来，機械学習コミュニティによって広く利用されてきた。

• [4] A. Dempster, N. Laird, and D. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,”J. Roy. Statis. Soc. A, vol. 39, no. 1, pp. 1–38, 1977.
• [5] W. Jones, L. Byars, and M. Casey, “Positron emission tomographic images and expectation maximization: A VLSI architecture for multiple iterations per second,” IEEE Trans. Nuclear Sci., vol. 35, no. 1, pp. 620–624, 1988.
• [6] Z. Liang and H. Hart, “Bayesian reconstruction in emission computerized tomography,” IEEE Trans. Nuclear Sci., vol. 35, no. 1, pp. 788–792, 1988.
• [7] S. Newcomb, “A generalized theory of the combination of observations so as to obtain the best result,” Amer. J. Math., vol. 8, pp. 343–366, 1886.
• [8] G. McLachlan and T. Krishnan, The EM Algorithm and Extensions. New York: Wiley, 1997.
• [9] T.K. Moon, “The EM algorithm in signal processing,” IEEE Signal Process. Mag., vol. 13, no. 6, pp. 47–60, 1996.
• [10] V. Smídl and A. Quinn, The Variational Bayes Method in Signal Processing. New York: Springer-Verlag, 2005.

2. ベイズ推論の基本

$x$ を観測値，$\theta$ を $x$ を生成したモデルの未知のパラメータと仮定する。 本稿では，パラメータを指す場合は推定，確率変数を指す場合は推論という用語を厳密に使用する ことにする。

推定という用語は，不完全で不確かで雑音の多いデータから，パラメータの値を計算で近似することを意味する。

対照的に，推論という用語はベイズ推論を意味するものとして使われ，観察値 $x$ を与えられた確率変数 $\theta$ の事後確率 $p(\theta\vert x)$ を推論するために，事前証拠と観測が使われる過程を意味する。 パラメータ推定のための最も一般的なアプローチの 1 つは ML である。 このアプローチによると，ML 推定値は次のように得られる: \(\hat{\theta}_{ML}=\arg\max_{\theta} p(x;\theta)\tag{tz1}\)

ここで，$p(x;\theta)$ は，観測値 $x$ を生成した仮定モデルに基づく観測値とパラメータとの間の確率的関係を記述する。 このとき，$p(x;\theta)$ と $p(x\vert\theta)$ の表記の違いを明確にしたい。 $p(x;\theta)$ と書くと $\theta$ がパラメータであることを意味し，$\theta$ の関数として尤度関数と呼ばれる。 これに対して $p(x\vert\theta)$ と書くと，$\theta$ が確率変数であることを意味する。

尤度関数 $p(x;\theta)$ の直接評価は複雑で，直接計算することも最適化することも難しいか不可能な場合が多く，興味深い。 このような場合，隠れ変数 $z$ を導入することで，この尤度の計算が非常に容易になる。 これらの確率変数は，Bayes 則によって観測値を未知のパラメータに接続するリンクとして機能する。 隠れ変数の選択は問題に依存する。 しかし，その名前が示すように，これらの変数は観測されず，条件付き確率 $p(x|z)$ が計算しやすいように，観測に関する十分な情報を提供する。 この役割とは別に，隠れ変数は統計的モデリングにおいて別の役割を果たす。 それは，観測値を生成したと仮定される確率的機構の重要な部分であり，「グラフィカルモデル」と呼ばれるグラフによって簡潔に記述することができる。 グラフィカルモデルの詳細については 「グラフィカル・モデル」の項で説明する。

隠れ変数とその事前確率 $p(z;\theta)$ が導入されると，以下のように隠れ変数を積分 (周辺化) することによって，尤度または周辺尤度 (時折呼ばれる) を得ることができる。

\[\tag{tz2} p(x;\theta) = \int p(x;z,\theta)\,dz = \int p(x\vert z,\theta) p(z;\theta)\,dz.\]

この一見単純な積分が Bayes 的手法の肝であり，この方法で尤度関数とベイズの定理による隠れ変数の事後分布を得ることができるからである。

\[p(z\vert x,\theta)=\frac{p(x\vert z,\theta)\, p(z;\theta)}{p(x;\theta)} \tag{tz3}\]

事後情報が得られれば，隠れ変数について上述のごとき推論も可能である。 上記の定式化は簡単であるが，ほとんどの場合 式 (2) の積分は閉形式で計算することが不可能か非常に困難である。 したがって，Bayes 推定では，この積分を迂回または近似的に評価する技術に主な取り組みが集中している。

このような方法は，大きく 2 つに分類される。 1 つはモンテカルロ法とも呼ばれる数値サンプリング法，もう 1 つは決定論的近似法である。 この記事では，モンテカルロ法については一切触れない。 このような手法に興味のある読者は，例えば [11] や [12] のようなこのトピックに関する多くの書籍やサーベイ記事を参照されたい。 さらに ML 法を拡張した最大事後推定 (MAP) は，非常に粗いベイズ近似と考えることができる。 Maximum A Posteriori: Poor Man's Bayesian Inference 参照。

以下に示すように，EM アルゴリズムは，事後確率 $p(z\vert x;\theta)$ の知識を前提とし，明示的に計算することなく尤度関数を反復的に最大化するベイズ推論手法である。 この方法論の重大な欠点は，多くの興味深い事例においてこの事後情報が得られないことである。 しかし，最近の Bayes 推論の発展により，事後分布を近似することでこの難点を回避することができるようになった。 これらは「変分 Bayes」と呼ばれ，このチュートリアルの焦点となるものである。

3. グラフィカルモデル (p4L)

グラフィカルモデルは，統計モデリング問題の確率変数間の依存関係を表現するための枠組みを提供し，確率系に関わる項目間の相互作用をグラフィカルに表現する包括的でエレガントな方法を構成する。 グラフィカルモデルは，問題の確率変数に対応するノードと，変数間の依存関係を表すエッジからなるグラフである。 グラフのノード A からノード B への有向エッジは，変数 B が変数 A の値に確率的に依存することを表す。 グラフィカルモデルには有向グラフと無向グラフとがある。 後者の場合，マルコフ確率場 (13,14,15) とも呼ばれる。 このチュートリアルの残りの部分では，ベイジアンネットワークとも呼ばれる有向グラフィカルモデルに焦点を当てる。 さらに，有向グラフは無サイクル (すなわちサイクルを含まない) であると仮定する。

$G=(V, E)$ を有向無サイクルグラフとし，$V$ をノードの集合， $E$ を有向エッジの集合とする。 また $x_ {s}$ はノード $s$ に関連する確率変数，$\pi(s)$ はノード $s$ の親の集合を表すとする。 また，各ノード $s$ には，その親変数の値が与えられたときの $x_ {s}$ の分布を定義する条件付き確率密度 $p(x_ {s}\vert x_ {\pi(s)})$ が関連付けられている。 したがって，グラフモデルを完全に定義するためには，グラフ構造とは別に，各ノードにおける条件付き確率分布も指定する必要がある。 これらの分布がわかれば，全変数の集合に対する結合分布は積として計算することができる。

\[p(x)=\prod_{s} p(x_{s}\vert x_{\pi(s)}).\tag{tz4}\]

上式は，有向グラフモデル [13] を，上式で指定された方法で因数分解する確率分布の集まりとして正式に定義したものである (もちろん，これは基礎となるグラフの構造に依存する)。

図 2 に有向グラフモデルの一例を示す。 ノードに描かれた確率変数は $a,b,c,d$ である。 各ノードは条件付き確率密度を表し，その親からのノードの依存性を定量化する。 ノードの密度は正確にはわからないかもしれないが，パラメータ $\theta_ {i}$ の集合によってパラメータ化することができる。 確率の連鎖則を用いて，結合分布は以下のように書ける: \(p(a,b,c,d;\theta) = P(a;\theta_{1}) p(b\vert a;\theta_{2}) p(c\vert a,b;\theta_{3}) p(d\vert a,b,c;\theta_{4})\tag{tz5}\)

しかし，グラフ構造が意味する独立性を考慮することで，この式を簡略化することができる。 一般に，グラフィカルモデルでは，各ノードはその親を与えられた先祖から独立している。 つまり，ノード $c$ はノード $b$ があればノード $a$ に依存しないし，ノード $d$ はノード $b$ と $c$ があれば $a$ に依存しない。 したがって (tz4) 式からは次のように書くことができる:

\[p(a,b,c,d;\theta) = p(a;\theta_{1}) p(b\vert a;\theta_{2}) p(c\vert a;\theta_{3}) p(d\vert b,c;\theta_{4}).\tag{tz6}\]

グラフィカルモデリングで生じるもう一つの有用な特徴は，通常データセットと呼ばれるいくつかの観測値が存在する場合，確率変数は観測値が存在する観測型 (または可視型) と直接観測値が利用できない非観測型に区別されることである。 観測されたデータは，グラフィカルモデル構造によって記述される生成機構によって生成されると仮定することが有用であり，それは中間サンプリングおよび計算ステップとしての隠れ変数を含む。 また，グラフィカルモデルにはパラメトリックとノンパラメトリックがあることに注意しなければならない。 パラメトリックモデルであれば，グラフの一部のノードで条件付確率分布にパラメータが現れる，すなわち，これらの分布はパラメータ化された確率モデルである。

グラフィカルモデルが完全に決定されると (すなわち，すべてのパラメータが指定されると)，確率変数の部分集合の周辺分布の計算，残りの変数の値を与えられた変数の部分集合の条件付き分布の計算，いくつかの以前の密度における最大点の計算など，いくつかの推論問題を定義することができる。 グラフィカルモデルがパラメトリックである場合，ある観測データセットが与えられたときに，パラメータの適切な値を学習する問題がある。 通常，パラメータ学習の過程では，いくつかの推論段階が含まれる。

4. EM アルゴリズムの別見解 (p5L)

本稿では [16] と [13] の EM の解説に従うことにする。 対数尤度は次のように書けることは簡単である。

• [13] C. Bishop, Pattern Recognition and Machine Learning. New York: Springer-Verlag, 2006.
• [16] R.M. Neal and G.E. Hinton, “A view of the EM algorithm that justifies incremental, sparse and other variants,” in Learning in Graphical Models, M.I. Jordan, Ed. Cambridge, MA: MIT Press, 1998, pp. 355–368.

\(\ln p(x;\theta) = F(q,\theta)+ D_ {\text{KL}}(q\vert\vert p)\tag{tz7}\) ここで \(F(q,\theta) = \int q(z) \ln\left(\frac{p(x,z;\theta)}{q(z)}\right)\,dz\tag{tz8}\) であり，かつ， \(D_{\text{KL}}(q\vert\vert p) = − \int q(z)\ln\left(\frac{p(z|x;\theta)}{q(z)}\right)\,dz\tag{tz9}\)

ここで $q(z)$ は任意の確率密度関数である。 $D_{\text{KL}}(q||p)$ は $p(z|x;\theta)$ と $q(z)$ の間の カルバック・ライブラー 情報量 (Kullback-Leibler divergence; あるいは相互情報量) である。 $D_{\text{KL}} (q||p) \ge 0$ なので $\ln p(x;\theta)\ge F(q,\theta)$ が成立する。 つまり $F(q,\theta)$ は対数尤度の下界である。 等式が成り立つのは $D_{KL}(q||p) = 0$ のときだけである。 これは $p(z|x;\theta)=q(z)$ を意味する。 EM アルゴリズムやベイズ推定の決定論的近似における最近のいくつかの進歩は，密度 $q$ とパラメータ $\theta$ に関する下界 $F(q,\theta)$ の最大化として (7) 式の分解に照らして見ることができる。

特に EM は下界 $F(q,\theta)$ を最大化する 2 段階の反復アルゴリズムである。 したがって EM アルゴリズムは，モデルの対数尤度を最大化する。 パラメータの現在値を $θ^{\text{OLD}}$ とする。 E-ステップでは，下界 $F(q, \theta^{\text{OLD}}$ は $q(z)$ に関して最大化される。 これは $D_ {\text{KL}}(q\vert\vert p) = 0$ のとき，換言すれば $q(z) = p(z\vert x;\theta^{\text{OLD}})$ のときに起こることは容易に理解できる。 この場合，下界は対数尤度に等しくなる。 続く M-step では $q(z)$ は固定され，下界 $F(q,\theta)$ は $\theta$ に関して最大化されて，ある新しい値 $\theta^{\text{NEW}}$ を与える。 これにより，下界は増加し，その結果，対応する対数尤度も増加する。 $q(z)$ は $\theta^{\text{OLD}}$ を用いて決定され，M ステップで固定されるため，新しい事後確率 $p(z|x;\theta^{\text{NEW}})$ と等しくはならず，したがって KL ダイバージェンスはゼロにはならない。 従って，対数尤度の増加は下界の増加より大きい。 $q(z) = p(z|x;\theta^{OLD})$ を下界に代入し，式 (8) を展開すると次のようになる。

\[\begin{aligned} F(q,\theta) &= \int p(z\vert x;\theta^{\text{OLD}}) \ln p(x,z;\theta) dz - \int p(z\vert x;\theta^{\text{OLD}}) dz\\ &= Q(\theta,\theta^{\text{OLD}})+\text{定数}\\ \end{aligned}\tag{tz10}\]

ここで定数は $\theta$ に依存しない $p(z\vert x;\theta^{\text{OLD}})$ のエントロピーを単純に表したものである。 また，以下の関数は M-ステップで最大化される完全データ (観測値＋隠れ変数) の対数尤度の期待値である:

\[\begin{aligned} Q(\theta,\theta^{\text{OLD}}) &= \int p(z\vert x;\theta^{\text{OLD}})\ln p(x,z;\theta)\,dz\\ &= \left<\ln p(x,z;\theta\right>_{p(z\vert x;\theta^{\text{OLD}})}\\ \end{aligned}\tag{tz:11}\]

信号処理の文献で EM アルゴリズムを示す通常の方法は $Q(\theta,\theta^{\text{OLD}})$ 関数を直接使用することである。 ([9] と [17] を参照)。

要約すると EM アルゴリズムは，以下の 2 段階を含む反復算法である。

\[\begin{aligned} \text{E-step: Compute } & P(z\vert x;\theta^{\text{OLD}})\\ \text{M-step: Evaluate } & \theta^{\text{NEW}}=\arg\max Q(\theta,\theta^{\text{OLD}})\\ \end{aligned}\tag{tz12,tz13}\]

さらに EM アルゴリズムでは，$p(z|x;\theta)$ が明示的に分かっているか，少なくともその十分統計量 $\ln p(z|x;\theta) p(z|x;\theta^{\text{OLD}})$ の条件付き期待値を計算できる必要があることを指摘したい (11 参照)。 換言すれば，EM アルゴリズムを使用するためには，オブザベーションが与えられたときの隠れ変数の条件付き pdf を知らなければならない． $p(z|x;\theta)$ は一般的に $p(x;\theta)$ よりもずっと推論しやすいが，多くの興味深い問題ではこれは不可能で，したがって EM アルゴリズムが適用できない。

5. 変分 EM の枠組み (p5R)

式 (eq:tz7) の分解において適切な $q(z)$ を仮定することで $p(z|x;\theta)$ を正確に知るという要件を回避することができる。 E ステップでは $\theta$ を固定したまま $F(q,\theta)$ を最大化するような $q(z)$ を見つける。 この最大化を行うためには $q(z)$ の特定の形式を仮定しなければならない。 特定のケースでは $q(z;\omega)$ の形式の知識を仮定することが可能です ($\omega$ はパラメータのセット)。 したがって，下界 $F(\omega,\theta)$ はこれらのパラメータの関数となり，E-step では $\omega$ に関して，M-step では $\theta$ に関して最大化される。

しかし，一般的な形式では，下界 $F(q,\theta)$ は $q$ についての 汎関数 である。 換言すれば $q(z)$ の関数を入力として受け取り，関数の値を出力として返す写像である。 これは自然に汎関数の微分の概念につながり，関数の微分と同様，入力関数の無限小の変化に対する関数の変化を与える。 この分野の数学は 変分計算 と呼ばれ，流体力学，熱伝導，制御理論など，数学，物理科学，工学の多くの分野に応用されている。

変分理論には近似値はない。 だが，変分法はベイズ推論問題の近似解を求めるのに使用できる。 これは，最適化が行われる関数が特定の形を持っていると仮定することによって行われる。 例えば 2 次関数や，固定基底関数の線形結合である関数のみを仮定することができる。 ベイジアン推論において，大きな成功を収めている特定の形式は因子化された形式である ([19] と [20] を参照)。 この因子化近似のアイデアは，理論物理学に由来しており，それは 平均場理論 mean field theory と呼ばれている。

この近似によれば，隠れ変数 $z$ は $i=1,\ldots,M$ で $M$ 個の分割 $z_{i}$ に分割されると仮定される。 また $q(z)$ はこれらの分割に関して次のように因数分解すると仮定する。

したがって，下界 $F(q,\theta)$ を最大化する (tz14) の形の $q(z)$ を求めたい。 (tz14) を用いて，簡単のために $q_{j} (z_{j})=q_{j}$ とすると，次のようになる。

\[q(z)=\prod_ {i=1}^{M} q_ {i}(z_ {i}).\tag{tz14}\]

したがって，下界 $F(q,\theta)$ を最大化する 式(14) の形の $q(z)$ を求めたい。 式 (14) を用いて，簡単のために $q_{j}(z_{j})=q_{j}$ とすると，次のようになる。

\(\begin{aligned} F(q,\theta) &= \int \prod_{i} q_ {i}\left[\ln q(x,z;\theta) -\sum_{i}\ln q_{i}\right]dz\\ &=\int \prod_{i} q_ {i}\ln p(x,z;\theta)\prod_{i} dz_{i} - \sum_{i}\int\prod_{j} q_{j}\ln q_{i} dz_{i}\\ &=\int q_{j}\left[\ln p(x, z;\theta) \prod_{i\ne j}(q_{i}, dz_{i})\right] dz_{j} -\int q_{j}\ln q_{j} dz_{j} -\sum_{i\ne j}\int q_{i} \ln q_{i} dz_{i}\\ &= \int q_ {j}\ln \hat{p}(x,z_ {j};\theta) dz_ {j} - \int q_{j} \ln q_ {j}\, dz_ {j} -\sum\int q_ {i}\ln q_ {i} dz_ {i}\\ &=D_{\text{KL}}(q_{i}\vert\vert \hat{p})-\sum_{i\ne j} \int q_{i}\ln q_{i}dz_{i}\\ \end{aligned}\tag{tz15}\) ここで， \(\ln \hat{p}(x,z_ {j};\theta)=\left<\ln p(x,z;\theta)\right>_ {i\ne j}=\int \ln p(x,z;\theta)\prod_ {i\ne j}(q_ {i},dz_ {i}).\)

明らかに 式 (15) の境界は，Kullback-Leibler 距離が 0 になったときに最大となり，それは $q_{j}(z_ {j})=\hat{p}(x,z_{j}; \theta)$ の場合である。 換言すれば最適分布 $q^ {∗}{j}(z{j})$ の式は以下の通りである:

\[\ln q_{j}^ {\star} (z_{j})=\left< p(x,z;\theta)\right>_{i\ne j}+\text{const}.\tag{tz16}\]

式 (16) の加法定数は正規化によって求めることができるので

\[q^{\star}_{z_{j}}=\frac{\exp^{\left< \ln p(x,z;\theta)\right>_{i\ne j}}}{\int \exp^{\left< \ln p(x,z;\theta)\right>_{i\ne j}}dz_{j}}.\tag{tz17}\]

$j=1,\ldots,M$ に対する上式は，式 (14) の因数分解を前提とした下界の最大値に対する整合性条件の集合である。 これらは $i\ne j$ の他の因子 $q_{i}(z_{i})$ に依存するため，明示的な解を提供しない。 したがって，これらの因子を循環させ，それぞれを修正された推定値で順番に置き換えることにより，一貫した解を求めることができる。

要約すると，変分 EM アルゴリズムは，以下の 2 つのステップで与えられる。

• 変分 E-ステップ: 式 (17) を解いて，$F(q,\theta^{\text{OLD}})$ を最大化するために，$q^{\text{NEW}}(z)$ を評価する
• 変分 M-ステップ: $\theta^{\text{NEW}} = \arg\max F(q^{\text{NEW}},\theta)$ を探索する

このとき，ベイズモデルが隠れ変数のみを含み，パラメータを含まない場合があることに注目する必要がある。 そのような場合，変分 EM アルゴリズムは，式 (17) を用いて $q(z)$ を得る E-ステップのみを有する。 この 関数 $q(z)$ は，隠れ変数の推論に使用できる $p(z|x)$ の近似を構成する。

6. 線形回帰 Linear Regression (p6L)

本節では，線形回帰問題を例にして，前節までのベイズ推論の手法を示す。 線形回帰が選ばれた理由は，単純であり，入門的な例として優れているためである。 さらに，線形回帰は，デコンボリューション，チャネル推定，音声認識，周波数推定，時系列予測，システム同定に至るまで，多くの信号処理アプリケーションで発生する。

この問題では，未知信号 $y(x) \in\mathbb{R}, x\in\Omega\subseteq \mathbb{R}^{N}$ を考え，任意の位置 $x_{\star}\in\Omega$ におけるその値 $t_{\star}=y(x_{\star})$ を予測したいとする。 $n=1,\ldots,N$, $x=(x_{1},\ldots,t_{N})$,$x_{n}\in\Omega$, $n=1,\ldots,N$ の位置で$t_{n}=y(x_{n}) +\epsilon_{n}$ の雑音付き観測値を $\mathbf{t}=(T_{1},\ldots,T_{N})$ ベクトルとすると，$n$ の雑音がある。 加法性雑音 $\epsilon_{n}$ は一般に独立，平均ゼロのガウス分布と仮定する。

\[p(\epsilon) = \mathcal{N}(\epsilon\vert 0,\beta^{-1}\mathbf{I}),\tag{tz18}\]

ここで，$\beta$ は，分散行列の逆行列 $\epsilon=(\epsilon_{1},\ldots,\epsilon_{N})^{\top}$ である。

信号 $y$ は一般に $M$ 個の基底関数 $\phi_ {m}(x)$ の線形結合としてモデル化される。

\[y(x)=\sum_{m=1}^{M} w_{m}\phi_{m}(x),\tag{tz19}\]

ここで $w=(w_{1},\ldots, w_{M})^{\top}$ は線形結合の重みを表す。 設計行列 $\mathbf{\Phi}=( \phi_{1},\ldots, \phi_{M}$) で，$\phi_{m}=(\phi_{m}(x_ {1}),\ldots,\phi_{m}(x_{N})^{\top}$ を定義し，観測値 $t$ は次のようにモデル化される:

\[\mathbf{t}=\mathbf{\Phi}\mathbf{w}+\epsilon\tag{tz20}\]

そして，尤度は次式で与えられる:

\[p(t;w,\beta) = \mathcal{N}(t\vert \Phi W,\beta^{-1}\mathbf{I}). \tag{tz21}\]

以下では，前節までの理論を線形回帰問題に適用し，この線形モデルの未知の重み $w$ を計算する 3 つの方法論を示す。 まず，パラメータと仮定された重みの典型的な ML 推定を適用する。 これから示されるように，パラメータの数は観測値の数と同じであるため，ML 推定はモデルの雑音に非常に敏感で，観測値に過剰に適合してしまう。 そこで，この問題を改善するために，確率変数と仮定される重みに事前分布を課した。 まず，重みの定常ガウス事前分布に基づく単純なベイズモデルが使用される。 このモデルでは，EM アルゴリズムを用いてベイズ推定が行われ，結果として得られる解は雑音に対して頑健である。 しかしながら，このベイズモデルは非常に単純であり，局所的な信号の特性を捕らえる能力を持たない。 この目的のために，重みの非定常ガウス事前分布とハイパー事前確率に基づく，より洗練された空間的に変化する階層的モデルを導入することが可能である。 このモデルは，EM アルゴリズムで解くには複雑すぎる。 この目的のために 変分推論 EM の枠組み 節で説明した変分ベイズ法が，このモデルの未知数の値を推論するために使用されている。 最後に，ベイズモデルの複雑性が解の質を向上させることを実証するために，3 つの方法を用いて簡単な人工信号の推定値を得る。 図 3 a, b, c は，線形回帰の 3 つのアプローチのグラフィカルモデルである。

線形モデルの重み $w$ の最も単純な推定値は，モデルの尤度を最大化することによって得られる。 この ML 推定は，図 3a のグラフィカルモデルに示すように，重み $w$ をパラメータと仮定している。 ML 推定値は尤度関数を最大化することにより得られる \(p(t;w,\beta) =\left(2\pi\right)^{-\frac{N}{s}}\beta^{\frac{N}{2}}\exp\left(-\frac{\beta}{2}\left\|\mathbf{t}-\mathbf{\phi w}\right\|\right)\) これは，$E_{LS(w)}=\left|t-\mathbf{\phi}w\right|^{2}$ を最小化することと等価である。 したがって，この場合，ML は最小二乗 (LS) 推定値と等価である \(w_{LS}={\arg\max}_{w} p(t;w,\beta)={\arg\min}_{w}\mathbb{E}_{LS(w)}=\left(\mathbf{\phi}^{\top}\mathbf{\phi}\right)^{-1}\mathbf{\phi}^{\top}t\tag{tz22}\)

多くの状況で (そして使用される基底関数に依存して)，行列 $\mathbf{\phi}^{\top}\mathbf{\phi}$ は条件不一致で反転が困難な場合がある。 これは，雑音 $\epsilon$ が信号の観測値に含まれる場合，重みの推定値 $w_{LS}$ に大きく影響することを意味する。 したがって，ML 線形回帰を使用する場合，行列 $\mathbf{\phi}^{\top}\mathbf{\phi}$ が反転できるように基底関数を慎重に選択する必要がある。 これは一般に，基底関数が少ないスパースモデルを使用することで達成され，また推定しなければならないパラメータが少ないという利点もある。

6.1 EM アルゴリズムに基づく，ベイジアン線形回帰 EM-Based Bayesian Linear Regression (p7L)

線形モデルのベイズ的な取り扱いは，まずモデルの重みに 事前分布 を割り当てることから始まる。 これは推定にバイアスをもたらすが，ML 推定の大きな問題であるその分散を大幅に減少させる。 ここでは，線形モデルの重みに独立で，かつ，平均ゼロのガウス型事前分布をよく選択することを考える。

\[p(w;\alpha)= \prod_{m=1}^{M} \mathcal{N}(w_{m}\vert 0,\alpha^{-1}).\tag{tz23}\]

これは定常的な事前分布であり，すべての重みの分布が同一であることを意味する。 この問題のグラフィカルモデルを 図 [@2008TzikasFig3(b)] に示す。 ここで，重み $w$ は隠れ確率変数であり，モデルパラメータは $w$ の事前分布のパラメータ $\alpha$ と加法性ノイズの逆分散 $\eta$ であることに注意されたい。

ベイズ推論は，隠れ変数の事後分布を計算することで進められる。 \(p(w\vert t;\alpha,\beta) = \frac{p(t\vert w;\beta) p(w;\alpha)}{p(t;\alpha,\beta)}.\tag{tz24}\)

分母に現れる周辺尤度 $p(t;\alpha,\beta)$ は解析的に計算できることに注意。 \(p(t;\alpha,\beta) =\int p(t\vert w;\beta)\;p(w;\alpha)\;dw = \mathcal{N}\left(t\vert 0,\beta^{-1}\mathbf{I}+\alpha^{-1}\mathbf{\Phi\Phi}^{\top}\right).\tag{tz25}\)

このとき，隠れ変数の事後確率は: \(p(w\vert t;\alpha,\beta)= \mathcal{N}(w\vert \mu,\sigma),\tag{tz26}\) ここで \(\begin{aligned} \mathbf{\mu} &=\beta\,\mathbf{\Sigma\Phi}^{\top} t.\\ \mathbf{\Sigma} &=\beta\,\mathbf{\Phi}^{\top} \mathbf{\Phi} + \alpha\mathbf{I}^{-1}.\\ \end{aligned} \tag{tz27,tz28}\)

モデルのパラメータは、周辺尤度 $p(t;\alpha,\beta)$ の対数を最大化することにより推定できる:

\[(\alpha_{ML},\beta_{ML}) =\arg\min_ {\alpha,\beta} \left[\log\left|\beta^{-1}\mathbf{I} + \alpha^{-1}\mathbf{\Phi\Phi}^{\top}\right|\right] + \mathbf{t}^{\top} \left[ \beta^{-1} \mathbf{I} + \alpha^{-1} \mathbf{\Phi\Phi}^{\top-1}t \right].\tag{tz29}\]

式 (29) はパラメータ $(\alpha,\beta)$ に関する導関数が計算しにくいため，直接最適化を行うにはいくつかの計算上の困難が伴う。 また $(\alpha,\beta)$ の推定値は逆分散を表すため正でなければならず，制約付き最適化アルゴリズムが必要である。 その代わりに，先に述べた EM アルゴリズムにより $(\alpha,\beta)$ の推定値の取得と $w$ の値の推定を同時に行う効率的な枠組みを提供する。 EM アルゴリズムでは，周辺尤度 式(tz:25) の計算を行わないが，その局所最大値に収束することに注意すること。 パラメータ $\left{\alpha^ {0},\beta^ {0}\right}$ を適当な値に初期化した後，アルゴリズムは以下のステップを繰り返し実行することで進められる。

6.1.1 E step:

完全尤度対数の期待値を計算:

\[\begin{aligned} Q^{(t)}(t,w;\alpha,\beta) &= \left<\log p(t,w;\alpha,\beta)\right>_{p(w\vert t;\alpha^{(t)},\beta^{(t)})}\\ & =\left<\log p(t\vert w;\alpha,\beta) p(w;\alpha,\beta)\right>_{p(w\vert t;\alpha^{(t)},\beta^{(t)})}.\\ \end{aligned}\tag{tz30}\]

上式は，(21)，(23) 式を用いて計算され:

\[\begin{aligned} Q^{(t)}(t,w;\alpha,\beta) &= \left<\frac{N}{2} \log\beta - \frac{\beta}{2}\left(\left\|\mathbf{t}-\mathbf{\Phi\mu}^{(t)}\right\|^{2}\right)\right> + \text{const}\\ &= \frac{N}{2} \log\beta - \frac{\beta}{2}\left(\left\|\mathbf{t}-\mathbf{\Phi\mu}^{(t)}\right\|^{2}\right) + \frac{M}{2}\log\alpha - \frac{\alpha}{2}\left(\left<\left\|\mathbf{w}\right\|^{2}\right>\right)+\text{const}\\ \end{aligned} \tag{tz31}\]

これらの期待値は $p(w\vert t;\alpha^{t},\beta^{(t)})$ に対するもので，式 (26) から計算すると次のようになる:

\[Q^{(t)}(t,w;\alpha,\beta) = \frac{N}{2} \log\beta - \frac{\beta}{2}\left(\left\|\mathbf{t}-\mathbf{\Phi\mu}^{(t)}\right\|^{2} + \text{tr}\left| \mathbf{\Phi}^{\top}\mathbf{\Sigma}^{(t)} \mathbf{\Phi} \right|\right) + \frac{M}{2} \log\alpha - \frac{\alpha}{2}\left(\left|\mathbf{\mu}^{(t)}\right|^{2} + \text{tr}\left|\mathbf{\Sigma}^{(t)}\right|\right) + \text{const} \tag{tz32}\]

ここで $\mathbf{\mu}^{(t)}$ と $\mathbf{\Sigma}^{(t)}$ は，パラメータ $\alpha^{(t)}$ と $\beta^{(t)}$ の現在の推定値を使用して計算される:

\[\begin{aligned} \mathbf{\mu}^{(t)} &= \beta^{(t)} \mathbf{\Sigma}^{(t)} \mathbf{\Phi}^{\top} \mathbf{t},\\ \mathbf{\Sigma}^{(t)} &= \left(\beta^{(t)} \mathbf{\Phi}^{\top} \mathbf{\Phi} + \alpha^{(t) \mathbf{I}}\right)^{-1}.\\ \end{aligned}\tag{tz33,tz34}\]

6.1.2 M step:

$Q^{(t)}(t,w;\alpha,\beta)$ をパラメータ $\alpha$ と $\beta$ に関して最大化する。

\[\left(\alpha^{(t+1)},\beta^{(t+1)}\right)=\arg\max_{(\alpha,\beta)}Q^{(t)}\left(t,w;\alpha,\beta\right).\tag{tz35}\]

$Q^{(t)}(t,w;\alpha,\beta)$ のパラメータに関する導関数は次の通りである。

\[\begin{aligned} \frac{\partial Q^{(t)}}{\partial\alpha} &=\frac{M}{2\alpha}-\frac{1}{2}\left(\left\|\mathbf{\mu}^{(t)}\right\|^{2}+\text{tr}\left|\Sigma^{(t)}\right|\right)\\ \frac{\partial Q^{(t)}}{\partial\beta} &=\frac{N}{2\beta}-\frac{1}{2}\left(\left\|t-\mathbf{\phi\mu}^{(t)}\right\|^{2} +\text{tr}\left|\mathbf{\phi}^{\top}\mathbf{\Sigma}^{(t)}\mathbf{\phi}\right|\right). \end{aligned}\tag{tz36,tz37}\]

これらを 0 とすると，パラメータ $\alpha$, $\beta$ を更新する以下の式が得られる:

\[\begin{aligned} \alpha^{(t+1)} &= \frac{M}{\left\|\mathbf{\mu}^{(t)}\right\|^{2}+\text{tr}\left|\mathbf{\Sigma}^{(t)}\right|}\\ \beta^{(t+1)} &= \frac{N}{\left\|t-\mathbf{\phi\mu}^{(t)}\right\|^{2}+\text{tr}\left|\mathbf{\Phi}^{\top}\mathbf{\Sigma}^{(t)}\mathbf{\Phi}\right|} \end{aligned}\tag{tz38,tz39}\]

数値最適化を必要とする (25) の周辺尤度の直接最大化とは対照的に，最大化ステップは解析的に実行できることに注意されたい。 さらに (38) と (39) はパラメータ $\alpha$ と $\beta$ の正の推定を保証しており，これらは逆分散パラメータを表しているので，これは必要条件である。 ただし，パラメータの初期化によっては，異なる局所最大値になる可能性があるため，初期化には注意が必要である。 E-step では事後統計量 $p(w|t;\alpha,\beta)$ が計算されているので，$w$ に対する推論が直接得られる。 この事後統計量の平均値 (33) は $w$ のベイズ線形最小平均二乗誤差 (LMMSE) 推論として用いることができる。

6.2 変分 EM に基づく Bayesian 線形回帰 (P9L)

前節で述べたベイズ的アプローチでは，線形モデルの重みに定常ガウス事前分布を用いるため，限界尤度の厳密な計算が可能で，ベイズ推定は解析的に行われる。 しかし，多くの場面で，信号の局所的な特性を柔軟にモデル化することが重要であり，これは単純な定常ガウス事前分布では不可能である。 このため，各重みに対して明確な逆分散 $\alpha_ {m}$ を持つ非定常ガウス事前分布が検討される。

\[p\left(w\vert\mathbf{\alpha}\right) = \prod_ {m=1} ^ {M}\mathcal{N}\left(w_ {m}\vert 0, \alpha_ {m}^{-1}\right)\tag{tz40}\]

しかし，このようなモデルは，推定すべきパラメータとほぼ同数の観測値があるため，パラメータが過剰になる。 このため 精度パラメータ $\mathbf{\alpha}=\left(\alpha_ {1},\ldots,\alpha_ {M}\right)^{\top}$ を確率変数として扱い，事前ガンマ分布を課すことにより制約をかける。

\[p\left(\mathbf{\alpha};a,b\right) = \prod_ {m=1}^ {M}\text{Gamma}(\alpha_ {m}\vert a,b)\tag{tz41}\]

この事前分布は Gaussian 分布と共役であることから選択された[13]。 さらに，ノイズの逆分散 $\beta$ の事前分布としてガンマ分布を仮定する

\[p(\beta;c,d)=\text{Gamma}(\beta\vert c, d).\tag{tz42}\]

この Bayes プローチのグラフィカルモデルは 図 3(c) のようになり，隠れ変数 $w$ の隠れ変数 $\alpha$ への依存性が明らかになる。 また，このモデルのパラメータ $a, b, c, d$ とそれらに依存する隠れ変数も明らかにされている。

Bayes 推論では，事後分布の計算が必要であり:

\[p(w,\alpha\beta\vert t) = \frac{p(t\vert w, \beta)p(w\vert\alpha)p(\alpha)p(\beta)}{p(t)}.\tag{tz43}\]

しかし，周辺尤度 $\displaystyle p(t)=\int p(t\vert w,\beta)p(w\vert \alpha)p(\alpha)dw\;d\alpha\;d\beta$ は解析的に計算できないため，式(43) の正規化定数を計算することができない。 そこで，近似的な Bayes 推論手法，特に変分推論手法に頼ることになる。 重み $w$ と分散パラメータ $\alpha,\beta$ の間に事後独立性があると仮定する。

\[p(w,\alpha,\beta\vert t; a, b, c, d)\approx q(w,\alpha,\beta) = q(w)q(\alpha)q(\beta),\tag{tz44}\]

は (16) 式から次のように近似的な事後分布 $q$ を計算することができる。 $\ln q(w)$ のうち $w$ に依存する項だけを残して，次のようになる。

$$\begin{aligned} \ln q(w) &=\left<\ln p(t,w,\mathbf{\alpha},\beta){q(\mathbf{\alpha})q(\beta)}\right>+\text{const}
&= \left<\ln p(t|w,\beta) p(w|\mathbf{\alpha}){q(\mathbf{\alpha})q(\beta)}\right>+\text{const}
&= \left<\ln p(t|w,\beta) + \ln p(w|\mathbf{\alpha}){q(\mathbf{\alpha})q(\beta)}\right>+\text{const}
&= \left<-\frac{\beta}{2}\left(t-\Phi w\right)^{\top}\left(t-\Phi w\right)-\frac{1}{2}\sum{m=1}^{M}\alpha_{m}w_{m}^{2}\right>+\text{const}
&=\frac{\left<\beta\right>}{2}\left[t^{\top}t-2t^{\top}\Phi w+w^{\top}\Phi^{\top}\Phi w\right]

• \frac{1}{2}\sum_{m=1}^{M}\left<\alpha_{m}\right>w_{m}^{2}+\text{const}
  &=-\frac{1}{2}w^{\top}\left(\left<\beta\right>\Phi^{\top}\Phi+\left<\mathbf{A}\right>\right)w-\left<\beta\right>w^{\top}\Phi^{\top}t+\text{const}
  &=-\frac{1}{2}w^{\top}\Sigma^{-1}w - w^{\top}\Sigma^{-1}\mathbf{\mu}+\text{const} \end{aligned}\tag{tz45}$$ ここで，$\mathbf{A}=\text{diag}\left(\alpha_{1},\ldots,\alpha_{M}\right)$.

これは，平均 $\mu$ 分散行列 $\Sigma$ に従う指数 Gauss 分布であり，それぞれ次式に従う:

\[\begin{aligned} \mathbf{\mu} &=\left<\beta\right>\mathbf{\Sigma\Phi}^{\top}t.\\ \mathbf{\Sigma} &= \left(\left<\beta\right>\mathbf{\Phi}^{\top}\mathbf{\Phi}+\left<A\right>\right)^{-1}. \end{aligned}\tag{taz46,tz47}\]

それゆえ，$q(w)$ は，次式で与えられる:

\[q(w)=\mathcal{N}\left(w|\mathbf{\mu},\mathbf{\Sigma}\right).\tag{tz48}\]

事後確率 $q(\mathbf{\alpha})$ は同様に $\mathbf{\alpha}$ に依存する $\ln q(\mathbf{\alpha})$ を計算することで得られる。

\[\begin{aligned} \ln q(\mathbf{\alpha}) &= \left<\ln p(t,w,\mathbf{\alpha},\beta\right>_{q(w) q(\beta)}\\ \end{aligned}\tag{tz49}\]

これは，パラメータ $\tilde{a},\tilde{b}_m$ を持つ $M$ 個の独立なガンマ分布の積の指数で，次式で与えられる:

\(\begin{aligned} \tilde{a} &= a+\frac{1}{2}\\ \tilde{b}_{m} &= b + \frac{1}{2}\left<w_{m}^{2}\right>. \end{aligned}\tag{tz50,tz51}\) それゆえ $q(\mathbf{\alpha})$ は，次式で与えられ: \(q(\alpha)=\prod_{m=1}^{M}\Gamma(\alpha_{m}|\tilde{a},\tilde{b}_{m}).\tag{tz52}\)

雑音の逆分散の事後分布も同様に次のように計算される:

\(\tag{tz53} q(\beta)=\Gamma(\beta|\tilde{c},\tilde{d}_m),\) with \(\tag{tz54,tz55} \begin{aligned} \tilde{c} &= c+\frac{N}{2},\\ \tilde{d} &= d+\frac{1}{2}\left<\left\|\mathbf{t}-\mathbf{\Phi w}\right\|^{2}\right>.\\ \end{aligned}\)

そして，(48), (52), (53) の近似事後分布は，互いの統計量に依存するため，収束するまで繰り返し更新される，詳細は [22] を参照のこと。

ここで，重みの真の事前分布は，ハイパーパラメータ $\alpha$ を周辺化することで計算できることに注意。

\[\tag{tz56} \begin{aligned} p(w,a,b) &= \int p(w|\alpha) p(\alpha;a,b)d\alpha\\ &= \int \prod_{m=1}^{M} N\left(w_m|0,\alpha_m^{-1}\right)\Gamma (a_m|a,b)d\alpha_{m}\\ &= \prod_{m=1}^{m} S t(w_m|\lambda,\nu) \end{aligned}\]

and is a student-t pdf, \(St(x|\mu,\lambda,\nu)=\frac{\Gamma((\nu+1)/2}{\Gamma(\nu/2)}\left(\frac{\lambda}{\pi\nu}\right)^{1/2}\times \left[1+\frac{\lambda(x-\mu)^{2}}{\nu}\right]^{-(\nu+1)/2},\)

で 平均 $\mu=0$, パラメータ $\lambda= a/b$, 自由度 $\nu=2a$。 この分布は自由度 $\nu$ が小さい場合，重尾分布となる。 したがって，この分布は，少数の基底関数だけを含み，残りの基底関数を対応する重みを非常に小さな値に設定することによって刈り込む，疎な解に有利である。 最終的なモデルで実際に使用される基底関数を関連性基底関数と呼ぶ。

簡単のために，Student-t 分布のパラメータ a, b, c, d を固定とした。 実際には，これらのパラメータを非常に小さな値，すなわち，$a=b=c=d= 10^{-6}$ に設定することによって得られる非情報的な分布を仮定することによって，しばしば良い結果を得ることができる。 あるいは，変分 EM アルゴリズムを用いて，これらのパラメータを推定することもできる。 このようなアルゴリズムは，これらのパラメータに関して変分境界を最大化する M-ステップを，説明した方法に追加することになる。 しかし，ベイズモデリングにおける典型的なアプローチは，ハイパーパラメータを固定し，モデルの最上位レベルで非情報的なハイパー事前分布を定義することである。

6.3 線形回帰の例 (p10R)

次に，先に説明した線形回帰モデルの特性を示す数値例を示す。 また，変分 Bayes 推論を用いることで得られる利点を示す。 人工的に生成された信号 $y(x)$ を用いるので，観測値を生成した元の信号は既知であり，したがって推定の品質を評価することができる。 信号の $N=50$ サンプルを得て，分散 $\sigma^{2}=4\times10^{-2}$ のガウス雑音を加え，SN 比 $\text{SNR}=6.6$ dB に相当させた。 $N$ 個の基底関数を用い，特に各信号観測位置を中心とした基底関数を 1 個用いた。 基底関数は以下の形のガウシアンカーネルである。

\[\phi_{i}(x)=K(x,x_{i})=\exp\left(-\frac{1}{2\sigma^{2}_{\phi}}\left\|x-x_{i}\right\|^{2}\right).\tag{tz57}\]

この例では，非常に重い尾を持つ，情報量の少ない Student-t 分布を得るために $a=b=0$ と設定した。

次に，この観測結果を用いて i) 最尤推定 (ML 推定) (22), ii) EM-に基づくベイズ推論 (33), iii) 変分 EM-に基づく ベイズ推論 (46) により、信号の出力を予測することができた。 結果は 図4(a) に示されている。 ML 推定は雑音の多い観測値に忠実に従うことに注目されたい。 したがって，平均二乗誤差の点では最悪である。 この定式化では，観測値と同じ数の基底関数を使用し，重みに制約がないため，これは予想されることである。 ベイズ法は，重みが事前分布によって制約されるため，この問題を克服している。 しかし，この信号には大きな分散を持つ領域と非常に小さな分散を持つ領域があるため，定常事前分布ではその局所的な振る舞いを正確にモデル化できないことは明らかである。 これに対して，階層的非定常事前分布はより柔軟で，より良い局所的な適合が得られる。 実際，後者の事前分布に対応する解は，基底関数の小さな部分集合しか使用しておらず，その位置は図 4 の丸で囲んだ観測値として示されている。 これは，我々は $a=b=0$ を設定したためで，情報量の少ない Student-t 分布を定義している。 したがって，ほとんどの重みは正確に 0 と推定され，信号推定に使用される基底関数はごくわずかである。 最終的なモデルで実際に使用されるこれらの基底関数は関連性基底関数と呼ばれ，それらが中心化されたベクトルは関連性ベクトル (RV: revelence vector) と呼ばれ，図 4 に示されている。

この例と同じ精神で，画像復元問題，画像超解像問題，画像ブラインドデコンボリューション問題に対して，それぞれ [23], [24], [25] で階層的な非定常事前分布が提案されている。 画像再構成問題において，このような事前分布は，画像のエッジを保持すると同時に，画像の平坦部における雑音を抑制する能力を実証した。 さらに，このような事前分布は，下絵が未知の場合の電子透かし検出器の設計にも用いられ，成功を収めている[26]。

7. Gauss 混合モデル (p10R)

Gauss 混合モデル (GMM) は，密度をモデル化するための貴重な統計ツールである。 任意の密度を高い精度で近似できる柔軟性があり，さらに，ソフトクラスタリング解として解釈することができる。 GMM は音声理解，画像モデリング，追跡，セグメンテーション，認識，電子透かし，ノイズ除去など様々な信号処理問題で広く用いられている。

GMM はガウス分布の凸組み合わせとして定義され，単一の分布では十分でない場合にデータセットの密度を記述するために広く用いられている。 M 個の成分を持つ混合モデルを定義するには，各成分 $j$ の確率密度 $p_j(x)$ と，成分の混合重み $\pi_{j}(\pi_{j}\ge0$ と $\sum_{j=1}^{M}\pi_{j}=1$) を含む確率ベクトル $(\pi_{1},\ldots,\pi_{M})$ を指定しなければならない。

このような混合物を用いてデータセット X の密度をモデル化する場合の重要な仮定は，各データムが以下の手順で生成されていることである：

1. 成分 k を確率ベクトル $\pi_1,\ldots,pi_M$ を用いてランダムにサンプリングする。
2. 成分 k の密度 $p_k(x)$ からサンプリングして観測を生成する。

ここで，離散隠れ確率変数 $Z$ は，観測サンプル $x$ を生成するために，すなわち，観測確率変数 $X$ に値 $X=x$ を割り当てるために選択された成分を表す。 このグラフモデルでは，ノード分布は $P(Z=j)=\pi_j$ と $P(X=x|Z=j)=p_j(x)$ である。 X と Z との同時確率密度関数については，次式が成り立つ \(p(\mathbf{X},Z) = p(\mathbf{X}|Z) p(Z)\tag{tz58}\) そして，$Z$ を周辺化することにより，混合モデルのよく知られた公式が得られる。 \(p(\mathbf{X}=x)=\sum_{j=1}^{M}p(\mathbf{X}=x|Z=j)p(Z=j)=\sum_{j=1}^{M}\pi_{j}p_{j}(x).\tag{tz59}\)

GMM の場合，各成分 $j$ の密度は $p_j(x)=N(x;\mu_j,\Sigma_j)$となり，$\mu_j\in\mathcal{R}^d$ は平均 $\Sigma_j$ は $d\times d$ の共分散行列を表す。 それゆえ， \(p(x)=\sum_{j=1}^{M}\pi_{j}N\left(x;\mathbf{\mu}_{j},\mathbf{\Sigma}_{j}\right).\tag{tz60}\)

図 4. グラフィカルモデル (a) 単一ガウス成分，(b) ガウス金剛モデル

混合モデルにおける注目すべき利便性は，ベイズの定理を用いて，ある観察 $x$ が混合成分 $j$ の分布からサンプリングして生成されたという事後確率 $P(j|x)=p(Z=j|x)$ を計算するのが簡単であることである。 \(P(j|X)=\frac{p(x|Z=j)p(Z=j)}{p(x)}=\frac{\pi_j N(x|\mu,\Sigma_j)}{\sum_{l=1}^{M}\pi_l N(x|\mu_l,\Sigma_l)}\tag{tz61}\) この確率は，観察 $x$ を生成した成分 $j$ の責任と呼ばれることもある。 さらに，データ点 $x$ を最大事後確率を持つ成分に割り当てることで，データ集合 $X$ を$M$ 個のクラスタにクラスタリングすることが簡単にできる。

7.1 Gauss 混合モデルの訓練のための EM

$X=\left{x_n|x_n\in\mathcal{R}^{d},n=1,\ldots,N\right}$ を，$M$ 個の成分を持つ GMM を使ってモデル化されるデータ点の集合とする： $p(x)=\sum_{j=1}^{M}\pi_j N(x_n|\mu_j,\Sigma_j)$. 成分数 $M$ はあらかじめ指定されているものとする。 推定する混合パラメータのベクトル $\theta$ は，混合重みと各成分のパラメータから構成される，すなわち $\theta=\left{\pi_j,\mu_j,\Sigma_j|j=1,\ldots,M\right}$.

パラメータ推定は，対数尤度の最大化によって達成できる。 \(\theta_{ML}=\arg\max_{\theta} \log p(\mathbf{X};\mathbf{\theta}),\tag{tz62}\) ここで，独立で同一分布の観察を仮定すると，尤度は次のように書ける。 \(p(X;\theta) =\prod_{n=1}^{N}p(X_n;\theta)=\prod_{n=1}^{N}\prod_{j=1}^{M}\pi_{j}N\left(x_n;\mu_j,\Sigma_{j}\right).\tag{tz63}\)

図 5(b) のグラフィカルモデルから，各観測変数 $x_{n}\in X$ は，$x_n$ を生成するために使用された成分を表す隠れ変数 $z_n$ に対応することが明らかである。 この隠れ変数は，$x_n$ が混合成分 $j$ から生成された場合に $z_{jn}=1$，そうでない場合に $z_{jn}=0$ となるような，$M$ 個の要素を持つ二値ベクトル $z_{jn}$ を用いて表現することができる。 $z_{jn}=1$ は確率 $\pi_{j}$ と $\sum_{j=1}^{M}\pi_{j}=1$ であるので，$z_n$ は多項分布に従う。 $Z=\left{z_n,n=1,\ldots,N\right}$ を隠れ変数の集合とする。 すると $(Z|\theta)$ は次のように書かける: \(p(Z;\theta)=\prod_{n=1}^{N}\prod_{j=1}^{M}\pi_j^{z_{jn}}\tag{tz64}\) かつ \(p(X|z;\theta)=\prod_{n=1}^{N}\prod_{m=1}^{M}N\left(x_n;\mu_j,\Sigma_j\right)^{z_{jn}}\tag{tz65}\)

前述したように，混合モデルの便利な問題は，(61) 式を用いて観察を与えられた隠れ変数の正確な事後値 $p\left(z_{jn}=1|x_n;\theta\right)$ を簡単に計算できることである。 したがって，厳密 EM アルゴリズムの適用が可能である。

より具体的には，$\theta^{(t)}$ が EM 反復 $t$ におけるパラメータベクトルを表すとすると，隠れ変数 $z_{jn}$ の事後 $p(z|x;\theta^{(t)})$ の期待値は次式で与えられる: \(\tag{tz66} \left<z_{jn}^{(t)}\right>=\frac{\pi_{j}^{(t)} N\left(\right)x_{n};\mu_{j}^{(t),\sigma_{j}^{(t)}}}{\sum_{j=1}^{M}\pi_{j}^{(t)} N\left(\right)x_{n};\mu_{j}^{(t),\sigma_{j}^{(t)}}}\)

上式は，$j=1,\ldots,M$ および $n=1,\ldots N$ に対して E ステップで実行されるべき計算を指定する。

事後分布 $p(Z|X;\theta(t))$ に対する完全対数尤度 $\log P(X,Z)$ の期待値は次式で与えられる。 \(\tag{tz67}\begin{aligned} Q\left(\theta,\theta^{(t)}\right) &= \left<\log p\left(X,Z;\theta\right)\right>\\ &=\left<\log p\left(X|Z;\theta\right)\right>+\left<\log p\left(Z;\theta\right)\right>_{p\left(z|x;\theta^{(t)}\right)}\\ &=\sum_{n=1}^{N}\sum_{j=1}^{M}\left<z_{jn}^{(t)}\right>\log\pi_{j} +\sum_{n=1}^{N}\sum_{j=1}^{M}\log N\left(x_{n};\mu_j,\Sigma_{j}\right)\\ \end{aligned}\)

M ステップでは，期待完全対数尤度 Q はパラメータ $\theta$ に関して最大化される。 対応する偏導関数をゼロに等しくし，制約条件 $\sum_{j=1}^{M}\pi_{j}=1$ に対してラLagrange 乗数を用いると，M ステップの更新について以下の式が導ける： \(\tag{68} \pi_{j}^{(t+1)}=\frac{1}{N}\sum_{n=1}^{N}\left<z_{jn}^{(t)}\right>\) \(\tag{69} \mu_{j}^{(t+1)}=\frac{\sum_{n=1}^{N}\left<z_{jn}^{(t)}\right>x_n}{\sum_{n=1}^{N}\left<z_{jn}^{(t)}\right> }\) \(\tag{70} \Sigma_{j}^{(t+1)}=\frac{\sum_{n=1}^{N}\left<x_{jn}^{(t)}\right>\left(x_{n}-\mu_{j}^{(t)}\right)\left(x_{n}-\sum_{j}^{(t)}\right)^{\top}}{\sum_{n=1}^{N}\left<z_{jn}^{(t)}\right>}\)

GMM 訓練のための上記の更新式は非常に簡潔で実装が簡単である。 これらは，EM の採用が尤度最大化問題の解決をどのように促進するかについての注目すべき例である。

上記のアプローチで起こりうる問題は，図 6 に示す例のように，共分散行列が特異になる可能性があるという事実に関連している。 この図は，2 次元 (2-D) データ集合に対して，EM を適用して 20 個の成分を持つ GMM を訓練したときに得られた解の等高線プロットである。 GMM 成分の一部が特異であること，すなわち，それらの密度がデータ点の周りに集中し，いくつかの主軸に沿った分散が 0 になる傾向があることは明らかである。 GMM 学習のための典型的な ML アプローチのもう 1 つの欠点は，モデル選択，すなわち成分数の決定に使えないことである。 これらの問題に対する解決策は，Bayes GMM を使うことで得られるかもしれない。

7.2 変分 GMM 訓練

7.2.1 完全 Bayesian GMM

$X=\left{x_{n}\right}$ を N 個の観察集合とし，各 $x_{n}\in\mathcal{R}^{d}$ を特徴ベクトルとする。 また，$p(x)$ を $M$ 個のガウス成分を持つ混合とする \(\tag{71} p(x)=\sum_{j=1}^{M}\pi_{j} N(x;\mu_{j},T_{j})\) ここで $\pi=\left{\pi_{j}\right}$ は混合係数 (重み)，$\mu=\left{\mu_j\right}$ は成分の平均 (中心)，$T=\left{T_j\right}$ は精度 (逆共分散) 行列である (Bayes GMM では共分散行列の代わりに精度行列を使用する方が便利であることに注意)。

Bayes 型 Gauss 混合モデルは，モデルパラメータ $\pi,\mu,T$ に事前分布を与えることで得られる。 典型的には共役事前分布が用いられ，$\pi$ に対しては Dirichlet, $(\mu,T)$ に対しては Gauss-Wishart 分布である。 パラメータ $\left{\alpha_{j}\right}$ を持つ $\pi$ の Dirichlet 事前分布は次式で与えられる。 \(\text{Dir}\left(\pi|\alpha_1,\ldots,\alpha_{M}\right)=\frac{\Gamma\left(\sum_{j=1}^{M}\alpha_{j}\right)}{\prod_{j=1}^{M}\gamma\left(\alpha_{j}\right)}\prod_{j=1}^{M}\pi_{j}^{\alpha_{j}-1},\) ここで $\Gamma(x)$ はガンマ関数を表す。 通常，$\alpha_j$ はすべて等しい，すなわち $\alpha_j=\alpha_0, j=1,\ldots,M$.

$(\mu,T)$ の Gauss-Wishart 事前分布は $p(\mu,T)=\prod_{j=1}^{M}p(\mu_j,T_j)=p(\mu_j |T_j)p(T_j)$ であり，$p(\mu_j|T_j=N(\mu_j;\mu_0,\beta_0T_j)$ (パラメータ $\mu_0$ と $\beta_0$) であり，$p(T_j )$ は Wishart 分布である。 \(W(T_j|\nu,V)=\frac{\left|T_j\right|^{(v-d-1)/2}\exp\text{tr}\left\{-\frac{1}{2}VT_j\right\}}{2^{vd/2}\pi^{d(d-1)/4}\left|V\right|^{-n/2}\prod_{i=1}^{d}\Gamma((v+1-i)/2)},\) パラメータ $\nu$ と $V$ はそれぞれ自由度とスケール行列を表す。 Wishart 分布は，Gamma 分布の多次元一般化であることに注意。 線形回帰では，Gauss-Gamma 事前分布を用い，独立の事前分布 $\alpha_i$ を仮定し，したがって，それらに独立の Gamma 事前分布を割り当てる。 しかし，ここでは，データ間に有意な相関があるかもしれないので，これらの相関を捕捉するために Wishart 事前分布を使うことができる。

この Bayes GMM に対応するグラフィカルモデルを 図 7(a) に示す。 これは完全 Bayes GMM であり，すべてのハイパーパラメータ (すなわち，事前分布のパラメータ $\alpha,\mu_0,\beta_0,\nu,V$) が事前に指定されている場合，モデルは推定されるべきパラメータを含まず，データが与えられた事後値 $p(h|x)$ を計算しなければならない隠れ確率変数 $h=(Z,\pi,\mu,T)$ だけを含む。 この場合，事後値が解析的に計算できないことは明らかであり，したがって，変分平均場 (16) を特定の Bayes モデルに適用することによって近似 $q(h)$ が計算される[27]。

平均場近似では，$q$ は次の形式の積であると仮定する。 \(\tag{72} q(h)=q_{z}(Z) q_{\pi}(\pi) q_{\mu T}(\mu, T).\) となり，解は (16) で与えられる。 必要な計算を行った結果，次のような密度の集合が得られる： \(\tag{73} q_z(Z)=\prod_{n=1}^{N}\prod_{j=1}^{M}r_{jn}^{z_{jn}}\)

\[\tag{74} q_{\pi}(\pi)=\text{Dir}\left(\pi|\left\{\lambda_{j}\right\}\right)\] \[\tag{75} q_{\mu T}(\mu,T)=\prod_{j=1}^{M}q_{\mu}(q_{\mu}\left(\mu_j|T_{j}\right) q_{T}\left(T_{j}\right)\] \[\tag{76} q_{\mu}(\mu_j|T)=\prod_{j=1}^{M}N\left(\mu_j;m_j,\beta_jT_j\right)\]

\(\tag{77} q_{T}=\prod_{j=1}^{M}W\left(T_j;n_j,U_j\right)\) そして，密度のパラメータ $(r_{jn},\lambda_j,m_j,\beta_j,\nu_j,U_j)$ を更新するための詳細な公式は [27] にある。 簡単な反復更新手順を用いて上記の方程式系を解くことにより，平均場制約の下での真の事後 $p(h|x)$ に対する最適な近似 $q(h)$ が得られる。

Bayes モデリングにおける典型的なアプローチは，モデルのハイパーパラメータ $\alpha,\nu,V,\mu_0,\beta_0$ を指定し，情報量のない事前分布が定義されるようにすることである。 これらのパラメータを調整するために，アルゴリズムに M ステップを組み込むことは可能であるが，我々はこのアプローチに従う。 しかし，これは通常行われない。

完全 Bayes GMM の利点の 1 つは，事前分布を持たない GMM と比較して，Gauss 成分が 1 つのデータ点を担当するようになる ML アプローチでしばしば生じる特異解を許さないことである。 2 番目の利点は，交差検証法などの手法に頼ることなく，最適な成分数を直接決定するためにベイズ GMM を使用できることである。 しかしながら，Direchlet 事前分布は，成分の混合重みがゼロになって混合から除去されることを許さないので，この問題では完全 Bayes 混合分布の有効性は制限される。 また，この場合，最終的な結果は，事前に指定されなければならない事前分布のハイパーパラメータ (特に Direchlet 事前分布のパラメータ) に大きく依存する[13]。 特定のハイパーパラメータの集合に対して，混合成分の数 $M$ のいくつかの値に対して変分アルゴリズムを実行し，変分下界の最良の値に対応する解を保持することが可能である。

7.2.2 混合重みからの事前分布除去

[28] では，Bayes GMM モデルのもう 1 つの例が提案されており，これは混合重み $\left{\pi_j\right}$ の事前分布を仮定せず，混合重みは確率変数ではなくパラメータとして扱われる。 このアプローチのグラフモデルは図 7(b) に描かれている。

このベイズ GMM (本研究の 2 人の著者の頭文字をとって CB モデルと呼ぶ) では，$\mu$ と $T$ に対してそれぞれ Gauss 事前分布と Wishart 事前分布が仮定される

\[\tag{78} p(\mu)=\prod_{j=1}^{M}N(\mu_j|0,\beta\mathbf{I})\] \[\tag{79} p(T)=\prod_{j=1}^{M}W(T_j|\nu,\mathbf{V})\]

この Bayes モデルは，(ある程度) 最適な成分数を推定することができる。 これは，隠れ変数 $h=\left{Z,\mu,T\right}$ を積分することによって得られる周辺尤度 $p(X;\pi)$ の最大化によって達成される。 \(\tag{80} p(X;\pi)=\int p\left(X,h;\pi\right)\,dh\) パラメータとして扱う混合重み $\pi$ に関して。 変分近似は，対数周辺尤度 \(\tag{81}\) の下界の最大化を提案する。 ここで $q(h)$ は事後 $p(h|X)$ を近似する任意の分布である。 注目すべき性質は，$F$ を最大化するとき，いくつかの成分がデータ空間の同じ領域に入る場合，この領域のデータが少ない成分で十分に説明できるようになると，モデルには冗長な成分を排除する (つまり，それらの $\pi_j$ を 0 にする) 強い傾向があることである。 その結果，混合成分間の競合は，モデル選択問題に対処するための自然なアプローチを示唆する： 多数の成分で初期化された混合を適合させ，競合によって冗長な成分を除去させる。

変分法に従って，我々の目的は対数周辺尤度 $\log p(X;\pi)$ の下界 $F$ を最大化することである。 \(\tag{82} F\left[q,\pi\right]=\sum_{z}\int q(Z,\mu,T)\log\frac{p(X,Z,\mu,T)}{q(Z,\mu,T)}du\,dT.\) ここで $q$ は事後分布 $p(Z,\mu,T|X;\pi)$ を近似する任意の分布である。 $F の最大化は，変分 EM アルゴリズムを用いて反復的に実行される。 各反復で 2 つのステップが行われる：まず $q$ に関する境界の最大化，次に $\pi$ に関する境界の最大化。

$q$ に関する最大化を実行するために，$q$ が次の形式の積であると仮定する平均場近似が採用されている(14)。 \(q(h)=q_z(Z)\,q_\mu(\mu)\,q_T(T).\tag{83}\). (16)で必要な計算を行った結果，以下のような密度が得られた： \(q_z(Z)=\prod_{n=1}^{N}\prod_{j=1}^{M}r_{jn}^{z_{jn}}\tag{84}\) \(q_\mu(\mu)=\prod_{j=1}^{M}N(\mu_j|m_j,S_j)\tag{85}\) \(q_T(T)=\prod_{j=1}^{M}W(T_j|n_j,U_j)\tag{86}\) ここで，密度のパラメータは次のように計算できる。 \(\begin{aligned} r_{jn} &=\frac{\tilde{r}_{jn}}{\sum_{k=1}^{M}\tilde{r}_{jn}} &&& (tz:87)\\ \tilde{r}_{jn} &=\pi_{j}\exp\left\{\frac{1}{2}\left<\log\left|T_j\right|\right>-\frac{1}{2}\text{tr}\left\{\left<T_j\right>\left(x_nx_n^{\top}-x_n\left<u_j\right>^{\top}+ \left<u_j\right>x_n^{\top}+\left<\mu_j\mu_j^{\top}\right>\right)\right\} \right\} &&& (tz:88)\\ m_J &= S_j^{-1}\left<T_j\right>\sum_{n=1}^{N}\left<z_{jn}\right>x_n &&& (tz:89)\\ S_j & =\beta\mathbf{I}+\left<T_j\right>\sum_{n=1}^{N}\left<z_{jn}\right> &&& (tz:90)\\ n_J & =\nu+\sum_{n=1}^{N}\left<z_{jn}\right> &&& (tz:91)\\ U_j&=V+\sum_{n=1}^{N}\left<z_{jn}\right>\left(x_nx_n^{\top}-x_n\left<u_j\right>^{\top}+\left<u_j\right>x_n^{\top}+\left<u_ju_j^{\top}\right>\right)&&&(tz:92)\\ \end{aligned}\)

上式で用いた $q(h)$ に関する期待値は，次式を満たす：   $\left<T_j\right>=\eta_jU^{-1}j,\left<\log\left|T_j\right|\right>=\sum{i=1}^{d}\Psi(0.5(\eta_j+1-i)) +d\ln 2 - \ln\left|U_j\right|,\left<\mu_j\right>=m_j,\left<\mu_j\mu^{\top}j\right>=S^{-1}{j}+m_jm^{\top}j$ [ここで，$\psi$ は次式で定義されるディガンマ関数を示す。 $d/dx\ln\Gamma(x)=\Gamma^{\prime}(x)/\Gamma(x)$] かつ，$z{jn}=r_{jn}$。 密度が期待値を通して結合していることがわかるので，パラメータの反復推定が必要である。 しかし実際には，変分 E ステップでは 1 回のパスで十分である。

$q$ に関して $F$ を最大化した後，訓練法の各反復の第 2 ステップでは，$\pi$ に関して $F$ を最大化する必要があり，次のような変分 M ステップの簡単な更新式が導かれる： \(\tag{93} \pi_j=\frac{\sum_{n=1}^{N}r_{jn}}{\sum_{k=1}^{M}\sum_{n=1}^{N}r_{kn}}\) 上記の変分 EM 更新方程式は反復的に適用され，変分境界の局所最大値に収束する。 最適化中に混合係数の一部がゼロに収束し，対応する成分が混合物から除去される。 このようにして複雑さの制御が達成される。 これは，$\mu$ と $T$ に関する事前分布が重複する成分に罰則を課しているからである。 定性的には，変分境界は 2 つの項の和として書くことができる。 1 つは尤度項 (データ適合の質に依存する) で，もう 1 つは複雑なモデルに罰則を与える事前分布による罰則項である。

図 8 は，すでに図 6 で示した 2 次元データセットを用いたこの手法の性能の例示である。 この手法は 20 成分から始まり，反復回数が増えるにつれて，成分数は徐々に減少し(いくつかの $\pi_j$ はゼロになる)，最終的にこのデータセットに対する良好な GMM モデルが達成される。 また，共分散行列に事前分布が存在することで，図 6 に示された事前分布なしの GMM 解とは対照的に，特異解に到達しないことが観察される。

一般に，CB は，成分がよく分離されている場合に良好な性能を示す効果的な手法である。 しかし，その性能は，精度行列に課される Wishart 事前分布のスケール行列 V の指定に敏感である。 上記の混合モデルを構築するための漸進的手法が [29] で提案されている。 各ステップにおいて，学習は特定の混合成分jが占めるデータ領域に制限されるため，精度行列 $T_j$ に基づいて局所的な精度事前分布を指定することができる。 この動作を実現するために，図 7 の生成モデルに対して，成分の下位集合でのみ競合を制限する修正を加えた。

8. まとめ

EM アルゴリズムは，ML 推定に多くの利点をもたらす反復法である。 尤度関数の直接最適化が困難な問題に対して，局所収束が保証された単純な反復解を提供する。 多くの場合，共分散行列が正定値，確率ベクトルが正で和が 1 など，推定パラメータに対するいくつかの制約を満たす解を提供する。 さらに，EM の適用には尤度関数を明示的に評価する必要はない。

しかし EM アルゴリズムを適用するためには，観察が与えられたときの隠れ変数の事後知識が必要である。 これは，複雑なベイズモデルには EM を適用できないので重大な欠点である。 しかし，複雑なベイズモデルは，適切に構築されれば，データ生成機構の顕著な特性をモデル化する能力を持ち，困難な問題に対して非常に良い解を提供するため，非常に有用である。

変分法は EM アルゴリズムのこの欠点を改善するために，信号処理界で人気を集めている反復アプローチである。 この方法論によれば，観測値が与えられたときの隠れ変数の事後的近似値が使用される。 この近似に基づき，尤度関数の下界を最大化することでベイズ推定が可能となり，局所収束が保証される。 この方法論は，複雑なグラフモデルの推論を可能にし，ある場合には EM で解ける単純なモデルと比較して，著しい改善をもたらす。

この問題は，信号処理アプリケーションの 2 つの基本的な問題である線形回帰とガウス混合モデリングの文脈で，本論文で実証された。 具体的には，線形回帰の文脈では，変分法によって解かれた複雑なベイズモデルが，局所的な信号特性をよりよく捉え，信号の不連続領域でのリンギング ringing を回避できることを実証した。 ガウス混合モデリングの文脈では，変分法によって解かれたモデルは，特異点を回避し，モデルの構成要素の数を推定することができた。 これらの結果は，変分法が信号処理アプリケーションを長い間悩ませてきた難問に解を提供する力を持つことを示すものである。 この方法の主な欠点は (少なくとも今のところ) 使用される境界の厳密さを評価できる結果がないことである。

B.1 ELBO: Evidence Lower BOund

from https://www.cs.princeton.edu/archive/fall11/cos597C/lectures/variational-inference-i.pdf

変分推論では潜在変数空間 $\mathcal{Z}$ 上での $z$ の確率密度を考える。 ここでの目標は，KL ダイバージェンスの意味で最良の以下の式を最適化することである:

\[q^{\star}(z)=\operatorname{argmin}_{g(x)\in\mathcal{Z}} \text{KL}(q(z)\vert\vert p(z\vert x))\tag{eq:blei10}.\]

この解が得られれば，$q^{\star}(\cdot)$ は，条件付近似となる。

変分推論では，潜在変数に対する密度の系列 $\mathcal{Z}$ を指定する。 各 $q(z\in\mathcal{Z})$ は，正確な条件に近似する候補である。 我々の目的は，最良の候補を見つけることであり，正確な条件に KL ダイバージェンスで最も近いものを見つけることである。 推論は，以下の最適化問題を解くことになる。

\[q^{\star}(z)=\arg\min_ {g(x)\in\mathcal{Z}}\text{KL}(q(z)\vert\vert p(z\vert x))\tag{eq:blei11}.\]

見つかった $q^{\star}(\cdot)$ は，$\mathcal{Z}$ 族の中で，条件式の最適な近似値となる。 分布族の複雑さは，この最適化の複雑さを決定する。

しかし，この目的は $p(x)=\int p(z,x)\,dz$ という式で $\log p(x)$ の証拠を計算する必要があるため，計算できない。 (証拠を計算するのが難しいからこそ，そもそも近似推論に訴えるのである)。 なぜかというと KL ダイバージェンスを思い出せば:

\[\text{KL}(q(z)\vert\vert p(z\vert x)) =\mathbb{E}(\log q(z)) - \mathbb{E}(\log p(z\vert x))\]

ここで，すべての期待値は $q(z)$ を基準にしています。 ここで，条件式を展開する。

\[\text{KL}(q(z)\vert\vert p(z\vert x)) =\mathbb{E}(\log p(z)) - \mathbb{E}(\log p(z,x)) +\log p(x).\tag{eq:blei12}\]

これにより $\log p(x)$ への依存性が明らかになった。

KL を計算することができないので，定数を追加した上で KL と同等の代替目的を最適化する。 \(\text{ELBO}(q)=\mathbb{E}(\log p(z,x)) - \mathbb{E}(\log p(z)).\tag{eq:blei13}\)

この関数は evidence lower bound (ELBO) と呼ばれる。 ELBO は 式(\ref{eq:blei12}) の 負の KL ダイバージェンス に $q(z)$ に対する定数である $\log p(x)$ を加えたものである。 ELBO を最大化することは KL ダイバージェンスを最小化することと同等である。

ELBO を調べることで， 最適な変分密度を直感的に理解することができる。 ELBO は，データの期待対数尤度と，事前 の $p(z)$ と $q(z)$ の間の KL ダイバージェンス の和と書き換えられる。 \(\begin{aligned} \text{ELBO}(q) &= \mathbb{E}(\log p(z)) + \mathbb{E}(\log p(x\vert z)) - \mathbb{E}(\log q(z))\\ &= \mathbb{E}(\log p(x\vert z)) -\text{KL}(q(z)\vert\vert p(z)).\\ \end{aligned}\)

この目的語は $q(z)$ がどのような値に質量を置くことを促すのだろうか？ 第 1 項は，期待尤度で， 観測データを説明する潜在変数の構成に質量を置く密度を促す。 第 2 項は，変分密度と事前分布の間の負の発散で，事前分布に近い密度を奨励する。 このように，変分目標は，尤度と事前分布の通常のバランスを反映している。

ELBO のもう一つの特性は，任意の $q(z)$ に対して (対数) 証拠 である $\log p(x)\ge\text{ELBO}(q)$ を下界にすることである。 これが ELBO の名前の由来である。 これを見るには 式(eq:blei12) と(eq:blei13) が次のような証拠の表現を与えることに注意。 \(\log p(x)=\text{KL}(q(z)\|p(z|x)) + \text{ELBO}(q).\tag{eq:blei14}\)

そして，この境界は $\text{KL}(\cdot)\ge0$ [@1951KullbackLeibler] という事実から導かれる。 変分推論のオリジナル文献では，これは Jensen’s inequality [@Jordan1999] によって導かれた。

\[\lambda^{\star}=\operatorname{argmin}_{\lambda} D_{\text{KL}}{q(h|e,\lambda)\|p(h\vert e}),\tag{eq:jordan41}\]

ここで、任意の確率分布 $q(s)$ と $p(s)$ に対して，カルバック・ライブラーダイバージェンスは以下のように定義される:

\[D_ {\text{KL}} (q\vert\vert p) = \sum_ {s} q(s)\log\frac{q(s)}{p(s)}.\tag{eq:jordan42}\]

変分パラメータ $\lambda^{\star}$ の最小値は，特定の分布 $q(h\vert e,\lambda^{\star})$ を定義し， これを $q(h\vert e,\lambda)$ 族における $p(h\vert e)$ の最良の近似値として扱う。

KL ダイバージェンスを近似精度の指標として用いる単純な理由の一つは，KL ダイバージェンスが近似 $q(h\vert e,\lambda)$ 族における証拠 $p(e)$ の確率 (すなわち尤度) の最良の 下界 をもたらすからである。 実際 $p(e)$ の対数を イェンセンの不等式 を用いて以下のように束縛する。

\[\begin{aligned} \log p(e) &= \log\sum_ {H} p(h,e) \\ &= \log\sum_ {h} q(h\vert e) \frac{p(h,e)}{q(h\vert e)}\\ &\ge \sum_{h} q(h\vert e) \log\left[\frac{p(h,e)}{q(h\vert e)}\right].\\ \end{aligned}\tag{eq:jordan43}\]

ELBO と $\log p(x)$ の関係から，モデル選択基準として変分境界を用いることになった。 これは，混合モデル (Ueda and Ghahramani, 2002; McGrory and Titterington, 2007) やより一般的なモデル (Beal and Ghahramani, 2003) で検討されている。 その前提は，境界が周辺尤度の良い近似であり，モデルを選択するための根拠となることである。 これは，実際にはうまくいくこともあるが，境界に基づいて選択することは，理論的には正当化されない。 他の研究では，対数予測密度の変分近似を用いて，交差妥当性に基づくモデル選択に VI を使用している (Nott et al.2012)。

最後に，多くの読者は 式(eq:blei13) の ELBO の第 1 項が EM アルゴリズム [@Dempster1977] によって最適化された期待完全対数尤度であることに気づくだろう。 EM アルゴリズムは潜在変数を持つモデルの最尤推定値を求めるために設計されたものである。 このアルゴリズムは，ELBO が対数尤度 $\log p(x)$ (すなわち 対数尤度 $\log p(x)$) に等しいという事実を利用している。 ($q(z)=p(z\vert x)$ のとき ELBO は対数尤度 $\log p(x)$ (すなわち対数証拠) に等しいという事実を利用している。 EM は $p(z\vert x0$ に従った期待される完全な対数尤度の計算 (E ステップ) と，モデルパラメータに関する最適化（M ステップ) を交互に行う。 変分推論とは異なり，EM は $p(z\vert x)$ の下での期待値が計算可能であることを前提とし，それ以外の難しいパラメータ推定問題に使用する。 EM とは異なり，変分推論は固定モデルのパラメータを推定するものではなく，古典的なパラメータが潜在変数として扱われるベイズの設定でよく使われる。 変分推論は，潜在変数の正確な条件を計算できないモデルに適用される。

C. Probabilistic Machine Learning

• from 2016Blei NIPS, slide 13

• 確率モデルは，潜在変数 $z$ と観測変数 $x$ の結合分布 $p(z,x)$ である。
• 未知数についての推論は，観測値が与えられたときの潜在変数の条件付き分布である 事後分布 を通して行われる。

\[p(z\vert x)=\frac{p(z,x)}{p(x)}.\]

• ほとんどの興味深いモデルにおいて，分母は扱いにくい。

近似事後推論 に訴える。

Schematic graph of VI

• VIは 推論を最適化に変える。
• 潜在変数に対する分布の変数族を仮定する。

\[q(z;\nu)\]

• 変動パラメータ $\nu$ を正確な事後処理に (KL で) 近づけるように適合させる。 (EP, BP などのアルゴリズムにつながる代替ダイバージェンスもある)。

C.1 歴史

• 変分推論は，統計物理学の考え方を確率推論に応用したものである。 おそらく 80年代後半に Peterson と Anderson (1987) が平均場近似法を用いてニューラルネットワークを適合させたのが始まりと思われる。
• このアイデアは 1990 年代初頭に Jordan 研究室 (Tommi Jaakkola, Lawrence Saul, Zoubin Gharamani) によって取り上げられ，多くの確率的モデルに一般化された (総説論文は Jordan et al.1999)。
• 並行して Hinton and Van Camp (1993) はニューラルネットワークのための平均ば近似を開発した。 Neal and Hinton (1993) はこのアイデアを EM アルゴリズムに結びつけ，さらに混合エキスパートモデルに対する変分法 (Waterhouse et al., 1996) や HMM (MacKay,1997) につなげた。

C.2 今日

• 現在，変分推論に関する新しい研究が盛んに行われており，スケーラブルにし，導出を容易にし，より速く，より正確にし，より複雑なモデルやアプリケーションに適用している。
• 現代の変分推論は，確率的プログラミング，強化学習，ニューラルネットワーク，凸最適化，ベイズ統計学，そして無数のアプリケーションなど，多くの重要な領域に関連している。
• 今日の我々の目標は，基本を学び，新しいアイデアを説明し，新しい研究のオープンな領域を提案することである。