机器学习-白板推导系列(十三)- 马尔可夫链&蒙特卡洛方法(MCMC, Markov Chain&Monte Carlo)之Gibbs Sampling和平稳分布

13. 马尔可夫链&蒙特卡洛方法(MCMC)

13.4 Gibbs Sampling

13.4.1 概念

1. 思想
   假设有一随机向量$X=(x_1,x_2,...,x_d)$，其中$d$表示有$d$维，每一维是一随机变量，且并不是常见的相互独立前提。那么，如果已知这个随机向量的概率分布，如何从这个分布中进行采样呢？我们的思想就是一维一维的来，在对每一维进行采样的时候固定住其他的维度，这就是Gibbs Sampling。
2. 数学表示
   从一个随机的初始化状态$x^{(0)}=(x_1|x_2^{(0)},x_3^{(0)},\cdots ,x_d^{(0)})$开始，对每个维度单独进行采样，其采样顺序大致如下：
   $$\begin{gathered} x_1^{(1)}\sim p(x_1|x_2^{(0)},x_3^{(0)},\cdots ,x_d^{(0)}) \\ {x}_2^{(1)}\sim p(x_2|x_1^{(0)},x_3^{(0)},\cdots ,x_d^{(0)}) \\ ⋮ \\ {x}_d^{(1)}\sim p(x_d|x_1^{(0)},x_2^{(0)},\cdots ,x_{d-1}^{(0)}) \\ ⋮ \\ {x}_1^{(t)}\sim p(x_1|x_2^{(t-1)},x_3^{(t-1)},\cdots ,x_d^{(t-1)}) \\ ⋮ \\ {x}_d^{(t)}\sim p(x_d|x_1^{(t-1)},x_2^{(t-1)},\cdots ,x_{d-1}^{(t-1)})\text{\hspace{1em}(13.4.1)} \end{gathered}$$
   遵从上面的采样步骤，我们最终能够采样得到所需要的高维分布的样本。Gibbs Sampling的过程更像是一个单步迭代的过程，实际上Gibbs是一种特殊的MH采样，为什么呢？我们来证明一下。

   迭代的最开始采样得到的样本并不是完全满足所需要的分布的样本，因为$采样之初采样的分布是提议分布，一般是均匀分布$。

3. 迭代过程
   
   上图所示，右图是我们需要的分布，左边是迭代的过程，最开始抽样的点0和1都是均匀分布抽样得到的，而越到后面，抽样的点都越满足我们右边的分布，所以这个过程可以说明Gibbs Sampling抽样的过程是可行的。

13.4.2 Gibbs是特殊的MH采样

