【机器学习】Tensorflow概率编程: 贝叶斯线性回归、变分贝叶斯与黑盒变分推断

　　我们首先看看概率编程最简单的实例：贝叶斯线性回归。
线性回归的基本公式为$y=wx+b$。贝叶斯线性回归之前已经有相关的博客说明。从该博客中看出，贝叶斯方法在这里的作用是定义了w的先验：$w \backsim N(0,\sigma^2)$。从另外一个角度看，等价于添加L2正则化项: $argmin{\frac{1}{2}(y-wx-b)^2+|w|^2}$。求解上述最优化方程比较简单。如果b的先验也服从正太分布$b\backsim N(0,\sigma^2)$，那么最优化方程就变成了$argmin{\frac{1}{2}(y-wx-b)^2+|w|^2+|b|^2}$。tensorflow求解起来还是比较简单的。

传统贝叶斯线性回归的局限

如果w、b的先验不服从正太分布，我们就不容易写出类似上面的最优化方程，也就无从求解。而概率编程则可以部分解决这些问题。
下面的例子还是使用最简单的正太分布，效果最显著。
　　首先写出概率分布方程组：

$$p(w)\backsim N(0,\sigma_w^2) \tag{1}$$ $$p(b)\backsim N(0,\sigma_b^2) \tag2$$ $$p(y|w,b)\backsim \prod_1^N N(wx+b,\sigma_y^2) \tag{3}$$ 综上，w、b是所谓的隐变量；先验分布需要已知，即 $\sigma_w,\sigma_b，\sigma_y$均已知。
根据概率分布方程组，我们可以看出方程组可以归为一个泛函方程。概率编程上求解这类方程的方法主要有两大类：变分推断和蒙特卡洛方法。变分推断计算复杂，不过较为精确。我们来看看如何求解。
　　泛函方程有解析解的情况并不多，一般求解数值近似解。对于一个复杂的泛函分布，我们可以用多个简单的独立分布来拟合，进行近似推断。请注意独立这个要求，可以类比为重新构建一类基空间。这一大类方法称为 平均场。分布是函数，用简单独立分布去拟合真实分布，其实就是用简单基函数去构建真实复杂函数。衡量函数近似度的一个基准就是大名鼎鼎的KL散度/KL距离/相对熵。其定义为 $KL(q||p)=-\int{q(Z,\lambda)\ln{\frac{p(Z|X)}{q(Z,\lambda)}}dZ}$, $\lambda$是q分布的参数,确定 $\lambda$就知道q的表达式，在贝叶斯学派中被视为随机变量。即用分布q逼近分布p的程度。请注意， 分布q的形式是已知的。下面数学家的脑洞来了。

变分推断

我们不知道样本Y的真实分布，但是客观真理是不会改变的，所以$p(Y)$是一个未知的常量。对$p(Y)$求对数，加入KL散度，于是得到：

