简介

综合转载以下文章：

前言

在概率模型的应⽤中，⼀个中⼼任务是在给定观测（可见）数据变量 $X$ 的条件下，计算潜在变量 $Z$ 的后验概率分布 $p(Z|X)$，以及计算关于这个概率分布的期望。模型可能也包含某些确定性参数，我们现在不考虑它。模型也可能是⼀个纯粹的贝叶斯模型，其中任何未知的参数都有⼀个先验概率分布，并且被整合到了潜在变量集合中，记作向量 $Z$。

例如，在EM算法中，我们需要计算完整数据对数似然函数关于潜在变量后验概率分布的期望。对于实际应⽤中的许多模型来说，计算后验概率分布或者计算关于这个后验概率分布的期望是不可⾏的。这可能是由于潜在空间的维度太⾼，以⾄于⽆法直接计算，或者由于后验概率分布的形式特别复杂，从⽽期望⽆法解析地计算。在连续变量的情形中，需要求解的积分可能没有解析解，⽽空间的维度和被积函数的复杂度可能使得数值积分变得不可⾏。对于离散变量，求边缘概率的过程涉及到对隐含变量的所有可能的配置进⾏求和。这个过程虽然原则上总是可以计算的，但是我们在实际应⽤中经常发现，隐含状态的数量可能有指数多个，从⽽精确的计算所需的代价过⾼。

举个例子：高斯函数的贝叶斯混合，

有 $\mu_k\sim \cal N(0,\tau^2)$，其中 $k=1,\dots,K$；
对于 $i=1,\dots,n$：
1. 有 $z_i\sim \text{Mult}(\pi)$；
2. 有 $x_i\sim \cal N(\mu_{z_i},\sigma^2)$。

固定其他参数，后验分布为：

$p(\mu_{1:K},z_{1:n}|x_{1:n})=\frac{\prod_{k=1}^Kp(\mu_k)\prod_{i=1}^np(z_i)p(x_i|z_i,\mu_{1:K})}{\int_{\mu_{1:K}}\sum_{z_{1:n}}\prod_{k=1}^Kp(\mu_k)\prod_{i=1}^np(z_i)p(x_i|z_i,\mu_{1:K})}$

对于隐变量的任何结构，分子都很容易计算，问题是分母。

我们试着计算一下：首先，给定簇中心，我们可以利用 $z_i$ 的条件独立性，

$p(x_{1:n})=\int_{\mu_{1:K}}\prod_{k=1}^Kp(\mu_k)\prod_{i=1}^n\sum_{z_i}p(z_i)p(x_i|z_i,\mu_{1:K})$

这出现了一个积分，计算困难。

或者，我们可以把隐变量的和移到外面，

$p(x_{1:n})=\sum_{z_i}p(z_i)\int_{\mu_{1:K}}\prod_{k=1}^Kp(\mu_k)\prod_{i=1}^np(x_i|z_i,\mu_{1:K})$

结果是我们可以在这个和式中计算每一项。但是，也有 $K^n$ 项。当 $n$ 相当大时，这是很难计算的。这种情况出现在大多数有趣的模型中。这就是为什么近似后验推理是贝叶斯统计的核心问题之一。

在这种情况下，我们需要借助近似⽅法。根据近似⽅法依赖于随机近似还是确定近似，⽅法⼤体分为两⼤类。随机⽅法，例如马尔可夫链蒙特卡罗⽅法，使得贝叶斯⽅法能够在许多领域中⼴泛使⽤。这些⽅法通常具有这样的性质：给定⽆限多的计算资源，它们可以⽣成精确的结果，近似的来源是使⽤了有限的处理时间。在实际应⽤中，取样⽅法需要的计算量会相当⼤，经常将这些⽅法的应⽤限制在了⼩规模的问题中。并且，判断⼀种取样⽅法是否
⽣成了服从所需的概率分布的独⽴样本是很困难的。

变分推断（variational inference）或者变分贝叶斯（variational Bayes），它使⽤了更加全局的准则，并且被⼴泛应⽤于实际问题中。上世纪90年代，变分推断在概率模型上得到迅速发展，在贝叶斯框架下一般的变分法由Attias的两篇文章给出。Matthew J.Beal的博士论文《Variational Algorithms for Approximate Bayesian Inference》中有比较充分地论述，作者将其应用于隐马尔科夫模型，混合因子分析，线性动力学，图模型等。它主要应用于复杂的统计模型中，这种模型一般包括三类变量：观测变量（observed variables，data），未知参数（parameters）和隐变量（latent variables）。在贝叶斯推断中，参数和潜变量统称为不可观测变量（unobserved variables）。变分贝叶斯方法主要是两个目的：

近似不可观测变量的后验概率，以便通过这些变量作出统计推断。
对一个特定的模型，给出观测变量的边缘似然函数（或称为证据，evidence）的下界。主要用于模型的选择，认为模型的边缘似然值越高，则模型对数据拟合程度越好，该模型产生Data的概率也越高。

对于第一个目的，随机方法，例如蒙特卡罗模拟，特别是用Gibbs取样的MCMC方法，可以近似计算复杂的后验分布，能很好地应用到贝叶斯统计推断。此方法通过大量的样本估计真实的后验，因而近似结果带有一定的随机性。与此不同的是，变分贝叶斯方法提供一种局部最优，但具有确定解的近似后验方法。

从某种角度看，变分贝叶斯可以看做是EM算法的扩展，因为它也是采用极大后验估计（MAP），即用单个最有可能的参数值来代替完全贝叶斯估计。另外，变分贝叶斯也通过一组相互依赖（mutually dependent）的等式进行不断的迭代来获得最优解。

变分贝叶斯和EM的区别：

EM算法的E-step计算的 $q(Z)$ 是精确的，即 $q(Z)=\frac{p(X,Z)}{\sum_Zp(X,Z)}$，不是近似的。只是因为每一步的E-step受限于M-step，M-step又受限于E-step，所以才需要迭代优化，最终得到的是近似后验。而变分贝叶斯的动机就是用简单的 $q(Z)$ 来近似 $p(Z|X)$。这是两者的区别，所以才有花书中所说的：EM并不是一个近似推断算法，而是一种能够学到近似后验的算法。
EM算法中的E-step和M-step分别是针对隐变量和未知参数来做迭代的，而变分贝叶斯中，未知参数被划分到隐变量中，之后通过划分隐变量，来对每一次优化做迭代。

变分推断

问题描述

重新考虑一个问题：有一组观测数据 $\pmb X$，并且已知模型的形式，求参数与潜变量（或不可观测变量）$\pmb Z = \{ {Z_1},…,{Z_n}\}$ 的后验分布：$p(\pmb Z|\pmb X)$。

正如上文所描述的后验概率的形式通常是很复杂（Intractable）的，对于一种算法如果不能在多项式时间内求解，往往不是我们所考虑的。因而我们想能不能在误差允许的范围内，用更简单、容易理解（tractable）的数学形式 $q(\pmb Z)$ 来近似 $p(\pmb Z \vert\pmb X)$，即 $p(\pmb Z \vert\pmb X) \approx q(\pmb Z)$。

由此引出如下两个问题：

假设存在这样的 $q(\pmb Z)$，那么如何度量 $q(\pmb Z)$ 与 $p(\pmb Z|\pmb X)$ 之间的差异性（dissimilarity）。
如何得到简单的 $q(\pmb Z)$？

对于问题一，幸运的是，我们不需要重新定义一个度量指标。在信息论中，已经存在描述两个随机分布之间距离的度量，即相对熵，或者称为Kullback-Leibler散度，即下一节。

