ELBO（変分下界）とは？対数尤度の下から支える不等式を導出する

手元に猫の写真が10万枚あるとして、「この画像たちを生み出した確率モデル」を学習したいとします。素直な方針は最尤推定です。モデルのパラメータ $\theta$ を、データの対数尤度 $\log p_\theta(x)$ が最大になるように選べばいい。ところが、画像のような複雑なデータをまともにモデル化しようとすると、必ず「潜在変数」$z$ — 目に見えない、画像の背後にある要因（猫の姿勢、毛色、照明…）— を導入したくなります。そして潜在変数を導入した瞬間、対数尤度は

$$ \log p_\theta(x) = \log \int p_\theta(x, z)\, dz $$

という積分の中身の対数になり、$z$ が高次元だと事実上計算できなくなります。100次元の潜在空間の積分を数値的に潰すのは、宇宙の年齢をかけても終わりません。

では諦めるしかないのでしょうか。ここで登場するのが ELBO（Evidence Lower Bound、証拠下界・変分下界） です。ELBO は「計算できない $\log p_\theta(x)$ を、下から支える計算可能な量」です。しかも、$\log p_\theta(x)$ と ELBO の差がちょうど「近似事後分布と真の事後分布の KL ダイバージェンス」になる、という美しい等式が成り立ちます。この等式ひとつで、「尤度を上げること」と「事後分布をうまく近似すること」が同じ最適化の表と裏であることがわかります。

この図が本記事の全体像です。赤い天井が「計算したいけれど計算できない対数尤度」、青い床が「計算できる ELBO」で、床の高さは自分で用意した近似分布 $q$ の出来によって決まります。雑な $q$ なら床は低い位置にとどまり、良い $q$ を選ぶほど床は天井のすぐ下まで持ち上がります。そして床と天井のすき間の正体こそ、この記事で正面から扱う「$q$ と真の事後分布の KL ダイバージェンス」です。

ELBO は、現代の機械学習・統計のあちこちで顔を出します。

• VAE（変分オートエンコーダ） — 学習の目的関数はそのまま ELBO です。「再構成誤差 + KL 正則化」という見慣れた損失関数は、ELBO を符号反転しただけのものです。
• 変分ベイズ / トピックモデル — LDA のようなベイズモデルの事後分布を、MCMC より遥かに速く近似する平均場変分推論は、ELBO を座標上昇法で最大化しています。
• EM アルゴリズム — 混合ガウスモデルの EM は、実は ELBO の交互最大化として統一的に理解できます。
• 拡散モデル — 拡散モデルの学習損失も、変分下界から導かれるノイズ予測誤差です。

この記事では、ELBO を「公式として覚える」のではなく、2通りの導出（イェンセンの不等式による下界化と、完全な等式による分解）を最後まで手を動かして追いかけ、ギャップの正体を突き止めるところまでやります。

本記事の内容

• 潜在変数モデルで対数尤度が計算できなくなる理由
• イェンセンの不等式による ELBO の導出（不等式の側から）
• $\log p(x) = \mathcal{L} + D_{\mathrm{KL}}(q \| p(z|x))$ という完全な等式の導出（等式の側から）
• ELBO を「再構成項 − 正則化項」に分解する
• ガウス事前・ガウス近似事後のときの KL の閉じた式の導出
• 再パラメータ化トリックによる勾配推定
• Python で ELBO とギャップ、平均場変分推論の単調増加、小さな VAE の学習を確認する

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• KLダイバージェンス完全ガイド — 確率分布間の情報量距離を直感で理解する
• ベイズ推論完全ガイド — 事前分布から事後分布へ、確率で学ぶ機械学習
• 【深層学習】VAEとは？VAEの理論を徹底的に分かりやすく解説する
• 機械学習で登場するイェンゼンの不等式を学ぶ

なぜ対数尤度は直接計算できないのか

まず「何が困っているのか」をはっきりさせましょう。困っている中身が曖昧なままだと、ELBO のありがたみが伝わりません。

袋の中に3種類のコイン（表が出る確率がそれぞれ違う）が入っていて、目をつぶって1枚引いて10回投げた、という状況を考えます。あなたが観測できるのは「表が7回、裏が3回」という結果 $x$ だけで、「どのコインを引いたか」$z$ は見えません。このとき $x$ の確率は、3通りの可能性をすべて足し合わせたもの

$$ p(x) = \sum_{k=1}^{3} p(z = k)\, p(x \mid z = k) $$

になります。$z$ が3通りしかないなら、足し算3回で終わりです。何の問題もありません。

問題は $z$ が連続かつ高次元になったときです。VAE の潜在変数は典型的に $z \in \mathbb{R}^{d}$（$d = 2 \sim 512$ 程度）で、和は積分に置き換わります。

$$ \begin{equation} p_\theta(x) = \int p_\theta(x \mid z)\, p(z)\, dz \end{equation} $$

この積分を周辺化（marginalization） と呼び、$p_\theta(x)$ をエビデンス（evidence、証拠）、あるいは周辺尤度と呼びます。ELBO の “E” はこの Evidence の E です。

なぜこれが計算できないのでしょうか。理由は2つあります。

理由1: 次元の呪い。$d$ 次元の積分を格子で近似しようとすると、1次元あたり $M$ 点刻んだだけで $M^d$ 回の評価が必要です。$M = 10$、$d = 20$ でも $10^{20}$ 回。これは現実的ではありません。

理由2: 素朴なモンテカルロが効かない。「じゃあ事前分布 $p(z)$ から $z$ をサンプリングして平均を取ればいいのでは」と思うかもしれません。実際、

$$ p_\theta(x) \approx \frac{1}{S}\sum_{s=1}^{S} p_\theta(x \mid z_s), \quad z_s \sim p(z) $$

は不偏推定量です。しかしこれは実用になりません。事前分布からランダムに引いた $z_s$ のほとんどは、目の前の $x$ とまるで無関係な潜在変数です。猫の画像 $x$ に対して、ランダムに引いた $z_s$ が「その猫を生成しうる潜在コード」である確率はほぼゼロで、$p_\theta(x \mid z_s)$ は天文学的に小さい値ばかりになります。ごくまれに当たりを引いたときだけ巨大な値が出るので、推定量の分散が爆発します。当たりくじが $10^{-30}$ の確率でしか出ないくじ引きの期待値を、100回引いて推定するようなものです。

実際に $p(z) = \mathcal{N}(0, I_d)$、$p(x\mid z) = \mathcal{N}(x; z, 0.5 I_d)$ というモデルで、事前分布から1万サンプル引いて $\log p(x)$ を推定した結果です。左のグラフでは $d \le 10$ ならほぼ正しく推定できているのに、$d = 20$ で約 $-3.6$、$d = 40$ では約 $-19$ も過小評価しています。右のグラフがその理由で、実際に推定に効いているサンプルの割合（有効サンプル率）は $d = 1$ で 41% あったものが $d = 40$ では $1.7\times10^{-4}$、つまり1万サンプル引いても実質2本弱しか使えていません。残りの9998本は「その $x$ とまるで関係のない $z$」だったわけです。

だとすれば、やるべきことははっきりしています。「$x$ を説明できそうな $z$ の周辺」を狙って重点的にサンプリングしたい。その「狙いを定める分布」こそが、これから導入する 変分分布 $q(z \mid x)$ です。理想的には $q(z\mid x)$ を真の事後分布 $p_\theta(z \mid x)$ に一致させたいのですが、それも計算できません。なぜなら

$$ p_\theta(z \mid x) = \frac{p_\theta(x, z)}{p_\theta(x)} $$

の分母に、まさに求めたかった $p_\theta(x)$ がいるからです。「尤度が計算できないのは事後分布が計算できないから、事後分布が計算できないのは尤度が計算できないから」という循環にはまっています。

ELBO は、この循環をきれいに断ち切ります。次の節で、まず結論から見ていきましょう。

ELBOとは — 30秒でわかる要点

先に答えを書きます。ELBO とは、計算できない対数尤度 $\log p_\theta(x)$ を下から押さえる「床」であり、その床の高さは、あなたが用意した近似分布 $q(z\mid x)$ の出来次第で上下する量です。

$$ \begin{equation} \log p_\theta(x) \;\ge\; \mathcal{L}(q, \theta) \;=\; \mathbb{E}_{q(z\mid x)}\!\left[\log \frac{p_\theta(x, z)}{q(z \mid x)}\right] \end{equation} $$

この床がありがたいのは、次の3点が同時に成り立つからです。

3番目が肝心です。ELBO をひとつ最大化するだけで、「モデルの学習」と「事後分布の推論」が同時に進みます。下界を持ち上げれば天井（対数尤度）も一緒に上がり、下界を天井にぴったり貼り付けるほど事後分布の近似が良くなる。この二重の意味があるからこそ、ELBO は変分推論と深層生成モデルの共通言語になっています。

ちなみに $-\mathcal{L}$ は変分自由エネルギーとも呼ばれます。後で見るように $-\mathcal{L} = \mathbb{E}_q[-\log p(x,z)] – H[q]$ と書けて、これは統計力学の自由エネルギー $F = U – TS$（エネルギー引く温度かけるエントロピー）とまったく同じ形をしています。ELBO 最大化は自由エネルギー最小化であり、「エネルギーの低い状態に寄せつつ、エントロピー（分布の広がり）は保て」という物理でおなじみの綱引きなのです。

では、この不等式がどこから出てくるのか。まずは古典的な、イェンセンの不等式による導き方から見ていきます。

導出その1: イェンセンの不等式で下から支える