$$\ln{p(Y)} = -\int{q(Z,\lambda)\ln{\frac{p(Z|Y)}{q(Z,\lambda)}}dZ}+ELBO(\lambda)\\ ELBO(\lambda)= \int{q(Z,\lambda)\ln{\frac{p(Z,Y)/p(Y)}{q(Z,\lambda)}}dZ}+\ln{p(Y)}\\=\int{q(Z,\lambda)\ln{\frac{p(Z,Y)}{q(Z,\lambda)}}dZ}-\int{q(Z,\lambda)\ln{p(Y)}dZ}+{\ln{p(Y)}} \\ =\int{q(Z,\lambda)\ln{\frac{p(Z,Y)}{q(Z,\lambda)}}dZ}-\ln{p(Y)}+{\ln{p(Y)}}=\int{q(Z,\lambda)\ln{\frac{p(Z,Y)}{q(Z,\lambda)}}dZ} \tag{4}$$
$$\rightarrow \ln{p(Y)} = C = KL(q||p)+ELBO(\lambda) \tag{5}$$
其中C是常数，ELBO表示Evidence Lower Bound Objective。我们的目标是最小化KL散度，其中 $p(Z|Y)$有点麻烦，因为： $$p(Z|Y)=\frac{p(Z,Y)}{\int{p(Z,Y)dZ}}=\frac{p(Z,Y)}{p(Y)} \tag{6}$$ $p(Y)$虽然是常量，但是不能直接约去，因为函数逼近不仅要比较形状还要比较大小， $p(Y)$的大小未知，只能在形状上逼近，达不到我们的要求；分母这个积分通常不好求。这两种表达式都没有办法求解，所以直接处理KL散度是不好办的。
好在 $ELBO(\lambda)$和KL是此消彼长的关系，其和为常量，所以转而最大化 $ELBO(\lambda)$，进而求取 $\lambda$，由 $\lambda$得到q(Z)。又有 $q(Z)=\prod{q(Z_i)}$,各自独立，
$$\lambda = argmax{\int{q(Z,\lambda)\ln{\frac{p(Z,Y)}{q(Z,\lambda)}}dZ}} \\ = argmax{\int{\prod{q(Z_i,\lambda)}\ln{p(Z_i,Y)}dZ}-\int{\prod{q(Z_i,\lambda)}\ln{\prod{q(Z_i,\lambda)}}dZ}} \tag{7}$$ 其中 $q(Z_i,\lambda)$其实就是PRML第10章说的 $q(Z_i)$，这本书没有考虑 $\lambda$，硬性求解，我们来看看推导有多麻烦。
　　这么多未知量，一个求解方法就是先优化一类变量 $Z_i,q(Z_i)$，固定其他变量为常量，依次优化。上述方程变为： $$\int{q(Z_i)[\int{\ln{p(Z,Y)\prod_{i\neq j}{q(Z_j)}}dZ_j}]dZ_i}-\sum_i\int{q(Z_i)\prod_{i\neq j}{q(Z_j,\lambda)}\ln{q(Z_i,\lambda)}dZ} \tag{8}$$ 左边有 $p(Z)$所以 $dZ_j$项去不掉。右边项，把 $dZ_i$有关的项抽取出来，根据 $q_i$分布的独立性、概率总和为1的特性得到得到 $$\int{q(Z_i)\ln{q(Z_i)}dZ_i}\int{q(Z_{j1})q(Z_{j2})...q(Z_{jn})dZ_{j1}...dZ_{jn}} \\ = \int{q(Z_i)\ln{q(Z_i)}dZ_i}\int{q(Z_{j1})dZ_{j1}...\int{q(Z_{jn})dZ_{jn}}} \\ =\int{q(Z_i)\ln{q(Z_i)}dZ_i} \tag{9}$$
合并在一起，我们看出来这长得很像KL散度的负数形式： $$\int{q(Z_i)[\int{\ln{p(Z,Y)\prod_{i\neq j}{q(Z_j)}}dZ_j}]dZ_i}-\int{q(Z_i)\ln{q(Z_i)}dZ_i} \\ =\int{q(Z_i)\ln{exp^{\int{\ln{p(Z)\prod_{i\neq j}{q(Z_j)}}dZ_j}}}dZ_i}-\int{q(Z_i)\ln{q(Z_i)}dZ_i}\tag{10}$$ 接下来，最大化ELBO的思路就是把上述方程当做负数KL散度来求解，那么强制令 $q(Z_i)$逼近exp函数，得到 $$q(Z_i)=\frac{1}{N}exp^{\int{\ln{p(Z,Y)\prod_{i\neq j}{q(Z_j)}}dZ_j}}\tag{11}$$ 其中N是归一化数值。
　　这个计算能得到较为精确的结果，不过解析式太复杂了，尤其是面对非正太分布的时候，计算简直了不得。所以实际上是提供了ELBO可解的证明。

Tensorflow Edward与黑盒变分推断

Tensorflow Edward加入了自由参数$\lambda$，对ELBO的最优化求解有两种方式，将蒙特卡洛和梯度下降法结合起来，这正是多核CPU和GPU擅长干的事情。
对$\lambda$求导，得到$\nabla$梯度：

