变分推断

为什么需要变分推断？

在贝叶斯推断问题中，对于复杂高维的后验分布难以进行解析计算，前面小节介绍了如何使用马尔可夫-蒙特卡洛算法（MCMC）避免复杂的解析计算，通过马尔可夫链和采样的方式获得后验分布的样本。然而，MCMC面对高维参数空间，存在收敛速度慢、计算成本高、收敛判断难等问题。因此，在机器学习研究与实际应用中，变分推断成为更高效的替代方案。从理论上讲，很多时候我们永远无法知道真实的生成模型是什么，使用近似分布作为代替是合理的方案。

变分推断是什么？

变分推断是一种近似推断技术，通过求解特定的优化问题寻找一个近似的后验分布。不同于解析求解与MCMC对后验分布 $p(\theta|d)$ 进行准确估计，变分推断的核心思想是引入一个简单的参数分布 $q(\theta)$，来近似复杂的真实后验分布。这种方法将求解计算后验分布问题转化成了优化问题。

现在，让我们具体介绍变分优化的概念是如何应用于推断问题中。假设我们有一个贝叶斯模型，每个参数都有给定的先验分布。这个模型具有一些潜变量及参数，我们将模型的潜变量与参数的集合记为 $θ$， 将所有观测变量记为$d$，概率模型指定了联合分布 $p(\theta,d)$。我们做贝叶斯推断的目标是求解后验分布 $p(\theta|d)$，由于后验分布有时难以计算，我们就会使用变分推断，目标是找到后验分布的近似分布。

$$p(θ|d) = \frac{p(θ,d)}{p(d)} \qquad \qquad \tag{1}$$

该公式也能写成另一个形式

$$p(d)= \frac{p(θ,d)}{p(θ|d) } \qquad \qquad \tag{2}$$

此时在公式右边分子分母同时引入一个分布 $q(θ)$，无论 $q(θ)$是什么，只要不为0，等式依然成立

$$p(d)= \frac{p(θ,d)/q(θ))}{p(θ|d) /q(θ)} \qquad \qquad \tag{3}$$

两边同时求对数，

$$log(p(d))= log(\frac{p(θ,d)}{q(θ)})-log(\frac{p(θ|d)} {q(θ)}) \qquad \qquad \tag{4}$$

在公式两边，我们对 $q(θ)$求期望

$$\int log(p(d)) q(θ)dθ =\int log(\frac{p(θ,d)}{q(θ)}) q(θ)dθ-\int log(\frac{p(θ|d)}{q(θ)})q(θ)dθ$$

$$log(p(d)) =\int log(\frac{p(θ,d)}{q(θ)}) q(θ)dθ+\int log(\frac{q(θ)}{p(θ|d)})q(θ)dθ \qquad \qquad \tag{5}$$

其中，公式（5）右边第二项为 $q(θ)$与 $p(θ|d)$两个分布之间的KL散度

$$D_{KL}(q(θ)||p(θ|d))=\int log(\frac{q(θ)}{p(θ|d)})q(θ)dθ \qquad \qquad \tag{6}$$

同时，定义公式右边第一项为

$$L=\int log(\frac{p(θ,d)}{q(θ)}) q(θ)dθ \qquad \qquad \tag{7}$$

由此可得，

$$log(p(d)) = L + D_{KL}(q(θ)||p(θ|d)) \qquad \qquad \tag{8}$$

我们的目标是求解后验概率分布 $p(θ|d)$，由于这个分布很难计算，我们引入了另一个关于 $θ$ 的分布 $q(θ)$，那么我们只要让 $q(θ)$尽可能地去近似 $p(θ|d)$，就可以使用 $q(θ)$代替 $p(θ|d)$作为我们要求的解。让 $q(θ)$尽可能近似 $p(θ|d)$，就是要最小化KL散度$D_{KL}(q(θ)||p(θ|d))$。不过直接从公式（6）求解最小值是不可行的，其中$p(θ|d)$正是需要求解的部分。考虑到 $log(p(d))$是一个定值，根据公式（8）最小化KL散度的求解可以等价转化为最大化 $L$的优化问题，这里的 $L$就叫做证据下界 （Evidence Lower Bound, ELBO）。

