在上一篇文章中,我们介绍了“受限文本生成”这个概念,指出可以通过量化目标并从中采样的方式来无监督地完成某些带条件的文本生成任务。同时,上一篇文章还介绍了“重要性采样”和“拒绝采样”两个方法,并且指出对于高维空间而言,它们所依赖的易于采样的分布往往难以设计,导致它们难以满足我们的采样需求。

此时,我们就需要引入采样界最重要的算法之一“Markov Chain Monte Carlo(MCMC)”方法了,它将马尔可夫链和蒙特卡洛方法结合起来,使得(至少理论上是这样)我们从很多高维分布中进行采样成为可能,也是后面我们介绍的受限文本生成应用的重要基础算法之一。本文试图对它做一个基本的介绍。

马尔可夫链 #

马尔可夫链实际上就是一种“无记忆”的随机游走过程,它以转移概率p(yx)为基础,从一个初始状态x0出发,每一步均通过该转移概率随机选择下一个状态,从而构成随机状态列x0,x1,x2,,xt,,我们希望考察对于足够大的步数txt所服从的分布,也就是该马尔可夫链的“平稳分布”。

假设xt所服从的分布为p(t)(xt),那么xt+1所服从的分布就是
p(t+1)(xt+1)=xtp(xt+1xt)p(t)(xt)
如果已经达到平衡,那么应该有p(t+1)=p(t),所以如果p(x)是平稳分布,那么它满足
p(y)=xp(yx)p(x)
或者说,平稳分布就是上述关于p(x)的方程的非零解。

那么,有一个需要回答的问题是:上述方程的非零解唯一吗?这个问题的答案关系到马尔可夫链的平稳分布是否依赖于初始状态x0。很遗憾,答案是不一定。简单来说,非零解唯一的条件是该马尔可夫链的任意两个状态都是连通的,说具体点就是对于任意两个状态x,y,存在状态链x=x0,x1,,xn1,xn=y,使得
p(x1x0)>0p(x2x1)>0p(xnxn1)>0
说白了,就是从状态x跳转到y(不要求一步到位,可以多步)的概率大于0。其实这个判据也不难理解,形象比喻就是指平稳分布唯一的条件是马尔可夫链不能存在“孤岛”,要不然如果x0在孤岛内,那么就永远在孤岛内游走了,如果x0在孤岛外,那么也就永远在孤岛外游走了,从而不同的初始状态导致了不同的平稳分布。

细致平稳 #

由于马尔可夫链存在一个平稳分布,而构建马尔可夫链则只需要一个转移矩阵,如果我们可以构建一个易于采样的转移矩阵,并且它的平稳分布就是我们要采样的分布p(x),那么迭代足够多步后,不就相当于从p(x)中采样了?这就是所有MCMC方法的思想。当然,马尔可夫链往往需要迭代多步才能达到平稳,因此采样往往比较费时,但即便如此,终究能在有限成本内完成了采样过程,因此在很多情况下依然是实用的。

假设MCMC是可行的,那么很显然关键之处就是如何构造出一个转移概率p(yx),使得它的平稳分布是我们给定的分布p(x),同时,根据上一篇文章的介绍,我们往往只知道跟p(x)成正比的ρ(x),而不知道归一化因子,因此构建p(yx)的过程中,不能依赖于p(x)的归一化因子。

为此,我们需要用到一个“细致平稳条件”:

细致平稳条件 如果分布p(x)和转移概率p(yx)满足恒等式p(yx)p(x)=p(xy)p(y) 那么p(x)就是p(yx)的平稳分布。

名字看上去高端,实际上它非常容易证明,只要两端对x求和,马上就得到式(13)了,满足式(13)就说明p(x)p(yx)的平稳分布。“细致平稳条件”是判断平稳分布的一个充分条件,但不是必要的,它的价值在于提供了一种比较方便的为任意分布构建转移概率的方式。

转移概率 #

给定平稳分布p(x)以及任意的参考转移概率q(yx),如果它满足细致平稳条件,那么就万事大吉,直接用参考转移概率q(yx)作为最终的概率就好,如果它不满足,那么记α(xy)=q(yx)p(x),那么有下述显然成立的恒等式:
α(yx)q(yx)˜q(yx)p(x)=α(xy)q(xy)˜q(xy)p(y)
事实上,这个恒等式近乎废话,它相当于说“虽然ab,但ab=ba始终成立”。但就是这么个“废话”般的恒等式,给了我们启示:如果用˜q(yx)作为转移概率,那么对应的平稳分布不就是p(x)了?