計算を始める前に、なぜ「不等式」が登場するのかを直感で押さえておきます。

私たちが計算したいのは $\log(\text{和})$ の形、つまり「たくさん足してから対数を取る」です。一方、計算しやすいのは $\text{和}(\log)$ の形、「対数を取ってから足す」です。この2つは一般に一致しませんが、$\log$ が上に凸（凹関数）であるおかげで、大小関係だけは常に決まっています。凹関数では「平均を取ってから関数に入れた値」のほうが「関数に入れてから平均を取った値」より大きい。これがイェンセンの不等式です。

$$ \begin{equation} \log \mathbb{E}[Y] \;\ge\; \mathbb{E}[\log Y] \end{equation} $$

グラフで言えば、凹関数の曲線は、その上の2点を結ぶ弦より必ず上にあります。「先に平均してから曲線に乗せる」のが弦の話ではなく曲線の話になるので、有利になる、というわけです。

$Y$ が確率 0.55 で 0.35、確率 0.45 で 3.2 を取る場合を描いたものです。平均 $\mathbb{E}[Y] = 1.633$ の位置で縦に見ると、曲線上の赤い点（$\log \mathbb{E}[Y] = 0.490$）が弦上のオレンジの点（$\mathbb{E}[\log Y] = -0.054$）より 0.544 だけ高い位置にあります。この高さの差がそのまま ELBO のギャップになるので、「$Y$ のばらつきが大きいほど弦は曲線から離れ、下界は緩くなる」という関係が見て取れます。逆に $Y$ が一点に集中していれば2点は重なり、すき間はゼロになります。

この道具を使って、$\log p_\theta(x)$ を変形します。ゴールは「計算できない $\log \int$ を、期待値の形に持ち込む」ことです。

まず周辺化の定義から出発します。

$$ \log p_\theta(x) = \log \int p_\theta(x, z)\, dz $$

ここで恒等的な仕掛けを入れます。任意の確率密度 $q(z \mid x) > 0$ を持ってきて、被積分関数に $\frac{q(z\mid x)}{q(z\mid x)} = 1$ を掛けます。何も変えていませんが、「$q$ で割ったものを $q$ で期待値を取る」形が現れます。

$$ \log p_\theta(x) = \log \int q(z \mid x)\, \frac{p_\theta(x, z)}{q(z \mid x)}\, dz $$

$\int q(z\mid x)\, f(z)\, dz = \mathbb{E}_{q}[f(z)]$ という期待値の定義を使うと、右辺は期待値そのものです。

$$ \log p_\theta(x) = \log \mathbb{E}_{q(z\mid x)}\!\left[\frac{p_\theta(x, z)}{q(z \mid x)}\right] $$

ここで $Y = \dfrac{p_\theta(x,z)}{q(z\mid x)}$ と置いてイェンセンの不等式 $\log \mathbb{E}[Y] \ge \mathbb{E}[\log Y]$ を適用すると、$\log$ が期待値の内側に潜り込みます。

$$ \begin{equation} \log p_\theta(x) \;\ge\; \mathbb{E}_{q(z\mid x)}\!\left[\log \frac{p_\theta(x, z)}{q(z \mid x)}\right] \;=\; \mathcal{L}(q, \theta) \end{equation} $$

これで ELBO が出ました。導出自体は3行ですが、中身は深いです。ポイントを2つ確認しましょう。

ポイント1: $q$ は何でもよい。上の導出で $q(z\mid x)$ に課した条件は「$p_\theta(x,z) > 0$ となる領域で $q(z\mid x) > 0$」（ゼロ割りを避けるため）だけです。$q$ をどんな分布に選んでも不等式は成り立ちます。だから $q$ を自由に動かして、下界をできるだけ高く持ち上げられるわけです。

ポイント2: 不偏推定量に $\log$ をかぶせると下向きにバイアスが乗る。$\hat{p} = \frac{p_\theta(x,z)}{q(z\mid x)}$（$z \sim q$）は $\mathbb{E}[\hat p] = p_\theta(x)$ を満たす、エビデンスの不偏推定量です。ところが $\log \hat p$ の期待値は $\log p_\theta(x)$ より小さい。「不偏推定量の対数は、対数の不偏推定量ではない」——ELBO のギャップは、この $\log$ の非線形性が生む系統的なズレなのです。

この見方を推し進めると、$S$ 個のサンプルで平均してから対数を取る

$$ \mathcal{L}_S = \mathbb{E}\!\left[\log \frac{1}{S}\sum_{s=1}^{S}\frac{p_\theta(x, z_s)}{q(z_s\mid x)}\right], \quad z_s \stackrel{\text{iid}}{\sim} q(z\mid x) $$

という一般化が得られます。$S$ を増やすほど平均の分散が減り、イェンセンのギャップが縮むので $\mathcal{L}_1 \le \mathcal{L}_2 \le \cdots \le \log p_\theta(x)$ と単調に締まっていきます。$S=1$ が通常の ELBO で、$S>1$ にしたものが IWAE（Importance Weighted Autoencoder） の目的関数です。「ELBO は多サンプル下界の $S=1$ の場合にすぎない」という位置づけが見えると、視界が広がります。

後の節で使う1次元モデル（$p(z)=\mathcal{N}(0,1)$、$p(x\mid z)=\mathcal{N}(x;z,0.5)$、$x=1.5$）で、$q$ をわざと事前分布 $\mathcal{N}(0,1)$ のまま固定して $\mathcal{L}_S$ をモンテカルロ推定した結果です。$S=1$ では $-3.813$ と天井 $-1.872$ から大きく離れているのに、$S=8$ で $-1.978$、$S=128$ では $-1.878$ とほぼ天井に貼り付きます。$q$ の質を1ミリも改善していないのに、サンプルを増やすだけで下界が締まるのがポイントで、これは「イェンセンのすき間が平均の分散とともに縮む」ことの直接的な現れです。

ここまでで下界が得られました。しかし、この導出ではギャップがどれだけあるのかがわかりません。「$\ge$」の記号は差の大きさを教えてくれないからです。次の節では、不等式ではなく完全な等式から出発して、ギャップの正体を突き止めます。

導出その2: ギャップの正体を暴く完全な等式

イェンセンの不等式は便利ですが、「どれだけ損をしているか」を隠してしまいます。同じ ELBO を、今度は一切の不等式を使わずに導いてみましょう。驚くべきことに、ギャップは完全に明示的な形で書けます。

出発点は、拍子抜けするほど単純な事実です。$\log p_\theta(x)$ は $z$ を含みません。$z$ を含まない量を、$z$ についての確率分布 $q(z\mid x)$ で期待値を取っても、値は変わりません。

$$ \log p_\theta(x) = \int q(z\mid x)\, \log p_\theta(x)\, dz = \mathbb{E}_{q(z\mid x)}\big[\log p_\theta(x)\big] $$

「定数の期待値は定数」という当たり前の操作ですが、これで $\log p_\theta(x)$ を $\mathbb{E}_q[\cdot]$ の中に閉じ込められました。ここから中身をいじっていきます。

ベイズの定理 $p_\theta(z \mid x) = \dfrac{p_\theta(x,z)}{p_\theta(x)}$ を $p_\theta(x)$ について解くと $p_\theta(x) = \dfrac{p_\theta(x,z)}{p_\theta(z\mid x)}$ です。これを期待値の中に代入します。

$$ \log p_\theta(x) = \mathbb{E}_{q}\!\left[\log \frac{p_\theta(x, z)}{p_\theta(z \mid x)}\right] $$

次に、先ほどと同じ仕掛けで $\frac{q(z\mid x)}{q(z\mid x)} = 1$ を掛け込みます。今度は分子と分母それぞれに $q$ を配るのがコツです。

$$ \log p_\theta(x) = \mathbb{E}_{q}\!\left[\log \left( \frac{p_\theta(x, z)}{q(z\mid x)} \cdot \frac{q(z\mid x)}{p_\theta(z \mid x)} \right)\right] $$

積の対数は対数の和なので、期待値の線形性を使って2つに分けます。

$$ \log p_\theta(x) = \underbrace{\mathbb{E}_{q}\!\left[\log \frac{p_\theta(x, z)}{q(z\mid x)}\right]}_{= \;\mathcal{L}(q,\theta)\;（\text{ELBO}）} + \underbrace{\mathbb{E}_{q}\!\left[\log \frac{q(z\mid x)}{p_\theta(z \mid x)}\right]}_{= \;D_{\mathrm{KL}}(q(z\mid x)\,\|\,p_\theta(z\mid x))} $$

第1項は前節で得た ELBO そのもの、第2項は KL ダイバージェンスの定義そのものです。つまり、

$$ \begin{equation} \log p_\theta(x) = \mathcal{L}(q, \theta) + D_{\mathrm{KL}}\big(q(z \mid x)\,\big\|\,p_\theta(z \mid x)\big) \end{equation} $$

これが変分推論の中心等式です。不等号ではなく等号であることに注目してください。KL ダイバージェンスは常に $D_{\mathrm{KL}} \ge 0$ で、等号成立は $q(z\mid x) = p_\theta(z\mid x)$ のときだけ（ギブスの不等式）。したがって自動的に $\log p_\theta(x) \ge \mathcal{L}$ が従い、前節のイェンセンの結果を含んでいます。