$$\hat{q(θ)} =\argmax_{q(θ)}L \\ = \argmax_{q(θ)} \int log(\frac{p(θ,d)}{q(θ)}) q(θ)dθ \\=\argmax_{q(θ)}\int log(\frac{p(θ)*p(d|θ)}{q(θ)}) q(θ)dθ \\=\argmax_{q(θ)} \int (logp(θ)+logp(d|θ)-logq(θ))q(θ)dθ \qquad \qquad \tag{9}$$

当ELBO取得最大值时，相应的 $q(θ)$便是最接近 $p(θ|d)$的分布。如果我们允许 $q(\theta)$的任意取值（不加任何约束），那么证据下界的最大值将会发生在 $q(θ)$与后验分布 $p(θ|d)$相等的时候，也就是KL散度为0时。然而，我们在实际运算中面对难以解析的后验分布，如果不对 $q(\theta)$的分布进行一定约束，将难以用参数化的方法对 $q(\theta)$进行优化。

因此，我们需要考虑一类有限的分布族 $q(\theta)$，接着从中寻找使KL散度最小化的分布。我们添加约束的目标是充分地限制分布族，使它们仅包含可解析的分布，同时要保证分布族足够丰富与灵活，以此很好地近似于真实后验分布。必须要强调的是，施加约束纯粹是为了保证分布达到可解析的程度，在这要求之下我们应该使用尽可能丰富的近似分布族。尤其是，使用高度灵活的分布并不存在过拟合的现象，使用更多的灵活分布近似只会让我们更接近真实的后验分布。

其中一种限制近似分布族的方法是，采用一种参数化分布 $q(\theta|\omega)$，由一系列参数 $\omega$控制。证据下界L将变成 $\omega$的函数 $L(\omega)$，我们便可采用标准的非线性优化技术来求解这些参数的最优取值。这个方法的一个例子是，采用高斯分布作为近似，我们可以对高斯分布的均值和方差这两个参数进行优化。

变分推断实例：高斯分布近似

接下来，我们介绍一个变分推断的编程实例，以高斯分布作为近似分布。假设目标采样分布p(x)为一个混合高斯分布，由于目标采样分布作为后验分布是已知的，我们可以利用公式6，直接对KL散度进行最小化的优化求解。

def p(x):
return ((0.3 * np.exp(-(x-0.3)**2))+0.7*np.exp(-(x-2.)**2/0.3))/1.2113

$$p(x)=\frac{0.3 * e^{-(x-0.3)^2}+0.7*e^{\frac{-(x-2)^2}{0.3}}}{1.2113}$$

我们尝试采用一个高斯分布q(x)去近似p(x)，假设分布的均值为u，标准差为 $\sigma$

$$q(x)=\frac{e^{-(x-u)^2}}{\sigma^2}$$

对于给定的参数组合（u和 $\sigma$），我们可以使用蒙特卡洛方法对KL散度进行计算。我们从 $q(x|u,\sigma)$进行采样，假设抽取500个符合分布的样本 $X$，即可计算每个样本 $x_i$在两个分布中概率密度的对数差，对所有样本的概率密度对数差求平均，即为KL散度的离散形式。

$$D_{KL}=\frac{\sum_i{logq(x_i)-logp(x_i)}}{samplesize}$$

eps = np.finfo(float).eps #一个大于0的极小值ε，保证在进行对数计算时大于0.

def lossfun(params):
    samplesize = 500 #蒙特卡洛法的采样数
    # q-高斯分布的均值与标准差
    u = params[0]
    sigma = params[1]
    # 从q分布中采样
    ss = norm.rvs(loc=u,scale=sigma,size = samplesize)
    # q中所采样本在p分布的概率密度
    pp1 = p(ss) * 0.999 + eps 
    # q中所采样本在q分布的概率密度
    pp2 = norm.pdf(ss,loc=u,scale=sigma) * 0.999 + eps 
    
    dd = np.log(pp2)-np.log(pp1) 
    return dd.sum()/samplesize

接下来便可使用一般参数优化的方式，在参数空间 $(u,\sigma)$中求解KL散度最小对应的参数。