其实没那么简单,因为这样构造出来的˜q(yx)通常是没有归一化的,而没有归一化即意味着它本身就不是一个合理的转移概率,也就没有什么平稳分布之说了。注意,如果要进行采样,我们或许可以不需要知道归一化因子,但是理论分布的时候,我们是需要保证归一化成立的。有的读者可能会想,那我手动除以个归一化因子行不?不行,因为除以归一化因子之后,细致平稳条件就未必成立了,甚至式(13)都可能不成立了。笔者查阅了很多介绍资料,都没有澄清这个问题,都是说直接用˜q(yx)作为转移概率的,困惑了很长时间。

事实上,MCMC方法根本就没有去归一化˜q(yx),而是将剩余部分的概率概率全部转移到自身去了。用数学的话来说,就是它最后使用的转移概率实际上是:
p(yx)=˜q(yx)+(1y˜q(yx))δ(yx)
其中δ(yx)={1,y=x0,yx,代表了只能转移到自身(即“永远不变”)的转移概率。所以,p(yx)的定义其实就是把多出来的概率统统叠加到“不变”这一操作上,所以它显然满足归一化,至于细致平稳条件可以代入检验,同样是满足的,其中主要是因为对于任意的p(x)f(x,y),下式总是恒成立的。
f(x,y)δ(yx)p(x)=f(y,x)δ(xy)p(y)
所以p(yx)才是我们所要寻求的真正的转移概率。如果状态空间有限的话,那么转移概率就对应一个有限的矩阵,上述结果其实就是说“通过调整转移矩阵的对角线元素可以完成归一化,并且不影响细致平稳条件的成立”。

MCMC方法 #

怎么实现从上述p(yx)从采样呢?很简单,先实现从˜q(yx)采样,然后如果不满足从˜q(yx)采样的条件,那么就保持不变;而˜q(yx)=α(yx)q(yx),我们假定从q(yx)中采样是容易的,那么回忆上一篇文章介绍的拒绝采样,我们可以通过进一步采样一个εU[0,1],通过α(yx)<ε来决定是否接受采样出来的y。综上所述,我们就得到了如下的MCMC采样流程:

Metropolis采样

初始状态为x0t时刻状态为xt

通过如下流程采样出xt+1
  1、采样yq(yxt)
  2、采样εU[0,1]
  3、计算α(yxt)=q(xty)p(y)
  4、如果εα(yxt),那么xt+1=y,否则xt+1=xt

这就是由Metropolis在1953年提出来的采样算法,一般称为Metropolis算法,或者干脆称为MCMC方法。它需要知道p(x)的精确表达式,并且需要知道一个易于采样的转移概率q(yx),然后根据上述流程迭代采样足够多步后,所得的xt我们就可以认为是从p(x)中采样出来的。这里的足够多步取决于p(x)q(yx)的具体选择,一般来说只能根据实验结果确定,没有统一的标准。

可能有读者觉得,这不就是一个普通的拒绝采样吗?刚才说的“把多出来的概率统统叠加到‘不变’这一操作上”体现在哪?事实上,跟普通的拒绝采样还是有点不一样的,在普通的拒绝采样中,如果ε>α(yx),那就重复1、2、3步采样过程,直到采样出εα(yx)y为止,然后让xt+1=y;而上述采样过程中,如果ε>α(yx),那么直接让xt+1=xt了,再下一步采样出来的就是xt+2了。假如我们认为从T时刻开始,分布就平稳了,那么{xt}t=T都是p(x)中采样出来的样本(同分布但不独立),显然两种不同的拒绝方式,会影响{xt}t=T最终的分布情况。

此外,值得指出的是,尽管本系列文章中,我们所处理的x是离散对象,但事实上同样的结论也适用于连续型对象的采样,只需要在前面的推导过程中,将求和改为对概率密度的积分就行了,没有实质的区别。

MH采样 #

上述Metropolis采样在很多场景下已经可用了,但它还不是特别完美。首先,刚才说了,Metropolis采样需要知道p(x)的精确表达式,这在不少任务下是难以做到的;其次,接受率α(yx)往往过小,导致长时间“原地踏步”,从而达到平稳状态的时间过长。幸运的是,有一个简单的技巧可以同时解决这两个问题。

该技巧就是使用如下式子作为接受率
A(yx)=min
该接受率其实也是源于恒等式\eqref{eq:feihua},对于接受率而言,我们唯一的要求就是它在0~1之间,此外是越接近于1越好,为此,我们可以让恒等式\eqref{eq:feihua}两边的接受率都除以\max(\alpha(\boldsymbol{y}\leftarrow\boldsymbol{x}),\alpha(\boldsymbol{x}\leftarrow\boldsymbol{y})),这样一来,它们俩其中一个就变成了1,另一个也不大于1,并且得到了放大,化简之后就得到了\mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x})。而且相当巧妙的是,此时\mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x})只依赖于p(\boldsymbol{y})的相对值,所以我们可以不用算出它的归一化因子了。