• 我们首先回顾一下，MH采样的方法。假设有一随机向量$X=(x_1,x_2,...,x_d)$，其中$d$表示有$d$维，$t$表示时间序列中的某一时刻。我们的目的是从${\mathbf{Q}}_{X^{(t)},X}(X^{(t)}\mapsto X)$中采样获得$X^{(next)}$，然后计算接受率
  $$\alpha (X^{(t)},X^{(next)})=min(\frac{P(X^{(next)}){\mathbf{Q}}_{next,t}}{P(X^{(t)}){\mathbf{Q}}_{t,next}},1)\text{\hspace{1em}(13.4.2)}$$
  1. 首先看$P_{(next)}(X)$，其可以写成：
     $$P(X^{(next)})=P(X_i^{(next)}|X_{-i}^{(next)})P(X_{-i}^{(next)})\text{\hspace{1em}(13.4.3)}$$
     其中$i$表示对第$i$维进行采样，$-i$表示除了第$-i$剩下的维度。
  2. 接着看${\mathbf{Q}}_{next,t}$($X^{(next)}\mapsto X^{(t)}$）：
     $${\mathbf{Q}}_{next,t}=P(X^{(t)}|X^{(next)})\text{\hspace{1em}(13.4.4)}$$
     则公式(13.4.2)可以写成：
     $$\begin{array}{ll}\alpha (X^{(t)},X^{(next)})& =min(\frac{P(X^{(next)}){\mathbf{Q}}_{next,t}}{P(X^{(t)}){\mathbf{Q}}_{t,next}},1)\\ & =min(\frac{P(X_i^{(next)}|X_{-i}^{(next)})P(X_{-i}^{(next)})P(X^{(t)}|X^{(next)})}{P(X_i^{(t)}|X_{-i}^{(t)})P(X_{-i}^{(t)})P(X^{(next)}|X^{(t)})},1)\end{array}\text{\hspace{1em}(13.4.5)}$$
  3. 然后优化${\mathbf{Q}}_{next,t}$
     • 根据Gibbs Sampling的方法，由公式(13.4.1)可知，每次迭代是固定$-i$,迭代$i$，那么对于${\mathbf{Q}}_{next,t}(P(X^{(t)}|X^{(next)}))$可以写成：
       $$\begin{gathered} {\mathbf{Q}}_{next,t}=P(X^{(t)}|X^{(next)})=P(X_i^{(t)}|X_{-i}^{(next)}) \\ {\mathbf{Q}}_{t,next}=P(X^{(next)}|X^{(t)})=P(X_i^{(next)}|X_{-i}^{(t)})\text{\hspace{1em}(13.4.6)} \end{gathered}$$
     • 每次迭代，$X_{-i}^{(t)}$和$X_{-i}^{(next)}$是一样的，即：
       $$P(X_{-i}^{(next)})=P(X_{-i}^{(t)})\text{\hspace{1em}(13.4.7)}$$
       因此${\mathbf{Q}}_{next,t}$还能写成：
       $$\begin{array}{ll}{\mathbf{Q}}_{next,t}& =P(X^{(t)}|X^{(next)})\\ & =P(X_i^{(t)}|X_{-i}^{(next)})\\ & =P(X_i^{(t)}|X_{-i}^{(t)})\\ {\mathbf{Q}}_{t,next}& =P(X_i^{(next)}|X_{-i}^{(next)})\end{array}\text{\hspace{1em}(13.4.8)}$$
  4. 整理
     因此公式(13.4.5)~(13.4.8)可知：
     $$\begin{array}{ll}\alpha (X^{(t)},X^{(next)})& =min(\frac{P(X^{(next)}){\mathbf{Q}}_{next,t}}{P(X^{(t)}){\mathbf{Q}}_{t,next}},1)\\ & =min(\frac{P(X_i^{(next)}|X_{-i}^{(next)})P(X_{-i}^{(next)})P(X^{(t)}|X^{(next)})}{P(X_i^{(t)}|X_{-i}^{(t)})P(X_{-i}^{(t)}))P(X^{(next)}|X^{(t)})},1)\\ & =min(\frac{P(X_i^{(next)}|X_{-i}^{(next)})P(X_{-i}^{(t)})P(X_i^{(t)}|X_{-i}^{(t)})}{P(X_i^{(t)}|X_{-i}^{(t)})P(X_{-i}^{(t)})P(X_i^{(next)}|X_{-i}^{(next)})},1)\\ & =1\end{array}\text{\hspace{1em}(13.4.9)}$$
     因此Gibbs Samplings是$\alpha =1$的MH Sampling的意义了。Gibbs Sampling就是把目标分布$P$对应的条件概率当作状态转移分布$Q$。
• 使用Gibbs Sampling是有使用前提的，即：固定其他维度后的一维分布时方便进行采样的，如果固定其他维度的时候得到的一维分布仍然是非常难进行采样的，就很难使用Gibbs Sampling。

---

13.5 回顾

这一小节将主要来介绍：什么是采样？为什么而采样？什么样的样本是好的样本？以及采样中主要遇到的困难？

1. 采样的动机
   1. 完成任务
      我们机器学习中经常需要进行采样来完成各种各样的任务。如：从一个$P(X)$中采出一堆样本。
   2. 求和求积分
      包括大名鼎鼎的Monte Carlo算法。求$P(X)$主要是为了求在$P(X)$概率分布下的一个相关函数的期望，也就是：
      $$\int P(x)f(x)dx={\mathbb{E}}_{P(X)}[f(X)]\approx \frac{1}{N}\sum _{i=1}^{N}f(x^{(i)})\text{\hspace{1em}(13.5.1)}$$
      通过采样来得到 P ( X ) ∼ { x ( 1 ) , x ( 2 ) , ⋯   , x ( N ) } P(X) \sim \{ x^{(1)},x^{(2)},\cdots, x^{(N)} \} P(X)∼{x(1),x(2),⋯,x(N)}样本点。