res = minimize(fun=lossfun,x0=(1,1),method='Powell',bounds=((-3,3),eps,3))
print(res.x)

参数优化的结果为 $(u=0.817,\sigma=1.034)$，两个分布的示意图如下。

x = np.arange(-3., 5., 0.01)
plt.plot(x, p(x), color='r', label='p(x)')
plt.plot(x, norm.pdf(x,loc=res.x[0], scale=res.x[1]), color='b', label='q(x)')
plt.legend()
plt.ylabel('pdf(x)')

变分推断在参数估计中的应用

如何使用变分推断的方法估计Weibull心理物理曲线的阈值

（此处简单讲解Weibull心理物理曲线）

Weibull函数为

$$y=1-(1-\gamma)exp^{-(\frac{k*x}{\alpha})^\beta}$$

其中

$$k=-log(\frac{1-g}{1-\gamma})^\frac{1}{\beta}$$

$\alpha$为该正确率下的阈值，$\beta$为函数的斜率，$\gamma$为基线概率水平的概率值， $g$为阈值对应的正确率。

假设我们当前的参数组合 $(\alpha,\beta,\gamma,g)$的真值为（0.25，3，0.5，0.82），在实验设置的一致性水平中对应的心理物理曲线如下，

我们根据实验设置与Weibull函数，利用伯努利随机过程生成了5000个模拟数据，即5种一致性水平下，各生成1000个试次数据，每个试次数据代表被试做对（1）和做错（0）。接下来，我们尝试对这批模拟数据进行参数估计。为了简化例子，这里我们只假定阈值 $\alpha$是需要估计的自由参数，其它参数皆为已知的固定参数。

根据二项选择任务常用的对数似然公式，计算出模拟数据在给定参数的Weibull函数中的似然值

$$log(p(d|\theta))=\sum_ix_i*lnp+(1-x_i)*ln(1-p)$$

eps = np.finfo(float).eps #大于0的极小数，防止计算对数时出现0
def loglikeli(alpha): # 计算模拟数据在给定参数的Weibull函数中的对数似然值
    g = 0.82
    beta = 3
    gamma = 0.5

    prob = np.zeros(cohTrial.size)
    
    prob = weibullfun(cohTrial, alpha, beta, gamma, g)
    prob = prob*data + (1-prob)*(1-data)
    prob = prob * 0.999 + eps #防止后面计算log时出现0值
   
    return np.log(prob).sum()

接着我们从假定 $q(x)$分布为高斯分布，并从中对参数 $\alpha$的值进行采样。接着，分别计算出 $logp(\alpha), logp(data|\alpha),logq(\alpha)$，根据公式9汇总并对所有采样进行平均得到 $L$值。

samplesize = 500 # 从q(x)里面取多少
from scipy.stats import norm,uniform
eps = np.finfo(float).eps

def lossfun(params):
    # q分布的参数
    u = params[0] # q高斯函数的均值
    sigma = params[1] # q高斯函数的std
    
    # 对q分布进行采样
    ss = norm.rvs(loc=u, scale=sigma, size=samplesize) # ss是(-inf, inf) # 从q分布中采500个样本
    # 计算这些样本的 logq
    llq = np.log(norm.pdf(ss, loc=u, scale=sigma)*0.999 + eps) li
    
    # 将采样的范围(-inf, inf)转换为系数的阈限（0，1）内
    ss = 1./(1+np.exp(-ss)) 
    llp = [loglikeli(s) for s in ss] # 求这些log(p(data|alpha))
    llp = np.array(llp) # 对于所有的ss，计算针对data的likelihood
    
    #计算先验值 log(p(alpha))
    llprior = np.log(uniform.pdf(ss,loc=0, scale=1)*0.999+eps) 

    # 构建优化目标，根据公式9汇总得到L值，并取相反数
    dd = llq-llp - llprior 
    
    return dd.sum()/samplesize/cohTrial.size # 我们再除以trial数量，以保证返回数量不会太大

接下来优化这个负对数似然函数，求解最小化负对数似然（等价于最大化L值）对应的参数 $\alpha$的分布，输出计算所得参数 $\alpha$的均值并画出分布