对于问题二，显然我们可以自主决定 $q(\pmb Z)$ 的分布，只要它足够简单，且与 $p(\pmb Z \vert\pmb X)$ 接近。然而不可能每次都手工给出一个与 $p(\pmb Z \vert\pmb X)$ 接近且简单的 $q(\pmb Z)$，其方法本身已经不具备可操作性。所以需要一种通用的形式帮助简化问题。那么数学形式复杂的原因是什么？奥卡姆剃刀，认为一个模型的参数个数越多，那么模型复杂的概率越大；此外，如果参数之间具有相互依赖关系（mutually dependent），那么通常很难对参数的边缘概率精确求解。

幸运的是，统计物理学界很早就关注了高维概率函数与它的简单形式，并发展了平均场理论。简单讲就是：系统中个体的局部相互作用可以产生宏观层面较为稳定的行为。于是我们可以作出后验条件独立（posterior independence）的假设。即，$\forall i,p(\pmb Z \vert\pmb X) = p(\pmb Z_i \vert\pmb X)p({\pmb Z_{-i}} \vert\pmb X)$。（$\pmb Z_{-i}$ 表示互斥集）

相对熵（KL散度）

熵的本质是香农信息量（$\text{log}\frac{1}{p}$）的期望。

现有关于样本集的2个概率分布 $p$ 和 $q$ ，其中 $p$ 为真实分布，$q$ 非真实分布。按照真实分布 $p$ 来衡量识别一个样本的所需要的编码长度的期望（即平均编码长度）为：

$H(p)=\sum_{i}^{} p(i)\ast\text{log}\frac{1}{p(i)}$

如果使用错误分布 $q$ 来表示来自真实分布 $p$ 的平均编码长度，则应该是：

$H(p,q)=\sum_{i}^{} p(i)*log\frac{1}{q(i)}$

因为用 $q$ 来编码的样本来自分布 $p$，所以期望 $H(p,q)$ 中概率是 $p(i)$。$H(p,q)$ 我们称之为“交叉熵”。

比如含有4个字母 $(A,B,C,D)$ 的数据集中，真实分布 $p=(1/2, 1/2, 0, 0)$，即 $A$ 和 $B$ 出现的概率均为 $1/2$，$C$ 和 $D$ 出现的概率都为0。计算 $H(p)$ 为 $1$ ，即只需要1位编码即可识别 $A$ 和 $B$ 。如果使用分布 $q=(1/4, 1/4, 1/4, 1/4)$ 来编码则得到 $H(p,q)=2$ ，即需要2位编码来识别 $A$ 和 $B$ （当然还有 $C$ 和 $D$ ，尽管$C$ 和 $D$ 并不会出现，因为真实分布 $p$ 中 $C$ 和 $D$ 出现的概率为0，这里就钦定概率为0的事件不会发生了）。

可以看到上例中根据非真实分布 $q$ 得到的平均编码长度 $H(p,q)$ 大于根据真实分布 $p$ 得到的平均编码长度 $H(p)$。事实上，根据Gibbs’ inequality可知，$H(p,q)\ge H(p)$ 恒成立，当 $q$ 为真实分布 $p$ 时取等号。我们将由 $q$ 得到的平均编码长度比由 $p$ 得到的平均编码长度多出的bit数称为“相对熵”：

$D_{KL}(p||q)=H(p,q)-H(p)=\sum_{i}^{} p(i)\ast \text{log}\frac{p(i)}{q(i)}$

其又被称为KL散度（Kullback–Leibler divergence，KLD）。它表示2个函数或概率分布的差异性：差异越大则相对熵越大，差异越小则相对熵越小，特别地，若2者相同则熵为0。

KL散度有如下性质：

${D_{KL}}(p \vert \vert q) \ne {D_{KL}}(q \vert \vert p)$；
${D_{KL}}(p \vert \vert q) \ge 0$，当且仅当 $p=q$ 时为零；
不满足三角不等式。

$q$ 分布与 $p$ 分布的KL散度为：

$D_{KL}(q \vert \vert p) = \sum\limits_{\pmb Z} q(\pmb Z)\ln \frac{q(\pmb Z)}{p(\pmb Z \vert\pmb X)} = \sum\limits_{\pmb Z} q(\pmb Z)\ln \frac{q(\pmb Z)}{p(\pmb Z,\pmb X)} + \ln p(\pmb X)$

则：

$\ln p(\pmb X) = {D_{KL}}(q \vert \vert p) - \sum\limits_{\pmb Z} {q(\pmb Z)\ln \frac{q(\pmb Z)}{p(\pmb Z,\pmb X)}} ={D_{KL}}(q \vert \vert p) + \cal L(q).$

上述是在离散变量的情况，所以使用了求和。在连续变量的情况中为，

$\text{ln}p(\pmb X)=\mathcal L(q)+D_{KL}(q||p)$

其中，

$\mathcal L(q)=\int q(\pmb Z)\text{ln}\left\{\frac{p(\pmb X,\pmb Z)}{q(\pmb Z)}\right\}d\pmb Z\\ D_{KL}(q||p)=-\int q(\pmb Z)\text{ln}\left\{\frac{p(\pmb Z|\pmb X)}{q(\pmb Z)}\right\}d\pmb Z$

这与我们关于EM的讨论的唯⼀的区别是参数向量 $θ$ 不再出现，因为参数现在是随机变量，被整合到了 $\pmb Z$ 中。

与之前⼀样，由于对数证据 $\ln p(\pmb X)$ 被相应的 $q$ 所固定，我们可以通过关于概率分布 $q(\pmb Z)$ 的最优化来使下界 $\cal L(q)$ 达到最⼤值，这等价于最⼩化KL散度，这样就可以得到后验 $p(\pmb Z|\pmb X)$ 的近似解析表达式和证据（log evidence）的下界 $\cal L(q)$ ，又称为变分自由能（variational free energy）：

$\mathcal L(q)=\sum\limits_{\pmb Z} {q(\pmb Z)\ln p(\pmb Z,\pmb X)}-\sum\limits_{\pmb Z} {q(\pmb Z)\ln q(\pmb Z)}={E_q}[\ln p(\pmb Z,\pmb X)]+H(q)$

上式是离散变量情况，连续变量为：

$\mathcal L(q)=\int {q(\pmb Z)\ln p(\pmb Z,\pmb X)}d\pmb Z-\int {q(\pmb Z)\ln q(\pmb Z)}d\pmb Z={E_q}[\ln p(\pmb Z,\pmb X)]+H(q)$

如果我们允许任意选择 $q(\pmb Z)$，那么下界的最⼤值出现在KL散度等于零的时刻，此时 $q(\pmb Z)$ 等于后验概率分布 $p(\pmb Z |\pmb X)$。然⽽，我们假定在需要处理的模型中，对真实的概率分布进⾏操作是不可⾏的。

于是，我们转⽽考虑概率分布 $q(\pmb Z)$ 的⼀个受限制的类别，然后寻找这个类别中使得KL散度达到最⼩值的概率分布。我们的⽬标是充分限制 $q(\pmb Z)$ 可以取得的概率分布的类别范围，使得这个范围中的所有概率分布都是可以处理的概率分布。同时，我们还要使得这个范围充分⼤、充分灵活，从⽽它能够提供对真实后验概率分布的⼀个⾜够好的近似。需要强调的是，施加限制条件的唯⼀⽬的是为了计算⽅便，并且在这个限制条件下，我们应该使⽤尽可能丰富的近似概率分布。特别地，对于⾼度灵活的概率分布来说，没有“过拟合”现象。使⽤灵活的近似仅仅使得我们更好地近似真实的后验概率分布。

限制近似概率分布的范围的⼀种⽅法是使⽤参数概率分布 $q(\pmb Z | \pmb \omega )$，它由参数集合 $\pmb \omega$ 控制。这样，下界 $\cal L(q)$ 变成了 $\pmb \omega$ 的函数，我们可以利⽤标准的⾮线性最优化⽅法确定参数的最优值。下图给出了这种⽅法的⼀个例⼦，其中变分分布是⼀个⾼斯分布，并且我们已经关于均值和协⽅差进⾏了最优化。

平均场理论（Mean Field Method）

