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

来源：互联网发布：淘宝没有访客编辑：程序博客网时间：2024/04/29 03:37

　　我们首先看看概率编程最简单的实例：贝叶斯线性回归。
线性回归的基本公式为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的先验不服从正太分布，我们就不容易写出类似上面的最优化方程，也就无从求解。而概率编程则可以部分解决这些问题。
下面的例子还是使用最简单的正太分布，效果最显著。
　　首先写出概率分布方程组：

p (w) ∽ N (0, σ 2 w) (1)

p (b) ∽ N (0, σ 2 b) (2)

p (y | w, b) ∽ \prod 1 N N (w x + b, σ 2 y) (3)

综上，w、b是所谓的隐变量；先验分布需要已知，即

σw,σb，σy均已知。
根据概率分布方程组，我们可以看出方程组可以归为一个泛函方程。概率编程上求解这类方程的方法主要有两大类：变分推断和蒙特卡洛方法。变分推断计算复杂，不过较为精确。我们来看看如何求解。
　　泛函方程有解析解的情况并不多，一般求解数值近似解。对于一个复杂的泛函分布，我们可以用多个简单的独立分布来拟合，进行近似推断。请注意独立这个要求，可以类比为重新构建一类基空间。这一大类方法称为平均场。分布是函数，用简单独立分布去拟合真实分布，其实就是用简单基函数去构建真实复杂函数。衡量函数近似度的一个基准就是大名鼎鼎的KL散度/KL距离/相对熵。其定义为

KL(q||p)=−∫q(Z,λ)lnp(Z|X)q(Z,λ)dZ,

λ是q分布的参数,确定

λ就知道q的表达式，在贝叶斯学派中被视为随机变量。即用分布q逼近分布p的程度。请注意，分布q的形式是已知的。下面数学家的脑洞来了。

变分推断

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

ln p (Y) = - \int q (Z, λ) ln p ( Z | Y ) q ( Z , λ ) d Z + E L B O (λ) E L B O (λ) = \int q (Z, λ) ln p ( Z , Y ) / p ( Y ) q ( Z , λ ) d Z + ln p (Y) = \int q (Z, λ) ln p ( Z , Y ) q ( Z , λ ) d Z - \int q (Z, λ) ln p (Y) d Z + ln p (Y) = \int q (Z, λ) ln p ( Z , Y ) q ( Z , λ ) d Z - ln p (Y) + ln p (Y) = \int q (Z, λ) ln p ( Z , Y ) q ( Z , λ ) d Z (4)

\to ln p (Y) = C = K L (q | | p) + E L B O (λ) (5)

其中C是常数，ELBO表示Evidence Lower Bound Objective。我们的目标是最小化KL散度，其中

p(Z|Y)有点麻烦，因为：

p (Z | Y) = p ( Z , Y ) \int p ( Z , Y ) d Z = p ( Z , Y ) p ( Y ) (6)

p(Y)虽然是常量，但是不能直接约去，因为函数逼近不仅要比较形状还要比较大小，

p(Y)的大小未知，只能在形状上逼近，达不到我们的要求；分母这个积分通常不好求。这两种表达式都没有办法求解，所以直接处理KL散度是不好办的。
好在

ELBO(λ)和KL是此消彼长的关系，其和为常量，所以转而最大化

ELBO(λ)，进而求取

λ，由

λ得到q(Z)。又有

q(Z)=∏q(Zi),各自独立，

λ = a r g m a x \int q (Z, λ) ln p ( Z , Y ) q ( Z , λ ) d Z = a r g m a x \int \prod q (Z i, λ) ln p (Z i, Y) d Z - \int \prod q (Z i, λ) ln \prod q (Z i, λ) d Z (7)

其中

q(Zi,λ)其实就是PRML第10章说的

q(Zi)，这本书没有考虑

λ，硬性求解，我们来看看推导有多麻烦。
　　这么多未知量，一个求解方法就是先优化一类变量

Zi,q(Zi)，固定其他变量为常量，依次优化。上述方程变为：

\int q (Z i) [\int ln p (Z, Y) \prod i \neq j q (Z j) d Z j] d Z i - \sum i \int q (Z i) \prod i \neq j q (Z j, λ) ln q (Z i, λ) d Z (8)

