最近在看LightLDA的论文,发现自己当时看到采样相关的知识都忘的差不多了,重新翻翻记录一下。 在无法求出来具体的分布的时候,一种是采用变分法、mean field近似,另一种就是使用采样的方法了。虽然无法求出来具体的分布,我们可以利用一些方法生成服从相应分布的样本,然后利用这些样本再去拟合想要的分布。这一点在Gibbs Sampling中已经见到了,由于求不出来doc-topictoken-topic矩阵,就采样出每一个隐变量的值,然后用统计量来估算其分布。

Rejection Sampling

在已经知道分布的形式的情况下,要是想要获得相应分布的样本的话其实只要求个逆,然后就可以利用\(U\sim(0,1)\)来生成相应的分布了。

但是之前的很多例子里面,我们要不是后验分布的形式及其复杂,要不就是下面有个归一化的项\(Z\)还是带积分的,都不好求,所以,才会有更复杂的采样方法。

Rejection Sampling就是个稍微复杂一点的方法,其实想法很直观,在求不出来归一化项Z的时候,我们就取一个上界,然后在这个上界函数上采样,然后再以一个\(U \sim (0,1)\)中的随机值,来判断是不是可以接受。

具体过程,取一个proposal distribution\(q(x)\)满足存在\(Mq(x) \geq \hat{p(x)}\),其中\(\hat{p(x)}\)是为归一化的概率项。则首先依据\(x \sim q(x)\)采样\(x\),然后在从\(u \sim U(0,1)\)中采样出来一个\(u\),如果\(u \gt \frac{\hat{p(x)}}{Mq(x)}\),那么就拒绝这个样本,不然的话,我们才接受它。如下图

enter image description here

enter image description here

\(Mq(x)\)是一个上界,我们只有在\(uMq(x) \leq p(x)\)的时候才会接受它。下面是证明 定义\(S = \{(x,u): u \leq \hat{p}(x)/Mq(x)\}\)以及\(S_0 = \{(x,u): x \leq x_0 , u \leq \hat{p}(x)/Mq(x)\}\),那么

\[\begin{split} P(x \leq x_0 | x \ \ accepted) & = \frac{P(x \leq x_0, x \ \ accepted)}{P(x \ \ accepted)} \\ & = \frac{\int \int \mathbb{I}((x,u) \in S_0)q(x)dudx}{\int \int \mathbb{I}((x,u) \in S)q(x)dudx} \\ & = \frac{\int_{- \infty}^{x_0}\hat{pc}(x)dx}{\int_{-\infty}^{\infty}\hat{p}(x)dx} \end{split}\]

其实\(S\)可以看成所有满足条件的接受条件的点,而\(S_0\)则是同时保证\(x \lt x_0\)的那部分,最终得到的正好就是\(p(x)\)的cdf。

还有一个问题就是我们有多大的概率能得到一个可以接受的结果呢?

\[ p(accepted) = \int \frac{\hat{p}(x)}{Mq(x)}q(x)dx = \frac{1}{M} \int \hat{p}(x)dx \]

这个结果也是很符合直觉,符合条件的M越小越好,当M趋向于正无穷的时候,我们几乎不可能得到合适的采样,我们希望M尽可能的小,来提高接受率。

当我们想要对后验分布进行采样的时候,即\(p(\theta | D) = p(D|\theta)p(\theta)/p(D)\),采用\(p(\theta)\)作为我们的proposal,选取\(M = p(D | \hat{\theta})\)其中\(\hat{\theta} = arg\ max\ p(D|\theta)\),这样接受的接受率为

\[\frac{\hat{p}(\theta)}{Mq(\theta)} = \frac{p(D | \theta)}{p(D | \hat{\theta})}\]

但是在高维的情况下,Rejection sampling几乎是不可用的,比如我们希望从一个\(\mathcal{N}(0, \sigma^2_p I)\)中采样,为此我们使用一个分布\(\mathcal{N}(0, \sigma^2_q I)\),这样我们选择的最优化的M应该为\((\sigma_q/\sigma_p)^D\),之前已经看到了接受率与M成反比,所以显然接受率会随着维度的提高指数级的下降,这显然是不可以接受的。

MCMC

