はじめに

お詫び：はてなブログのはてなTeXの仕組みのせいで、一部数式が正しく表示できない問題が発生しています...色々試行錯誤して見たものの、修正ができませんでした。読みにくい部分が一部あります...
　　
ベイズ統計学の世界では、複雑なモデルの事後分布を正確に計算することが長年の課題でした。高次元パラメータ空間における積分計算は解析的に解くことができず、統計的推論の大きな障壁となっています。この困難に立ち向かうために、様々な近似手法が開発されてきましたが、その中でも特に注目されているのが「変分ベイズ法」です。
変分ベイズ法は、事後分布を直接計算する代わりに、取り扱いやすい分布族の中から最も事後分布に近い分布を見つけるというアプローチをとります。このアイデアは物理学の変分原理に由来し、計算効率と近似精度のバランスを追求する方法として機械学習や統計学で広く応用されています。
しかし、変分ベイズ法について「KLダイバージェンスを最小化すればいい」「平均場近似を使えばよい」といった断片的な知識は広まっているものの、その深い理論的背景や実践上の注意点、そしてMCMCとの比較における長所・短所について、包括的に理解している人は意外と少ないのではないでしょうか。
本稿では、変分ベイズ法の基礎理論から実装手法、MCMCとの比較、そして結果をMCMCに近づけるための最新テクニックまで、徹底的に解説します。複雑な数式も直感的な例えと共に説明することで、初学者から実務者まで、すべての読者に深い理解をもたらすことを目指します。
変分ベイズ法の世界に足を踏み入れ、その美しさと実用性を一緒に探求していきましょう。

1. ベイズ推論と計算の課題

1.1 ベイズ推論の基本

ベイズ推論は、データから学習する際の不確実性を確率的に表現し、新しい証拠が得られるたびに信念を更新する体系的なフレームワークです。その中心にあるのが、18世紀の数学者トーマス・ベイズにちなんで名付けられたベイズの定理です：

$P(\theta|D) = \frac{P(D|\theta)P(\theta)}{P(D)}$

この式の各要素には明確な解釈があります：

$P(\theta|D)$ は事後分布（posterior distribution）と呼ばれ、データ $D$ を観測した後のパラメータ $\theta$ に関する確率分布を表します。これがベイズ推論の最終的な目標物です。

$P(D|\theta)$ は尤度（likelihood）と呼ばれ、パラメータ $\theta$ が与えられたときにデータ $D$ が観測される確率を表します。この項はデータがモデルにどれだけ適合しているかを示します。

$P(\theta)$ は事前分布（prior distribution）と呼ばれ、データを観測する前のパラメータ $\theta$ に関する確率分布、つまり事前知識や信念を表します。

$P(D)$ は周辺尤度（marginal likelihood）または証拠（evidence）と呼ばれ、すべての可能なパラメータ値について尤度と事前分布を積分したものです： $P(D) = \int P(D|\theta)P(\theta)d\theta$ 。これは正規化定数としての役割を果たします。

ベイズ推論のプロセスを直感的に説明すると、「事前の信念（事前分布）を持って観測を始め、データ（尤度）を通じて学習することで、より洗練された信念（事後分布）に更新する」ということになります。例えば、コインの公平さについての信念を考えてみましょう。最初は公平だと思っていても（事前分布）、何度も表が出続ければ（データと尤度）、そのコインは偏りがあるという信念（事後分布）に更新されます。
この枠組みは、「学習」と「不確実性の定量化」を自然に統合できる点で非常に強力です。パラメータの点推定だけでなく、その不確実性の程度も含めた完全な確率分布として推論結果を得られるのです。

1.2 計算の課題：積分の壁

ベイズ推論の理論的な美しさとは裏腹に、実際の計算は非常に難しい場合が多いです。その最大の原因は、事後分布を求める際に必要となる周辺尤度（エビデンス）の計算です：

$P(D) = \int P(D|\theta)P(\theta)d\theta$

この積分は、多くの実用的なモデルでは解析的に計算できません。特に、パラメータ空間の次元が高い場合や、尤度関数と事前分布の組み合わせが複雑な場合には、この積分は「計算不可能」な領域に入ってしまいます。これは「積分の壁」とも呼ばれる問題です。
例えば、ごく単純なケースとして、ガウス分布の平均と精度（分散の逆数）をベイズ推論する場合を考えましょう。事前分布としてガウス分布とガンマ分布を使用する「共役事前分布」を選べば、事後分布は解析的に求められます。しかし、ニューラルネットワークのような数万、数百万のパラメータを持つモデルでは、各パラメータの複雑な相互作用により、解析的な解を得ることはほぼ不可能になります。
この計算上の課題を解決するために、大きく分けて二つのアプローチが発展してきました：

サンプリング法：事後分布からサンプルを生成し、そのサンプルを用いて事後分布の特性（期待値や分散など）を推定するアプローチ。最も代表的なのがマルコフ連鎖モンテカルロ法（MCMC）です。
近似法：事後分布そのものを、扱いやすい分布族で近似するアプローチ。変分ベイズ法（VB）がこれに当たります。

これらのアプローチはそれぞれ長所と短所を持ちます。MCMCは理論上、十分な時間をかければ真の事後分布に収束するという保証がありますが、計算コストが高く、特に大規模データや複雑なモデルでは実用的でない場合があります。一方、変分ベイズ法は計算効率が良い反面、選択した分布族による表現力の制限から、真の事後分布を完全に捉えられない可能性があります。
これから詳しく見ていくように、変分ベイズ法の本質は「計算不可能な事後分布を、計算可能な分布で最もよく近似する」という点にあります。それは「精密な解を得ることを諦め、近似解を効率的に求める」という妥協の産物ですが、多くの実用的な問題では、この妥協が非常に価値のある選択肢となるのです。

2. 変分ベイズ法の基礎理論

2.1 変分ベイズ法の基本的考え方

変分ベイズ法の根本的なアイデアは、「計算が困難な事後分布 $p(\theta|D)$ を、扱いやすい近似分布 $q(\theta)$ で代用しよう」というものです。これは「最適な解を求めるのが難しい問題を、制約された条件の中での最適化問題に置き換える」という変分法の一般的な考え方に基づいています。

具体的には、パラメータ空間上で定義される真の事後分布 $p(\theta|D)$ を、ある関数族（変分族と呼ばれる） $\mathcal{Q}$ の中から選ばれた分布 $q(\theta) \in \mathcal{Q}$ で近似します。

この近似の「良さ」を測る指標として、変分ベイズ法では通常、カルバック・ライブラー（KL）ダイバージェンスが使用されます：

$KL(q(\theta) \| p(\theta|D)) = \int q(\theta) \log \frac{q(\theta)}{p(\theta|D)} d\theta$

KLダイバージェンスは二つの確率分布間の「距離」を測る尺度です。ただし、通常の距離とは異なり非対称であり、 $KL(p | q) \neq KL(q | p)$ です。KLダイバージェンスは常に非負で、二つの分布が同じときのみゼロになります。

変分ベイズ法の目標は、このKLダイバージェンスを最小化する $q(\theta)$ を見つけることです：

$q^*(\theta) = \arg\min_{q \in \mathcal{Q}} KL(q(\theta) \| p(\theta|D))$

これを概念的に説明すると、変分ベイズ法は「複雑な地形図（事後分布）を、単純な形状（近似分布）でどれだけうまく近似できるか」という問題に取り組んでいるわけです。例えば、複雑な山岳地帯を、いくつかのなだらかな丘の組み合わせで表現しようとしているイメージです。

ここで重要なのは選択する関数族 $\mathcal{Q}$ です。この選択によって、計算の難しさと近似の精度のトレードオフが決まります。例えば：

ガウス分布族：計算が比較的容易で、単峰性の事後分布に対しては良い近似を与えますが、多峰性や歪んだ分布の表現能力は限られます。
より複雑な分布族（混合ガウス分布など）：より広範な形状の事後分布を近似できますが、計算コストが増加します。

変分ベイズ法の美しさは、「複雑な積分計算」を「最適化問題」に置き換えることで、計算不可能に見えた問題を扱える形に変換する点にあります。しかし、選んだ関数族に真の事後分布が含まれていない場合、完全な近似は原理的に不可能であることに注意が必要です。これが変分ベイズ法の根本的な限界でもあります。

2.2 変分下界（ELBO）の導出

KLダイバージェンス $KL(q(\theta) | p(\theta|D))$ を直接最小化するのは依然として困難です。なぜなら、真の事後分布 $p(\theta|D)$ を計算するには、最初からの課題であった周辺尤度 $p(D)$ が必要だからです。この問題を回避するために、変分ベイズ法では「変分下界」（Evidence Lower Bound、ELBO）と呼ばれる等価な目的関数を最大化します。