2. 什么样的是好样本
   什么样的样本就是好样本呢？或者说是采样的效率更高一些。
   1. 样本的分布要趋向于原始的目标分布
      也就是说样本要趋向于高概率选择区域。或者是说，采出来的样本出现的概率和实际的目标分布的概率保持一致。
   2. 样本和样本之间是相互独立的
      如果采出来的一堆样本之间都差不多，那么就算采出来了趋向于高概率选择区域的样本，那采样效率太低了，样本中反映的信息量太有限了。
3. 实际采样中的困难
   1. $\text{Partation function is intractable.}$
      后验分布往往被写成$P(X)=\frac{1}{Z}\hat{P}(X)$，上面这个$\hat{P}(X)$都比较好求，就是等于 Likelihood $\times$ Prior。而$Z$是归一化常数，它非常的难以计算，$Z=\int \hat{P}(X)dX$，这几乎就是不可计算的。所以，有很多采样方法就是想要跳过求$P(X)$的过程，来从一个近似的分布中进行采样，当然这个近似的分布采样要比原分布简单。比如：Rejection Sampling和Importance Sampling。
   2. $\text{The curse of high dimension}.$
      如果样本空间$\mathcal{X}\in {\mathbb{R}}^p$，每个维度都有$K$维。那么总的样本空间就有$K^p$的状态。要知道那个状态的概率高，就必须要遍历整个样本空间，不然就不知道哪个样本的概率高，如果状态的数量是这样指数型增长的话，全看一遍之后进行采样时不可能的。所以，直接采样的方法是不可行的。
4. 采样方法
   • 借助$Q(x)$逼近目标分布$P(x)$
     Rejection Sampling和Importance Sampling，都是借助$Q(x)$逼近目标分布$P(x)$，通过从$Q(x)$中进行采样来达到在$P(x)$中采样的目的，而且在$Q(x)$中采样比较简单。如果$Q(x)$和$P(x)$直接的差距太大的话，采样效率会变得很低。
   • 利用马氏链
     MCMC方法，我们主要介绍了MH Sampling和Gibbs Sampling，主要通过构建一个马氏链去逼近目标分布。

---

13.6 平稳分布

MCMC采样借助马氏链，经过若干步以后会收敛到一个平稳分布。马尔可夫链的组成可以大致分成两个部分：

1. 状态空间： { 1 , 2 , 3 , ⋯   , k } \{ 1,2,3,\cdots,k \} {1,2,3,⋯,k}；
2. 状态转移空间$Q=[Q_{ij}{]}_{k\times k}$。

13.6.1 基本概念

马尔可夫链的模型可以被我们表达为：

• 每一个时间点有一个状态分布，$q^{(t)}(x)$表示当前时间点位于某个状态的概率分布情况。假设在$t=1$的时间节点有$K$个状态，状态的概率分布为$q^{(1)}(x)$表示为：

  $x$	$1$	$2$	…	$K$
  $q^{(1)}(x)$	$q^{(1)}(1)$	$q^{(1)}(2)$	…	$q^{(1)}(K)$

• 相邻时间节点之间的状态转移矩阵为：
  $$\begin{array}{c}Q={\left[\begin{array}{cccc}{Q}_{11}& Q_{12}& \cdots & Q_{1k}\\ Q_{21}& Q_{22}& \cdots & Q_{2k}\\ ⋮& ⋮& \ddots & ⋮\\ Q_{k1}& Q_{k2}& \cdots & Q_{kk}\end{array}\right]}_{k\times k}\end{array}\text{\hspace{1em}(13.6.1)}$$
  • 状态转移矩阵描述的是:$Q_{ij}=Q(X^{(t+1)}=j|X^{(t)}=i)$。描述的是从一个状态$i$转移到另外一个状态$j$的概率。
  • 状态转移矩阵的每一行$i$表示为目前状态是$i$时，到其他状态的概率，那么必然有${\sum }_{k=1}^k{Q}_{ik}=1$。