この等式が本当に成り立つことを、先取りして数値で見ておきます。左は5つの近似分布 $q$ について「青い ELBO」と「オレンジのすき間」を積み上げたもので、どの $q$ でも合計が赤い線（対数尤度 $-1.8717$）にぴったり届いています。右はさらに強い確認で、横軸が「$\log p(x)$ から ELBO を引いた値」、縦軸が「$q$ と真の事後分布の KL をガウスの公式で直接計算した値」です。5点すべてが傾き1の直線に乗っており、まったく別経路で計算した2つの量が一致することがわかります。

この等式を3つの角度から読み解きましょう。

読み方1: ギャップの正体は事後分布の近似誤差。$\log p_\theta(x)$ は $q$ を含みません。つまり $q$ をどう動かしても左辺は動かない。すると右辺の2項はシーソーの関係になります。ELBO を上げれば KL は必ず同じだけ下がる。逆も然り。だから、

$$ \arg\max_{q} \mathcal{L}(q,\theta) = \arg\min_{q} D_{\mathrm{KL}}\big(q(z\mid x)\,\|\,p_\theta(z\mid x)\big) $$

が成り立ちます。計算できない KL を最小化する問題が、計算できる ELBO を最大化する問題にすり替わった——これが変分推論という手法の全体です。前節で述べた「循環」は、こうして断ち切られます。

読み方2: 下界の締まり具合が近似の良さの通信簿。ギャップが小さければ、$q$ が真の事後分布をよく捉えているということです。逆にギャップが大きいまま ELBO を最大化しても、$\log p_\theta(x)$ が本当に上がっているかはわかりません。ELBO を対数尤度の代理として報告する論文が多いですが、それは常に過小評価であり、しかも過小評価の量はモデルごとに違う、という点は覚えておく価値があります。

読み方3: EM アルゴリズムとの関係。$\theta$ と $q$ を交互に最適化することを考えます。$q$ について最大化するステップで、もし $q$ に制約がなければ最適解は $q = p_\theta(z\mid x)$ となり、ギャップがゼロになって下界が対数尤度に接します。次に $\theta$ について ELBO を最大化する。この2ステップこそ、EMアルゴリズム の E ステップと M ステップに他なりません。「EM は ELBO の交互最大化」という統一的な見方は、混合ガウスの EM を丸暗記から解放してくれます。

ここまでで、ELBO が何であり、何を測っているかがわかりました。次は、この ELBO を実際に計算・実装できる形に書き換えます。

ELBOの2つの顔: 再構成項と正則化項

ELBO の定義式 $\mathcal{L} = \mathbb{E}_q\left[\log \frac{p(x,z)}{q(z\mid x)}\right]$ は美しいのですが、このままではコードに落としにくい。同時分布 $p(x,z)$ を、モデルの部品である尤度 $p_\theta(x\mid z)$ と事前分布 $p(z)$ に分解しましょう。

生成モデルは通常「まず潜在変数を引いて、それから観測を生成する」という因数分解を持ちます。

$$ p_\theta(x, z) = p_\theta(x \mid z)\, p(z) $$

これを ELBO に代入し、対数を分配します。

$$ \mathcal{L} = \mathbb{E}_{q}\big[\log p_\theta(x\mid z) + \log p(z) – \log q(z\mid x)\big] $$

期待値の線形性で第1項を切り離し、残った2項をまとめると、KL の定義がふたたび現れます。

$$ \mathcal{L} = \mathbb{E}_{q}\big[\log p_\theta(x\mid z)\big] – \mathbb{E}_{q}\!\left[\log \frac{q(z\mid x)}{p(z)}\right] $$

$$ \begin{equation} \mathcal{L}(q,\theta) = \underbrace{\mathbb{E}_{q(z\mid x)}\big[\log p_\theta(x \mid z)\big]}_{\text{再構成項}} – \underbrace{D_{\mathrm{KL}}\big(q(z\mid x)\,\|\,p(z)\big)}_{\text{正則化項}} \end{equation} $$

これが実装で最もよく使われる形です。2つの項の意味を、VAE の文脈で読み解いてみましょう。

再構成項 $\mathbb{E}_q[\log p_\theta(x\mid z)]$ は、「$q$ が提案した潜在コード $z$ から、元の $x$ をどれだけうまく復元できるか」を測ります。$p_\theta(x\mid z)$ をガウス分布にすれば $-\|x – \hat{x}(z)\|^2$ に比例する項（平均二乗誤差）、ベルヌーイ分布にすればバイナリクロスエントロピーになります。「エンコードしてデコードして、元に戻るか」というオートエンコーダの発想はここから出てくるわけです。

正則化項 $D_{\mathrm{KL}}(q(z\mid x)\|p(z))$ は、「エンコーダが出す潜在分布が、事前分布からどれだけ離れているか」のペナルティです。これがないと、エンコーダは各データ点をバラバラの遠い場所に、限りなく尖ったデルタ関数として配置してしまいます。それでは再構成は完璧でも、潜在空間はスカスカの穴だらけになり、$p(z)$ から $z$ を引いて画像を生成する、という使い方ができません。KL 項は「潜在空間を事前分布の形に押し込め、隙間を埋めさせる圧力」として働きます。

この2項の綱引きが VAE の本質です。再構成項だけなら普通のオートエンコーダ、KL 項が勝ちすぎると $q(z\mid x)$ が $x$ に依存しなくなって（＝事後分布の崩壊、posterior collapse）、デコーダは平均的な画像しか吐かなくなります。ちなみに KL 項に係数 $\beta$ を掛けて綱引きの強さを人為的に調整したものが $\beta$-VAE で、$\beta$ を大きくすると潜在次元が独立な意味を持ちやすくなる（disentanglement）ことが知られています。

もうひとつ、覚えておくと便利な変形があります。$-\log q$ の期待値は $q$ の微分エントロピー $H[q]$ なので、

$$ \begin{equation} \mathcal{L} = \underbrace{\mathbb{E}_{q}\big[\log p_\theta(x, z)\big]}_{\text{エネルギー項}} + \underbrace{H\big[q(z\mid x)\big]}_{\text{エントロピー項}} \end{equation} $$

とも書けます。第1項は「同時確率の高い場所に $q$ を集めろ」、第2項は「でも一点に潰れるな、広がっていろ」という要求です。尖りすぎるとエントロピーで罰せられ、散らばりすぎるとエネルギーで罰せられる。冒頭で触れた自由エネルギーの形 $F = U – TS$ そのものであり、変分推論が物理から名前を借りている理由がここにあります。

3つの表現をまとめておきます。どれも同じ量で、場面によって使い分けます。

さて、再構成 − 正則化の形が実用的なのは、KL 項を解析的に計算できる場合があるからです。次の節で、その代表例であるガウス同士の KL を、最後まで手で導いてみます。

ガウス同士のKLを閉じた式にする

サンプリングで期待値を近似すると、どうしてもノイズが乗ります。もし ELBO の一部でも解析的に計算できるなら、勾配の分散が減って学習が安定します。VAE で事前分布を標準正規、近似事後を対角ガウスに取るのは、まさに KL 項が閉じた式で書けるからです。

まずは1次元で、$q = \mathcal{N}(\mu, \sigma^2)$、$p = \mathcal{N}(0, 1)$ の KL を計算します。定義は

$$ D_{\mathrm{KL}}(q\|p) = \mathbb{E}_{q}\big[\log q(z) – \log p(z)\big] $$

なので、まず被期待値の中身を書き下します。

$$ \log q(z) = -\frac{1}{2}\log(2\pi\sigma^2) – \frac{(z-\mu)^2}{2\sigma^2}, \qquad \log p(z) = -\frac{1}{2}\log(2\pi) – \frac{z^2}{2} $$

差を取ると $\log(2\pi)$ が打ち消し合い、$\log \sigma^2$ だけが残ります。

$$ \log q(z) – \log p(z) = -\frac{1}{2}\log \sigma^2 – \frac{(z-\mu)^2}{2\sigma^2} + \frac{z^2}{2} $$

ここで $z \sim q = \mathcal{N}(\mu,\sigma^2)$ での期待値を取ります。必要なのは次の2つのモーメントだけです。第1に、分散の定義から $\mathbb{E}_q[(z-\mu)^2] = \sigma^2$。第2に、$\mathbb{E}_q[z^2] = \mathrm{Var}[z] + (\mathbb{E}[z])^2 = \sigma^2 + \mu^2$。これらを代入すると、

$$ D_{\mathrm{KL}}(q\|p) = -\frac{1}{2}\log\sigma^2 – \frac{\sigma^2}{2\sigma^2} + \frac{\sigma^2 + \mu^2}{2} $$

第2項の $\sigma^2/\sigma^2$ が 1 に約分されるので、

$$ \begin{equation} D_{\mathrm{KL}}\big(\mathcal{N}(\mu,\sigma^2)\,\|\,\mathcal{N}(0,1)\big) = \frac{1}{2}\left(\mu^2 + \sigma^2 – 1 – \log \sigma^2\right) \end{equation} $$

きれいな形になりました。式の妥当性を目で確かめておきましょう。$\mu = 0,\ \sigma^2 = 1$ を入れると $\frac12(0 + 1 – 1 – 0) = 0$。分布が一致するとき KL がゼロになる、という当然の性質が満たされています。また $\sigma^2 \to 0$（分布が一点に潰れる）で $-\log\sigma^2 \to +\infty$、$\sigma^2 \to \infty$（広がりすぎ）でも $\sigma^2 \to \infty$ で発散します。尖りすぎも広がりすぎも罰せられるという、前節のエントロピー項の議論と同じ構造がここに現れています。