ELBOの導出は以下の通りです。まず、周辺尤度（エビデンス）の対数をとり、観察します：

$\log p(D) = \log \int p(D, \theta) d\theta$

ここに変分分布 $q(\theta)$ を導入します：

$\log p(D) = \log \int q(\theta) \frac{p(D, \theta)}{q(\theta)} d\theta$

ここで、ジェンセンの不等式を適用します。ジェンセンの不等式は、凹関数（対数関数はその一例）と確率分布に対して、 $f(\mathbb{E}[X]) \geq \mathbb{E}[f(X)]$ が成り立つという定理です。この不等式を適用すると：

$\log p(D) = \log \mathbb{E}_{q(\theta)}\left[\frac{p(D, \theta)}{q(\theta)}\right] \geq \mathbb{E}_{q(\theta)}\left[\log \frac{p(D, \theta)}{q(\theta)}\right]$

右辺を展開すると：

$\mathbb{E}_{q(\theta)}\left[\log \frac{p(D, \theta)}{q(\theta)}\right] = \int q(\theta) \log \frac{p(D, \theta)}{q(\theta)} d\theta \$

$= \int q(\theta) \log p(D, \theta) d\theta - \int q(\theta) \log q(\theta) d\theta$

$= \mathbb{E}_{q(\theta)}[\log p(D, \theta)] + H[q(\theta)]$

ここで $H[q(\theta)] = -\int q(\theta) \log q(\theta) d\theta$ は $q(\theta)$ のエントロピーを表します。この右辺が変分下界（ELBO）です：

$\text{ELBO}(q) = \mathbb{E}_{q(\theta)}[\log p(D, \theta)] + H[q(\theta)]$

さらに、 $p(D, \theta) = p(D|\theta) p(\theta)$ を使うと：

$\text{ELBO}(q) = \mathbb{E}_{q(\theta)}[\log p(D|\theta)] + \mathbb{E}_{q(\theta)}[\log p(\theta)] + H[q(\theta)]$

この式から、ELBOは以下の3つの要素から構成されていることがわかります：

尤度の期待値 $\mathbb{E}_{q(\theta)}[\log p(D|\theta)$ ]：データがモデルにどれだけ適合しているかを表す項
事前分布と変分分布の類似度 $\mathbb{E}_{q(\theta)}[\log p(\theta)$ ]：変分分布が事前分布からどれだけ離れているかを表す項
変分分布のエントロピー $H[q(\theta)$ ]：変分分布の不確実性や広がりを表す項

さらに、KLダイバージェンスとELBOの関係を導くと：

$\log p(D) - \text{ELBO}(q) = KL(q(\theta) \| p(\theta|D))$

つまり、ELBOを最大化することは、KLダイバージェンスを最小化することと等価です。周辺尤度 $\log p(D)$ は $q(\theta)$ に依存しない定数なので、 $\text{ELBO}(q)$ を最大化することで $KL(q(\theta) | p(\theta|D))$ が最小化されるのです。

概念的に説明すると、ELBOは「データの説明のしやすさ」と「近似分布の単純さ」のトレードオフを表しています。第一項と第二項は「データと事前分布をどれだけ説明できるか」を表し、第三項は「近似分布がどれだけ幅広く（不確実性を保持して）柔軟か」を表します。これは「データに適合しながらも、過度に自信を持ちすぎない（不確実性を保持する）」バランスを取ろうとしていると解釈できます。

2.3 平均場近似

変分ベイズ法において、近似分布 $q(\theta)$ の構造をどう選ぶかは重要な問題です。最も一般的なアプローチの一つが「平均場近似」（Mean Field Approximation）です。

平均場近似の核心は、パラメータベクトル $\theta$ を複数の成分 $\theta = (\theta_1, \theta_2, \ldots, \theta_M)$ に分解し、これらの成分が互いに独立であると仮定することです：

$q(\theta) = \prod_{i=1}^{M} q_i(\theta_i)$

この仮定により、最適化問題が大幅に簡単になります。各成分の最適な分布 $q_i^*(\theta_i)$ は以下の式で与えられます：

$q_i^*(\theta_i) \propto \exp\left( \mathbb{E}_{j \neq i}\left[\log p(D, \theta)\right] \right)$

ここで $\mathbb{E}_{j \neq i}$ は $\theta_i$ 以外の全てのパラメータに関する期待値を表します。この期待値は、他のすべての成分 $\theta_j \quad (j \neq i)$ に関して計算されます。

この更新式は、座標降下法（coordinate descent）のように、各成分を順番に最適化していくことを示唆しています。つまり、一つの成分について最適化を行い、次に別の成分について最適化を行う、というプロセスを収束するまで繰り返します。

平均場近似の導出を詳しく見てみましょう。まず、ELBOを $\theta_i$ に関して表現します：

$\text{ELBO}(q) = \int \prod_{j=1}^M q_j(\theta_j) \log p(D, \theta) \, d\theta - \sum_{j=1}^M \int q_j(\theta_j) \log q_j(\theta_j) \, d\theta_j$

これを $\theta_i$ について最適化するために、他のすべての $\theta_j \quad (j \neq i)$ を固定して考えます。ラグランジュ乗数法を使い、 $\int q_i(\theta_i) d\theta_i = 1$ という制約条件の下でELBOを最大化します。

最適化の結果、最適な $q_i^*(\theta_i)$ は以下の形になります：

$q_i^*(\theta_i) \propto \exp\left( \mathbb{E}_{j \neq i}\left[\log p(D, \theta)\right] \right)$

$= \exp\left( \int \prod_{j \neq i} q_j(\theta_j) \log p(D, \theta) \prod_{j \neq i} d\theta_j \right)$

この更新式は直感的に理解できます。 $q_i^*(\theta_i)$ は、他のパラメータの期待値を考慮した上での、 $\theta_i$ に関する「有効な」事後分布と解釈できます。

平均場近似の大きな利点は、共役事前分布を使用した場合に更新式が解析的に解けることが多い点です。例えば、線形回帰モデルでガウス・ガンマ共役事前分布を使うと、更新式が閉形式で表現できます。
しかし、平均場近似にはいくつかの重要な限界があります：

パラメータ間の相関を無視する
最大の欠点は、パラメータ間の依存関係を無視することです。実際の事後分布では強い相関があるパラメータでも、独立として扱われるため、事後分布の構造が正確に捉えられません。
分散の過小評価
平均場近似は、事後分布の分散（不確実性）を過小評価する傾向があります。これは、パラメータ間の相関を無視することに起因します。
モードの集中
複数のモード（峰）を持つ事後分布に対して、平均場近似は通常、単一のモード周辺に集中します。

概念的に説明すると、平均場近似は「複雑な山岳地帯を、互いに独立した単純な丘の集まりとして近似する」ようなものです。実際の山が互いに連なっていて相関があるのに対し、平均場近似ではそれぞれの山を別々に扱ってしまいます。これはちょうど、複雑な多変量確率分布を、各次元の周辺分布の積として近似するようなものです。

この制約にもかかわらず、平均場近似は多くの実用的な問題で効果的に機能し、計算効率と近似精度のバランスが取れた方法として広く使用されています。

3. MCMCと変分ベイズ法の比較

3.1 MCMCの基本的考え方

MCMC（マルコフ連鎖モンテカルロ法）は、直接サンプリングが難しい確率分布からサンプルを生成するための強力な手法です。MCMCの中心的なアイデアは、目標とする確率分布（この場合は事後分布 $p(\theta|D)$ ）を定常分布とするマルコフ連鎖を構築し、そのマルコフ連鎖を長時間実行することでサンプルを得るというものです。

MCMCのプロセスを詳しく見てみましょう。マルコフ連鎖は、現在の状態 $\theta^{(t)}$ から次の状態 $\theta^{(t+1)}$ への遷移が、過去の履歴ではなく現在の状態のみに依存するような確率過程です。MCMCでは、この遷移確率 $T(\theta^{(t+1)}|\theta^{(t)})$ を、連鎖の定常分布が目標分布 $p(\theta|D)$ になるように設計します。

具体的には、「詳細釣り合い条件」（detailed balance）を満たすように遷移確率を設計します：