• 假设在$t=m$时刻之后到达了平稳分布状态，那么就可以得到：
  $$q^{(m)}(X)=q^{(m+1)}(X)=q^{(m+2)}(X)\text{\hspace{1em}(13.6.2)}$$

13.6.2 Markov Chain收敛性

1. Markov Chain状态转移表示
   假设在$t+1$时刻，状态是$x=j$。那么$q^{(t+1)}(x=j)$的分布为：$所有可能转移到这个状态的概率i$ 乘以$这个状态的分布q^{(t)}(x=i)$，我们用公式表达就是：
   $$q^{(t+1)}(x=j)=\sum _{i=1}^k{q}^{(t)}(x=i)Q_{ij}\text{\hspace{1em}(13.6.3)}$$
   • 当$x=j$时概率，即：在$t+1$时刻，可能出现的状态有$K$个。那么$q^{t+1}(X)$的分布可以表示为：
     $$\begin{array}{l}{q}^{(t+1)}(X)={\left[\begin{array}{ccccc}{q}^{(t+1)}(x=1)& q^{(t+1)}(x=2)& q^{(t+1)}(x=3)& \cdots & q^{(t+1)}(x=k)\end{array}\right]}_{1\times k}\end{array}\text{\hspace{1em}(13.6.4)}$$
     而，
     $$q^{(t+1)}(x=j)=\sum _{i=1}^K{q}^{(t)}(x=i)Q_{ij}\text{\hspace{1em}(13.6.5)}$$
   2. $q^{(t+1)}$可以被我们表示为：
      $$\begin{array}{ll}{q}^{(t+1)}(X)& ={\left[\begin{array}{cccc}{\sum }_{i=1}^k{q}^{(t)}(x=i)Q_{i1}& {\sum }_{i=1}^k{q}^{(t)}(x=i)Q_{i2}& \cdots & {\sum }_{i=1}^k{q}^{(t)}(x=i)Q_{ik}\end{array}\right]}_{1\times k}\\ & =q^{(t)}(X)\cdot Q\end{array}\text{\hspace{1em}(13.6.6)}$$
      其中，$q^{(t)}(X)={\left[\begin{array}{ccccc}{q}^{(t)}(x=1)& q^{(t)}(x=2)& q^{(t)}(x=3)& \cdots & q^{(t)}(x=k)\end{array}\right]}_{1\times k}$。
   3. 通过这个递推公式，我们可以得到：
      $$q^{(t+1)}(X)=q^{(t)}(X)Q=q^{(t-1)}(X)Q^2=\cdots =q^{(1)}(X)Q^t\text{\hspace{1em}(13.6.7)}$$
      通过上述的描述，已知：在Markov Chain中，每个时刻点的状态的分布$q^{(t)}(X)$的计算方法。
2. Markov Chain 收敛性
   由于$Q$是一个随机概率矩阵，那么我们可以得到，$Q每个值都是小于或等于1$，所以也必然有特征值的绝对值$\le 1$。所以，我们可以对概率转移矩阵做特征值分解，分解成对角矩阵：
   $$\begin{array}{l}Q=A\Lambda A^{-1}\Lambda =\left[\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_k\end{array}\right],|{\lambda }_i|\le 1(i=1,2,\cdots ,k)\end{array}\text{\hspace{1em}(13.6.8)}$$

   为什么特征值的绝对值$\le 1$？我们可以从特征值的几何意义上看，特征值代表变换中方向不变的向量的变化尺度。随机矩阵的变化尺度必然是小于1的。

   1. 假设只有一个${\lambda }_i=1$，则：
      $$q^{(t+1)}(X)=q^{(1)}(X)(A\Lambda A^{-1}{)}^t=q^{(1)}(X)A{\Lambda }^t{A}^{-1}\text{\hspace{1em}(13.6.9)}$$
      当$t\to \infty$时，必然有：
      $$\begin{array}{l}{\Lambda }^t=\left[\begin{array}{cccc}0& & & \\ & 1& & \\ & & \ddots & \\ & & & 0\end{array}\right]\end{array}\text{\hspace{1em}(13.6.10)}$$
      当$M$足够大，依然有：
      $$\begin{array}{l}{\Lambda }^M=\left[\begin{array}{cccc}0& & & \\ & 1& & \\ & & \ddots & \\ & & & 0\end{array}\right]\end{array}\text{\hspace{1em}(13.6.11)}$$
      所以，
      $$\begin{array}{ll}{q}^{(m+1)}& =q^{(1)}A{\Lambda }^m{A}^{-1}\\ q^{(m+2)}(X)& =q^{(m+1)}(X)A\Lambda A^{-1}\\ & =q^{(1)}(X)A{\Lambda }^m{A}^{-1}A\Lambda A^{-1}\\ & =q^{(1)}(X)A{\Lambda }^{(m+1)}{A}^{-1}\\ & =q^{(m+1)}(X)\end{array}\text{\hspace{1em}(13.6.12)}$$
   2. 当特征值小于1时，依然可证明：综上所述：
      $$当t>m，q^{(m+1)}(X)=q^{(m+2)}(X)=\cdots =q^{(\infty )}(X)\text{\hspace{1em}(13.6.13)}$$
      这就是平稳分布，当Markov Chain经过足够大的步数$m$之后，一定会收敛到一个平稳分布。