左は導いた式 $\frac12(\mu^2+\sigma^2-1-\log\sigma^2)$ を $(\mu, \sigma)$ 平面に描いたもので、谷底（白い部分）はただ一点 $(0, 1)$、すなわち事前分布と完全に一致する場所だけです。VAE の KL 項が「潜在分布を標準正規に引き寄せる力」として働く様子が、そのまま地形の形になっています。右は $\mu = 0$ に固定した断面で、KL がオレンジの実線＝緑の破線（尖りを罰する $-\frac12\log\sigma^2$）＋青の点線（広がりを罰する $\frac12\sigma^2 – \frac12$）に分解できることを示しています。片方は $\sigma \to 0$ で、もう片方は $\sigma \to \infty$ で発散し、和がちょうど $\sigma = 1$ でゼロになるという綱引きの構造が見て取れます。

$d$ 次元で共分散が対角行列 $\mathrm{diag}(\sigma_1^2,\dots,\sigma_d^2)$ のときは、各次元が独立なので KL は単純に足し算になります（独立な分布同士の KL は各成分の KL の和）。

$$ \begin{equation} D_{\mathrm{KL}}\big(\mathcal{N}(\bm{\mu}, \mathrm{diag}(\bm{\sigma}^2))\,\|\,\mathcal{N}(\bm{0}, \bm{I})\big) = \frac{1}{2}\sum_{j=1}^{d}\left(\mu_j^2 + \sigma_j^2 – 1 – \log \sigma_j^2\right) \end{equation} $$

VAE の実装で見かける

import torch
mu, logvar = torch.zeros(4, 2), torch.zeros(4, 2)   # エンコーダの出力（例: 事前分布と一致する場合）
kl = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())
print(kl.item())   # -> -0.0（符号付きゼロ。q が事前分布と一致するので KL はゼロ）

という1行は、まさにこの式です。$\log \sigma_j^2$ をネットワークに直接出力させる（logvar）のは、$\sigma_j^2 > 0$ の制約を自動で満たしつつ、$\log$ を取り直す数値誤差を避けるためです。

一般の多変量ガウス同士でも閉じた式があります。$q = \mathcal{N}(\bm{\mu}_1, \bm{\Sigma}_1)$、$p = \mathcal{N}(\bm{\mu}_2, \bm{\Sigma}_2)$ のとき、

$$ D_{\mathrm{KL}}(q\|p) = \frac{1}{2}\left[\operatorname{tr}(\bm{\Sigma}_2^{-1}\bm{\Sigma}_1) + (\bm{\mu}_2 – \bm{\mu}_1)^{\top}\bm{\Sigma}_2^{-1}(\bm{\mu}_2 – \bm{\mu}_1) – d + \ln\frac{|\bm{\Sigma}_2|}{|\bm{\Sigma}_1|}\right] $$

となります。$\bm{\mu}_2 = \bm{0}$、$\bm{\Sigma}_2 = \bm{I}$、$\bm{\Sigma}_1$ が対角、と置けば先ほどの式に一致することを確かめてみてください。

KL が閉じた式になったので、ELBO のうち残る不確定要素は再構成項の期待値だけです。それを具体的な数値で追いかけるために、次はすべてが解析的に計算できる1次元モデルを持ち出します。

具体例: 1次元共役ガウスで数値を追う

抽象的な等式を、数字で殴って納得しましょう。ここでは対数尤度も、真の事後分布も、ELBO も、ギャップも、すべて紙の上で計算できるモデルを使います。

設定は次の通りです。潜在変数 $z$ を標準正規から引き、それをノイズ付きで観測する、というだけのモデルです。

$$ p(z) = \mathcal{N}(z; 0, 1), \qquad p(x \mid z) = \mathcal{N}(x; z, \sigma_n^2) $$

観測ノイズの分散を $\sigma_n^2 = 0.5$、観測値を $x = 1.5$ とします。

まず真の答え（対数尤度）。$x = z + \varepsilon$（$z \sim \mathcal{N}(0,1)$、$\varepsilon\sim\mathcal{N}(0,\sigma_n^2)$、独立）なので、$x$ 自身は分散が足し算になったガウス分布 $\mathcal{N}(0, 1 + \sigma_n^2)$ に従います。

$$ \log p(x) = -\frac{1}{2}\log\big(2\pi(1+\sigma_n^2)\big) – \frac{x^2}{2(1+\sigma_n^2)} $$

$\sigma_n^2 = 0.5$、$x = 1.5$ を代入すると $\log p(x) = -1.8717$ です。

次に真の事後分布。ガウス同士の共役性から $p(z\mid x)$ もガウスです。精度（分散の逆数）が足し算になることを使うと、事後精度は $1 + 1/\sigma_n^2 = 1 + 2 = 3$、つまり事後分散は $s_*^2 = \frac{\sigma_n^2}{1+\sigma_n^2} = \frac{0.5}{1.5} = 0.3333$。事後平均は $\mu_* = \frac{x}{1+\sigma_n^2} = \frac{1.5}{1.5} = 1.0$ です。事前分布（平均0）に引っ張られて、観測値 $1.5$ より原点寄りの $1.0$ に落ち着いているのが、ベイズらしい振る舞いです。

左の図が、いま計算した3つの分布です。尤度（緑の破線）のピークは観測値と同じ $z = 1.5$ にありますが、事前分布（灰色）に引き戻された結果、真の事後分布（赤）のピークは $z = 1.0$ まで戻され、幅も両者より狭くなっています。右の図はこれから比べる4つの近似分布 $q$ を、目標である真の事後分布（太い赤）と重ねたもので、形が真の事後分布に近い $q$ ほど凡例の ELBO が高く、すき間が小さいという対応が読み取れます。事前分布のまま（灰色の破線）は位置も幅も外しているので、すき間が最大の 1.951 になっています。

そして ELBO。近似分布を $q = \mathcal{N}(m, s^2)$ と置きます。再構成項は

$$ \mathbb{E}_q[\log p(x\mid z)] = -\frac{1}{2}\log(2\pi\sigma_n^2) – \frac{\mathbb{E}_q[(x-z)^2]}{2\sigma_n^2} = -\frac{1}{2}\log(2\pi\sigma_n^2) – \frac{(x-m)^2 + s^2}{2\sigma_n^2} $$

です（$\mathbb{E}_q[(x-z)^2] = (x-m)^2 + s^2$ は、平均まわりに展開して交差項が消えることから出ます）。正則化項は前節の公式そのままで $\frac12(m^2 + s^2 – 1 – \log s^2)$。

いくつかの $q$ で数値を並べてみます。

3列目（ELBO と対数尤度の差）と4列目（$q$ と真の事後分布の KL を直接計算した値）が、小数点以下4桁まで完全に一致しています。これが前々節で導いた等式 $\log p(x) = \mathcal{L} + D_{\mathrm{KL}}(q\|p(z\mid x))$ の、数値による確認です。

表からもうひとつ読み取れることがあります。$q$ を真の事後 $\mathcal{N}(1.0, 0.3333)$ に近づけるほど ELBO が上がり、ぴったり一致したところで ELBO $= \log p(x) = -1.8717$ となって下界が天井に接触します。逆に $q$ を事前分布 $\mathcal{N}(0,1)$ のまま放置すると、ギャップは 1.95 も開いてしまいます。「$q$ の設計が下界の締まり具合を決める」という主張が、具体的な数字で確認できました。

表の5行を、$(m, s)$ 平面全体に広げて等高線で描いたのがこの図です。地形はただ一つの山になっていて、山頂（赤い星）は $(m, s) = (1.00, 0.577)$、まさに真の事後分布の位置です。そして山頂の高さはカラーバー上端の $-1.872$、つまり対数尤度そのもので、どこをどう探してもこれより高い場所はありません。表に並べた4つの $q$（白い丸）は山の斜面に散らばっており、山頂から遠いほどすき間が大きい、という関係がそのまま距離として見えています。

この地形を断面で切って、ELBO の中身の綱引きまで見てみましょう。

Python実装1: ELBO・対数尤度・KLギャップを同時に描く

上の1次元モデルで、近似事後の平均 $m$ と標準偏差 $s$ を動かしたときに、ELBO・対数尤度（定数）・ギャップがどう動くかをプロットします。狙いは「ELBO が下から天井を目指して持ち上がる」様子を目で見ることです。

import numpy as np
import matplotlib, matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

# --- モデル: p(z)=N(0,1), p(x|z)=N(x; z, sn2) ---
sn2, x = 0.5, 1.5

# 真の対数尤度（x ~ N(0, 1+sn2)）
log_px = -0.5*np.log(2*np.pi*(1+sn2)) - x**2/(2*(1+sn2))
# 真の事後分布 N(mu_star, s_star2)
s_star2 = sn2/(1+sn2)
mu_star = x/(1+sn2)

def elbo(m, s):
    """再構成項 - 正則化項 の形でELBOを計算"""
    rec = -0.5*np.log(2*np.pi*sn2) - ((x-m)**2 + s**2)/(2*sn2)  # E_q[log p(x|z)]
    kl  = 0.5*(m**2 + s**2 - 1 - np.log(s**2))                  # KL(q || p(z))
    return rec - kl, rec, kl

def kl_to_posterior(m, s):
    """KL(q || p(z|x)) をガウス同士の閉じた式で直接計算"""
    return np.log(np.sqrt(s_star2)/s) + (s**2 + (m-mu_star)**2)/(2*s_star2) - 0.5