这样一来,我们就得到了改进版的Metropolis采样算法,这就是大名鼎鼎的Metropolis-Hastings Sampling(MH采样)

Metropolis-Hastings采样

初始状态为\boldsymbol{x}_0t时刻状态为\boldsymbol{x}_t

通过如下流程采样出\boldsymbol{x}_{t+1}
  1、采样\boldsymbol{y}\sim q(\boldsymbol{y}\leftarrow\boldsymbol{x}_t)
  2、采样\varepsilon\sim U[0,1]
  3、计算\mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}_t) = \min\left(1, \frac{q(\boldsymbol{x}_t\leftarrow\boldsymbol{y})p(\boldsymbol{y})}{q(\boldsymbol{y}\leftarrow\boldsymbol{x}_t)p(\boldsymbol{x}_t)}\right)
  4、如果\varepsilon \leq \mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}_t),那么\boldsymbol{x}_{t+1} = \boldsymbol{y},否则\boldsymbol{x}_{t+1}=\boldsymbol{x}_t

分析思考 #

上一篇文章中,我们指出直接从p(\boldsymbol{x})中采样是很困难的,就算是拒绝采样,如果\boldsymbol{x}的维度过高,那么从近似分布q(\boldsymbol{x})采样后其接受率通常也会很低,导致拒绝采样的效率极低而无法使用。那么很自然的一个问题是:

在MCMC方法中,同样要从转移概率q(\boldsymbol{y}\leftarrow \boldsymbol{x})中采样,同样有接受概率\mathcal{A}(\boldsymbol{y}\leftarrow \boldsymbol{x}),那为什么MCMC就是实用的,而普通的拒绝采样就是不实用的?

其实,对于直接拒绝采样来说,它要从近似分布q(\boldsymbol{x})中直接采样的是整个高维序列,如果q(\boldsymbol{x})与精确的p(\boldsymbol{x})差得比较大,那么接受概率往往是指数衰减的,因此接受概率极低;而对于MCMC方法来说,我们对q(\boldsymbol{y}\leftarrow \boldsymbol{x})的形式没有过多的限制,那么我们可以适当地设计q(\boldsymbol{y}\leftarrow \boldsymbol{x}),使得概率分布只集中在跟\boldsymbol{x}相似的那些\boldsymbol{y}中,换句话说,只有当\boldsymbol{y}\boldsymbol{x}很相似时,q(\boldsymbol{y}\leftarrow \boldsymbol{x})才不为0,否则就为0。这样一来,从q(\boldsymbol{y}\leftarrow \boldsymbol{x})采样的结果跟输入\boldsymbol{x}只有小的差别,变化不大的情况下,接受率往往更高些,所以使得拒绝采样成为可能。

所以,说白了,MCMC方法就是将“直接生成\boldsymbol{x}”,转换为了“\boldsymbol{x}_0出发,反复对它微调、润色,直到生成符合条件的\boldsymbol{x}”这么一个渐进式过程,从而使得它行之有效,这跟上一篇文章提到的“一步步来”的思想是一致的。

接下来,我们将介绍两个MH采样的例子:Gibbs采样和模拟退火,它们都体现了MCMC的这种微调、渐进的思想。

Gibbs采样 #