这⾥，我们考虑另⼀种⽅法，这种⽅法⾥，我们限制概率分布 $q(\pmb Z)$ 的范围，具体地说，用平均场理论来分解 $\pmb Z$。

数学上说，平均场的适用范围只能是完全图，或者说系统结构是well-mixed，在这种情况下，系统中的任何一个个体以等可能接触其他个体。反观物理，平均场与其说是一种方法，不如说是一种思想。其实统计物理的研究目的就是期望对宏观的热力学现象给予合理的微观理论。物理学家坚信，即便不满足完全图的假设，但既然这种“局部”到“整体”的作用得以实现，那么个体之间的局部作用相较于“全局”的作用是可以忽略不计的。

根据平均场理论，变分分布 $q(\pmb Z)$ 可以通过参数和潜在变量的划分（partition）因式分解，比如将 $\pmb Z$ 划分为 ${\pmb Z_1} \ldots {\pmb Z_M}$，

$q(\pmb Z) = \prod\limits_{i = 1}^M {q_i(\pmb Z_i)}$

需要强调的是，我们关于概率分布没有做更多的假设。特别地，我们没有限制各个因⼦ $q_i(Z_i)$ 的函数形式。

注意：这里并非一个不可观测变量一个划分，而应该根据实际情况做决定。当然你也可以这么做，但是有时候，将几个潜变量放在一起会更容易处理。

泛函极值

上文已经提到我们要找到一个更加简单的函数 $q(\pmb Z)$ 来近似 $p(\pmb Z \vert\pmb X)$，同时问题转化为求解证据 $\text{log} p(D)$ 的下界 $\cal L(q)$，或者 $\mathcal L(q(\pmb Z))$。应该注意到 $\cal L(q)$ 并非普通的函数，而是以整个函数为自变量的函数，这便是泛函。我们先介绍一下什么是泛函，以及泛函取得极值的必要条件。

泛函

设对于（某一函数集合内的）任意一个函数 $y(x)$，有另一个数 $J[y]$ 与之对应，则称 $J[y]$ 为 $y(x)$ 的泛函。泛函可以看成是函数概念的推广。这里的函数集合，即泛函的定义域，通常要求 $y(x)$ 满足一定的边界条件，并且具有连续的二阶导数。这样的 $y(x)$ 称为可取函数。

泛函不同于复合函数

例如 $g=g(f(x))$；对于后者，给定一个 $x$ 值，仍然是有一个 $g$ 值与之对应；对于前者，则必须给出某一区间上的函数 $y(x)$，才能得到一个泛函值 $J[y]$。（定义在同一区间上的）函数不同，泛函值当然不同，为了强调泛函值 $J[y]$ 与函数 $y(x)$ 之间的依赖关系，常常又把函数 $y(x)$ 称为变量函数。

泛函的形式多种多样，大部分以积分形式表示：$J[y] = \int_{x_1}^{x_2} {L(x,y,y’)} dx$。

泛函取极值的必要条件

泛函的极值

“当变量函数为 $y(x)$ 时，泛函 $J[y]$ 取极大值”的含义就是：对于极值函数 $y(x)$ 及其“附近”的变量函数 $y(x)+δy(x)$，恒有 $J[y+δy]≤J[y]$；

所谓函数 $y(x)+δy(x)$ 在另一个函数 $y(x)$ 的“附近”，指的是：

$|δy(x)|<ε$；
有时还要求 $|(δy)′(x)|<ε$。

这里的 $δy(x)$ 称为函数 $y(x)$ 的变分。

Euler–Lagrange方程

可以仿造函数极值必要条件的导出办法，导出泛函取极值的必要条件，这里不做严格的证明，直接给出。泛函 $J[y]$ 取到极大值的必要条件是一级变分 $δJ[y]$ 为0，其微分形式一般为二阶常微分方程，即Euler-Largange方程：

$\frac{\partial F}\partial y - \frac{d}dx\frac{\partial F}\partial y' = 0$

泛函的条件极值

在约束条件下求函数 $J[y]$ 的极值，可以引入Largange乘子 $λ$，从而定义一个新的泛函， $\tilde J[y] = J[y] - \lambda {J_0}[y]$。仍将 $\delta y$ 看成是独立的，则泛函 $\tilde J[y]$ 在边界条件下取极值的必要条件就是，