$p(\theta|D)T(\theta'|\theta) = p(\theta'|D)T(\theta|\theta')$

この条件が満たされると、マルコフ連鎖の定常分布は $p(\theta|D)$ になることが保証されます。

MCMCの代表的なアルゴリズムを見てみましょう：

メトロポリス・ヘイスティングス（Metropolis-Hastings）アルゴリズム

現在の状態 $\theta^{(t)}$ から提案分布 $q(\theta'|\theta^{(t)})$ を使って新しい状態 $\theta'$ を提案する
受容率 $\alpha = \min\left(1, \frac{p(\theta'|D)q(\theta^{(t)}|\theta')}{p(\theta^{(t)}|D)q(\theta'|\theta^{(t)})}\right)$ を計算する
確率 $\alpha$ で新しい状態 $\theta'$ を受け入れ、 $\theta^{(t+1)} = \theta'$ とする。確率 $1-\alpha$ で棄却し、 $\theta^{(t+1)} = \theta^{(t)}$ とする

ギブスサンプリング（Gibbs Sampling）
ギブスサンプリングは、パラメータを成分ごとに更新するMCMCの特殊ケースです：

現在の状態を $\theta^{(t)} = (\theta_1^{(t)}, \theta_2^{(t)}, \ldots, \theta_d^{(t)})$ とする
各成分について、他の全ての成分を固定した条件付き分布からサンプリングする：

$\theta_i^{(t+1)} \sim p(\theta_i | \theta_1^{(t+1)}, \ldots, \theta_{i-1}^{(t+1)}, \theta_{i+1}^{(t)}, \ldots, \theta_d^{(t)}, D)$

ギブスサンプリングの利点は、提案分布が条件付き事後分布そのものであるため、常に受け入れられる（受容率が1）点です。しかし、これらの条件付き分布から効率的にサンプリングできることが前提となります。

ハミルトニアンモンテカルロ（Hamiltonian Monte Carlo, HMC）
HMCは、物理学のハミルトン動力学を利用して、より効率的なサンプリングを実現します：

各パラメータ $\theta_i$ に対応する「運動量」変数を導入する
ハミルトニアン $H(\theta, p) = U(\theta) + K(p)$ を定義する

（ $U(\theta) = -\log p(\theta|D)$ は「ポテンシャルエネルギー」、 $K(p) = \frac{1}{2}p^Tp$ は「運動エネルギー」）

ハミルトン方程式に従って $(\theta, p)$ の値を更新する
メトロポリス・ヘイスティングスステップで新しい状態を受け入れるか決定する

HMCの大きな利点は、勾配情報を利用して効率的に状態空間を探索できる点です。特に高次元のパラメータ空間で効果を発揮します。
　　
MCMCは強力な方法ですが、以下のような課題があります：

計算コストが高い：収束するまでに多数のサンプルが必要で、特に高次元パラメータ空間や大規模データセットでは計算負荷が大きい
収束判定が難しい：マルコフ連鎖が目標分布に収束したかどうかを厳密に判断するのは困難
連続的なサンプル間の自己相関：生成されたサンプルは独立ではなく、効果的なサンプル数は実際のサンプル数よりも少なくなる

これらの課題にもかかわらず、MCMCは理論的に十分な反復回数を重ねれば真の事後分布に収束することが保証されている点で、非常に価値のある手法です。

3.2 変分ベイズ法とMCMCの違い

変分ベイズ法とMCMCは、ベイズ推論における計算の困難さに対処するための二つの主要なアプローチですが、その根本的な発想と特性は大きく異なります。ここではその違いを多角的に分析します。

1. 基本的なアプローチ

MCMC：
MCMCは「サンプリング法」に分類されます。その目標は、事後分布からの代表的なサンプルを生成することです。これらのサンプルを使って、期待値や分位数などの事後分布の特性を推定します。MCMCは事後分布の形状に関する仮定をほとんど置かず、理論上は任意の複雑な分布をサンプリングできます。
　　
変分ベイズ法：
変分ベイズ法は「近似法」に分類されます。その目標は、事後分布を解析的に扱いやすい分布族で近似することです。変分ベイズ法は、KLダイバージェンスを最小化する（またはELBOを最大化する）という最適化問題に帰着させて解を求めます。この方法は計算効率が良い反面、選択した分布族に真の事後分布が含まれていない場合、完全な近似は原理的に不可能です。

2. 計算効率と精度のトレードオフ

MCMC：
MCMCは一般に計算コストが高く、特に高次元パラメータ空間や大規模データセットでは収束に時間がかかります。しかし、十分な時間をかければ理論上は真の事後分布に収束することが保証されています。MCMCは「正確だが遅い」方法と言えます。
　　
変分ベイズ法：
変分ベイズ法は計算効率が非常に良く、大規模データや複雑なモデルにも適用可能です。しかし、選択した近似分布族の制約内でしか最適化できないため、真の事後分布を完全に捉えられない場合があります。変分ベイズ法は「近似的だが速い」方法と言えます。

3. 出力形式

MCMC：
MCMCの出力は事後分布からのサンプル集合です。これらのサンプルを使って、ヒストグラムや密度推定によって事後分布の形状を可視化したり、様々な統計量（平均、分散、分位数など）を計算したりできます。
　　
変分ベイズ法：
変分ベイズ法の出力は、事後分布の近似として解析的な関数形（例：ガウス分布のパラメータ）です。この関数形から直接、期待値や分散などの統計量を計算できます。また、近似分布からサンプリングすることも可能です。

4. 不確実性の捉え方

MCMC：
MCMCは事後分布の形状に制約をほとんど課さないため、多峰性（マルチモーダル）分布や複雑な相関構造を持つ分布も正確に捉えることができます。これにより、パラメータの不確実性をより正確に表現できます。
　　
変分ベイズ法：
変分ベイズ法は、選択した近似分布族（例：平均場近似を使った独立分布の積）の制約により、複雑な分布構造を完全に捉えられないことがあります。特に、パラメータ間の相関や多峰性分布の表現が難しいという制約があります。

5. 適用シナリオ

MCMC：
MCMCは以下のようなシナリオに適しています：

事後分布の正確な特性化が重要な場合
計算時間よりも精度が優先される場合
モデルが複雑で、事後分布が多峰性や強い相関構造を持つ可能性がある場合
データセットのサイズが比較的小さい、または並列計算リソースが豊富にある場合

　　
変分ベイズ法：
変分ベイズ法は以下のようなシナリオに適しています：

大規模データセットや複雑なモデル（深層ニューラルネットワークなど）に対するベイズ推論
リアルタイムに近い応答が必要な場合
近似解でも十分な場合（例：予測性能が主な関心事の場合）
オンライン学習や逐次的更新が必要な場合

概念的に説明すると、MCMCは「実際に地形をあちこち歩き回ってサンプルを集める」アプローチであり、変分ベイズ法は「地形を単純な数学的関数で近似する」アプローチです。MCMCは時間をかければより正確な描写ができますが、変分ベイズ法は素早く全体像の「スケッチ」を提供します。
両者は相互に排他的ではなく、相補的な手法として、問題の性質や計算リソースに応じて適切に選択することが重要です。また、両者を組み合わせたハイブリッドアプローチも研究されています。

3.3 変分ベイズの推定値がMCMCと異なる原因

変分ベイズ法による推定値がMCMCの結果と異なる場合がしばしばあります。ここでは、その主な原因を詳細に分析します。

1. 近似分布の表現力不足

ガウス分布族の制約：
変分ベイズ法でよく使用されるガウス分布は、単峰性で対称な分布です。しかし、実際の事後分布は：

多峰性（複数のモードを持つ）かもしれません
歪んでいるかもしれません
特異点や不連続点を持つかもしれません
重い裾（heavy tail）を持つかもしれません

これらの特性は単一のガウス分布では捉えられません。例えば、混合モデルのパラメータは典型的に多峰性の事後分布を持ちますが、ガウス変分分布はこれを単一のモードとして近似してしまいます。
　　
平均場近似の制約：

平均場近似（ $\theta = \prod_i q_i(\theta_i)$ ）を用いると、パラメータ間の依存関係が完全に無視されます。しかし実際の事後分布では、パラメータ間に強い相関がある場合が多いです。例えば、線形回帰モデルでの係数間の相関や、ニューラルネットワークの重みパラメータ間の複雑な依存関係などです。

具体例として、二つのパラメータ $\theta_1$ と $\theta_2$ が強い負の相関を持つ事後分布を考えてみましょう。MCMCサンプルはこの相関構造を正確に捉えますが、平均場近似を用いた変分ベイズ法では $q(\theta_1)$ と $q(\theta_2)$ が独立となるため、この相関を完全に無視してしまいます。

2. KLダイバージェンスの非対称性

変分ベイズ法で最小化するKLダイバージェンス $KL(q(\theta) || p(\theta|D))$ （フォワードKL）は非対称です。この順序のKLダイバージェンスの特性を理解することが重要です：

$KL(q || p) = \int q(\theta) \log \frac{q(\theta)}{p(\theta|D)} d\theta$

この式からわかるように、 $q(\theta)$ がゼロの領域では、たとえ $p(\theta|D)$ が大きな値を持っていても、積分への寄与はゼロになります。しかし、 $p(\theta|D)$ がゼロに近いのに $q(\theta)$ が大きな値を持つ領域では、 $\log \frac{q(\theta)}{p(\theta|D)}$ は非常に大きくなり、KLダイバージェンスに大きく寄与します。

これは、フォワードKLが「ゼロ回避型」（zero-avoiding）と呼ばれる性質を持つことを意味します。変分分布 $q(\theta)$ は、事後分布 $p(\theta|D)$ がゼロに近い領域をカバーしないようにする（ゼロを避ける）傾向があります。

その結果、変分ベイズ法は事後分布のモードを捉えつつも、その広がり（分散）を過小評価する傾向があります。つまり、変分分布は事後分布よりも「尖った」（より確信的な）形状になりがちです。

3. 平均場近似によるパラメータ間の独立性の仮定

平均場近似はパラメータ間の独立性を仮定しますが、この仮定が実際の事後分布と大きく異なる場合、推定値に重大な影響を与えます：

分散の過小評価：パラメータ間の正の相関がある場合、平均場近似はこれを無視するため、各パラメータの分散を過小評価します。
条件付き平均のシフト：パラメータ間に相関がある場合、一方のパラメータの条件付き期待値は他方の値に依存します。平均場近似ではこれを捉えられず、推定された平均値が真の事後平均からシフトすることがあります。

4. 最適化の局所解

変分ベイズ法はELBOの最大化という最適化問題を解きますが、これが局所解に陥る場合があります：

非凸最適化問題：複雑なモデルでは、ELBO最大化が非凸最適化問題になります。
確率的最適化の課題：大規模データセットでは、確率的勾配法などの使用によって収束解にバラツキが出ます。

一方、MCMCは（理論上は）局所解問題に悩まされず、十分な実行時間があれば真の事後分布全体を探索できます。

5. 具体例：バイモーダル分布

バイモーダル（二峰性）の事後分布を考えてみましょう。例えば、次のような混合ガウスモデルの事後分布：

$p(\theta|D) \approx 0.5 \times \mathcal{N}(\theta|-3, 1) + 0.5 \times \mathcal{N}(\theta|3, 1)$

MCMCはこの二つのモード両方からサンプリングでき、ヒストグラムは二つのピークを持ちます。結果として事後平均は約0になります。
しかし、単一のガウス分布で近似する変分ベイズ法の場合：

$q(\theta) = \mathcal{N}(\theta|\mu_q, \sigma_q^2)$ を仮定
KLダイバージェンス最小化の結果、一方のモードに集中するか、または両モードの間に大きな分散で位置する可能性

どちらの場合も、真の事後分布の重要な構造を捉えられていません。

4. 変分ベイズ法の実装手法

4.1 平均場変分ベイズ（MFVB）

平均場変分ベイズ（Mean Field Variational Bayes, MFVB）は、前述の平均場近似に基づく最も基本的な変分ベイズの実装です。MFVBの核心は、パラメータ $\theta$ を複数の成分 $\theta_1, \theta_2, \ldots, \theta_M$ に分解し、これらのパラメータが互いに独立であると仮定することです：

$q(\theta) = \prod_{i=1}^{M} q_i(\theta_i)$

この仮定の下で、各成分の最適な変分分布 $q_i^(\theta_i)$ は以下の更新式で与えられます：

$q_i^(\theta_i) \propto \exp(\mathbb{E}{j \neq i}[\log p(D, \theta)])$

ここで $\mathbb{E}{j \neq i}$ は $\theta_i$ 以外の全てのパラメータに関する期待値を表します。

MFVBのアルゴリズムは以下の通りです：

各 $\theta_i$ を適当な分布で初期化します。
収束するまで以下を繰り返します：
　a. 各 $\theta_i$ について、以下の更新を行います：
　 $q_i^{new}(\theta_i) \propto \exp(\mathbb{E}_{j \neq i}[\log p(D, \theta)$ )]
　b. 収束判定を行います（例：ELBO増加率が閾値以下になるか、あるいは最大反復回数に達するまで）。最終的な近似事後分布

$q^(\theta) = \prod_i q_i^(\theta_i)$

を出力します。

共役事前分布の場合

MFVBの最大の利点は、共役事前分布を使用した場合に更新式が解析的に解けることが多い点です。例えば、一般化線形モデルの多くは、適切な共役事前分布を選ぶことで閉形式の更新式が導出できます。
具体的な例として、単純な線形回帰モデルを考えてみましょう：

$y = X\beta + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2 I)$

ここで $\beta$ は回帰係数、 $\sigma^2$ は誤差分散です。これらのパラメータに対して以下の共役事前分布を設定します：

$\beta \sim \mathcal{N}(\beta_0, \Sigma_0)$

$\sigma^2 \sim \text{InvGamma}(a_0, b_0)$

平均場近似では、 $\beta, \sigma^2 = q(\beta)q(\sigma^2)$ と仮定します。このとき、最適な変分分布は以下の形をとります：

$q^(\beta) = \mathcal{N}(\beta | \mu_\beta, \Sigma_\beta)$

$q^(\sigma^2) = \text{InvGamma}(\sigma^2 | a, b)$

これらのパラメータは以下の閉形式の更新式で計算できます：

$\Sigma_\beta = (\Sigma_0^{-1} + \mathbb{E}[1/\sigma^2]X^TX)^{-1}$

$\mu_\beta = \Sigma_\beta(\Sigma_0^{-1}\beta_0 + \mathbb{E}[1/\sigma^2]X^Ty)$

$a = a_0 + n/2$

$b = b_0 + \frac{1}{2}\mathbb{E}[(y - X\beta)^T(y - X\beta)]$

$\mathbb{E}[1/\sigma^2] = a/b$

$\mathbb{E}[(y - X\beta)^T(y - X\beta)] = ||y - X\mu_\beta||^2 + \text{Tr}(X^TX\Sigma_\beta)$

これらの更新式を収束するまで繰り返すことで、近似事後分布を得ることができます。
　　

MFVBの長所と短所

長所：

共役事前分布を使用した場合、解析的な更新式が得られる
計算効率が良く、低次元から中程度の次元のモデルに適している
実装が比較的容易
各パラメータ成分の更新を並列化できる可能性がある

短所：

パラメータ間の相関を無視するため、事後分布の構造を正確に捉えられない
事後分布の分散（不確実性）を過小評価する傾向がある
非共役モデルでは更新式が解析的に解けないことが多い
確率変数が多い複雑なモデルでは、全ての変数に対して更新式を導出する必要があり、実装が煩雑になる

MFVBは、比較的単純なモデルや、パラメータ間の相関があまり重要でない問題に適しています。また、初期の探索的分析や、より複雑な方法の初期値を提供するためにも有用です。

4.2 固定形変分ベイズ（FFVB）

固定形変分ベイズ（Fixed-Form Variational Bayes, FFVB）では、近似分布 $q(\theta)$ の関数形を予め指定し、そのパラメータを最適化します。最も一般的なのがガウス分布を使用する方法です：

$q(\theta) = \mathcal{N}(\theta | \mu, \Sigma)$

ここで $\mu$ は平均ベクトル、 $\Sigma$ は共分散行列です。ELBOを $\mu$ と $\Sigma$ に関して最大化することで、最適な近似分布を求めます。

共分散行列の表現

完全な共分散行列を使うとパラメータ数が $d^2$ （ $d$ はパラメータの次元数）と増加するため、計算コストが高くなります。そこで、以下のような様々な共分散構造が考えられます：

　　
対角共分散行列：

$\Sigma = \text{diag}(\sigma_1^2, \sigma_2^2, \ldots, \sigma_d^2)$

この場合、パラメータ数は $d$ に削減されますが、パラメータ間の相関を全く捉えられません。
　　
因子分析型共分散：

$\Sigma = \text{diag}(\sigma_1^2, \sigma_2^2, \ldots, \sigma_d^2) + VV^T$

ここで $V$ は $d \times k$ 行列（ $k \ll d$ ）です。これにより、 $dk$ のパラメータで主要な相関構造を捉えることができます。

　　
コレスキー分解：

$\Sigma = LL^T$

ここで $L$ は下三角行列です。この表現は共分散行列の正定値性を保証しますが、パラメータ数は $d^2$ です。

最適化手法

固定形変分ベイズでは、通常の勾配ベースの最適化アルゴリズム（勾配降下法、L-BFGS、Adamなど）を使ってELBOを最大化します。ELBOの勾配は以下のように計算できます：

$\nabla_\lambda \text{ELBO} = \nabla_\lambda \mathbb{E}{q(\theta|\lambda)}[\log p(D, \theta) - \log q(\theta|\lambda)]$

ここで $\lambda$ は変分分布のパラメータ（例： $\mu$ と $\Sigma$ の要素）を表します。

この期待値の勾配を直接計算するのは難しいですが、以下の二つの主要なアプローチがあります：　　
スコア関数勾配（Score Function Gradient）：

$\nabla\lambda \text{ELBO} = \mathbb{E}{q(\theta|\lambda)}[(\log p(D, \theta) - \log q(\theta|\lambda)) \nabla\lambda \log q(\theta|\lambda)]$

この方法はブラックボックス変分推論（後述）の基礎にもなっています。　　
リパラメトリゼーション勾配（Reparameterization Gradient）：

変数変換 $\theta = g(\epsilon, \lambda)$ （ $\epsilon$ は標準分布からのサンプル）を使って：

$\nabla_\lambda \text{ELBO} = \mathbb{E}{p(\epsilon)}[\nabla\lambda \log p(D, g(\epsilon, \lambda)) - \nabla_\lambda \log q(g(\epsilon, \lambda)|\lambda)]$

例えば、ガウス分布の場合は $\theta = \mu + L\epsilon$ （ $\epsilon \sim \mathcal{N}(0, I)$ ）とできます。

FFVBの長所と短所

長所：

適切な共分散構造を使えばパラメータ間の相関を捉えられる
一般的な最適化アルゴリズムを使用できる
平均場近似よりも事後分布の構造をより正確に近似できる可能性がある
モデル固有の更新式を導出する必要がなく、汎用的なアプローチとして使える

短所：

選択した関数形（例：ガウス分布）が真の事後分布と大きく異なる場合、近似精度が低下する
高次元パラメータの場合、共分散行列の表現と最適化が計算上の課題になる
勾配計算のためのサンプリングが必要で、計算コストが増加する可能性がある

FFVBは、パラメータ間の相関が重要で、なおかつ真の事後分布がガウス分布などの標準的な分布で良く近似できる問題に適しています。特に、事後分布が単峰性で比較的対称な場合に効果的です。

4.3 ブラックボックス変分推論（BBVI）

上記の手法はモデルごとに導出や実装が必要でしたが、ブラックボックス変分推論（Black-Box Variational Inference, BBVI）は汎用的なアプローチを提供します。BBVIはモンテカルロサンプリングと確率的勾配降下法を用いて、モデルの詳細に関わらずELBOを最適化します。
BBVIのキーアイデアは、ELBOの勾配をサンプリングによって近似することです：

$\nabla_\lambda \text{ELBO} = \mathbb{E}{q(\theta|\lambda)}[(\log p(D, \theta) - \log q(\theta|\lambda)) \nabla\lambda \log q(\theta|\lambda)]$

モンテカルロ近似を用いると：

$\nabla_\lambda \text{ELBO} \approx \frac{1}{S} \sum_{s=1}^{S} (\log p(D, \theta^s) - \log q(\theta^s | \lambda)) \nabla_\lambda \log q(\theta^s | \lambda)$

ここで $\theta^s \sim q(\theta | \lambda)$ は変分分布からのサンプルです。 $\nabla_\lambda \log q(\theta | \lambda)$ はスコア関数（score function）と呼ばれ、多くの標準的な分布で解析的に計算できます。

勾配分散低減

上記の勾配推定は分散が大きく、収束が遅い原因になります。そこで、分散を低減するための様々なテクニックが提案されています：

制御変量（Control Variates）：期待値は変わらずに分散を減らす関数 $h(\theta)$ を導入します：

$\nabla_\lambda \text{ELBO} \approx \frac{1}{S} \sum_{s=1}^{S} (f(\theta^s) - a \cdot h(\theta^s)) \nabla_\lambda \log q(\theta^s | \lambda)$

ここで $f(\theta) = \log p(D, \theta) - \log q(\theta | \lambda)$ 、 $a$ は分散を最小化する係数です。よく使われる制御変量としては、 $h(\theta) = \nabla_\lambda \log q(\theta | \lambda)$ などがあります。

Rao-Blackwellization：
条件付き期待値を取ることで分散を減らします：

$\nabla_{\lambda_i} \text{ELBO} \approx \frac{1}{S} \sum_{s=1}^{S} \mathbb{E}{q(\theta{-i}|\theta_i^s, \lambda)}[f(\theta)] \nabla_{\lambda_i} \log q(\theta_i^s | \lambda_i)$

ここで $\theta_{-i}$ は $\theta_i$ 以外の全てのパラメータを表します。

BBVIのアルゴリズム

BBVIのアルゴリズムは以下の通りです：

変分パラメータ $\lambda$ を初期化します。
収束するまで以下を繰り返します：

　a. $q(\theta | \lambda)$ から $S$ 個のサンプル $\theta^1, \ldots, \theta^S$ を生成します。

　b. 各サンプルについて $f(\theta^s) = \log p(D, \theta^s) - \log q(\theta^s | \lambda)$ とスコア関数 $\nabla_\lambda \log q(\theta^s | \lambda)$ を計算します。

　c. ELBOの勾配を近似します：

$\nabla_\lambda \text{ELBO} \approx \frac{1}{S} \sum_{s=1}^{S} f(\theta^s) \nabla_\lambda \log q(\theta^s | \lambda)$

　d. 確率的勾配降下法で $\lambda$ を更新します：

$\lambda \leftarrow \lambda + \rho \nabla_\lambda \text{ELBO}$

（ $\rho$ は学習率）

　e. 必要に応じて学習率 $\rho$ を調整します。

最終的な近似事後分布 $q(\theta | \lambda^*)$ を出力します。

BBVIの長所と短所

長所：

モデルに対して汎用的：尤度関数と事前分布が計算できれば適用可能
モデル固有の導出が不要で、「ブラックボックス」として使用可能
様々な変分分布（ガウス分布、混合分布など）に適用可能
自動微分と組み合わせることで、実装が容易になる

短所：

サンプリングに基づく勾配推定の分散が大きく、収束が遅い場合がある
分散低減テクニックが必要で、実装が複雑になる可能性がある
勾配計算のためのサンプリング数 $S$ の選択が結果に影響する
リパラメトリゼーショントリックを使わない場合、連続パラメータに対する勾配推定の効率が悪い

BBVIは、特に既製のソフトウェアパッケージを使用する場合や、カスタムモデルを迅速に開発する必要がある場合に、非常に魅力的なアプローチです。モデルに特化した導出をせずに変分ベイズ法を適用できる柔軟性は、実践的な観点から大きな利点です。

4.4 自動微分変分推論（ADVI）

自動微分変分推論（Automatic Differentiation Variational Inference, ADVI）は、BBVIをさらに発展させ、自動微分を活用して効率的に勾配を計算する手法です。ADVIの主な特徴は以下の2点です：

制約のある変数の変換：制約のあるパラメータ空間（例：正の値のみ、[0,1]の範囲など）を、制約のない実数空間に変換します。
リパラメトリゼーショントリック：サンプリングの分散を低減するために、リパラメトリゼーションを用います。

変数変換

ADVIでは、以下のような変換を使って制約のあるパラメータを制約のない空間に変換します：

正の値（ $\theta > 0$ ）：対数変換 $\zeta = \log(\theta)$
区間（ $\theta \in (a,b)$ ）：ロジット変換 $\zeta = \log\left(\frac{\theta-a}{b-\theta}\right)$
シンプレックス（ $\theta \in \Delta^{K-1}$ ）：スティック破り変換とロジット変換の組み合わせ
相関行列：コレスキー因子と対角成分の変換

これらの変換により、あらゆるパラメータを制約のない空間で扱えるようになり、単一のガウス変分分布を適用できます。

変換されたパラメータ $\zeta$ に対して、変分分布を以下のように設定します：

$q(\zeta | \lambda) = \mathcal{N}(\zeta | \mu, \Sigma)$

ここで $\lambda = {\mu, \Sigma}$ は変分パラメータです。

リパラメトリゼーショントリック

ADVIでは、リパラメトリゼーショントリックを用いて勾配の低分散推定を行います：

$\zeta = \mu + L\epsilon, \quad \epsilon \sim \mathcal{N}(0, I)$

ここで $L$ は $\Sigma = LL^T$ を満たす行列（例：コレスキー因子）です。このトリックにより、期待値の勾配が以下のように表現できます：

$\nabla_\mu \text{ELBO} = \mathbb{E}{\epsilon \sim \mathcal{N}(0, I)}[\nabla\zeta \log p(D, T^{-1}(\zeta)) - \nabla_\zeta \log q(\zeta | \lambda)]$

$\nabla_L \text{ELBO} = \mathbb{E}{\epsilon \sim \mathcal{N}(0, I)}[(\nabla\zeta \log p(D, T^{-1}(\zeta)) - \nabla_\zeta \log q(\zeta | \lambda)) \epsilon^T + \nabla_L H[q(\zeta | \lambda)]]$

ここで $T^{-1}$ は変換の逆関数です。

ADVIのアルゴリズム

ADVIのアルゴリズムは以下の通りです：

制約のあるパラメータ $\theta$ を制約のない変数 $\zeta$ に変換します。
変分パラメータ $\lambda = {\mu, L}$ を初期化します。
収束するまで以下を繰り返します：

　a. 標準正規分布から $S$ 個のサンプル $\epsilon^1, \ldots, \epsilon^S$ を生成します。

　b. 各サンプルについて $\zeta^s = \mu + L\epsilon^s$ を計算し、さらに $\theta^s = T^{-1}(\zeta^s)$ を計算します。

　c. 各サンプルについて $\nabla_\zeta \log p(D, \theta^s)$ と $\nabla_\zeta \log q(\zeta^s | \lambda)$ を計算します。

　d. ELBOの勾配を近似し、 $\lambda$ を更新します。

4. 最終的な近似事後分布を $q(\theta | \lambda^*)$ （元のパラメータ空間）として出力します。

ADVIの長所と短所

長所：

制約のあるパラメータも扱える汎用的な手法
リパラメトリゼーショントリックによる効率的な勾配推定
自動微分と組み合わせることで、実装が容易になる
様々なモデルに対して統一的なアプローチが可能

短所：

単一のガウス分布を使用するため、複雑な事後分布（多峰性や歪んだ分布）の近似には限界がある
変数変換によるヤコビアン補正が必要で、計算コストが増加する
高次元パラメータの場合、共分散行列の表現が計算上の課題になる

ADVIは特に、制約のあるパラメータを含むモデル（例：階層モデル、ベイズ非負行列分解など）に対して効果的です。また、Stan、PyMC、TensorFlowなどのソフトウェアパッケージに実装されており、実用的なベイズ推論のデファクトスタンダードとなりつつあります。

4.5 確率的変分推論（SVI）

確率的変分推論（Stochastic Variational Inference, SVI）は、大規模データセットに対する変分ベイズ法の計算効率を向上させるために開発されました。SVIの核心は、データのミニバッチを使ってELBOの勾配を確率的に推定することです。

ELBOのデータ分解

データ $D = {x_1, \ldots, x_N}$ に対するELBOは以下のように分解できます：

$\text{ELBO} = \mathbb{E}{q(\theta)}[\log p(\theta)] - \mathbb{E}{q(\theta)}[\log q(\theta)] + \sum_{i=1}^{N} \mathbb{E}_{q(\theta)}[\log p(x_i | \theta)]$

データ項 $\sum_{i=1}^{N} \mathbb{E}_{q(\theta)}[\log p(x_i | \theta)$ ]は、各データポイントの寄与の合計です。大規模データセットでは、この項の計算が時間的にボトルネックになります。

ミニバッチ勾配

SVIでは、ミニバッチ $\mathcal{B} \subset {1, \ldots, N}$ を使って、データ項を近似します：

$\sum_{i=1}^{N} \mathbb{E}{q(\theta)}[\log p(x_i | \theta)] \approx \frac{N}{|\mathcal{B}|} \sum{i \in \mathcal{B}} \mathbb{E}_{q(\theta)}[\log p(x_i | \theta)]$

ここで $|\mathcal{B}|$ はミニバッチのサイズです。この近似により、計算コストを $O(N)$ から $O(|\mathcal{B}|)$ に削減できます。

ナチュラル勾配

SVIではしばしば、ELBOの通常の勾配ではなく、ナチュラル勾配を使用します。ナチュラル勾配は、変分分布のパラメータ空間における情報幾何学的な構造を考慮した勾配です。
指数型分布族（ガウス分布、ベルヌーイ分布など多くの標準的な分布はこれに含まれる）の場合、ナチュラル勾配は以下のように表現できます：

$\hat{\nabla}\lambda \text{ELBO} = G(\lambda)^{-1} \nabla\lambda \text{ELBO}$

ここで $G(\lambda)$ はフィッシャー情報行列です。例えば、平均場近似で各変分分布が指数型分布族に属する場合、ナチュラル勾配を用いた更新式は特に簡潔になります。

SVIのアルゴリズム

SVIのアルゴリズムは以下の通りです：

変分パラメータ $\lambda$ を初期化します。
収束するまで以下を繰り返します：

　a. データからミニバッチ $\mathcal{B}$ をランダムに選択します。

　b. ミニバッチを使ってELBOの勾配（またはナチュラル勾配）を近似します：

　c. 確率的勾配降下法で $\lambda$ を更新します：

$\lambda \leftarrow \lambda + \rho_t \nabla\lambda \text{ELBO}_{\mathcal{B}}$

（ $\rho_t$ はステップ $t$ での学習率）

最終的な近似事後分布 $q(\theta | \lambda^*)$ を出力します。

学習率スケジューリング

SVIの収束には適切な学習率スケジューリングが重要です。通常、以下のようなスケジュールが使われます：

$\rho_t = (t_0 + t)^{-\kappa}$

ここで $t_0 > 0$ と $\kappa \in (0.5, 1$ ]はハイパーパラメータです。このスケジュールは確率的近似の理論に基づいており、適切な条件下での収束を保証します。

SVIの長所と短所

長所：

大規模データセットに対する計算効率の良さ
オンライン学習や逐次的更新が自然に可能
様々な変分ベイズ法（MFVB、FFVB、BBVI、ADVIなど）と組み合わせて使用できる
現代の分散計算フレームワークで並列処理が容易

短所：

ミニバッチによる勾配推定の分散が大きく、収束に時間がかかる場合がある
学習率やミニバッチサイズなどのハイパーパラメータの選択が結果に影響する
限られたミニバッチでは、完全なデータセットの情報を捉えきれない場合がある

SVIは特に、大規模データセットを扱うベイズモデル（トピックモデル、ベイズニューラルネットワーク、大規模階層モデルなど）に適しています。現代の機械学習において、データサイズの増大とともにその重要性はますます高まっています。

5. 変分ベイズ法の結果をMCMCに近づける方法

変分ベイズ法の結果をMCMCの結果に近づけるために、様々な方法が提案されています。ここでは、主要なアプローチを詳しく解説します。

5.1 より表現力のある変分族の使用

変分ベイズ法の近似精度を向上させる最も直接的な方法は、より表現力のある変分分布族を使用することです。

5.1.1 フルランク共分散行列

最も基本的なアプローチは、ガウス変分分布で対角共分散行列ではなく、フルランク共分散行列を使用することです：

$q(\theta) = \mathcal{N}(\theta | \mu, \Sigma)$

ルランク共分散行列を使用すると、パラメータ間の相関関係を捉えることができます。これにより、平均場近似の主要な制約であるパラメータ独立性の仮定を緩和できます。
例えば、線形回帰モデルでは係数間に強い相関がありますが、フルランク共分散行列を使えばこの相関構造を適切に表現できます。

ただし、パラメータ数が $O(d^2)$ （ $d$ はパラメータの次元数）に増加するため、高次元の場合は計算コストが高くなるという欠点があります。このトレードオフに対処するために、以下のような中間的なアプローチも考えられます：

低ランク近似：

$\Sigma = \text{diag}(\sigma^2) + VV^T$

ここで $V$ は $d \times k$ 行列（ $k \ll d$ ）です。これにより、 $O(dk)$ のパラメータで主要な相関構造を捉えることができます。

スパース近似： $\Sigma$ の多くの要素をゼロと仮定し、重要な相関のみを表現する方法です。グラフィカルモデルの構造を利用して、条件付き独立関係を特定できる場合に特に有効です。

フルランク共分散行列の使用は、パラメータ次元が中程度（数十～数百）で、パラメータ間の相関が重要な問題に特に効果的です。

5.1.2 ノーマライジングフロー

ノーマライジングフローは、単純な分布（例：標準ガウス分布）を一連の可逆変換を通じて複雑な分布に変換する手法です：

$z_0 \sim q_0(z_0), \quad z_K = f_K \circ f_{K-1} \circ \cdots \circ f_1(z_0)$

ここで $f_1, \ldots, f_K$ は可逆な変換関数です。最終的な変分分布は以下のようになります：

$q_K(z_K) = q_0(z_0) \prod_{k=1}^{K} \left| \det \frac{\partial f_k^{-1}}{\partial z_k} \right|$

この式は変数変換の公式に基づいており、ヤコビアン行列式の絶対値の積が分布の変換に対する「補正項」となっています。

ノーマライジングフローの鍵は、以下の条件を満たす変換関数 $f_k$ の設計です：

可逆であること（逆関数 $f_k^{-1}$ が存在する）
ヤコビアン行列式が計算可能であること
表現力が十分であること

主要なノーマライジングフローの例としては以下があります：

プラナーフロー（Planar Flow）：

$f(z) = z + u h(w^T z + b)$

ここで $h$ は活性化関数（例：tanh）、 $u, w$ はベクトル、 $b$ はスカラーパラメータです。

放射状フロー（Radial Flow）：

$f(z) = z + \beta h(\alpha, r)(z - z_0)$

ここで $r = |z - z_0|$ 、 $h(\alpha, r) = \frac{1}{\alpha + r}$ 、 $\alpha, \beta, z_0$ はパラメータです。

アフィンカップリングフロー（Affine Coupling Flow）：

$z$ を二つの部分 $z_1, z_2$ に分割し、

$f(z_1, z_2) = (z_1, z_2 \odot \exp(s(z_1)) + t(z_1))$

ここで $s, t$ はニューラルネットワークで表現される関数、 $\odot$ は要素ごとの積です。

これらの変換を複数層重ねることで、非常に複雑な分布も表現可能になります。ノーマライジングフローは、特に多峰性や歪んだ事後分布の近似に有効です。ノーマライジングフローの利点は、十分な変換層を使えば非常に複雑な分布（マルチモーダルや非対称な分布）も表現できる点です。欠点は、計算コストが高く、適切な変換関数の設計が難しい点です。

5.1.3 混合ガウス分布

複数のガウス分布の混合を変分分布として使用する方法も効果的です：

$q(\theta) = \sum_{i=1}^{K} \pi_i \mathcal{N}(\theta | \mu_i, \Sigma_i)$

ここで $\pi_i$ は混合重み（ $\sum_i \pi_i = 1$ ）、 $\mu_i$ と $\Sigma_i$ はそれぞれ $i$ 番目の成分の平均と共分散です。混合ガウス分布を用いた変分ベイズ法では、各ガウス成分のパラメータに加えて混合重みも最適化パラメータとなります。ELBOの計算は複雑になりますが、以下のようなアプローチがあります：

モンテカルロ近似：サンプルを生成して期待値を近似する方法です。各サンプルは、まず混合成分 $i$ を確率 $\pi_i$ で選び、次にその成分からサンプリングします。
変分下界のさらなる近似：混合分布の対数の期待値は解析的に計算できないため、ジェンセンの不等式などを用いてさらに近似することもあります。

混合ガウス分布は、異なるモードを持つマルチモーダルな事後分布を捉えるのに適しています。例えば、混合モデルのパラメータや、複数の局所解を持つ非線形モデルの事後分布などに効果的です。
ただし、混合成分数 $K$ の選択や、各成分のパラメータ最適化が難しくなります。また、成分数に比例して計算コストが増加する点も考慮する必要があります。

5.2 改良されたKLダイバージェンス

5.2.1 リバースKLダイバージェンス

変分ベイズ法では通常、 $KL(q(\theta) || p(\theta|D))$ （フォワードKL）を最小化しますが、これを $KL(p(\theta|D) || q(\theta))$ （リバースKL）に変更する手法もあります。フォワードKLとリバースKLの振る舞いの違いを理解することが重要です：

フォワードKL： $KL(q || p)$

$q(\theta)$ が $p(\theta|D)$ を「カバー」しようとする（zero-avoiding）
$p(\theta|D)$ がゼロの領域で $q(\theta)$ もゼロになることを強制する
事後分布の単一のモードに集中しがちで、分散を過小評価する傾向がある

リバースKL： $KL(p || q)$

$q(\theta)$ が $p(\theta|D)$ の高密度領域に「一致」しようとする（zero-forcing）
$p(\theta|D)$ の複数のモードがある場合、 $q(\theta)$ はそのいずれかに集中しがち
事後分布の単一のモードをより正確に捉える傾向がある

リバースKLの課題は、事後分布 $p(\theta|D)$ を直接計算できない場合、これを最小化することが困難である点です。そこで、事後分布のサンプルを使って近似的に最小化する方法が提案されています。例えば、MCMCやその他のサンプリング法で事後分布からサンプル ${\theta^s}_{s=1}^S$ を生成し、以下の目的関数を最小化します：

$\hat{KL}(p || q) \approx \frac{1}{S} \sum_{s=1}^S (-\log q(\theta^s | \lambda)) + \text{const.}$

これは、サンプルの対数尤度を最大化することに相当します。MCMC-VAEなどのモデルでは、この考え方が応用されています。
リバースKLアプローチの利点は、事後分布のモードをすべて捉える傾向があることです。欠点は、事後分布からのサンプリングが必要であり、計算コストが高くなることです。

5.2.2 アルファダイバージェンス

KLダイバージェンスを一般化したアルファダイバージェンスを使用する方法もあります：

$D_\alpha(p || q) = \frac{4}{1-\alpha^2} \left(1 - \int p(x)^{\frac{1+\alpha}{2}} q(x)^{\frac{1-\alpha}{2}} dx \right)$

パラメータ $\alpha \in [-1, 1$ ]を変化させることで、様々な特性を持つダイバージェンスを得ることができます：

$\alpha \to 1$ のとき、フォワードKL： $KL(p || q)$
$\alpha \to -1$ のとき、リバースKL： $KL(q || p)$
$\alpha = 0$ のとき、ヘリンジャー距離（Hellinger distance）

アルファダイバージェンスを用いた変分推論（Variational Inference with α-divergence, VI-α）では、一般化された変分下界を最大化します：

$L_\alpha(\lambda) = -\frac{1}{\alpha} \log \mathbb{E}_{q(\theta | \lambda)} \left[ \left( \frac{p(D, \theta)}{q(\theta | \lambda)} \right)^\alpha \right]$

VI-αの実装には、重点サンプリング（importance sampling）などの技術が用いられます。アルファダイバージェンスの利点は、 $\alpha$ を調整することで、近似の特性（モードを捉える vs. 広がりを捉える）をコントロールできる点です。

5.2.3 最大平均乖離推定（MMD）

カーネル平均埋め込みに基づく最大平均乖離（Maximum Mean Discrepancy, MMD）を使用する方法もあります。MMDは二つの分布間の距離を、再生核ヒルベルト空間（RKHS）における平均埋め込みの距離として定義します：

$\text{MMD}(p, q) = ||\mu_p - \mu_q||_{\mathcal{H}}^2$

ここで $\mu_p = \mathbb{E}{p(x)}[k(x, \cdot)]$ 、 $\mu_q = \mathbb{E}{q(x)}[k(x, \cdot)]$ はそれぞれの分布のカーネル平均埋め込みで、 $k(\cdot, \cdot)$ は正定値カーネル（例：ガウスカーネル）です。サンプル ${x_i}{i=1}^m \sim p$ と ${y_j}{j=1}^n \sim q$ を用いて、MMDは以下のように推定できます：

$\widehat{\text{MMD}}(p, q) = \frac{1}{m^2} \sum_{i,i'} k(x_i, x_{i'}) - \frac{2}{mn} \sum_{i,j} k(x_i, y_j) + \frac{1}{n^2} \sum_{j,j'} k(y_j, y_{j'})$

MMDを用いた変分推論では、カーネルパラメータの選択が重要になります。適切なカーネルを使えば、KLダイバージェンスでは捉えられない分布の特性を考慮できる可能性があります。 MMDの利点は、パラメトリックな分布族の制約を受けず、より柔軟な分布の比較が可能な点です。欠点は、カーネルの選択や計算コストの高さが挙げられます。

5.3 変分推論の後処理

変分推論の結果を後処理することで、近似の質を向上させる方法もあります。

5.3.1 重点サンプリング（Importance Sampling）

変分分布 $q(\theta | \lambda)$ からサンプル ${\theta^s}_{s=1}^S$ を生成し、重み $w_s = \frac{p(D, \theta^s)}{q(\theta^s | \lambda)}$ を計算します。これらの重み付きサンプルを使って、期待値などの統計量をより正確に推定できます：

$\mathbb{E}{p(\theta | D)}[f(\theta)] \approx \frac{\sum{s=1}^S w_s f(\theta^s)}{\sum_{s=1}^S w_s}$

重点サンプリングは特に、予測分布の計算などに有効です。変分推論自体の目的である「近似事後分布 $q(\theta | \lambda)$ の最適化」はそのままに、その結果を用いた推論をより正確にするアプローチです。この方法は、変分近似で捉えられなかった事後分布の特性を、ある程度補正できる可能性があります。ただし、変分分布と真の事後分布の乖離が大きい場合、重みの分散が大きくなり効果が限定的になることがあります。

5.3.2 マルコフ連鎖変分推論（MCVI）

変分分布 $q(\theta | \lambda)$ を初期分布として、短いマルコフ連鎖を走らせる方法もあります。具体的には：

変分ベイズ法で近似事後分布 $q(\theta | \lambda)$ を得る
$q(\theta | \lambda)$ からサンプル $\theta^{(0)}$ を生成
このサンプルを初期値として、短いMCMCチェーン（例：メトロポリス・ヘイスティングス法）を実行し、 $\theta^{(1)}, \ldots, \theta^{(T)}$ を得る
最終的なサンプル $\theta^{(T)}$ （またはburn-in後のサンプル）を使用

このアプローチは、変分分布が真の事後分布の高密度領域に近い場合に特に効果的です。MCMCの収束を加速しつつ、真の事後分布からのサンプルを得られる可能性があります。

5.3.3 多変量デルタ法（Multivariate Delta Method）

変分分布の平均と分散を用いて、非線形変換された量の期待値や分散を近似する方法もあります。例えば、関数 $g(\theta)$ の期待値は以下のように近似できます：

$\mathbb{E}{p(\theta | D)}[g(\theta)] \approx g(\mathbb{E}{q(\theta | \lambda)}[\theta]) + \frac{1}{2} \text{Tr}(\nabla^2 g(\mathbb{E}{q(\theta | \lambda)}[\theta]) \cdot \text{Var}{q(\theta | \lambda)}[\theta])$

この方法は、特に予測分布のモーメント（平均、分散など）を計算する際に有用です。ただし、関数 $g(\theta)$ が非常に非線形で、変分分布と真の事後分布の乖離が大きい場合には、近似精度が低下する可能性があります。

5.4 実践的なヒント：変分ベイズ法とMCMCの併用

実践的な観点から、変分ベイズ法とMCMCを相補的に使用することが有効な場合があります：

5.4.1 変分ベイズ法による初期化

MCMCの収束を加速するために、変分ベイズ法の結果を初期値として使用する方法があります：

変分ベイズ法で近似事後分布 $q(\theta | \lambda)$ を得る
この分布からサンプル $\theta^{(0)}$ を生成し、MCMCの初期値として使用
MCMCを実行して事後分布からのサンプルを得る

この方法は、変分分布が真の事後分布の高密度領域に近い場合、MCMCの混合時間を大幅に短縮できる可能性があります。

5.4.2 モデル選択のための変分ベイズ法

複数のモデルを比較する場合、変分ベイズ法で各モデルの周辺尤度（の下界）を迅速に計算し、有望なモデルを選択した後で、そのモデルに対してのみMCMCを適用するアプローチも有効です：

候補モデル $M_1, M_2, \ldots, M_K$ に対して変分ベイズ法を適用し、各ELBOを計算
最も高いELBOを持つモデル（または上位数モデル）を選択
選択したモデルに対してMCMCを適用し、より正確な事後分布を得る

このアプローチにより、計算リソースを効率的に使いながら、モデル選択と事後推論の両方を行うことができます。

5.4.3 分割統治アプローチ

大規模モデルでは、パラメータを複数のブロックに分割し、一部のブロックには変分ベイズ法を、他のブロックにはMCMCを適用するハイブリッドアプローチも考えられます：

パラメータ $\theta$ を $\theta_1$ と $\theta_2$ に分割
$\theta_1$ に対しては変分近似 $q(\theta_1)$ を使用
$\theta_2$ に対しては、 $\theta_1$ を固定した条件付き分布 $p(\theta_2 | \theta_1, D)$ からのMCMCサンプリングを行う
両者を交互に更新

この方法は、モデルの一部が解析的に扱いやすく、別の部分が複雑な構造を持つ場合に特に有効です。

6. 変分ベイズ法の応用例

変分ベイズ法は多くの実用的な問題に応用されています。ここでは主要な応用例を紹介します。

6.1 トピックモデリング

潜在ディリクレ配分法（Latent Dirichlet Allocation, LDA）などのトピックモデルは、変分ベイズ法の初期の成功例の一つです。LDAでは、文書はトピックの混合として、トピックは単語の分布として表現されます。
変分ベイズ法を用いたLDAでは、潜在変数（各単語のトピック割り当てなど）に対して平均場近似を適用します：

$q(\beta, \theta, z) = q(\beta)q(\theta)\prod_{d,n} q(z_{d,n})$

ここで $\beta$ はトピック-単語分布、 $\theta$ は文書-トピック分布、 $z_{d,n}$ はトピック割り当てです。大規模テキストコーパスに対しても効率的に適用できるため、トピックモデリングでは変分ベイズ法が標準的なアプローチとなっています。確率的変分推論（SVI）を用いることで、オンライン学習も可能です。

6.2 ベイズ ニューラルネットワーク

ニューラルネットワークのパラメータ（重みとバイアス）に事前分布を設定し、事後分布を推論するベイズニューラルネットワークでも、変分ベイズ法が広く使われています。
一般的なアプローチでは、各パラメータに対してガウス変分分布を使用します：

$q(w_{ij}) = \mathcal{N}(w_{ij} | \mu_{ij}, \sigma_{ij}^2)$

これにより、パラメータの不確実性を考慮した予測が可能になります。
自動微分変分推論（ADVI）やブラックボックス変分推論（BBVI）は、ベイズニューラルネットワークの学習に特に適しています。また、ドロップアウトはガウス変分推論の特殊ケースとみなせることも示されています（MC Dropout）。

6.3 状態空間モデル

時系列データの解析に使われる状態空間モデル（隠れマルコフモデルやカルマンフィルタなど）にも、変分ベイズ法が適用されています。
例えば、線形ガウス状態空間モデルでは：

$z_t = Az_{t-1} + \epsilon_t, \quad \epsilon_t \sim \mathcal{N}(0, Q)$

$y_t = Bz_t + \delta_t, \quad \delta_t \sim \mathcal{N}(0, R)$

潜在状態 $z_t$ と遷移パラメータ $A, B, Q, R$ の事後分布を変分ベイズ法で近似します。変分状態空間モデルは、平均場近似の限界を克服するために構造化変分推論（SVI）やノーマライジングフローなどの高度な手法が使われることが多いです。

最後に

変分ベイズ法は近似ベイズ推論の主要なアプローチとして、理論・アルゴリズム・応用の各面で急速に発展しています。特に、深層学習と組み合わせた変分推論の可能性はまだ十分に探求されていない領域です。
一方で、変分ベイズ法にはまだ多くの課題があることも事実です。真の事後分布との乖離をより正確に評価する方法、近似分布族の選択の体系化、高次元モデルにおける計算効率のさらなる向上、KLダイバージェンス以外の目的関数の探求など、理論的・実践的な研究課題が残されています。
また、MCMCと変分ベイズ法のハイブリッドアプローチや、変分ベイズ法をベースとした新たな近似手法の開発など、計算ベイズ推論の分野には多くの可能性があります。
変分ベイズ法は、ベイズ推論の計算的課題と実用的応用のバランスを取る上で、今後も中心的な役割を果たしていくでしょう。その理論的基盤の美しさと実践的有用性の両面から、統計学と機械学習の交差点に位置する魅力的な研究分野であり続けます。
この記事が、変分ベイズ法の背後にある考え方と具体的な実装手法の理解の一助となれば幸いです。実際の問題に適用する際には、選択した近似手法の仮定と制約を常に念頭に置き、結果を慎重に解釈することが重要です。状況に応じて、変分ベイズ法、MCMC、あるいはその他の手法を適切に選択・組み合わせることで、より強力なベイズ推論が可能になるでしょう。
最後までお読みいただきありがとうございました。とにかく書き上げるのが本当に大変でしたし、何もすることが無かったので良かったのですが、ゴールデンウイークがいつの間にか終わっていました...
　　
本当に最後にお詫びなのですが、はてなブログのはてなTeXの仕組みのせいで、一部数式が正しく表示できない問題が発生しています...色々試行錯誤して見たものの、修正ができませんでした。読みにくい中お読みいただきありがとうございました。