print(f"log p(x) = {log_px:.4f},  真の事後 = N({mu_star:.4f}, {s_star2:.4f})")
for m, s in [(0.0,1.0), (1.2,1.3), (0.5,0.4), (1.0,0.3), (mu_star,np.sqrt(s_star2))]:
    L, rec, kl = elbo(m, s)
    print(f"q=N({m:5.3f},{s:5.3f}^2)  ELBO={L:8.4f}  ギャップ={log_px-L:7.4f}"
          f"  KL(q||事後)={kl_to_posterior(m,s):7.4f}")

出力は次の通りです。

log p(x) = -1.8717,  真の事後 = N(1.0000, 0.3333)
q=N(0.000,1.000^2)  ELBO= -3.8224  ギャップ= 1.9507  KL(q||事後)= 1.9507
q=N(1.200,1.300^2)  ELBO= -3.1550  ギャップ= 1.2833  KL(q||事後)= 1.2833
q=N(0.500,0.400^2)  ELBO= -2.3537  ギャップ= 0.4820  KL(q||事後)= 0.4820
q=N(1.000,0.300^2)  ELBO= -2.1613  ギャップ= 0.2897  KL(q||事後)= 0.2897
q=N(1.000,0.577^2)  ELBO= -1.8717  ギャップ= 0.0000  KL(q||事後)= 0.0000

「ELBO と対数尤度の差」を引き算で求めた値と、「$q$ と真の事後分布の KL」をガウスの公式で独立に計算した値が、すべての行で一致しています。2つはまったく別の経路で計算しているのに同じ数字になる——これは中心等式が正しいことの、かなり強い実験的証拠です。最終行では ELBO がちょうど $\log p(x)$ に到達しており、下界が天井に接することも確認できました。

次に、$m$ と $s$ をスキャンして地形を描きます。

# --- (左) 平均 m をスキャン、(右) 標準偏差 s をスキャン ---
fig, axes = plt.subplots(1, 3, figsize=(15, 4.2))

ms = np.linspace(-1.0, 3.0, 400)
L_m = np.array([elbo(m, np.sqrt(s_star2))[0] for m in ms])
axes[0].axhline(log_px, color="crimson", lw=2, label="対数尤度 log p(x)（天井）")
axes[0].plot(ms, L_m, color="royalblue", lw=2, label="ELBO（s は最適値に固定）")
axes[0].fill_between(ms, L_m, log_px, color="orange", alpha=0.25, label="ギャップ = KL(q||事後)")
axes[0].axvline(mu_star, ls="--", color="gray"); axes[0].set_ylim(-4.2, -1.5)
axes[0].set_xlabel("近似事後の平均 $m$"); axes[0].set_ylabel("値")
axes[0].set_title("平均を動かすとELBOは天井に接する"); axes[0].legend(fontsize=8)

ss = np.linspace(0.15, 2.0, 400)
L_s = np.array([elbo(mu_star, s)[0] for s in ss])
axes[1].axhline(log_px, color="crimson", lw=2, label="対数尤度 log p(x)")
axes[1].plot(ss, L_s, color="royalblue", lw=2, label="ELBO（m は最適値に固定）")
axes[1].fill_between(ss, L_s, log_px, color="orange", alpha=0.25, label="ギャップ")
axes[1].axvline(np.sqrt(s_star2), ls="--", color="gray"); axes[1].set_ylim(-4.2, -1.5)
axes[1].set_xlabel("近似事後の標準偏差 $s$"); axes[1].set_title("尖りすぎても広がりすぎても損をする")
axes[1].legend(fontsize=8)

# 再構成項と正則化項の綱引き
rec_s = np.array([elbo(mu_star, s)[1] for s in ss])
kl_s  = np.array([elbo(mu_star, s)[2] for s in ss])
axes[2].plot(ss, rec_s, color="seagreen", lw=2, label="再構成項 $E_q[\\log p(x|z)]$")
axes[2].plot(ss, -kl_s, color="darkorange", lw=2, label="$-$正則化項 $-KL(q||p(z))$")
axes[2].plot(ss, rec_s-kl_s, color="royalblue", lw=2.5, label="ELBO（両者の和）")
axes[2].axvline(np.sqrt(s_star2), ls="--", color="gray"); axes[2].set_ylim(-6, 1)
axes[2].set_xlabel("近似事後の標準偏差 $s$"); axes[2].set_title("再構成項と正則化項の綱引き")
axes[2].legend(fontsize=8)

for ax in axes: ax.grid(alpha=0.3)
plt.tight_layout(); plt.show()

3枚のグラフから読み取れることを整理します。左：平均 $m$ を動かすと ELBO は上に凸の放物線を描き、$m = \mu_* = 1.0$ でちょうど赤い天井（対数尤度）に接します。オレンジで塗った領域がギャップ＝事後分布との KL で、$m$ が真の事後平均からずれるほど厚くなります。中央：標準偏差 $s$ を動かすと、$s = 0.577$（真の事後の標準偏差）で最大になり、そこから小さくしても大きくしても ELBO は下がります。$s \to 0$ で急落するのは KL の $-\log s^2$ 項が発散するためで、「デルタ関数に潰れるのは禁止」という圧力が効いている証拠です。右：その綱引きの内訳です。緑の再構成項は $s$ が小さいほど良く（$z$ が定まるほど $x$ を当てやすい）、オレンジの $-$正則化項は $s$ が 1 に近いほど良い。両者の和である青い ELBO のピークが、その中間の 0.577 に落ち着いています。ELBO の最適解は「当てやすさ」と「事前分布への従順さ」の妥協点である、という VAE の核心が1枚の図に出ています。

ここまでは $q$ をガウスに固定して2つのパラメータを動かしただけでした。次は、複数の潜在変数がある実際のモデルで、ELBO を体系的に最大化するアルゴリズムを見ていきます。

平均場変分推論とELBOの単調増加

潜在変数が複数ある場合、$q(\bm z)$ を自由に動かすとまた計算不能になります。そこで実用的な制約を入れます。最も広く使われるのが平均場近似（mean-field approximation）です。

イメージとしては、「相互作用している粒子たちを、それぞれが平均的な場の中を独立に動いていると見なす」という統計力学の近似そのものです。相関を切り捨てる代わりに、扱いやすさを得ます。

$$ \begin{equation} q(\bm{z}) = \prod_{j=1}^{J} q_j(z_j) \end{equation} $$

この制約のもとで ELBO を最大化しましょう。ここで効くのが、前の節で導いた「エネルギー + エントロピー」形です。ある一つの因子 $q_j$ だけに注目し、他は固定します。

$$ \mathcal{L} = \int \prod_i q_i(z_i) \left[\log p(x,\bm z) – \sum_i \log q_i(z_i)\right] d\bm z $$

$q_j$ に関係する部分だけを取り出すため、$z_j$ 以外について先に積分してしまいます。第1項は $z_j$ を残して他を積分すると、$q_{-j}$ による条件付き期待値になります。

$$ \mathcal{L} = \int q_j(z_j)\, \mathbb{E}_{q_{-j}}\big[\log p(x, \bm z)\big]\, dz_j – \int q_j(z_j)\log q_j(z_j)\, dz_j + \text{const} $$

ここで「$\text{const}$」は $q_j$ を含まない項（他の因子のエントロピーなど）をまとめたものです。次に、期待値の部分を未規格化の対数密度と見なします。

$$ \log \tilde{p}(z_j) \equiv \mathbb{E}_{q_{-j}}\big[\log p(x, \bm z)\big] + \text{const} $$

と定義して定数を規格化定数に吸わせると、上式は

$$ \mathcal{L} = \int q_j \log \frac{\tilde{p}(z_j)}{q_j(z_j)}\, dz_j + \text{const} = -D_{\mathrm{KL}}\big(q_j \,\|\, \tilde{p}\big) + \text{const} $$

という形になります。ELBO を $q_j$ について最大化することは、$-\mathrm{KL}$ を最大化すること、つまり KL をゼロにすることです。よって最適解は一発で書けます。

$$ \begin{equation} \log q_j^{*}(z_j) = \mathbb{E}_{q_{-j}}\big[\log p(x, \bm{z})\big] + \text{const} \end{equation} $$

これが CAVI（Coordinate Ascent Variational Inference、座標上昇変分推論） の更新式です。「同時分布の対数を、自分以外の因子で期待値を取れ」という、覚えやすくて強力な処方箋になっています。

そして重要な性質があります。各更新は「他を固定して ELBO を厳密に最大化」しているので、ELBO は反復ごとに必ず増加（少なくとも減少しない）します。しかも ELBO には上限 $\log p(x)$ があるので、単調増加＋上に有界により収束が保証されます。実装のバグを見つける最良の方法が「ELBO をプロットして単調増加か確かめる」なのは、この定理があるからです。

Python実装2: CAVIでELBOが単調増加することを確認する

Blei らの変分推論レビューで紹介されている、1次元ガウス混合モデルの CAVI を実装します。モデルは次の通りです。

• 各成分の平均: $\mu_k \sim \mathcal{N}(0, \sigma^2)$、$k = 1,\dots,K$（$\sigma^2$ は大きく取った事前分散）
• 各データ点の所属: $c_i \sim \mathrm{Categorical}(1/K,\dots,1/K)$
• 観測: $x_i \mid c_i, \bm{\mu} \sim \mathcal{N}(\mu_{c_i}, 1)$（観測分散は既知で1）

平均場近似 $q(\bm{\mu}, \bm{c}) = \prod_k q(\mu_k; m_k, s_k^2)\prod_i q(c_i; \bm{\varphi}_i)$ を置いて CAVI 更新式を作ると、次の2つになります（どちらも先ほどの $\log q_j^* = \mathbb{E}_{-j}[\log p]$ から導けます）。