$(\frac{\partial}{\partial y} - \frac{d}{dx} \frac{\partial}{\partial y'}) (F - \lambda G) = 0$

需要注意的是，Euler–Lagrange方程只是函数算子取得极值的必要条件，而不是充分条件，但是我们总算是对于变分法有一定的初步认识，也知道其中一个求函数算子极值的方法。

所以我们也清楚了为什么“变分推断”和“变分贝叶斯”里含有“变分”两个字，因为变分法是用来求泛函的极值，而我们需要最优化 $q(\pmb Z)$ 来最大化 $\cal L(q)$，$\cal L(q)$ 为泛函。

问题求解

由以上，我们希望对 $\cal L(q)$ 关于所有的概率分布 $q_i(\pmb Z_i)$ 进⾏⼀个⾃由形式的（变分）最优化。我们通过关于每个因⼦进⾏最优化来完成整体的最优化过程。之前有：

$\mathcal L(q)=\int {q(\pmb Z)\ln p(\pmb Z,\pmb X)}d\pmb Z-\int {q(\pmb Z)\ln q(\pmb Z)}d\pmb Z$

将平均场理论引入，为了记号的简洁，我们简单地将 $q_j(\pmb Z_j)$ 记作 $q_j$，得：

$\begin{align*} \mathcal L(q)&=\int\prod_iq_i\left\{\ln p(\textbf X,\textbf Z)-\sum_i\ln q_i\right\}d\textbf Z\\ &=\int q_j\left\{\int\ln p(\textbf X,\textbf Z)\prod_{i\neq j}q_id\textbf{Z}_i\right\}d\textbf Z_j-\int q_j\ln q_j d\textbf Z_j+\text{const}\\ &=\int q_j\ln \tilde p(\textbf X,\textbf Z_j)d\textbf Z_j-\int q_j\ln q_jd\textbf Z_j+\text{const} \end{align*}$

其中，$\ln \tilde p(\textbf X,\textbf Z_j)=\mathbb E_{i\neq j}[\ln p(\textbf X,\textbf Z)]+\text{const}$。这里，记号 $\mathbb E_{i:i\neq j}[\dots]$ 表示关于定义在所有 $\pmb Z_i(i\neq j)$ 上的 $q$ 概率分布的期望，即

$\mathbb E_{i\neq j}[\ln p(\pmb X, \pmb Z)]=\int \prod_{i:i\neq j}q_i\ln p(\pmb X,\pmb Z)d\pmb Z_i$

现在假设我们保持 $\{q_{i:i\neq j}\}$ 固定，关于概率分布 $q_j(\pmb Z_j)$ 的所有可能的形式最⼤化 $\cal L(q)$。于是，我们得到了最优解 $q_j^\ast(\pmb Z_j)$ 的一般的表达式，形式为：

$\ln q_j^\ast(\pmb Z_j)=\mathbb E_{i\neq j}[\ln p(\pmb X,\pmb Z)]+常数 \tag 1$

这个解表明，为了得到 $q_j$ 的最优解的对数，我们只需考虑所有隐变量和可见变量上的联合概率分布的对数，然后序所有的其他的因子 $\{q_{i:i\neq j}\}$ 取期望。

上式的可加性常数通过对概率分布 $q_j^\ast(\pmb Z_j)$ 进行归一化的方式来设定。因此，如果我们取两侧的指数，然后归⼀化，我们有

$q_j^\ast(\pmb Z_j)=\frac{\exp(\mathbb E_{i\neq j}[\ln p(\pmb X,\pmb Z)])}{\int\exp(\mathbb E_{i\neq j}[\ln p(\pmb X,\pmb Z)])d\pmb Z_j}$

在实际应⽤中，我们会发现，更⽅便的做法是对公式(1)进⾏操作，然后在必要的时候，通过观察的⽅式恢复出归⼀化系数。这⼀点通过下⾯的例⼦就会变得逐渐清晰起来。

上述的最优解也可以直接用变分法求得，由条件极值的必要条件：

$\frac{\partial}{\partial q_j(\pmb Z_j)}\left\{\int q_j(\pmb Z_j)\ln q_j^\ast (\pmb Z_j)d\pmb Z_j-\int q_j(\pmb Z_j)\ln q_j(\pmb Z_j)d\pmb Z_j+\lambda_i\left(\int_i q_i(\pmb Z_i) d\pmb Z_i-1\right)\right\}$

经过一系列推导，最后也可以得到最优解。

由公式(1 给定的⽅程的集合（其中 $j=1,\dots,M$）表示在概率能够进⾏分解这⼀限制条件下，下界的最⼤值满⾜的⼀组相容的条件。然⽽，这些⽅程并没有给出⼀个显式的解，因为最优化 $q_j^\ast(\pmb Z_j)$ 的公式(1)的右侧表达式依赖于关于其他的因⼦ $q_{i:i\neq j}(\pmb Z_j)$ 计算的期望。于是，我们会⽤下⾯的⽅式寻找出⼀个相容的解：

恰当地初始化所有的因⼦ $q_i(\pmb Z_i)$；
在各个因⼦上进⾏循环，每⼀轮⽤⼀个修正后的估计来替换当前因⼦。这个修正后的估计由公式(1)的右侧给出，计算时使⽤了当前对于所有其他因⼦的估计。

算法保证收敛，因为下界关于每个因⼦ $q_i(\pmb Z_i)$ 是⼀个凸函数。

虽然从理论上推导了变分推断的框架算法，但是对于不同模型，我们必须手动推导 $q_j^\ast(\pmb Z_j)$，简要来说，推导变分贝叶斯模型一般分为四个步骤：

如果想做全贝叶斯模型，确定好研究模型各个参数的的共轭先验分布；
写出研究模型的联合分布 $p(\pmb Z,\pmb X)$；
根据联合分布确定变分分布的形式 $q(\pmb Z)$；
对于每个变分因子 $q_j(\pmb Z_j)$ 求出 $p(\pmb Z,\pmb X)$ 关于不包含变量 $\pmb Z_j$ 的数学期望，再规整化为概率分布。

分解近似的性质

我们的变分推断的⽅法基于的是真实后验概率分布的分解近似。让我们现在考虑⼀下使⽤分解概率分布的⽅式近似⼀个⼀般的概率分布的问题。⾸先，我们讨论使⽤分解的⾼斯分布近似⼀个⾼斯分布的问题，这会让我们认识到在使⽤分解近似时会引⼊的不准确性有哪些类型。考虑两个相关的变量 $\pmb z=(z_1, z_2)$ 上的高斯分布 $p(\pmb z)=\mathcal N(\pmb z|\pmb \mu, \pmb \Lambda^{-1})$，其中均值和精度的元素为：

$\pmb \mu=\left(\begin{aligned}\mu_1\\\mu_2\end{aligned}\right), \pmb\Lambda=\left(\begin{aligned}&\Lambda_{11}&\Lambda_{12}\\&\Lambda_{21}&\Lambda_{22}\end{aligned}\right)$

并且由于精度矩阵的对称性， $Λ_{21} = Λ_{12}$。现在，假设我们希望使⽤⼀个分解的⾼斯分布 $q(\pmb z) = q_1(z_1)q_2(z_2)$ 来近似这个分布。⾸先，我们使⽤公式(1)来寻找最优因⼦ $q_1^\ast(z_1)$ 的表达式。在寻找表达式的过程中，我们注意到，在等式右侧，我们只需要保留哪些与 $z_1$ 有函数依赖关系的项即可，因为所有其他的项都可以被整合到归⼀化常数中。因此我们有

$\begin{align*} \ln q_1^\ast(z_1)=&\mathbb E_{z_2}[\ln p(\pmb z)]+常数\\ =&\mathbb E_{z_2}\left[-\frac{1}{2}(z_1-\mu_1)^2\Lambda_{11}-(z_1-\mu_1)\Lambda_{12}(z_2-\mu_2)\right]+常数\\ =&-\frac{1}{2}z_1^2\Lambda_{11}+z_1\mu_1\Lambda_{11}-z_1\Lambda_{12}(\mathbb E[z_2]-\mu_2)+常数 \end{align*}$

接下来，我们观察到这个表达式的右侧是 $z_1$ 的⼀个⼆次函数，因此我们可以将 $q^\ast(z_1)$ 看成⼀个⾼斯分布。值得强调的是，我们不假设 $q(z_i)$ 是⾼斯分布，⽽是通过对所有可能的分布 $q(z_i)$ 上的KL散度的变分最优化推导出了这个结果。还要注意，我们不需要显式地考虑公式(1)中的可加性常数，因为它表⽰归⼀化常数。如果需要的话，这个常数可以在计算的最后阶段通过观察的⽅式得到。使⽤配平⽅的⽅法，我们可以得到这个⾼斯分布的均值和⽅差，有

$q_1^\ast(z_1)=\mathcal N(z_1|m_1,\Lambda_{11}^{-1})$

其中，

$m_1=\mu_1-\Lambda_{11}^{-1}\Lambda_{12}(\mathbb E[z_2]-\mu_2)$

根据对成型，$q_2^\ast(z_2)$ 也是一个高斯分布性质，可以写成

$q_2^\ast(z_2)=\mathcal N(z_2|m_2,\Lambda_{22}^{-1})$

其中，

$m_2=\mu_2-\Lambda_{22}^{-1}\Lambda_{21}(\mathbb E[z_1]-\mu_1)$

这些解是相互依赖的。通常，是按照上节所说的循环更新法。但是这里，可以找到一个解析解。特别地，由于 $\mathbb E[z_1]=m_1$ 且 $\mathbb E[z_2]=m_2$，因此只要取 $\mathbb E[z_1]=\mu_1$ 且 $\mathbb E[z_2]=\mu_2$，那么这两个方差会得到满足。并且很容易证明，只要概率分布⾮奇异，那么这个解是唯⼀解。这个结果如下图所⽰。我们看到，均值被正确地描述了，但
是 $q(\pmb z)$ 的⽅差由 $p(\pmb z)$ 的最⼩⽅差的⽅向所确定，沿着垂直⽅向的⽅差被强烈地低估了。这是⼀个⼀般的结果，即分解变分近似对后验概率分布的近似倾向于过于紧凑。

作为⽐较，假设我们最⼩化相反的Kullback-Leibler散度 $D_{KL}(p||q)$。正如我们将看到的那样，这种形式的KL散度被⽤于另⼀种近似推断的框架中，这种框架被称为期望传播（expectation propagation）。于是，我们考虑⼀般的最⼩化 $D_{KL}(p||q)$ 的问题，其中 $q(\pmb Z)$ 还是同样的分解近似。这样， KL散度可以写成

$D_{KL}(p||q)=-\int p(\pmb Z)\left[\sum_{i=1}^M\ln q_i(\pmb Z_i)\right]d\pmb Z+常数$

其中，常数项就是 $p(\pmb Z)$ 的熵，因此不依赖于 $q(\pmb Z)$。我们现在可以关于每个因⼦ $q_j(\pmb Z_j)$ 进⾏最优化。使⽤拉格朗⽇乘数法，很容易得到结果

$q_j^\ast(\pmb Z_j)=\int p(\pmb Z)\prod_{i:i\neq j}d\pmb Z_j=p(\pmb Z_j)$

在这种情况下，我们看到 $q_j(\pmb Z_j)$ 的最优解等于对应的边缘概率分布 $p(\pmb Z)$。注意，这是⼀个解析解，不需要迭代。将其用于这个例子上，如(b)，我们再⼀次看到，对均值的近似是正确的，但是它把相当多的概率权重放到了实际上具有很低的概率的变量空间区域中。

这两个结果的区别可以⽤下⾯的⽅式理解。注意到， $\pmb Z$ 空间中 $p(\pmb Z)$ 接近等于零的区域对于Kullback-Leibler散度，即下式

$D_{KL}(q||p)=-\int p(\pmb Z)\ln\left\{\frac{p(\pmb Z)}{q(\pmb Z)}\right\}d\pmb Z$

有⼀个⼤的正数的贡献，除⾮ $q(\pmb Z)$ 也接近等于零。因此最⼩化这种形式的KL散度会使得概率分布 $q(\pmb Z)$ 避开 $p(\pmb Z)$ 很⼩的区域。相反地，使得Kullback-Leibler散度 $D_{KL}(p||q)$ 的散度取得最⼩值的概率分布 $q(\pmb Z)$ 在 $p(\pmb Z)$ 很小的区域也会有较大的概率权重。

如果我们考虑⽤⼀个单峰分布近似多峰分布的问题，我们会更深刻地认识两个KL散度的不同⾏为，如下图所⽰。在实际应⽤中，真实的后验概率分布经常是多峰的，⼤部分后验概率质量集中在参数空间中的某⼏个相对较⼩的区域中。这些多个峰值可能是由于潜在空间的不可区分性所造成的，也可能是由于对参数的复杂的⾮线性依赖关系造成的。基于最⼩化 $D_{KL}(q||p)$ 的变分⽅法倾向于找到这些峰值中的⼀个。相反，如果我们最⼩化 $D_{KL}(p||q)$，那么得到的近似会在所有的均值上取平均。在混合模型问题中，这种⽅法会给出较差的预测分布（因为两个较好的参数值的平均值通常不是⼀个较好的参数值）。可以使⽤ $D_{KL}(p||q)$ 定义⼀个有⽤的推断步骤，在期望传播中讨论。

两种形式的Kullback-Leibler散度都是散度的alpha家族（alpha family）的成员，定义为

$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)$

其中，$-\infty<\alpha<\infty$ 是一个连续参数。Kullback-Leibler散度 $D_{KL}(p||q)$ 对应于极限 $α\rightarrow 1$，⽽ $D_{KL}(q||p)$ 对应于极限 $α\rightarrow -1$。对于所有 $α$ 的值，我们有 $D_α(p||q) ≥ 0$，当且仅当 $p(x) = q(x)$ 时等号成⽴。假设 $p(x)$ 是⼀个固定的分布，我们关于某个概率分布 $q(x)$ 的集合最⼩化 $D_α(p||q)$。那么对于 $\alpha ≤ -1$ 的情况，散度是零强制的（zero forcing），即对于使得 $p(x) = 0$ 成⽴的任意 $x$ 值，都有 $q(x) = 0$，通常 $q(x)$ 的权重比 $p(x)$ 的权重大，所以会低估 $p(x)$，因此倾向于寻找具有最⼤权重的峰值。相反，对于 $α ≥ 1$ 的情况，散度是零避免的（zero avoiding），即对于使得 $p(x) > 0$ 成⽴的任意 $x$ 值，都有 $q(x) > 0$ ，通常 $q(x)$ 会进⾏拉伸来覆盖到所有的 $p(x)$ 值，从⽽⾼估了p(x)。当 $α = 0$ 时，我们得到了⼀个对称的散度，它与Hellinger距离线性相关，定义
为

$D_H(p||q)=\int(p(x)^{\frac{1}{2}}-q(x)^{\frac{1}{2}})^2dx$

Hellinger距离的平⽅根是⼀个合法的距离度量。

一元高斯分布

我们现在使⽤⼀元变量 $x$ 上的⾼斯分布来说明分解变分近似。我们的⽬标是在给定 $x$ 的观测值的数据集 $\mathcal D =\{x_1,\dots,x_N\} $ 的情况下，推断均值 $\mu$ 和精度 $τ$ 的后验概率分布。其中，我们假设数据是独⽴地从⾼斯分布中抽取的。似然函数为

$p(\mathcal D|\mu,\tau^{-1})=\left(\frac{\tau}{2\pi}\right)^{\frac{N}{2}}\text{exp}\left\{-\frac{\tau}{2}\sum_{n=1}^N(x_n-\mu)^2\right\}$

现在引入 $\mu$ 和 $\tau$ 的共轭先验分布，形式为：

$p(\mu|\tau)=\mathcal N(\mu|\mu_0,(\lambda_0\tau)^{-1})\\ p(\tau)=\text{Gam}(\tau|a_0,b_0)$

其中 $\text{Gam}(\tau|a_0,b_0)$ 是Gamma分布，这些分布共同给出了一个高斯-Gamma共轭先验分布。

Gamma分布

$f(x;\alpha,\beta)=\frac{\beta^\alpha x^{\alpha-1}\exp(-\beta x)}{\Gamma(\alpha)},\text{for } x>0 \text{ and }\alpha,\beta>0.$

考虑对后验概率分布的一个分解变分近似，形式为

$q(\mu,\tau)=q_\mu(\mu)q_\tau(\tau)$

注意，真实的后验概率分布不可以按照这种形式进⾏分解。最优的因⼦ $q_\mu(\mu)$ 和 $q_\tau(\tau)$ 可以从公式(1)得出。对于 $q_\mu(\mu)$，有：

$\begin{align*} \ln q_\mu^\ast(\mu)&=\mathbb E_\tau[\ln p(\mathcal D|\mu,\tau)+\ln p(\mu|\tau)+\ln p(\tau)]+常数\\ &=\mathbb E_\tau[\ln p(\mathcal D|\mu,\tau)+\ln p(\mu|\tau)]+\mathbb{E}_\tau[\ln p(\tau)]+常数\\ &=\mathbb E_\tau[\ln p(\mathcal D|\mu,\tau)+\ln p(\mu|\tau)]+常数\\ &=-\frac{\mathbb E[\tau]}{2}\left\{\lambda_0(\mu-\mu_0)^2+\sum_{n=1}^N(x_n-\mu)^2\right\}+常数 \end{align*}$

对于 $\mu$ 配平方，我们看到 $q_\mu(\mu)$ 是一个高斯分布 $\mathcal N(\mu|\mu_N,\lambda_N^{-1})$，其中，均值和方差为：

$\mu_N=\frac{\lambda_0\mu_0+N\bar x}{\lambda_0+N}\tag 2$ $\lambda_N=(\lambda_0+N)\mathbb E[\tau]\tag 3$

注意，对于 $N\rightarrow \infty$，这给出了最大似然的结果，其中 $\mu_N=\bar x$，精度为无穷大。

类似地，因子 $q_\tau(\tau)$ 的最优解为：

$\begin{align*} \ln q_\tau^\ast(\tau)&=\mathbb E_\mu[\ln p(\mathcal D|\mu, \tau)+\ln p(\mu|\tau)]+\ln p(\tau)+常数\\ &=\mathbb E_\mu[\ln p(\mathcal D|\mu, \tau)+\ln p(\mu|\tau)]+(a_0-1)\ln \tau-b_0\tau+常数\\ &=-\frac{\tau}{2}\mathbb {E}_\mu\left[\sum_{n=1}^N(x_n-\mu)^2+\lambda_0(\mu-\mu_0)^2\right]+\frac{N+1}{2}\ln \tau+(a_0-1)\ln \tau-b_0\tau+常数 \end{align*}$

因此 $q_\tau(\tau)$ 是一个Gamma分布 $\text{Gam}(\tau|a_N,b_N)$，参数为：

$a_N=a_0+\frac{N+1}{2}\tag 4$ $b_N=b_0+\frac{1}{2}\mathbb E_\mu\left[\sum_{n=1}^N(x_n-\mu)^2+\lambda_0(\mu-\mu_0)^2\right]\tag 5$

与之前⼀样，当 $N\rightarrow \infty$ 时，它的⾏为与预期相符。

应该强调的是，我们不假设最优概率分布 $q_\mu(\mu)$ 和 $q_τ(τ)$ 的具体的函数形式。它们的函数形式从似然函数和对应的共轭先验分布中⾃然地得到。

因此，我们得到了最优概率分布 $q_\mu(\mu)$ 和 $q_τ(τ)$ 的表达式，每个表达式依赖于关于其他概率分布计算得到的矩。因此，⼀种寻找解的⽅法是对例如 $E[τ]$ 进⾏⼀个初始的猜测，然后使⽤这个猜测来重新计算概率分布 $q_\mu(\mu)$。给定这个修正的概率分布之后，我们接下来可以计算所需的矩 $E[\mu]$ 和 $E[\mu^2]$，并且使⽤这些矩来重新计算概率分布 $q_τ(τ)$，以此类推。由于这个例⼦中，隐含变量空间是⼆维的，因此我们可以⽤图形来说明后验概率分布的变分近似过程。我们画出了真实后验概率的轮廓线和分解近似的轮廓线，如下图所示。

通常，我们需要使⽤⼀种迭代的⽅法来得到最优分解后验概率分布的解。然⽽，对于我们这⾥讨论的⾮常简单的例⼦来说，我们可以通过求解最优因⼦ $q_\mu(\mu)$ 和 $q_τ(τ)$ 的⽅程，得到⼀个显式的解。在做这件事之前，我们可以通过考虑⽆信息先验来简化表达式。⽆信息先验分布中，$\mu_0=a_0=b_0=λ_0=0$。虽然这些参数设置对应于⼀个反常先验，但是我们看到后验概率分布仍然具有良好的定义。使⽤Gamma分布的均值的标准结果 $\mathbb E[\tau]=\frac{a_N}{b_N}$，以及公式(4)(5)，有：

$\frac{1}{\mathbb E[\tau]}=\mathbb E\left[\frac{1}{N+1}\sum_{n=1}^N(x_n-\mu)^2\right]=\frac{N}{N+1}(\bar {x^2}-2\bar x\mathbb E[\mu]+\mathbb E[\mu^2])$

之后，使用公式(2)(3)，得到了 $q_\mu(\mu)$ 的一阶矩和二阶矩，形式为：

$\mathbb E[\mu]=\bar x,\ \mathbb E[\mu^2]=\bar x^2+\frac{1}{N\mathbb E[\tau]}$

现在，将这些矩代入公式，然后解出 $\mathbb E[\tau]$，可得

$\frac{1}{\mathbb E[\tau]}=(\bar {x^2}-\bar x^2)=\frac{1}{N}\sum_{n=1}^N(x_n-\bar x)^2$

模型比较

除了在隐含变量Z上进⾏推断之外，我们可能还希望对⽐⼀组候选模型。索引为 $m$ 的模型的先验概率分布为 $p(m)$。这样，我们的⽬标是近似后验概率分布 $p(m|\pmb X)$，其中 $\pmb X$ 是观测数据。这⽐我们⽬前为⽌考虑的情况稍微复杂⼀些，因为不同的模型可能具有不同的结构，并且隐含变量 $\pmb Z$ 的维度实际上可能不同。因此我们不能简单地考虑考虑分解近似 $q(\pmb Z)q(m)$，⽽是必须意识到 $\pmb Z$ 的后验概率分布必须以 $m$ 为条件，所以我们必须考虑 $q(\pmb Z,m) =q(\pmb Z|m)q(m)$。我们已经可以验证下⾯的基于变分概率分布的分解⽅式

$\ln p(\pmb X)=\mathcal L-\sum_m\sum_{\pmb Z}q(\pmb Z|m)q(m)\ln\left\{\frac{p(\pmb Z,m|\pmb X)}{q(\pmb Z|m)q(m)}\right\}$

其中，$\cal L$ 是 $\ln p(\pmb X)$ 的下界，形式为：

$\mathcal L=\sum_m\sum_{\pmb Z}q(\pmb Z|m)q(m)\ln\left\{\frac{p(\pmb Z,\pmb X,m)}{q(\pmb Z|m)}\right\}\tag 6$

这⾥，我们假定 $\pmb Z$ 是离散变量，但是同样的分析也适⽤于连续潜在变量，只要我们把求和替换为积分即可。我们可以使⽤拉格朗⽇乘数法关于概率分布 $q(m)$ 最⼤化 $\cal L$，结果为？？？

$q(m)\propto p(m)\exp\{\cal L_m\}$

其中，

$\mathcal L_m=\sum_{\pmb Z}q(\pmb Z|m)\ln\left\{\frac{p(\pmb Z,\pmb X|m)}{q(\pmb Z|m)}\right\}$

然⽽，如果我们关于 $q(\pmb Z|m)$ 最⼤化 $\cal L$，那么我们发现对于不同的 $m$ 值，解是相互偶合的，这与我们预期相符，因为这些概率分布是以 $m$ 为条件的。我们接下来⾸先通过最优化公式(6)，或者等价地，最优化 $\cal L_m$，来独⽴地最优化每个 $q(\pmb Z | m)$，然后使⽤公式来确定 $q(m)$。在对求得的 $q(m)$ 值进⾏归⼀化之后，它的值可以⽤于模型选择或者模型平均。

高斯的变分混合

许多贝叶斯模型，对应于复杂得多的概率分布，可以通过对本节中的分析进⾏简单的扩展进⾏求解。

对于每个观测 $\pmb x_n$，有一个对应的隐变量 $\pmb z_n$，它是一个”1-of-K”的二值向量，元素为 $z_{nk}$，其中 $k=1,\dots,K$。将观测数据集记作 $\pmb X=\{\pmb x_1, \dots, \pmb x_N\}$，将隐变量记作 $\pmb Z=\{\pmb z_1,\dots,\pmb z_N\}$。给定混合系数 $\pmb \pi$，写出 $\pmb Z$ 的条件概率分布，形式为：

$p(\pmb Z|\pmb \pi)=\prod_{n=1}^N\prod_{k=1}^K\pi_k^{z_{nk}}$

类似地，给定隐变量和分量参数，写出观测数据向量的条件概率分布，形式为：

$p(\pmb X|\pmb Z,\pmb \mu, \pmb \Lambda)=\prod_{n=1}^N\prod_{k=1}^K\mathcal N(\pmb x_n|\pmb \mu_k,\pmb \Lambda_k^{-1})^{z_{nk}}$

其中，$\pmb \mu=\{\pmb \mu_k\}$ 且 $\pmb\Lambda=\{\pmb\Lambda_k\}$。注意，我们计算时使用的是精度绝阵而不是协方差矩阵，因为这在一定程度上简化了数学计算的复杂度。

接下来，引入参数 $\pmb \mu,\pmb \Lambda$ 和 $\pmb \pi$ 上的先验概率分布。如果使用共轭先验分布，那么分析过程会得到极大的简化。于是，选择混合系数 $\pmb \pi$ 上的狄力克雷分布。

$p(\pmb \pi)=\text{Dir}(\pmb \pi|\pmb \alpha_0)=C(\pmb\alpha_0)\prod_{k=1}^K\pi_k^{\alpha_0-1}$

其中，根据对称性，我们为每个分量选择了同样的参数 $\alpha_0$，$C(\pmb\alpha_0)$ 是狄力克雷分布的归一化常数。参数 $\alpha_0$ 可以看成与混合分布的每个分量关联的观测的有效先验数量。如果 $\alpha_0$ 的值很小，那么后验概率分布会主要被数据集影响，⽽受到先验概率的影响很⼩。

类似地，我们引⼊⼀个独⽴的⾼斯-Wishart先验分布，控制每个⾼斯分布的均值和精度，形式为

$\begin{align*} p(\pmb \mu,\pmb\Lambda)&=p(\pmb\mu|\pmb\Lambda)p(\pmb\Lambda)\\ &=\prod_{k=1}^K\mathcal N(\pmb \mu_k|\pmb m_0,(\beta_0\pmb\Lambda_k)^{-1})\mathcal W(\pmb\Lambda_k|\pmb W_0,\nu_0) \end{align*}$

这是由于当均值和精度均未知的时候，它表⽰共轭先验分布。通常根据对称性，我们选择 $\pmb m_0 = 0$。⽣成的模型可以表⽰为下图所⽰的有向图。

这个例⼦很好地说明了潜在变量和参数之间的区别。像 $\pmb z_n$ 这样出现在⽅框内部的变量被看做隐变量，因为这种变量的数量随着数据集规模的增⼤⽽增⼤。相反，像 $\pmb \mu$ 这样出现在⽅框外的变量的数量与数据集的规模⽆关，因此被当做参数。然⽽，从图模型的观点来看，它们之间没有本质的区别。

变分分布

为了形式化地描述这个模型的变分⽅法，接下来写出所有随机变量的联合概率分布，形式为：

$p(\pmb X,\pmb Z,\pmb \pi,\pmb \mu,\pmb\Lambda)=p(\pmb X|\pmb Z,\pmb\mu,\pmb\Lambda)p(\pmb Z|\pmb \pi)p(\pmb \pi)p(\pmb\mu|\pmb\Lambda)p(\pmb\Lambda)$

注意，只有变量 $\pmb X=\{\pmb x_1,\dots,\pmb x_N\}$ 是观测变量。现在考虑一个变分分布，它在隐变量与参数之间进行分解，即

$q(\pmb Z,\pmb \pi, \pmb \mu,\pmb \Lambda)=q(\pmb Z)q(\pmb\pi,\pmb\mu,\pmb\Lambda)$

需要注意的是，为了让我们的贝叶斯混合模型能够有⼀个合理的可以计算的解，这是我们需要做出的唯⼀的假设。特别地，因⼦ $q(\pmb Z)$ 和 $q(\pmb π, \pmb \mu, \pmb Λ)$ 的函数形式会在变分分布的最优化过程中⾃动确定。注意，我们省略了 $q$ 分布的下标，依赖参数来区分不同的分布。

通过使⽤⼀般的结果，考虑因⼦ $q(\pmb Z)$ 的更新⽅程的推导。最优因⼦的对数为

$\ln q^\ast(\pmb Z)=\mathbb E_{\pmb \pi,\pmb\mu,\pmb\Lambda}[\ln p(\pmb X,\pmb Z,\pmb \pi, \pmb \mu,\pmb \Lambda)]+常数$

对上式分解，任何与变量 $\pmb Z$ ⽆关的项都可以被整合到可加的归⼀化系数中，从⽽有

$\ln q^\ast(\pmb Z)=\mathbb E_{\pmb \pi}[\ln p(\pmb Z|\pmb \pi)]+\mathbb E_{\pmb \mu,\pmb\Lambda}[\ln p(\pmb X|\pmb Z,\pmb \mu,\pmb \Lambda)]+常数\tag 7$

代入其值，得：

$\ln q^\ast(\pmb Z)=\sum_{n=1}^N\sum_{k=1}^Kz_{nk}\ln \rho_{nk}+常数$

其中定义了

$\ln \rho_{nk}=\mathbb E[\ln \pi_k]+\frac{1}{2}\mathbb E[\ln |\pmb\Lambda_k|]-\frac{D}{2}\ln(2\pi)-\frac{1}{2}\mathbb E_{\pmb \mu_k,\pmb\Lambda_k}[(\pmb x_n-\pmb \mu_k)^T\pmb\Lambda_k(\pmb x_n-\pmb \mu_k)]\tag 8$

其中 $D$ 是数据变量 $x$ 的维度。公式(7)两侧取指数，我们有

$q^\ast(\pmb Z)\propto \prod_{n=1}^N\prod_{k=1}^K\rho_{nk}^{z_{nk}}$

其中，

$r_{nk}=\frac{\rho_{nk}}{\sum_{j=1}^K\rho _{nj}}$

我们看到因子 $q(\pmb Z)$ 最优解的函数形式与先验概率分布 $p(\pmb Z|\pmb \pi)$ 的函数形式相同。注意，由于 $\rho_{nk}$ 是一个实数值的指数，因此 $r_{nk}$ 是非负的，且加和等于1，满足要求。

对于离散概率分布 $q^\ast(\pmb Z)$，我们有标准的结果：

$\mathbb E[z_{nk}]=r_{nk}\tag 9$

从中我们看到 $r_{nk}$ 扮演着“责任”的⾓⾊。注意，$q^\ast(\pmb Z)$ 的最优解依赖于关于其他变量计算得到的矩，因此与之前⼀样，变分更新⽅程是互相依赖的，必须⽤迭代的⽅式求解。

现在，我们会发现定义观测数据关于“责任”的下⾯三个统计量会⽐较⽅便，即

$N_k=\sum_{n=1}^Nr_{nk}\\ \bar x_k=\frac{1}{N_k}\sum_{n=1}^Nr_{nk}\pmb x_n\\ \pmb S_k=\frac{1}{N_k}\sum_{n=1}^Nr_{nk}(\pmb x_n-\bar {\pmb x}_k)(\pmb x_n-\bar{\pmb x}_k)^T$

注意，这些类似于⾼斯混合模型的最⼤似然EM算法中计算的量。

现在让我们考虑变分后验概率分布中的因⼦ $q(\pmb π; \pmb \mu; \pmb Λ)$。与之前⼀样，使⽤公式(1)给出的⼀般的结果，我们有

$\ln q^\ast(\pmb \pi,\pmb \mu,\pmb \pi)=\ln p(\pmb \pi)+\sum_{k=1}^K\ln p(\pmb \mu_k,\pmb \Lambda_k)+\mathbb E_{\pmb Z}[\ln p(\pmb Z|\pmb \pi)]+\sum_{k=1}^K\sum_{n=1}^N\mathbb E[z_{nk}]\ln \mathcal N(\pmb x_n|\pmb \mu_k,\pmb \Lambda_k^{-1})+常数\tag {10}$

我们观察到，这个表达式的右侧分解成了若⼲项的和，⼀些项只与 $\pmb π$ 相关，⼀些项只与 $\pmb \mu$ 和 $\pmb Λ$ 相关，这表明变分后验概率 $q(\pmb π; \pmb \mu;\pmb Λ)$ 可以分解为 $q(\pmb π)q(\pmb \mu; \pmb Λ)$。此外，与 $\pmb \mu$ 和 $\pmb \Lambda$ 相关的项本⾝由 $k$ 个与 $\pmb µ_k$ 和 $\pmb Λ_k$ 相关的项有关，因此可以进⼀步分解，即

$q(\pmb \pi,\pmb \mu,\pmb \Lambda)=q(\pmb \pi)\prod_{k=1}^Kq(\pmb \mu_k,\pmb \Lambda_k)$

分离出右侧的与 $\pmb \pi$ 相关的项，我们有：

$\ln q^\ast(\pmb \pi)=(\alpha_0-1)\sum_{k=1}^K\ln \pi_k+\sum_{k=1}^K\sum_{n=1}^Nr_{nk}\ln \pi_k+常数$

上式我们使用了公式(9)。两侧取指数，我们将 $q^\ast(\pmb \pi)$ 看成狄力克雷分布，

$q^\ast(\pmb \pi)=\text{Dir}(\pmb \pi|\pmb \alpha)\tag {11}$

其中 $\pmb \alpha$ 的元素为 $\alpha_k$，形式为

$\alpha_k=\alpha_0+N_k$

最后，变分后验概率分布 $q^\ast(\pmb \mu_k,\pmb \Lambda_k)$ 无法分解成边缘概率分布的乘积，但是总可以使用概率的乘积规则，将其写成 $q^\ast(\pmb \mu_k,\pmb \Lambda_k)=q^\ast(\pmb \mu_k|\pmb \Lambda_K)q^\ast(\pmb \Lambda_k)$ 。两个银子可以通过观察公式(10)得到，并且可以读出 $\pmb \mu_k$ 和 $\pmb \Lambda_k$。与预期相符，结果是⼀个⾼斯-Wishart分布，形式为

$q^\ast(\pmb \mu_k,\pmb \Lambda_k) = \mathcal N(\pmb \mu_k|\pmb m_k,(\beta_k\pmb\Lambda_k)^{-1})\mathcal W(\pmb\Lambda_k|\pmb W_k,\nu_k)\tag {12}$

其中定义了

$\beta_k=\beta_0+N_k\\ \pmb m_k=\frac{1}{\beta_k}(\beta_0\pmb m_0+N_k\pmb {\bar x}_k)\\ \pmb W_k^{-1}=\pmb W_0^{-1}+N_k\pmb S_k+\frac{\beta_0N_k}{\beta_0+N_k}(\pmb {\bar x}_k-m_0)(\pmb {\bar x}_k-m_0)^T\\ \nu_k=\nu_0+N_k$

更新⽅程类似于混合⾼斯模型的最⼤似然解的EM算法的M步骤的⽅程。我们看到，为了更新模型参数上的变分后验概率分布，必须进⾏的计算涉及到的在数据集上的求和操作与最⼤似然⽅法中的求和操作相同。

为了进⾏这个变分M步骤，我们需要得到表⽰“责任”的期望 $\mathbb E[z_{nk}] = r_{nk}$。这些可以通过对公式(8)给出的 $ρ_{nk}$ 进⾏归⼀化的⽅式得到。我们看到，这个表达式涉及到关于变分分布的参数求期望，这些期望很容易求出，即（下面三个式子应该都是分布的性质？？？）

$\mathbb E_{\pmb \mu_k,\pmb \Lambda_k}[(\pmb x_n-\pmb \mu_k)^T\pmb\Lambda_k(\pmb x_n-\pmb \mu_k)]=D\beta_k^{-1}+\nu_k(\pmb x_n-\pmb m_k)^T\pmb W_k(\pmb x_n-\pmb m_k)\\ \ln \tilde \Lambda_k \equiv \mathbb E[\ln|\pmb \Lambda_k|]=\sum_{i=1}^D\psi(\frac{\nu_k+1-i}{2})+D\ln 2+\ln |\pmb W_k|\\ \ln \tilde \pi_k\equiv \mathbb E[\ln \pi_k]=\psi(\alpha_k)-\psi(\hat \alpha)$

其中我们引⼊了 $\tilde Λ_k$ 和 $\tilde π_k$ 的定义，$\psi(·)$ 是Digamma函数， $\hat α = \sum_k α_k$。上式是从Wishart分布和狄利克雷分布的标准性质中得到的。

将上式代入公式(8)，再归一化，得：

$r_{nk}\propto \tilde\pi_k\tilde \Lambda_k^{\frac{1}{2}}\exp\left\{-\frac{D}{2\beta_k}-\frac{\nu_k}{2}(\pmb x_n-\pmb m_k)^T\pmb W_k(\pmb x_n-\pmb m_k)\right\}$

注意这个结果与最⼤似然EM算法得到的“责任”的对应结果的相似性，后者可以写成，

$r_{nk}\propto \pi_k|\pmb\Lambda_k|^{\frac{1}{2}}\exp\left\{-\frac{1}{2}(\pmb x_n-\pmb \mu_k)^T\pmb\Lambda_k(\pmb x_n-\pmb \mu_k)\right\}$

其中我们使用精度代替了协方差，来强调两条公式之间的相似性。

因此变分后验概率分布的最优化涉及到在两个阶段之间进⾏循环，这两个阶段类似于最⼤似然EM算法的E步骤和M步骤。在变分推断的与E步骤等价的步骤中，我们使⽤当前状态下模型参数上的概率分布来计算上三式中的各阶矩，从⽽计算 $\mathbb E[z_{nk}] = r_{nk}$。然后，在接下来的与M步骤等价的步骤中，我们令这些“责任”保持不变，然后使⽤它们通过公式(11)和(12)重新计算参数上的变分分布。在任何⼀种情形下，我们看到变分后验概率的形式与联合概率分布中对应因⼦的函数形式相同。这是⼀个⼀般的结果，是由于选择了共轭先验所造成的。

下图给出了将这种⽅法应⽤于⽼忠实间歇喷泉数据集上的结果。使⽤的模型是⾼斯混合模型，有 $K = 6$ 个分量。

我们看到，在收敛之后，只有两个分量的混合系数的期望值可以与它们的先验值区分开。这种效果可以根据贝叶斯模型中数据拟合与模型复杂度之间的折中来定性地理解。这种模型中的复杂度惩罚的来源是参数被推离了它们的先验值。对于解释数据点没有作⽤的分量满⾜ $r_{nk}\simeq 0$，从而 $N_k\sim 0$。根据 $\alpha_k=\alpha_0+N_k$，可以看到 $\alpha_k\simeq\alpha_0$。所以其他的参数回到了\趋近于它们的先验值。原则上，这些分量会微⼩地适应于数据点，但是对于⼀⼤类先验分布来说，这种微⼩的调整的效果太⼩了，以⾄于⽆法在数值上看出来。对于⾼斯混合模型，后验概率分布中的混合系数的期望值为

$\mathbb E[\pi_k]=\frac{\alpha_0+N_k}{K\alpha_0+N}$

考虑一个分量，其中 $N_k\simeq 0$ 且 $\alpha \simeq \alpha_0$。如果先验概率分布很宽，从而 $\alpha_0\rightarrow 0$，那么 $\mathbb E[\pi_k]\rightarrow 0$，分量对模型不起作用。而如果先验概率与混合系数密切相关，即 $\alpha_0\rightarrow \infty$，那么 $\mathbb E[\pi_k]\rightarrow \frac{1}{K}$。

上图中，混合系数上的先验概率分布是⼀个狄利克雷分布，对于 $α_0 < 1$，先验概率分布倾向于选择某些混合系数趋近于零的解。上图是使⽤ $α_0 = 10^{-3}$ 得到的结果，产⽣了两个混合系数⾮零的分量。如果我们选择 $α_0 = 1$，那么我们得到三个混合系数⾮零的分量，对于 $α = 10$，所有六个分量的混合系数都不等于零。

正如我们已经看到的那样，⾼斯分布的贝叶斯混合的变分解与最⼤似然的EM算法的解很相似。事实上，如果我们考虑 $N \rightarrow \infty$ 的极限情况，那么贝叶斯⽅法就收敛于最⼤似然⽅法的EM解。对于不是特别⼩的数据集来说，⾼斯混合模型的变分算法的主要的计算代价来⾃于“责任”的计算，以及加权数据协⽅差矩阵的计算与求逆。这些计算与最⼤似然EM算法中产⽣的计算相对应，因此使⽤这种贝叶斯⽅法⼏乎没有更多的计算代价。然⽽，这种⽅法有⼀些重要的优点。⾸先，在最⼤似然⽅法中，当⼀个⾼斯分量“退化”到⼀个具体的数据点时，会产⽣奇异性，⽽这种奇异性在贝叶斯⽅法中不存在。实际上，如果我们简单地引⼊⼀个先验分布，然后使⽤MAP估计⽽不是最⼤似然估计，这种奇异性就会被消除。此外，当我们在混合分布中将混合分量的数量 $K$ 选得较⼤时，不会出现过拟合问题，正如我们在图中看到的那样。最后，变分⽅法使得我们可以在确定混合分布中分量的最优数量时不必借助于交叉验证的技术。

总结

一开始看这个是想搞懂变分推理的思想和原理，后面就钻着PRML看了。思想和原理现在懂了，但是PRML后面还有很多很多很多很多例子，后面有空再看看吧。