我们首先看看概率编程最简单的实例:贝叶斯线性回归。
线性回归的基本公式为
y=wx+b
。贝叶斯线性回归之前已经有相关的博客说明。从该博客中看出,贝叶斯方法在这里的作用是定义了w的先验:
w∽N(0,σ2)
。从另外一个角度看,等价于添加L2正则化项:
argmin12(y−wx−b)2+|w|2
。求解上述最优化方程比较简单。如果b的先验也服从正太分布
b∽N(0,σ2)
,那么最优化方程就变成了
argmin12(y−wx−b)2+|w|2+|b|2
。tensorflow求解起来还是比较简单的。
传统贝叶斯线性回归的局限
如果w、b的先验不服从正太分布,我们就不容易写出类似上面的最优化方程,也就无从求解。而概率编程则可以部分解决这些问题。
下面的例子还是使用最简单的正太分布,效果最显著。
首先写出概率分布方程组:
根据概率分布方程组,我们可以看出方程组可以归为一个泛函方程。概率编程上求解这类方程的方法主要有两大类:变分推断和蒙特卡洛方法。变分推断计算复杂,不过较为精确。我们来看看如何求解。
泛函方程有解析解的情况并不多,一般求解数值近似解。对于一个复杂的泛函分布,我们可以用多个简单的独立分布来拟合,进行近似推断。请注意独立这个要求,可以类比为重新构建一类基空间。这一大类方法称为 平均场。分布是函数,用简单独立分布去拟合真实分布,其实就是用简单基函数去构建真实复杂函数。衡量函数近似度的一个基准就是大名鼎鼎的KL散度/KL距离/相对熵。其定义为 KL(q||p)=−∫q(Z,λ)lnp(Z|X)q(Z,λ)dZ , λ 是q分布的参数,确定 λ 就知道q的表达式,在贝叶斯学派中被视为随机变量。即用分布q逼近分布p的程度。请注意, 分布q的形式是已知的。下面数学家的脑洞来了。
变分推断
我们不知道样本Y的真实分布,但是客观真理是不会改变的,所以
p(Y)
是一个未知的常量。对
p(Y)
求对数,加入KL散度,于是得到:
其中C是常数,ELBO表示Evidence Lower Bound Objective。我们的目标是最小化KL散度,其中 p(Z|Y) 有点麻烦,因为:
好在 ELBO(λ) 和KL是此消彼长的关系,其和为常量,所以转而最大化 ELBO(λ) ,进而求取 λ ,由 λ 得到q(Z)。又有 q(Z)=∏q(Zi) ,各自独立,
这么多未知量,一个求解方法就是先优化一类变量 Zi,q(Zi) ,固定其他变量为常量,依次优化。上述方程变为:
合并在一起,我们看出来这长得很像KL散度的负数形式:
这个计算能得到较为精确的结果,不过解析式太复杂了,尤其是面对非正太分布的时候,计算简直了不得。所以实际上是提供了ELBO可解的证明。
Tensorflow Edward与黑盒变分推断
Tensorflow Edward加入了自由参数
λ
,对ELBO的最优化求解有两种方式,将蒙特卡洛和梯度下降法结合起来,这正是多核CPU和GPU擅长干的事情。
对
λ
求导,得到
∇
梯度:
于是 ∇ELBO 变成一个求取期望的过程。在实际工程中,解析求解该过程是很难的;不过既然是求期望,那么我随机采样的话就能估计出一个期望值。期望看做加权平均,换到蒙特卡洛方法,那么就是根据 q(Z,λ) 的概率来进行采样,取出S个 Zs 值后求平均即为期望。那么有
接下来我们看看先验信息是怎么使用的。根据(3)式,我们把w,b看做 Z1,Z2 ,那么改写为 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.