左边有

p(Z)所以

dZj项去不掉。右边项，把

dZi有关的项抽取出来，根据

qi分布的独立性、概率总和为1的特性得到得到

\int q (Z i) ln q (Z i) d Z i \int q (Z j 1) q (Z j 2) . . . q (Z j n) d Z j 1 . . . d Z j n = \int q (Z i) ln q (Z i) d Z i \int q (Z j 1) d Z j 1 . . . \int q (Z j n) d Z j n = \int q (Z i) ln q (Z i) d Z i (9)

合并在一起，我们看出来这长得很像KL散度的负数形式：

\int q (Z i) [\int ln p (Z, Y) \prod i \neq j q (Z j) d Z j] d Z i - \int q (Z i) ln q (Z i) d Z i = \int q (Z i) ln e x p \int ln p (Z) \prod i \neq j q (Z j) d Z j d Z i - \int q (Z i) ln q (Z i) d Z i (10)

接下来，最大化ELBO的思路就是把上述方程当做负数KL散度来求解，那么强制令

q(Zi)逼近exp函数，得到

q (Z i) = 1 N e x p \int ln p (Z, Y) \prod i \neq j q (Z j) d Z j (11)

其中N是归一化数值。
　　这个计算能得到较为精确的结果，不过解析式太复杂了，尤其是面对非正太分布的时候，计算简直了不得。所以实际上是提供了ELBO可解的证明。

Tensorflow Edward与黑盒变分推断

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

\nabla E L B O = \nabla \int q (Z, λ) ln p ( Z , Y ) q ( Z , λ ) d Z = \int \nabla q (Z, λ) ln p ( Z , Y ) q ( Z , λ ) d Z + \int q (Z, λ) [\nabla ln p (Z, Y) - \nabla q ( Z , λ ) q ( Z , λ )] d Z = \int \nabla q (Z, λ) ln p ( Z , Y ) q ( Z , λ ) d Z - \nabla \int q (Z, λ) d Z = \int \nabla q (Z, λ) ln p ( Z , Y ) q ( Z , λ ) d Z - \nabla 1 = \int \nabla q (Z, λ) ln p ( Z , Y ) q ( Z , λ ) d Z (12)

数学家开了个脑洞，把方程变为一个期望方程。一个典型的期望方程为：

E = \int p (x) f (x) d x (13)

ELBO的梯度函数经过类似变换后得到：

\nabla E L B O = \int \nabla q ( Z , λ ) q ( Z , λ ) q (Z, λ) ln p ( Z , Y ) q ( Z , λ ) d Z = \int q (Z, λ) \nabla ln q (Z, λ) ln p ( Z , Y ) q ( Z , λ ) d Z (14)

其中

q(Z,λ)扮演

p(x)的角色。
　　于是

∇ELBO变成一个求取期望的过程。在实际工程中，解析求解该过程是很难的；不过既然是求期望，那么我随机采样的话就能估计出一个期望值。期望看做加权平均，换到蒙特卡洛方法，那么就是根据

q(Z,λ)的概率来进行采样，取出S个

Zs值后求平均即为期望。那么有

\nabla E L B O = 1 S \sum s \nabla ln q (Z, λ) ln p ( Z , Y ) q ( Z , λ ) (15)

接下来就是使用梯度下降法了，为了便于求解，依然是先求解i系列变量，固定其他系列变量的套路：

λ t + 1 i = λ t i + μ \nabla λ i E L B O (16)

。为了求解过程快速收敛、稳定，可以添加一些moment等值，使用ada下降法等方法。
　　接下来我们看看先验信息是怎么使用的。根据(3)式，我们把w,b看做

Z1,Z2，那么改写为

p(y|w,b)=p(Y|Z)，然后乘以

p(Z)得到

p(Y,Z)代入(16)式即可。
公式虽然多，但是Tensorflow Edward已经写好了，最后放上一个拟合cos函数的代码，非常简洁：

import tensorflow as tffrom edward.models import Normalimport edward as eded.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_trainy = 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.

阅读全文

0 0