假设\boldsymbol{x}=(x_1,x_2,\dots,x_l)是一个长度为l的序列,Gibbs采样每次只调整其中一个元素。具体来说,Gibbs采样将参考转移概率定义为(严格来说还要再除以个l,但常数不改变结果,就不写出来了):
\begin{equation}q(\boldsymbol{x}_{[x_i=y]}\leftarrow \boldsymbol{x}) = p(y|\boldsymbol{x}_{-i})\triangleq\frac{p(x_1,\dots,x_{i-1},y,x_{i+1},\cdots,x_l)}{\sum\limits_y p(x_1,\dots,x_{i-1},y,x_{i+1},\cdots,x_l)}\end{equation}
其中\boldsymbol{x}_{[x_i=y]}是将\boldsymbol{x}的第i个位置替换为y后的序列。这也就是说,它根据待采样分布p(\boldsymbol{x})自身来构建p(y|\boldsymbol{x}_{-i})作为转移概率,每次先从1,2,\cdots,l中随机选择一个位置i,然后将位置i的元素替换成从分布p(y|\boldsymbol{x}_{-i})采样出来的结果。更妙的是,我们可以证明这种情况下接受概率恒为1:
\begin{equation}\begin{aligned} \mathcal{A}(\boldsymbol{x}_{[x_i=y]}\leftarrow \boldsymbol{x}) =&\, \min\left(1, \frac{p(x_i|\boldsymbol{x}_{-i})p(\boldsymbol{x}_{[x_i=y]})}{p(y|\boldsymbol{x}_{-i})p(\boldsymbol{x})}\right) \\ =&\, \min\left(1, \frac{p(x_i|\boldsymbol{x}_{-i})p(y|\boldsymbol{x}_{-i})p(\boldsymbol{x}_{-i})}{p(y|\boldsymbol{x}_{-i})p(x_i|\boldsymbol{x}_{-i})p(\boldsymbol{x}_{-i})}\right)\\ =&\, \min\left(1, 1\right)\\ =&\,1 \end{aligned}\end{equation}
由此,我们可以构建如下的Gibbs采样过程:

Gibbs采样

初始状态为\boldsymbol{x}_0=(x_{0,1},x_{0,2},\cdots,x_{0,l})t时刻状态为\boldsymbol{x}_t=(x_{t,1},x_{t,2},\cdots,x_{t,l})

通过如下流程采样出\boldsymbol{x}_{t+1}
  1、均匀地从1,2,\cdots,l中采样一个i
  2、计算p(y|\boldsymbol{x}_{t,-i})=\frac{p(x_{t,1},\dots,x_{t,i-1},y,x_{t,i+1},\cdots,x_{t,l})}{\sum\limits_y p(x_{t,1},\dots,x_{t,i-1},y,x_{t,i+1},\cdots,x_{t,l})}
  3、采样y\sim p(y|\boldsymbol{x}_{t,-i})
  4、\boldsymbol{x}_{t+1} = {\boldsymbol{x}_t}_{[x_{t,i}=y]}(即将\boldsymbol{x}_t的第i个位置替换为y作为\boldsymbol{x}_{t+1})。

模拟退火 #

另一个例子是模拟退火(Simulated Annealing)算法。前一篇文章中,我们提到了最大化一个目标,其实也可以看成是从这个目标中进行随机采样,其结果便对应着模拟退火算法。

首先,设我们要最大化的函数是f(\boldsymbol{x}),并且假设存在某个常数T > 0,使得对于任意0 < \tau \leq T,都有
\begin{equation}\sum_{\boldsymbol{x}} e^{f(\boldsymbol{x})/\tau} < \infty\end{equation}
这里的\tau的物理意义就是温度。这个假设看起来是额外的约束,但事实上模拟退火也就只适用于满足这个条件的场景。由于满足该条件,那么我们就可以构造分布
\begin{equation}p_{\tau}(\boldsymbol{x}) = \frac{e^{f(\boldsymbol{x})/\tau}}{\sum\limits_{\boldsymbol{x}} e^{f(\boldsymbol{x})/\tau}}\end{equation}
假设最大值点是唯一的,那么当\tau\to 0时,p_{\tau}(\boldsymbol{x})就会退化为one hot分布,也就是只有最大值点的概率为1,如果从p_{\tau}(\boldsymbol{x})中采样,那么采样结果就是最大值点了。即使\tau > 0,那么依然是最大值点的概率最大,所以从p_{\tau}(\boldsymbol{x})中采样,也会有趋于最大值点的趋势。所以,我们可以先固定温度\tau,构造一个从p_{\tau}(\boldsymbol{x})进行MH采样的随机过程,然后再慢慢让\tau下降,那么最终的收敛结果就是最大值点附近了,这就是所谓的“模拟退火”了。