$$ \varphi_{ik} \propto \exp\!\left\{ m_k x_i – \frac{m_k^2 + s_k^2}{2} \right\}, \qquad m_k = \frac{\sum_i \varphi_{ik}x_i}{1/\sigma^2 + \sum_i \varphi_{ik}}, \quad s_k^2 = \frac{1}{1/\sigma^2 + \sum_i \varphi_{ik}} $$

$\varphi_{ik}$ の式に $\mathbb{E}[\mu_k^2] = m_k^2 + s_k^2$ が現れているのが平均場らしいところです。点推定ではなく「不確かさ込み」で相手を評価するため、まだ位置が定まっていない（$s_k^2$ が大きい）成分には所属が振られにくくなります。

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
K, N, sigma2 = 3, 300, 25.0                  # 成分数 / データ数 / 事前分散
mu_true = np.array([-4.0, 0.0, 4.0])
z_true = rng.integers(0, K, N)
X = mu_true[z_true] + rng.standard_normal(N) # 観測分散1

# --- 初期化 ---
m, s2 = rng.standard_normal(K), np.ones(K)
phi = rng.dirichlet(np.ones(K), size=N)

def compute_elbo(X, phi, m, s2):
    """ELBO = E[log p(mu)] + E[log p(c)] + E[log p(x|c,mu)] + H[q(mu)] + H[q(c)]"""
    N = len(X)
    t1 = np.sum(-0.5*np.log(2*np.pi*sigma2) - (m**2 + s2)/(2*sigma2))   # E[log p(mu_k)]
    t2 = -N*np.log(K)                                                   # E[log p(c_i)]
    Exz = -0.5*np.log(2*np.pi) - (X[:, None]**2 - 2*X[:, None]*m[None, :]
                                  + (m**2 + s2)[None, :])/2
    t3 = np.sum(phi * Exz)                                              # E[log p(x|c,mu)]
    t4 = np.sum(0.5*np.log(2*np.pi*np.e*s2))                            # H[q(mu)]
    t5 = -np.sum(phi * np.log(phi + 1e-300))                            # H[q(c)]
    return t1 + t2 + t3 + t4 + t5

hist = [compute_elbo(X, phi, m, s2)]
for it in range(50):
    # (1) 所属 phi の更新（対数領域で計算して数値安定化）
    lp = X[:, None]*m[None, :] - (m**2 + s2)[None, :]/2
    lp -= lp.max(axis=1, keepdims=True)
    phi = np.exp(lp); phi /= phi.sum(axis=1, keepdims=True)
    # (2) 成分平均 q(mu_k) の更新
    nk = phi.sum(axis=0)
    m  = (phi * X[:, None]).sum(axis=0) / (1/sigma2 + nk)
    s2 = 1.0 / (1/sigma2 + nk)
    hist.append(compute_elbo(X, phi, m, s2))

hist = np.array(hist)
print("ELBO 推移(最初の6点):", np.round(hist[:6], 4))
print("最終 ELBO:", round(hist[-1], 4))
print("単調増加か:", np.all(np.diff(hist) >= -1e-8), " 最小の増分:", np.diff(hist).min())
o = np.argsort(m)
print("q(mu_k) の平均:", np.array2string(m[o], precision=4, suppress_small=True),
      " 標準偏差:", np.round(np.sqrt(s2[o]), 4))
print("実効的なクラスタサイズ:", np.round(nk[o], 2))

実行結果は次の通りです。

ELBO 推移(最初の6点): [-2399.5108 -1497.5981  -830.7013  -742.624   -739.7597  -739.6958]
最終 ELBO: -739.692
単調増加か: True  最小の増分: -4.547473508864641e-13
q(mu_k) の平均: [-4.2608 -0.0021  4.0412]  標準偏差: [0.1055 0.1029 0.0929]
実効的なクラスタサイズ: [ 89.78  94.41 115.81]

3つの確認ができました。第1に、ELBO は $-2399.5$ から $-739.7$ まで一度も下がらずに上昇しています。「最小の増分」が $-4.5\times10^{-13}$ という浮動小数点の丸め誤差レベルの負値にとどまっており、理論通りの単調増加です。第2に、上昇のほとんどが最初の3〜4反復で起き（初回だけで $+901.9$）、10反復目には改善が $10^{-6}$ を切ります。CAVI が非常に速く収束することがわかります。第3に、推定された成分平均 $(-4.26,\ -0.00,\ 4.04)$ は真値 $(-4, 0, 4)$ をよく捉え、実効的なクラスタサイズ $(89.8,\ 94.4,\ 115.8)$ も真の所属数 $(90,\ 95,\ 115)$ とほぼ一致しています。

推移を図にすると、単調増加と急速な収束がより鮮明になります。

fig, axes = plt.subplots(1, 2, figsize=(12, 4.2))

axes[0].plot(hist, "o-", color="royalblue", ms=4)
axes[0].set_xlabel("CAVI の反復回数"); axes[0].set_ylabel("ELBO")
axes[0].set_title("ELBOは反復ごとに単調増加する"); axes[0].grid(alpha=0.3)

axes[1].semilogy(np.maximum(np.diff(hist), 1e-16), "o-", color="crimson", ms=4)
axes[1].set_xlabel("CAVI の反復回数"); axes[1].set_ylabel("ELBO の増分（対数軸）")
axes[1].set_title("増分は指数的に小さくなる（収束）"); axes[1].grid(alpha=0.3, which="both")

plt.tight_layout(); plt.show()

左のグラフは、ELBO が階段を一段も踏み外さずに登り切る様子です。中央の対数軸のグラフを見ると、増分が反復ごとにおよそ1桁ずつ落ちていき、10反復あたりで機械精度に張り付きます。「ELBO をプロットして単調増加を確かめる」というのが変分推論のデバッグの定石である理由がここにあります。もしどこかで ELBO が下がったら、それは更新式の導出ミスか実装ミスのサインです。右のグラフは推定結果で、3つの $q(\mu_k)$ の山が真の成分平均（赤い破線）のすぐ近くに、しかも標準偏差 0.09〜0.11 という細い幅で立っています。点推定ではなく「平均の不確かさ」まで一緒に得られるのが、変分ベイズが最尤法の EM と違うところです。

CAVI が使えるのは、$\mathbb{E}_{q_{-j}}[\log p(x,\bm z)]$ が解析的に計算できる（＝共役性がある）モデルに限られます。ニューラルネットワークをデコーダに使う VAE では当然そんな共役性はありません。そこで必要になるのが、勾配法で ELBO を最大化する技術です。

再パラメータ化トリック — 期待値の中を微分する

ELBO を勾配法で最大化したいとき、意外な障害にぶつかります。期待値を取る分布そのものが、最適化したいパラメータに依存しているという点です。

再構成項を見てください。

$$ \mathcal{L} = \mathbb{E}_{q_\phi(z\mid x)}\big[\log p_\theta(x\mid z)\big] – D_{\mathrm{KL}}(q_\phi \| p) $$

デコーダのパラメータ $\theta$ については何の問題もありません。$q_\phi$ は $\theta$ に依存しないので、

$$ \nabla_\theta \mathbb{E}_{q_\phi}\big[\log p_\theta(x\mid z)\big] = \mathbb{E}_{q_\phi}\big[\nabla_\theta \log p_\theta(x\mid z)\big] $$

と微分をそのまま中に入れられ、サンプル平均で置き換えれば済みます。

問題はエンコーダのパラメータ $\phi$ です。$\phi$ を変えるとサンプルを引いてくる分布そのものが変わるため、

$$ \nabla_\phi \mathbb{E}_{q_\phi(z)}\big[f(z)\big] \neq \mathbb{E}_{q_\phi(z)}\big[\nabla_\phi f(z)\big] $$

です（そもそも右辺の $f(z)$ は $\phi$ を含まないので $0$ になってしまいます）。サイコロの目の期待値を上げたいとき、目の値を微分しても意味がなく、サイコロの出目の確率を変えないといけない、という話です。

方法A: スコア関数推定量（REINFORCE）。素直に微分を積分の中に入れると、

$$ \nabla_\phi \int q_\phi(z) f(z)\, dz = \int f(z)\, \nabla_\phi q_\phi(z)\, dz $$

となります。ここで対数微分の恒等式 $\nabla_\phi q_\phi = q_\phi \nabla_\phi \log q_\phi$（$\nabla \log u = \nabla u / u$ を移項しただけ）を使うと、ふたたび $q_\phi$ による期待値の形に戻せます。

$$ \begin{equation} \nabla_\phi \mathbb{E}_{q_\phi}[f(z)] = \mathbb{E}_{q_\phi}\big[f(z)\, \nabla_\phi \log q_\phi(z)\big] \end{equation} $$

これは任意の $f$ で使える万能な推定量ですが、致命的な弱点があります。$f$ が掛け算のスカラーとしてしか入ってこないので、「$z$ をどちらの方向に動かせば $f$ が増えるか」という情報を一切使っていません。結果として分散が大きく、実用には工夫（ベースライン、制御変量）が必須です。

方法B: 再パラメータ化トリック。発想を変えます。「$\phi$ に依存する分布からサンプリングする」のをやめて、「$\phi$ に依存しない分布からノイズを引き、それを $\phi$ に依存する決定的な関数で変換する」ことにします。

$$ z = g_\phi(\varepsilon, x), \qquad \varepsilon \sim p(\varepsilon)\quad(\text{$\phi$ に無関係}) $$

ガウス分布の場合、これは中学生でも知っている変換です。$\mathcal{N}(\mu, \sigma^2)$ からのサンプルは、標準正規からのサンプル $\varepsilon$ を使って