MCMC(Markov chain Monte Carlo)方法通过构造一个状态空间上的马尔科夫链,然后利用随机游走以及细致平稳分布的性质,对目标分布进行采样。之前的Gibbs Sampling就是这样一个例子,每个隐变量都是具有多个状态,在模型burn-in之后, 每一次迭代其实就是整个模型在其后验分布的状态空间随机游走,我们通过采样的得到的经验分布来近似其真实分布。

MCMC也是经常用在PGM上的解法,而且最近越来越火吧,它有很多优点 。

  • 好求解。对比Gibbs Sampling和VBEM的推倒就能发现,Gibbs Sampling要好求的多,解法形式都比较固定。实现也比较方便

  • 在大数据集上速度更快,采样出来的部分基本只依赖其马尔科夫毯上的部分,传递的时候也可以只传递充分统计量。

  • 对于一些先验非共轭或者中间有个什么广义线性模型的情况下,采样的方法也更好求解。

  • 尤其是最近PS等各种分布式的机器学习框架比如Petuum上的lightLDA,还有腾讯的Peacock都是采用采样的方法求解的。

Metropolis Hasting algorithm

之前已经见过的Gibbs Samping是利用full conditionals来计算,利用共轭的性质来优化,但是有一些限制。

  • 后验分布的形式未知。很多情况下都是非共轭分布引起的。

  • 后验分布的参数很多,grid approximations的方法不可行。

  • 一些full conditionals并不是已知的形式。这种情况下Gibbs Sampling也不可用。