为了进行MH采样,我们需要构建转移矩阵,模拟退火的选择比较简单,就是根据当前状态\boldsymbol{x}设计固定数目的\boldsymbol{y}作为新的\boldsymbol{x}的候选值(这步通常称为“变异”,往往只是对\boldsymbol{x}进行简单修改来得到\boldsymbol{y}),然后q(\boldsymbol{y}\leftarrow \boldsymbol{x})就是从候选中均允随机地选一个。既然是均匀的,那么q(\boldsymbol{y}\leftarrow \boldsymbol{x})就是一个常数,并且q(\boldsymbol{y}\leftarrow \boldsymbol{x})=q(\boldsymbol{x}\leftarrow \boldsymbol{y}),因此就有
\begin{equation}\mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}) = \min\left(1, \frac{q(\boldsymbol{x}\leftarrow\boldsymbol{y})p_{\tau}(\boldsymbol{y})}{q(\boldsymbol{y}\leftarrow\boldsymbol{x})p_{\tau}(\boldsymbol{x})}\right) = \min\left(1, e^{[f(\boldsymbol{y}) - f(\boldsymbol{x})]/\tau}\right)\end{equation}
所以,模拟退火就是一种搜索策略,如果f(\boldsymbol{y}) \geq f(\boldsymbol{x}),那么就更新\boldsymbol{x}=\boldsymbol{y},就算不是,那么还有一定的概率更新\boldsymbol{x}=\boldsymbol{y}。在整个搜索过程中,我们慢慢进行退火(\tau慢慢趋于0),在适当的退火策略下,模拟退火几乎总是能找到最大值点的。当然,如何根据具体问题制定“适当的”退火策略,那又是一个值得思考的问题了。

模拟退火

初始状态为\boldsymbol{x}_0,初始温度为\tau_0,温度按照事先制定的策略下降,t时刻状态为\boldsymbol{x}_t,温度为\tau_t

通过如下流程采样出\boldsymbol{x}_{t+1}
  1、采样\boldsymbol{y}\sim q(\boldsymbol{y}\leftarrow\boldsymbol{x}_t)
  2、采样\varepsilon\sim U[0,1]
  3、计算\mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}_t) = \min\left(1, e^{[f(\boldsymbol{y}) - f(\boldsymbol{x}_t)]/\tau_t}\right)
  4、如果\varepsilon\leq \mathcal{A}(\boldsymbol{y}\leftarrow\boldsymbol{x}_t),那么\boldsymbol{x}_{t+1} = \boldsymbol{y},否则\boldsymbol{x}_{t+1}=\boldsymbol{x}_t

如果简单修改一下,将接受策略改为“如果f(\boldsymbol{y}) \geq f(\boldsymbol{x}_t)则接受,否则就拒绝”,那么就变成了更朴素的“爬山法(Hill Climbing)”了。显然,爬山法目标更明确,前期可能收敛更快,但是一旦陷入了局部最大值后,通常就跳不出来了,最终的收敛结果没有模拟退火那么好。当然,也可以通过选择不同的初始值来重复实验,提升爬山法的效果。具体选择什么算法,那就看具体问题而定了。

本文小结 #

本文依然介绍“受限文本生成”所需要的基础采样算法——MCMC方法。经过反复修改,本文总算写出了笔者觉得比较满意的MCMC方法介绍,其中主要的区别是回答了一些常见的MCMC方法教程中没有提及到的细节问题。对于基础算法的介绍,到此就告一段落。从下一篇文章开始,我们将逐步介绍如何将这些看似枯燥的采样算法应用到生动、具体的文本生成任务中去。

转载到请包括本文地址:https://kexue.fm/archives/8084

更详细的转载事宜请参考:《科学空间FAQ》

如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。

如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!

如果您需要引用本文,请参考:

苏剑林. (Jan. 14, 2021). 《【搜出来的文本】⋅(二)从MCMC到模拟退火 》[Blog post]. Retrieved from https://kexue.fm/archives/8084

@online{kexuefm-8084,
        title={【搜出来的文本】⋅(二)从MCMC到模拟退火},
        author={苏剑林},
        year={2021},
        month={Jan},
        url={\url{https://kexue.fm/archives/8084}},
}