$$\nabla ELBO =\nabla \int{q(Z,\lambda)\ln{\frac{p(Z,Y)}{q(Z,\lambda)}}dZ} \\ =\int{\nabla q(Z,\lambda)\ln{\frac{p(Z,Y)}{q(Z,\lambda)}}dZ} + \int{q(Z,\lambda)[\nabla \ln{p(Z,Y)-\frac{\nabla q(Z,\lambda)}{q(Z,\lambda)}}]dZ} \\ =\int{\nabla q(Z,\lambda)\ln{\frac{p(Z,Y)}{q(Z,\lambda)}}dZ} - \nabla \int{q(Z,\lambda)dZ} \\ = \int{\nabla q(Z,\lambda)\ln{\frac{p(Z,Y)}{q(Z,\lambda)}}dZ} - \nabla 1 \\ = \int{\nabla q(Z,\lambda)\ln{\frac{p(Z,Y)}{q(Z,\lambda)}}dZ} \tag{12}$$ 数学家开了个脑洞，把方程变为一个期望方程。一个典型的期望方程为： $$E=\int{p(x)f(x)dx}\tag{13}$$ ELBO的梯度函数经过类似变换后得到： $$\nabla ELBO = \int{\frac{\nabla q(Z,\lambda)}{q(Z,\lambda)}q(Z,\lambda)\ln{\frac{p(Z,Y)}{q(Z,\lambda)}dZ }} \\ = \int{q(Z,\lambda)\nabla \ln{q(Z,\lambda)}\ln{\frac{p(Z,Y)}{q(Z,\lambda)}dZ }} \tag{14}$$ 其中 $q(Z,\lambda)$扮演 $p(x)$的角色 。
　　于是 $\nabla ELBO$变成一个求取期望的过程。在实际工程中，解析求解该过程是很难的；不过既然是求期望，那么我随机采样的话就能估计出一个期望值。期望看做加权平均，换到蒙特卡洛方法，那么就是根据 $q(Z,\lambda)$的概率来进行采样，取出S个 $Z_s$值后求平均即为期望。那么有 $$\nabla ELBO = \frac{1}{S}\sum_s{\nabla \ln{q(Z,\lambda)}\ln{\frac{p(Z,Y)}{q(Z,\lambda)}}} \tag{15}$$ 接下来就是使用梯度下降法了，为了便于求解，依然是先求解i系列变量，固定其他系列变量的套路： $$\lambda^{t+1}_i=\lambda_i^t+\mu\nabla_{\lambda_i} ELBO \tag{16}$$ 。为了求解过程快速收敛、稳定，可以添加一些moment等值，使用ada下降法等方法。
　　接下来我们看看先验信息是怎么使用的。根据(3)式，我们把w,b看做 $Z_1,Z_2$，那么改写为 $p(y|w,b)=p(Y|Z)$，然后乘以 $p(Z)$得到 $p(Y,Z)$代入(16)式即可。
公式虽然多，但是Tensorflow Edward已经写好了，最后放上一个拟合cos函数的代码，非常简洁：

import tensorflow as tf
from edward.models import Normal
import edward as ed
ed.set_seed(42)
sess = ed.get_session()
tf.global_variables_initializer().run()

W_0 = Normal(mu=tf.zeros([1, 2]), sigma=tf.ones([1, 2]))
W_1 = Normal(mu=tf.zeros([2, 1]), sigma=tf.ones([2, 1]))
b_0 = Normal(mu=tf.zeros(2), sigma=tf.ones(2))
b_1 = Normal(mu=tf.zeros(1), sigma=tf.ones(1))

x = x_train
y = Normal(mu=tf.matmul(tf.tanh(tf.matmul(x, W_0) + b_0), W_1) + b_1,
           sigma=0.1)
print(y.get_shape())

qW_0 = Normal(mu=tf.Variable(tf.zeros([1, 2])),
              sigma=tf.nn.softplus(tf.Variable(tf.zeros([1, 2]))))
qW_1 = Normal(mu=tf.Variable(tf.zeros([2, 1])),
              sigma=tf.nn.softplus(tf.Variable(tf.zeros([2, 1]))))
qb_0 = Normal(mu=tf.Variable(tf.zeros(2)),
              sigma=tf.nn.softplus(tf.Variable(tf.zeros(2))))
qb_1 = Normal(mu=tf.Variable(tf.zeros(1)),
              sigma=tf.nn.softplus(tf.Variable(tf.zeros(1))))

inference = ed.KLqp({W_0: qW_0, b_0: qb_0,
                     W_1: qW_1, b_1: qb_1}, data={y: y_train})
inference.run(n_iter=1000, n_samples=5)

在核心部分拟合的还可以，边角没拟合上，说明非线性不够。我很少见到用深度学习拟合周期函数的。

Ranganath, R., Gerrish, S., & Blei, D. (2014). Black box variational inference. In Artificial intelligence and statistics.