$$ \begin{equation} z = \mu_\phi(x) + \sigma_\phi(x) \odot \varepsilon, \qquad \varepsilon \sim \mathcal{N}(\bm{0}, \bm{I}) \end{equation} $$

と書けます。$\odot$ は要素ごとの積です。この置き換えで、期待値を取る分布は $\phi$ を含まない $p(\varepsilon)$ になり、$\phi$ は被期待値関数の中に移動します。

$$ \mathbb{E}_{q_\phi(z\mid x)}\big[f(z)\big] = \mathbb{E}_{p(\varepsilon)}\big[f(g_\phi(\varepsilon, x))\big] $$

こうなればもう障害はありません。期待値の外側の分布が $\phi$ と無関係なので、微分をそのまま中に入れられます。

$$ \begin{equation} \nabla_\phi \mathbb{E}_{q_\phi}[f(z)] = \mathbb{E}_{p(\varepsilon)}\big[\nabla_z f(z)\big|_{z=g_\phi(\varepsilon)} \cdot \nabla_\phi g_\phi(\varepsilon, x)\big] \end{equation} $$

決定的なのは $\nabla_z f(z)$ が現れている点です。「$z$ を動かすと $f$ がどう変わるか」という勾配情報が、そのまま $\phi$ へ流れ込みます。スコア関数推定量が $f$ の値しか見なかったのに対し、こちらは $f$ の傾きを使う。これがパスワイズ勾配（pathwise derivative）と呼ばれる所以で、分散が劇的に小さくなります。実装上は「$\mu$ と $\sigma$ を出力し、そこにノイズを足す」だけなので、自動微分がそのまま通ります。

どれくらい違うのか、数値で見てみましょう。$q = \mathcal{N}(m, 0.6^2)$、$f(z) = -\frac{(x-z)^2}{2\sigma_n^2}$（$x=1.5$、$\sigma_n^2=0.5$）として、$m = 0.5$ における $\nabla_m \mathbb{E}_q[f]$ を、8サンプルで2万回推定します。真値は $\frac{x-m}{\sigma_n^2} = \frac{1.0}{0.5} = 2.0$ です。

import numpy as np
rng = np.random.default_rng(1)

sn2, x, s, m = 0.5, 1.5, 0.6, 0.5
f  = lambda z: -(x - z)**2/(2*sn2)     # 再構成項に相当
df = lambda z: (x - z)/sn2             # その z 微分
print("真の勾配:", (x - m)/sn2)

S, trials = 8, 20000                   # 1回の推定に使うサンプル数 / 試行回数
g_sf, g_rp = np.zeros(trials), np.zeros(trials)
for t in range(trials):
    eps = rng.standard_normal(S)
    z = m + s*eps                                   # 再パラメータ化されたサンプル
    g_sf[t] = np.mean(f(z) * (z - m)/s**2)          # スコア関数推定量
    g_rp[t] = np.mean(df(z))                        # 再パラメータ化推定量

print(f"スコア関数    : 平均={g_sf.mean():.4f}  標準偏差={g_sf.std():.4f}")
print(f"再パラメータ化: 平均={g_rp.mean():.4f}  標準偏差={g_rp.std():.4f}")
print(f"分散比(スコア/再パラメータ化): {(g_sf.std()/g_rp.std())**2:.1f} 倍")

真の勾配: 2.0
スコア関数    : 平均=2.0072  標準偏差=1.6798
再パラメータ化: 平均=2.0035  標準偏差=0.4220
分散比(スコア/再パラメータ化): 15.8 倍

左のヒストグラムでは、どちらの山も赤い線（真の勾配 2.0）を中心にしていますが、青い再パラメータ化推定量が鋭く尖っているのに対し、オレンジのスコア関数推定量は右に長く裾を引いて 6〜8 まで外れる試行が珍しくありません。右の両対数プロットは、サンプル数 $S$ を増やしたときの標準偏差で、2本の直線がどちらも傾き $-1/2$（$1/\sqrt{S}$ 則）で平行に下がっています。分散比が $S$ によらずおよそ15倍で一定ということは、サンプルをいくら増やしてもスコア関数推定量が追いつくことはなく、推定量の選択そのものが効いていることを意味します。

どちらも平均は 2.0 付近で、両方とも不偏推定量であることが確認できます。違うのはばらつきです。スコア関数推定量の標準偏差 1.68 に対し、再パラメータ化は 0.42。分散にして約 16 倍の差があります。同じ精度の勾配を得るのに 16 分の1のサンプルで済むわけで、ミニバッチあたり1サンプルしか引かない VAE の学習が成立するのは、この分散削減のおかげです。逆に言えば、離散潜在変数のように再パラメータ化が使えない場合は、Gumbel-Softmax のような連続緩和や、ベースライン付きの REINFORCE に頼ることになります。

これで、ELBO を勾配法で最大化する道具がそろいました。最後に、実際に小さな VAE を学習させて、ELBO の2つの項がどう動くかを見届けましょう。

Python実装3: 小さなVAEでELBOの2項の推移を見る

scikit-learn に同梱されている 8×8 の手書き数字データ（1797枚、ダウンロード不要）に、潜在次元 2 の小さな VAE を学習させます。目的は精度を出すことではなく、再構成項と KL 項がそれぞれどんな軌跡を描くかを観察することです。

尤度は $p_\theta(x\mid z) = \prod_j \mathrm{Bernoulli}(x_j; \hat{x}_j(z))$ とします（画素値を $[0,1]$ に正規化してベルヌーイのパラメータと見なす、VAE で標準的な扱いです）。このとき $-\log p_\theta(x\mid z)$ はバイナリクロスエントロピーの和になります。KL 項は先に導いた閉じた式をそのまま使います。

import numpy as np, torch, torch.nn as nn
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits

torch.manual_seed(0); np.random.seed(0)
X = torch.tensor(load_digits().data / 16.0, dtype=torch.float32)  # (1797, 64) を [0,1] へ
D, H, L = 64, 128, 2                                              # 入力 / 隠れ / 潜在次元

class VAE(nn.Module):
    def __init__(self):
        super().__init__()
        self.enc = nn.Sequential(nn.Linear(D, H), nn.ReLU())
        self.mu, self.logvar = nn.Linear(H, L), nn.Linear(H, L)
        self.dec = nn.Sequential(nn.Linear(L, H), nn.ReLU(), nn.Linear(H, D))

    def forward(self, x):
        h = self.enc(x)
        mu, logvar = self.mu(h), self.logvar(h)
        std = torch.exp(0.5*logvar)
        z = mu + std*torch.randn_like(std)      # ← 再パラメータ化トリック
        return self.dec(z), mu, logvar

model = VAE()
opt = torch.optim.Adam(model.parameters(), lr=1e-3)
bce = nn.BCEWithLogitsLoss(reduction="none")

エンコーダは分散そのものではなく $\log \sigma^2$（logvar）を出力しています。これで $\sigma^2 > 0$ が自動的に保証され、exp(0.5*logvar) で標準偏差に戻せます。z = mu + std*randn_like(std) の1行が、前節で導いた $z = \mu + \sigma \odot \varepsilon$ そのものです。

hist = []
for epoch in range(200):
    perm = torch.randperm(len(X))
    r_sum = k_sum = 0.0
    for i in range(0, len(X), 128):
        xb = X[perm[i:i+128]]
        logits, mu, logvar = model(xb)
        rec = bce(logits, xb).sum(1).mean()                            # -E_q[log p(x|z)]
        kl  = (-0.5*(1 + logvar - mu**2 - logvar.exp()).sum(1)).mean() # KL(q||p(z))
        loss = rec + kl                                                # = -ELBO
        opt.zero_grad(); loss.backward(); opt.step()
        r_sum += rec.item()*len(xb); k_sum += kl.item()*len(xb)
    hist.append((r_sum/len(X), k_sum/len(X)))

hist = np.array(hist)
for e in [0, 40, 80, 120, 160, 199]:
    r, k = hist[e]
    print(f"epoch {e:3d}  再構成項 {r:7.3f}  KL項 {k:6.3f}  ELBO {-(r+k):8.3f}")

epoch   0  再構成項  42.237  KL項  0.034  ELBO  -42.271
epoch  40  再構成項  23.662  KL項  1.781  ELBO  -25.443
epoch  80  再構成項  23.088  KL項  2.033  ELBO  -25.121
epoch 120  再構成項  22.656  KL項  2.192  ELBO  -24.848
epoch 160  再構成項  22.436  KL項  2.276  ELBO  -24.712
epoch 199  再構成項  22.279  KL項  2.316  ELBO  -24.595

ここで注目すべきは、2つの項が正反対の方向に動いていることです。再構成項（負の対数尤度なので小さいほど良い）は 42.24 から 22.28 へ半減しています。一方 KL 項は 0.034 から 2.32 へ、70倍近く増えています。KL 項だけ見れば「悪化」ですが、その分を上回る再構成の改善があるので、和である ELBO は $-42.27$ から $-24.60$ へ着実に上昇しています。

epoch 0 で KL 項がほぼゼロだったことには意味があります。学習の初期、エンコーダはまだ何も学んでおらず、どの入力に対してもほぼ $\mathcal{N}(0,I)$ を返しています。つまり $q(z\mid x)$ が $x$ の情報を一切運んでいない、事後分布の崩壊（posterior collapse）そのものの状態です。KL 項は「潜在変数が何ビットの情報を運んでいるか」の上限を測る指標でもあるので、KL が育つということは、潜在コードが実際に使われ始めたことを意味します。