3. 采样
   通过平稳分布的性质，启发我们设计一个Markov Chain，收敛到我们想要采样的分布$p(x)$。那么。怎么才能让它收敛呢？实际上就是由状态转移矩阵$Q$所决定的。我们的核心问题就是设计一个合适的状态转移矩阵$Q$。
   即我们要做的就是：
   $设计一个MCMC，利用MarkovChain收敛到一个平稳分布q(x)，使得平稳分布\approx 目标分布p(x)$。当$m$足够大时，
   $$q^{(m)}(X)=q^{(m+1)}(X)=q^{(m+2)}(X)=\cdots =q(X)$$

   • 当Markov Chain解决了当维度很高的时候，$q(x)\approx p(x)$找不到的情况，在MCMC中不要显示的去找，而是构建一个Markov Chain去近似，跳过了直接去寻找的过程。
   • 定义： $从开始到收敛的这段时期称为bum-in$，也称这个时间$t$为$Mix-time$。

---

13.7 MCMC采样时遇到的困难

1. bum-in不确定
   虽然可以证明出MCMC最终可以收敛到一个平稳分布。但是并没有理论来判断Markov Chain是否进入了平稳分布，也就是不知道Markov Chain什么时候会收敛。

2. Mix-Time过长
   这就是有高维所造成的，维度和维度之间的相关性太强了，$p(x)$太过于复杂了。理论上MCMC是可以收敛的，但是如果$m$如果实在是太大的话，基本就是认为它是不收敛的。实际上，现在有各种各样的MCMC算法的变种都是在想办法解决这个Mix-Time过长的问题。

3. 独立性不确定
   我们希望采到的样本之间的样本相互独立，即：采到的样本之间的相关性越小越好。

   实际在高维分布中我们采用MCMC来进行采样很有可能造成样本单一，相关性太强的问题。举一个Mixture Gaussian Distribution的例子。下图所示是一个Mixture Gaussian Distribution的例子：
   
   用MCMC采样，样本都趋向于一个峰值附近，很有可能会过不了低谷，导致样本都聚集在一个峰值附近。

   这个问题出现的原因我们可以从能量的角度来解释这个问题。在无向图中，我们常用下列公式来进行表示：
   $$P(X)=\frac{1}{Z}\hat{P}(X)=\frac{1}{Z}ex{p}^{-\mathbb{E}(X)}\text{\hspace{1em}(13.7.1)}$$

   1. 实际上这里的$\mathbb{E}(X)$指的就是$能量函数$，能量和概率是成反比的，概率越大意味着能量越低，能量越低，越难发生跳跃的现象。所以，采样很容易陷入到一个峰值附近。
   2. 并且，多峰还可以分为均匀和陡峭，陡峭的情况中，能量差实在是太大了，就很难发生跳跃。就像孙悟空翻出如来佛祖的五指山一样，佛祖的维度很好，孙悟空在翻跟头的时候，一直在一个低维里面不同的打转，根本就跳不出来，就是来自佛祖的降维打击。

   所以，在高维情况下，很容易发生在一个峰值附近不停的采样，根本就跳不出来，导致采到的样本的多样性低，样本之间的关联性大，独立性低。