这种时候就需要Metropolis Hastings来救场了。大概过程如下:

  1. 选择一个初始值\(\theta^{(0)}\)

  2. 在第\(t\)次迭代,从Jumping distribution\(J_t(\theta^*|\theta^{(t - 1)})\)生成一个候选的\(\theta^*\)

  3. 计算接受率: \[ r = \frac{p(\theta^*|y)/J_t(\theta^* | \theta^{t-1}}{p(\theta^{t-1}|y)/J_t(\theta^{t - 1}|\theta^*)} \]

  4. 以min(r, 1)的概率接受\(\theta^*\)\(\theta^t\),如果不接受,就令\(\theta^t = \theta^{t - 1}\)

  5. 重复2-4,获取相应的分布\(p(\theta|y)\),要考虑burn-in的时间。

如果要求转移概率分布J为对称的,就称为random walk Metropolis sampling。如果每一步并不依赖前一步,即\(J_t(\theta^* | \theta^{t - 1}) = J_t(\theta^*)\)就称为independent Metropolis-Hasting sampling

在计算接受率的时候,当J对称时,可以化简为\(r = \frac{p(\theta^*|y)}{p(\theta^{t-1}|y)}\)。这里就可以看出来MH的好处,由于直接计算比例,所以不需要知道归一化项。其实这个加上下一步的以概率接受,可以看出如果我们这次得到的候选的\(\theta\)比之前的要好,那就可以直接接受,不然就以一个概率的比例接受,这里只要生成一个(0-1)之间的随机数,在其比\(r\)要小的时候接受就行了。也就是我们始终接受一个比之前可能性大的\(\theta\)。这个地方就是和rejection sampling不同的地方了,MH保证每次都会产生一个样本,而rejection sampling并不能保证这一点,所以采样率更低。

同时,采样的acceptance rate可以用来评判Metropolis-Hastings算法效果。比如,接受率太低,那采样可能不够充分,接受率太高,可能马尔科夫链拟合的并不好。一般来说random walk应该处于0.25到0.5之间,而independent应该尽可能接近1。

Gibbs Sampling is a special case of MH

对于Gibbs Sampling可以看成proposal distribution\(p(\theta^* | \theta) = p(\theta^*_i| \theta^{t-1}_{-i}) \mathbb{I}(\theta^*_{-i} = \theta_{-1})\)的MH采样。即每次只根据条件分布更新\(\theta_i\)而其他的\(\theta\)保留不变。

其接受率为

\[\begin{split} \alpha &= \frac{p(\theta^*)p(\theta|\theta^*)}{p(\theta)p(\theta|\theta^*)} = \frac{p(\theta^*_i|\theta^*_{-i})p(\theta^*_{-i})p(\theta_i|\theta_{-1}^*)}{p(\theta_i|\theta_{-i})p(\theta_{-i})p(\theta^*_i|\theta_{-1})} \\ & = \frac{p(\theta_i^*|\theta_{-i})p(\theta_{-i})p(\theta_i|\theta_{-1})}{p(\theta_i|\theta_{-i})p(\theta_{-i})p(\theta^*_i|\theta_{-1})} = 1 \end{split}\]

所以每次,Gibbs Sampling的结果都是接受的。

Proposal Distributions

对于给定一个目标分布\(p^*\),一个可行的proposal distribution是指可以以非0的概率转移到目标分布的任何非零的状态。比如,高斯分布的random walk没有概率为0的位置,所以对于任何continuous state space都是合法的。

如果实在是不确定使用哪个proposal,也可以采用多个凸的proposal加权一起使用。即

\(p(\theta^*|\theta) = \sum_{k = 1}^K w_kp_k(\theta^*|\theta)\)

LightLDA

lightLDA中提出了加速采样的方法。

回顾平时Gibbs Sampling采样的的式子

\[p(z_{di=k}| rest) \propto \frac{(n_{kd}^{-di} + \alpha_k)(n_{kw}^{-di}+\beta_w)}{n_k^{-di} + \bar{\beta}}\]

其中\(\bar{\beta} := \sum_w \beta_w\)。诸如SparseLDAAliasLDA都是对这个形式进行改进分解采样。LightLDA中将这个分解成

\[q(z_{di} = k | rest) \propto (n_{kd} + \alpha_k) * \frac{n_{kw} + \beta_w}{n_k + \bar{\beta}}\]

其中前面一部分只和文章相关可以看做doc-proposal,第二部分只和词有关系word-proposal部分。从直觉上来讲,合理的主题应该是在这两部分都有很大的概率。

Word-Proposal部分

重新定义word-proposal

\[p_w(k) \propto \frac{n_{kw}+\beta_w}{n_k + \bar{\beta}}\]

从状态s转义到t的概率可以写成

\[min\{1, \frac{p(t)p_w(s)}{p(s)p_w(t)}\}\]

定义\(\pi_w := \frac{p(t)p_w(s)}{p(s)p_w(t)}\)

\[ \pi_w = \frac{(n_{td}^{-di} + \alpha_t)(n_{tw}^{-di} + \beta_w)(n_s^{-di} + \bar{\beta})(n_{sw} + \beta_w)(n_t + \bar{\beta})}{(n_{sd}^{-di} + \alpha_s)(n_{sw}^{-di} + \beta_w)(n_t^{-di}+\bar{\beta})(n_{tw} + \beta_w)(n_s + \bar{\beta})} \]

只要保存了其所有的充分统计量\(n\),每一次完成采样之后在O(1)的时间里就可以计算出接受概率。为了采样\(p_w\)的复杂度保持在O(1),这里还采用了一个alias table的技巧,具体的可以看这里,挺有意思的技巧。简单来说,就是把多项分布平整为均匀分布,然后通过常数复杂度的采样得到。然后不过,由于使用了alias table的技巧,需要使用2K的空间来存储alias table,这种时候,稀疏性就又排上用场了。将\(p_w = \frac{n_{kw}}{n_k + \beta}+\frac{\beta_w}{n_k + \beta}\),这个同样也是稀疏的可以用SparseLDA里面的方法简化,这样空间复杂度也降到\(O(K_w)\)

Doc-Proposal部分

定义\(p_d\)

\[ p_d(k) \propto n_{kd} + \alpha_k \]

转义概率为

\[ min\{1,\frac{p(t)p_d(s)}{p(s)p_d(t)}\} \]

\(\pi_d := \frac{p(t)p_d(s)}{p(s)p_d(t)}\)。而且,其中\(n_{kd}\)本来就是利用\(z_{di}\)计算的,直接从\(1-n_d\)中取一个就行了。

总结

只看加速采样算法这一部分,采用Metropolis Hasting采样通过灵活的构造Proposal分布可以更好地利用模型的性质,混合使用两种分布也可以相互补充,更快的遍历所有的状态空间,更好地提高采样mixing效果。




X