推移を図で確認しましょう。

fig, axes = plt.subplots(1, 3, figsize=(15, 4.2))

axes[0].plot(-(hist[:,0] + hist[:,1]), color="royalblue", lw=2)
axes[0].set_xlabel("エポック"); axes[0].set_ylabel("ELBO")
axes[0].set_title("ELBO は学習とともに上昇する"); axes[0].grid(alpha=0.3)

axes[1].plot(hist[:,0], color="seagreen", lw=2, label="再構成項 $-E_q[\\log p(x|z)]$")
axes[1].set_xlabel("エポック"); axes[1].set_ylabel("再構成項（小さいほど良い）")
axes[1].set_title("再構成項は 42.2 → 22.3 へ半減"); axes[1].grid(alpha=0.3); axes[1].legend()

axes[2].plot(hist[:,1], color="darkorange", lw=2, label="KL項 $D_{KL}(q\\|p(z))$")
axes[2].set_xlabel("エポック"); axes[2].set_ylabel("KL項（潜在コードの情報量の目安）")
axes[2].set_title("KL項は 0.03 → 2.3 へ増加（潜在変数が使われ始める）")
axes[2].grid(alpha=0.3); axes[2].legend()

plt.tight_layout(); plt.show()

左の ELBO は最初の 20 エポックで急上昇し、その後ゆるやかに改善を続けます。中央の再構成項は同じタイミングで急落します。右の KL 項は逆に、ゼロ付近からじわじわと立ち上がっていきます。3枚を並べると「ELBO の上昇は、再構成の大幅な改善を KL の小幅な悪化で買っている」という取引の構造がはっきり見えます。潜在変数に情報を詰め込むほど KL のコストを払うことになるが、そのコストを上回るリターン（再構成の改善）がある限り、最適化はその方向に進む。ELBO は、この「情報を運ぶ対価」の会計をしている目的関数なのです。

学習の結果できあがった潜在空間と、再構成の実物も見ておきましょう。

上の散布図は、各データをエンコーダに通して得た $\mu_\phi(x)$ を2次元にプロットし、参考までに数字のラベルで色分けしたものです（ラベルは学習に一切使っていません）。ELBO を最大化しただけで、0 や 6 や 4 がそれぞれ固まり、しかも全体が原点まわりの狭い範囲（両成分とも $-2.8 \sim 2.6$）に収まっている——この「散らばりすぎない」性質こそ KL 項の仕事で、だからこそ $p(z)=\mathcal{N}(0,I)$ から $z$ を引いて新しい画像を作れます。下の比較は再構成項が測っているものの実物で、潜在次元をたった2まで絞ったため大まかな字形は戻るものの、線の細部はぼやけ、いくつかの数字は別の数字に見えるほど崩れています。このぼやけの量がそのまま再構成項の値（22.3）であり、KL 項を払ってでも潜在次元を増やせば減る、という取引の対象になっています。

ELBOの限界 — 下界であることの代償

ELBO は強力ですが、万能ではありません。実務で足をすくわれないよう、限界を3つ整理しておきます。

限界1: ELBO は対数尤度そのものではない。中心等式が示す通り、ELBO は常に $\log p_\theta(x)$ より $D_{\mathrm{KL}}(q\|p(z\mid x))$ だけ低い値です。しかもこのギャップは、モデルやデータによって大きく変わります。したがって「モデル A の ELBO のほうがモデル B より高いから、A のほうが尤度が高い」とは必ずしも言えません。A のギャップが小さくて B のギャップが大きいだけかもしれないからです。モデル比較で ELBO を使うときは、この危うさを意識する必要があります。より締まった評価が欲しければ、先に触れた多サンプル下界（IWAE bound）や、アニールド重点サンプリングを使います。

限界2: 逆向き KL のモード探索性。変分推論が最小化しているのは $D_{\mathrm{KL}}(q\|p)$ であって $D_{\mathrm{KL}}(p\|q)$ ではありません。この向きの KL は $q$ が正のところで $p$ が小さいと強く罰する（$q\log(q/p)$ が発散する）ため、$q$ は「$p$ が確実に大きい領域」に縮こまる性質を持ちます。真の事後分布が多峰性を持つとき、$q$ はそのうちひとつの峰にへばりつき、分散を過小評価する傾向があります。ベイズ的な不確かさを定量化したいのに、変分推論の信用区間が狭すぎる、という現象はここから来ます。MCMC が（十分回せば）分布全体を捉えるのと対照的です。詳しくは 平均場近似による変分推論をわかりやすく解説 で扱っています。

限界3: 近似族の表現力とアモチゼーションギャップ。ギャップは実は2段階に分けられます。ひとつは「近似族 $\mathcal{Q}$ の中で最良の $q^*$ を選んでも残るギャップ」（近似族ギャップ）。対角ガウスしか許さないのに真の事後が強く相関していれば、どうやっても埋まりません。もうひとつは VAE 特有で、「$x$ ごとに最適化する代わりに、エンコーダというひとつの関数に予測させる」ことで生じるアモチゼーションギャップです。エンコーダの容量が足りないと、$q^*$ に届かない $q$ しか出力できません。正規化流（Normalizing Flow）で $q$ の表現力を上げる、あるいはテスト時に $q$ のパラメータを追加で最適化する（semi-amortized VI）といった対策が研究されています。

これらの限界は ELBO を捨てる理由ではなく、「下界を使っている」という事実を忘れないための注意書きです。下界であることを自覚したうえで使えば、ELBO は依然として、計算不能だった問題を実際に解ける最適化問題へ翻訳してくれる、最も便利な道具であり続けます。

まとめ

本記事では、ELBO（変分下界）を2通りの導出から徹底的に解剖しました。

• 問題の出発点 — 潜在変数モデルの対数尤度 $\log p_\theta(x) = \log \int p_\theta(x,z)dz$ は高次元積分のせいで計算できず、真の事後分布 $p_\theta(z\mid x)$ もその尤度を分母に含むため計算できない、という循環がある
• 導出その1（イェンセン） — $\log$ の凹性から $\log \mathbb{E}_q[p/q] \ge \mathbb{E}_q[\log(p/q)]$ となり、3行で下界 $\mathcal{L}$ が得られる。ただしギャップの大きさはわからない
• 導出その2（完全な等式） — $\log p_\theta(x) = \mathcal{L}(q,\theta) + D_{\mathrm{KL}}(q(z\mid x)\|p_\theta(z\mid x))$。ギャップの正体は真の事後分布との KL であり、$q$ について ELBO を最大化することと事後分布との KL を最小化することは完全に同値
• 実装に落とす形 — $\mathcal{L} = \mathbb{E}_q[\log p_\theta(x\mid z)] – D_{\mathrm{KL}}(q(z\mid x)\|p(z))$。再構成項と正則化項の綱引きが VAE の損失関数そのもの
• ガウスのKL — $D_{\mathrm{KL}}(\mathcal{N}(\mu,\sigma^2)\|\mathcal{N}(0,1)) = \frac12(\mu^2+\sigma^2-1-\log\sigma^2)$ を1次元から導出し、対角多次元へ拡張した
• 再パラメータ化トリック — $z = \mu_\phi + \sigma_\phi \odot \varepsilon$ と書き換えることで、期待値の外の分布から $\phi$ を追い出し、勾配を中に入れられる。数値実験ではスコア関数推定量に対して分散が約16分の1になった
• 数値検証 — 1次元共役ガウスで「対数尤度 $-$ ELBO」と「$q$ と真の事後の KL」が小数点以下4桁まで一致すること、CAVI で ELBO が単調増加（最小増分 $-4.5\times10^{-13}$）すること、小さな VAE で再構成項 42.2→22.3・KL 項 0.03→2.32 と項が逆方向に動きながら ELBO が上昇することを確認した

ELBO を理解すると、一見バラバラに見えた手法群がひとつの視点でつながります。EM アルゴリズムは ELBO の交互最大化、変分ベイズは平均場制約つきの ELBO 最大化、VAE はニューラルネットで償却した ELBO 最大化、拡散モデルは多段階の変分下界。「計算できない対数尤度を、計算できる下界に置き換える」という一手が、現代の生成モデルのほぼすべてを支えています。

次のステップとして、以下の記事も参考にしてください。

• 【深層学習】VAEとは？VAEの理論を徹底的に分かりやすく解説する — ELBO を実際の生成モデルに実装する
• VAE（変分オートエンコーダ）のELBO導出と潜在空間の理解 — 潜在空間の構造と生成
• 平均場近似による変分推論をわかりやすく解説 — CAVI と平均場近似の詳細
• 【変分ベイズ】変分推論やELBOを理解する — 混合ガウスの変分ベイズ
• 【機械学習】EMアルゴリズムをゼロから理解する — ELBO の交互最大化としての EM
• モンテカルロ法と重点サンプリング — 高次元期待値計算の理論と実装 — 重点サンプリングと多サンプル下界の背景

参考文献

• D. P. Kingma, M. Welling, “Auto-Encoding Variational Bayes,” ICLR 2014.
• D. M. Blei, A. Kucukelbir, J. D. McAuliffe, “Variational Inference: A Review for Statisticians,” JASA 112(518), 2017.
• Y. Burda, R. Grosse, R. Salakhutdinov, “Importance Weighted Autoencoders,” ICLR 2016.
• C. M. Bishop, Pattern Recognition and Machine Learning, Springer, 2006（第10章 近似推論）.
• C. Cremer, X. Li, D. Duvenaud, “Inference Suboptimality in Variational Autoencoders,” ICML 2018.