from scipy.optimize import minimize
res = minimize(fun=lossfun, x0=(0, 1), method='Powell', bounds=((-5, 0), (eps, 5))) 

print('The estimated parameters is', res.x)
# Gaussian distribution的均值转化成(0,1)的阈限为
print('The mean of threshold is', 1./(1+np.exp(-res.x[0])))

The estimated parameters is [-1.0848973   0.02625092]
The mean of threshold is 0.25258036744019224

# 完成之后, 我们画出分布
x = np.arange(-5., 0., 0.01)
xx = 1/(1+np.exp(-x))
plt.plot(xx, norm.pdf(x, loc=res.x[0], scale=res.x[1]), color='b', label='q(x)')
plt.legend()
plt.ylabel('q(x)=p(x|data)')

---

变分推断在认知理论中的应用——自由能原理

上一节已经介绍了贝叶斯脑假说，即认为大脑是通过贝叶斯推断来完成各种任务的。我们大脑认识这个世界的方式就是根据观察到的感知状态 $s$ 来对环境的潜状态 $\phi$做推断，也就是求解后验概率 $p(\phi|s)$。这是认知科学中所有贝叶斯模型的基本思想。

如同本节开头所说，由于后验概率的求解具有一定难度，我们选择用变分分布 $q(\phi)$来近似它。尽管符号有所不同，但经过同样的推导过程，我们可以得到与公式（5）类似的公式。

$$log(p(s))=\int{log(\frac{p(\phi,s)}{q(\phi)})q(\phi)d\phi}-\int{log(\frac{p(\phi|s)}{q(\phi)})q(\phi)d\phi}$$

我们令 $F(s)$等于公式右边第一项的负数，也是证据下界ELBO的负数。右边第二项，也正是变分分布 $q(\phi)$和后验概率 $p(\phi|s)$的KL散度。

$$F(s)=-\int{log(\frac{p(\phi,s)}{q(\phi)})q(\phi)d\phi}$$

$$D_{KL}[q(\phi)||p(\phi|s)]=-\int{log(\frac{p(\phi|s)}{q(\phi)})q(\phi)d\phi}$$

那么公式可以重新表示成

$$log(p(s))=-F(s)+D_{KL}[q(\phi)||p(\phi|s)]\\F(s)=-log(p(s))+D_{KL}[q(\phi)||p(\phi|s)]$$

当 $q(\phi)$和后验概率 $p(\phi|s)$越接近，两者的KL散度也就越小（注意 $D_{KL}[q(\phi)||p(\phi|s)]>=0$）。同时，虽然我们不知道 $-log(p(s))$的值是多少，但肯定是个常数，因为给定一种感知状态 $s$，外部世界产生它的概率是恒定的，只不过我们无法获悉。所以当 $D_{KL}[q(\phi)||p(\phi|s)]$越小， $F(s)$也就越小；反之也成立，如果我们优化使得 $F(s)$越小，那么 $D_{KL}[q(\phi)||p(\phi|s)]$也就越小。这里的 $F(s)$就是自由能，最小化自由能就等价于最小化KL散度。

Karl Friston所提出的最小自由能原理，本质上就是认为我们大脑在做贝叶斯推断的时候，并非直接求解复杂的后验分布，而是用变分推断的形式用一个简单的近似分布去逼近后验分布。这也符合常理，我们每天看着太阳东升西落，但我们很难从这简单的观测中推断出地球以太阳为中心转动的精确轨迹方程，更可能会认为太阳在以地球为中心做圆周运动，这正是我们大脑在用一种有偏差但简单的近似分布去认识真实的世界。

（是否还需要说明完整的公式，即客观世界模型m，大脑内在模型u等内容）

$$F(s,u) = -log(s|m) + D_{KL}[q(\phi|u)||p(\phi|s,m)] \tag{12}$$

Karl Friston的自由能理论与预测性编码、最优控制、主动推断等多个概念存在复杂的关系，这里仅是简单介绍了自由能与变分推断相关的内容，读者有兴趣的话，可以阅读参考文献

Friston, K.(2010). The free-energy principle: a unified brain theory? Nature Review Neuroscience, 11, 127–138. https://doi.org/10.1038/nrn2787