Markov 链

Markov 链的定义

对于一个状态集合 \(S\)（\(S\) 有可能有限、无限但可数，其实也可以存在连续的 Markov 链，以下我们暂时不讨论无限且不可数的情况），概率分布向量 \(\mathbf p(t)\) 给出了 \(t\) 时刻处在某个状态的概率。（本文中的分布向量均为行向量）

如果存在概率状态转移矩阵 \(P\) 使得 \(\forall t\ge 0,\mathbf p(t)P=\mathbf p(t+1)\)，也就是每个时刻状态概率分布只和前一个时刻状态概率分布有关，则称该随机过程是一个 Markov 链。

本文的内容主要研究 Markov 链迭代无穷次或充分多次后的行为，并且浅谈其在 Monte Carlo 方法中的应用。

Markov 链的常返性

瞬态与常返态

对于 Markov 链状态空间无限的情况，有时它的行为不一定是非常良好的，可能从一个状态出发，走着走着，就再也不可能回到原来的状态了。为了刻画是否能返回某个状态，我们引入瞬态与常返态的概念。

对于某个 Markov 链，定义 \(f_i\) 表示从状态 \(i\) 出发至少一次返回 \(i\) 的概率。如果 \(f_i<1\)，则称状态 \(i\) 为瞬态 (transient)，否则为常返态 (recurrent)。

对于某个状态 \(i\)，再定义 \(N_i\) 表示随机游走过程中访问该状态的次数。那么 \(\mathbf E(N_i)=\frac{1}{1-f_i}\)，证明非常类似于计算几何分布的期望。这样我们同样可以得知如果 \(N_i\) 的期望收敛，那么状态 \(i\) 是瞬态。

网格对称随机游走

在一个 \(d\) 维的无限大的网格上随机游走，每一次等概率地选择一个相邻点（对于每个维度随机 \(\pm1\)），容易发现这是一个 Markov 过程。

可以发现，网格对称随机游走的每一个状态都是等价的。考虑计算任意一个 \(N_i\) 的期望，把这个过程看做 \(d\) 个一维网格对称随机游走，那么时刻 \(2t\) 回到原点的概率是 \(2^{-2td}{2t \choose t}^d\)

由斯特林公式，\(N_i\) 的期望可以看成：

\[\sum_{t=1}^{+\infty}\left\{ \frac{4^t}{\sqrt{\pi t}} \times [1+O(t^{-1})]\times 4^{-t}\right\}^d=\sum_{t=1}^{\infty} (\pi t)^{-\frac{d}{2}}\times [1+O(t^{-1})] \]

由于调和级数发散，而 \(\forall e>1,\sum_{t=1}^{\infty} t^{-e}\) 都收敛，所以我们得到了一个神奇的结论，对于 \(d\le 2\) 的网格对称随机游走，所有状态都是常返的，而 \(d>2\) 的高维网格对称随机游走所有状态都是瞬态。

正常返与零常返

即便如此，我们发现对于低维网格对称随机游走，它的性质还是非常奇怪。我们尽管知道期望上它会访问无数次的起点，然而这并不代表它返回原点的时间期望 \(\mathbf E(T_i)\) 有限（在无限的时间内，以期望无穷的时间间隔发生无穷次事件 QwQ）。

我们把常返但是期望返回时间 \(\mathbf E(T_i)\) 无穷的状态叫零常返态，反之叫做正常返态。注意到 \(\sum_{t=1}^{\infty} t(\pi t)^{-\frac{d}{2}}\) 在 \(d\le 2\) 时发散，我们知道了低维网格对称随机游走的所有状态均是零常返的。

Markov 链的可约性

不可约链

一个衡量 Markov 链的重要性质是不可约性，或者通俗点说就是连通性。如果对于任意一对状态 \(x,y\)，存在一个时间间隔 \(t\)，使得经过 \(t\) 个单位的时间后 \(x\) 到达 \(y\) 的概率为正，则称该链联通，即该链是不可约链。

链的可约性

那么为什么要起不可约链这个名字呢？对于不可约链的反面，任意可约链，状态集合 \(S\) 都可以分解成 \(T\cup C_1\cup C_2\cup C_3\dots\) 这种形式，其中 \(T\) 是瞬态集，该集合元素均为瞬态；\(C_i\) 都是不可约闭集，闭集即不存在从该集合指向外部的正概率转移。因此，我们研究任意 Markov 链在无穷次迭代之后的行为，都可以将其分解成不可约链再做研究。

一些定理

不考虑以下定理的证明。

不可约链所有状态的正常返性、零常返性、瞬态性都一致。

有限状态不可约链一定正常返。

所以常返性仅仅是对可数无穷状态的 Markov 链定义出来的概念。

Markov 链的收敛行为

由第一节的讨论，我们本节中可以只关心不可约链的情况。

平稳分布

对于概率分布向量 \(\pi\) 满足 \(\pi P=\pi\)，则称 \(\pi\) 为 Markov 链的平稳分布向量。可以发现，如果状态掉进了平稳分布中，那么无论经过多少次迭代概率分布都不再改变了。

我们有如下若干个定理，探讨常返性、不可约性和平稳分布的关系。

不可约链存在平稳分布等价于其正常返。

不可约链如果存在平稳分布，则平稳分布唯一。

正常返不可约链的平稳分布可以如此构造：\(\pi_i=\frac{1}{\mathbf E(T_i)}\)，\(T_i\) 为返回时间。

证明不考虑了。

当然，OIer 都熟稔的另一个计算平稳分布的方法是直接高斯消元。

细致平衡条件

判定一个向量 \(\pi\) 是不是平稳分布取决于平稳分布方程 \(\pi P=\pi\)。我们有一个比这个条件更强的判定条件 \(\forall i\ne j,\pi_i P_{i,j}=\pi_j P_{j,i}\)，被称为细致平衡条件。

并不是所有的平稳分布都满足细致平衡条件，然而在后文中构造指定平稳分布的 Markov 链时，细致平衡条件形式更简单，更容易检查。

更加重要的是，所有满足细致平衡条件的有限状态 Markov 链都可以看作带权无向图上的随机游走。对于一条边 \((i,j)\)，定义其权值 \(w_{i,j}=\pi_i P_{i,j}=\pi_j P_{j,i}=w_{j,i}\)。那么该 Markov 链可以看作这样的随机游走过程，对于所有与当前节点相邻的边，按照正比于边权的概率选择一条无向边转移出去。

因此，我们可以按照是否满足细致平衡条件对 Markov 链分类，一类叫对称链，一类叫非对称链。

对称链另一个描述方式是，它是“时间可逆”的。

无穷次迭代后的收敛判定

如果无穷次迭代之后 \(\mathbf p(t)\) 最终趋近于某个分布 \(\pi\)，那么称 Markov 链收敛。容易证明 \(\pi\) 一定是平稳分布。

欸，那么这是不是就意味着正常返不可约链就一定会收敛到它那个唯一存在的平稳分布上呢？

答案是否定的，还存在一种情况，即周期性 Markov 链。在一定的迭代次数过后，Markov 链近似地围绕一个周期打转。

不可约 Markov 链收敛的条件是正常返且非周期。

前缀平均分布的收敛性

周期性太烦了，如何刻画不可约正常返周期 Markov 链的性质呢？

虽然迭代分布 \(\mathbf p(t)\) 本身不收敛于平稳分布，但是其前缀平均值 \(\mathbf a(t)=\frac{1}{t}\sum_{i=0}^{t-1} \mathbf p(i)\) 一定收敛于平稳分布。也就是说，周期 Markov 链随着一遍遍迭代越来越接近一个概率分布的周期，这个周期极限的平均值同样是平稳分布。

\(\mathbf b(t)=\mathbf a(t)P-\mathbf a(t)=\frac{1}{t}\left[\sum_{i=0}^{t-1} \mathbf p(i+1)-\sum_{i=0}^{t-1} \mathbf p(i)\right]=\frac{1}{t}[\mathbf p(t)-\mathbf p(0)]\)。由于 \(|\mathbf p(t)-\mathbf p(0)|\le 2\)，因此当 \(t\to \infty\) 时，\(\mathbf a(t)P\) 和 \(\mathbf a(t)\) 的距离也趋近于零。这昭示着 \(\mathbf a(t)\) 的极限同样满足平稳方程 \(\pi P-\pi=0\)。

如何证明 \(\mathbf a(t)\) 收敛呢？这里贴一下课堂上的证明。

只考虑状态集合大小有限（设为 \(n\)）的简单情况，对于矩阵 \(P-I\)，在右边添加一行全 \(1\) 行向量形成一个 \(n\times (n+1)\) 的矩阵 \(B\)。

引理：矩阵 \(B\) 秩为 \(n\)。

具体地，该矩阵核中具有向量 \((1,1,\dots,1,0)\)。我们要证明该矩阵的核是一维的，这样由于秩-零化度定理，其像空间维数，即秩就是 \(n\)。

若该矩阵核维度更高，那么其中应该存在一个垂直于 \((1,1,\dots,1,0)\) 的非零向量 \((\mathbf x,\alpha)\)，满足 \(\mathbf x^T \mathbf 1=0,(P-I)\mathbf x+\alpha=0\)，即 \(\sum_i x_i=0,x_i=\alpha+\sum_j P_{i,j}x_j\)。

对于某个 \(i\),\(\forall P_{i,j}>0,x_i\le x_j\)，则称状态 \(i\) 是局部最小值，同理定义局部最大值。

由不可约性，一定存在某个局部最小值 \(x_i\)，满足 \(\exists j,P_{i,j}>0,x_i<x_j\)。否则的话考虑全局最小值，它一定是局部最小值，而它的邻接点的 \(x\) 值又全部等于它，依据连通性推导下去 \(\mathbf x\) 的所有分量相等，与 \(\mathbf x\) 的非零性矛盾。因此，对于该局部最小值 \(i\)，\(\alpha=x_i-\sum_j P_{i,j} x_j=\sum_j P_{i,j} (x_i-x_j)<0\)。

同理，对于局部最大值，\(\exists j,P_{i,j}>0,x_i>x_j\)，\(\alpha=x_i-\sum_j P_{i,j} x_j=\sum_j P_{i,j} (x_i-x_j)>0\)。

\(\alpha\) 既正又负矛盾！引理得证。

发现 \(\mathbf a(t)B=(\mathbf b(t),1)\)，由于 \(B\) 的秩为 \(n\)，所以像空间中的点与原像一一对应，存在线性映射 \(B'\) 把像空间中的点映射回原像，因此 \(\mathbf a(t)=B'(\mathbf b(t),1)\)。当 \(t\to \infty\) 的时候，\(\mathbf a(t)\) 就收敛于 \(B'(0,0,\dots,0,1)\)。

Markov 链与 Monte Carlo 方法

Ising 模型

学术之道学科史的展示正好涉及过 Ising 模型的介绍，结果马上就在数学课上碰到了，好巧。

这是一个统计物理学中描述磁场的模型。在规则的晶格（在图论中可以理解成一维、二维、三维的网格）上，每一个粒子 \(i\) 有两种自旋状态 \(s_i\in \{+1,-1\}\)。粒子只和在晶格中相邻的粒子发生相互作用。定义这个系统的总能量（哈密顿量）为：

\[H=-J\sum_{(u,v)\in E} s_us_v-\mu B\sum_{x\in V} s_x \]

（有些看不懂的常数就不用管了，我也不懂 QAQ）

在统计物理学中，还有一个概念：Gibbs 分布。意思是对于一个系统，假设有 \(|\Omega|\) 个微观态，能量为 \(E_1,E_2,\dots,E_{|\Omega|}\)，那么该系统处于某个微观态 \(i\) 的概率与 \(\exp(-\frac{E_i}{kT})\) 成正比（这个 \(T\) 是温度）。

将 Gibbs 分布应用到 Ising 模型上，我们发现当温度越高的时候系统处在每个状态的概率越均匀，整个系统表现出顺磁相；反之的话系统将会以很大概率处在能量低的状态，此时系统就会呈现铁磁相。

现在问题来了，我们想要对 Ising 模型进行采样，即生成一组样本，使其服从 Gibbs 分布，这样我们就可以数值计算很多关于该物理模型的问题了。

然而想要完全了解系统的概率分布，必须要遍历 \(2^{|V|}\) 个状态计算出整个分布，这样太慢了！如果只是想采采样观察统计规律，有没有什么更好的数值方法呢？这时候我们就需要我们的 MCMC 方法出场了。

Markov 链 Monte Carlo 方法（MCMC 方法）

在 Ising 模型的问题中，Gibbs 分布的归一化常数 \(Z=\frac{1}{\sum_{i=1}^{|\Omega|} \exp(-\frac{E_i}{kT})}\) 在没有枚举整个样本空间的情况下难以计算。我们又想在不枚举整个样本空间的情况下进行采样。

泛化这个问题，给定一个归一化常数未知的复杂分布 \(f(x)\)（假设这个分布定义在一个有限或者可数的集合上），我们要对这个分布生成一组服从该分布的采样点。

对于归一化常数未知的分布采样，我们已经学习过拒绝采样，即找到一个简单提议分布 \(g(x)\)，和一个常数 \(M\) 满足 \(\forall x,f(x)\le g(x)M\)，我们生成服从提议分布的样本点，以 \(\frac{f(x_i)}{g(x_i)M}\) 的概率接受一个样本点 \(x_i\)。

然而我们在对较为复杂的分布采样时，比如说像 Ising 模型这样定义在有限的高维空间上的分布，这种情况下合适的提议分布比较难找，简单的构造提议分布（如均匀分布）由于样本空间太大，导致 \(M\) 过大，接受率过低，效率低下。

这时候我们有一个想法，构造一个 Markov 链，使得其平稳分布正好是目标分布，我们就可以随机选取起始状态，每个时刻维护当前状态，模拟真正的概率转移过程，最后迭代足够多次之后，如果 Markov 链成功收敛于平稳分布，就可以得到一个服从目标分布的样本点。（注意这里并不是计算每个时刻的整个分布向量，而是只关注单个状态的概率转移。）

为了避免周期性导致的不收敛情况，我们可以计算所有迭代出来的样本点的平均值作为最终的样本。这样我们只需要保证 Markov 链不可约且正常返即可。事实上，像 Ising 模型这样的问题样本空间有限，所以只要不可约，正常返就可以自然满足。

当然，由于初始值随机，其可能对于平均值的影响非常大，如果将早期迭代的结果也纳入平均值会导致平均值的收敛变慢，所以好的方法是先迭代一会（burn-in 期）然后再将一个后缀的平均值当做最终结果。

关于如何具体设计合适的 Markov 链呢？这个 Markov 链不仅仅需要满足平稳分布是给定分布，还要使得这个 Markov 链联通、正常返。更重要的是，需要使得这个 Markov 链的状态转移容易计算，毕竟我们往往连遍历整个状态空间都做不到，只能依赖随机生成器模拟转移。

课堂上介绍了两种方法，主要思想都是构造满足细致平衡条件的 Markov 链。

Metropolis Hastings 算法

对于每一个状态 \(x\) 构造一个容易生成的提议分布 \(q_x(y)\)（实际上老师上课讲的版本直接就是构造某张连通图上的均匀分布，这里还是指出任意好生成的提议分布都是可以的）。

设 \(\alpha_{x,y}=\min\left(1,\frac{f(y)q_y(x)}{f(x)q_x(y)}\right)\)，构造 Markov 链：

\[P_{x,y}= \begin{cases} q_x(y)\alpha_{x,y} & (x\ne y)\\ 1-\sum_{x\ne z} q_x(z) \alpha_{x,z} & (x=y) \end{cases} \]

容易验证细致平衡条件成立。如果提议分布合适这个 Markov 链也是正常返的。

如何模拟这个 Markov 链呢？过程与拒绝采样极其类似。

对于当前状态 \(x\)，依照提议分布 \(q_x\) 生成一个 \(y\)，以 \(\alpha_{x,y}\) 的概率接受这个状态。如果接受，下一个时刻转移到 \(y\)，否则的话留在原地。

这个算法与拒绝采样拥有一些类似的问题，比如好的提议分布不好设计、接受率过低会导致效率低。

Gibbs 算法

注意到我们经常对于高维分布运行 MCMC 方法，所以我们可以针对高维分布的特性设计 Markov 链。

Gibbs 算法要求高维分布的在除了一个维度之外的其它维度的分量都已知的情况下，剩下一个维度的条件分布是易于生成的。比如说对于 Ising 模型，其它维度分量已知的情况下，当前维度只有两种情况，分别计算这两种情况的概率密度，按照比值生成 \(\pm 1\) 随机数就行了。

我们就有这样一个算法：随机生成一个起始点。不断地循环遍历所有维度，设当前维度为 \(x\)，计算固定所有其它维度上的分量后 \(x\) 的条件分布，更新维度 \(x\) 的分量为依据此分布生成的随机数。

如果我们将每次遍历所有维度的一轮迭代看作 Markov 链一次迭代，我们可以验证这个方法就是一种 MCMC 方法。

假设状态空间的维度为 \(d\)，设一轮迭代的概率转移矩阵为 \(P\)，设对于维度 \(i\)，单次迭代的概率转移矩阵为 \(P_i\)，我们有 \(P=\prod_{i=1}^d P_i\)。容易验证对于单个 \(P_i\)，细致平衡条件成立（运用条件概率的定义拆开易证），从而对于 \(P\) 细致平衡条件成立。

注意一下，这里为什么不直接把单个 \(P_i\) 当做转移矩阵构建 Markov 链呢？单个 \(P_i\) 也满足细致平衡条件啊？这是因为单个 \(P_i\) 构建的 Markov 链不联通！当然，这也不意味着 \(P\) 本身一定联通。如同作业题展示的那样，对于某些奇怪的目标分布，Gibbs 算法构建的 Markov 链也不一定联通。

MCMC 估算高维凸集体积

对于 \(d\) 维凸集，生成一系列同心高维球 \(S_1\subseteq S_2\subseteq \dots\subseteq S_k\)，使得 \(S_1\subseteq R\subseteq S_k\)。令相邻两个球的半径比值为 \(1+\frac{1}{d}\)，那么相邻两球的体积比值为 \((1+\frac{1}{d})^d\le e\)。一共有 \(O(\log_{(1+\frac{1}{d})} r)=O(d\ln r)\) 个球壳（这是因为 \(\ln(1+x)\sim x\)），其中 \(r\) 为 \(S_1\) 与 \(S_k\) 的半径之比。

我们考虑估算所有的 \(\frac{V(S_i\cap R)}{V(S_{i-1}\cap R)}\)，由凸集的性质，我们可以将 \(S_i\cap R\) 的所有点向着圆心按照比率 \(\frac{d}{d+1}\) 线性压缩进 \(S_{i-1}\cap R\) 中，因此该比值同样具有上界 \(e\)。

我们考虑在 \(S_i\cap R\) 中均匀采样，统计其出现在 \(S_{i-1} \cap R\) 中的频率，该频率不会显著低于 \(\frac{1}{e}\)。这样我们只要采样多项式级的点就可以保证误差了。

那么具体如何在 \(S_i\cap R\) 中均匀采样呢？直接构建 \(d\) 维网格，只保留网格中 \(S_i\cap R\) 中的点，构建对称随机游走的 Markov 链，运行 MCMC。可以证明（我不会证明）需要的迭代次数（即后面讲的 Mixing Time）是多项式级的。

这个算法要运行球体个数次，从而估算整个高维凸集体积在常数相对误差的情况下，只需要关于 \(d,r\) 的多项式级时间即可完成。

后面讲 Mixing Time 的时候有讨论网格的 Mixing Time，也许凸集网格的 Mixing Time 的上界道理是类似的。

MCMC 方法在其它问题上的应用

学习完 MCMC 算法，感觉这个算法跟 OI 中的某个算法思想很像？模拟退火！

事实上，注意到 Gibbs 分布的概率密度正比于 \(\exp(-\frac{E}{kT})\)，这不就是我们退火的接受概率吗？所以说 Gibbs 分布也就是模拟退火的物理意义了。

而且，对于 Metropolis Hastings，每次迭代的过程形如：如果对面的权值更高，一定转移状态，否则以一定概率接受新状态，这个过程也让我感受到退火的味道。

张老师在上课还举了一个例子，对于最小点覆盖问题，设计能量函数跟点覆盖的大小有关。然后同理对于这个能量，构造 Gibbs 分布跑 MCMC，迭代足够多次就可以求出一组很优秀的点覆盖。这个过程在我看来似乎与对于点覆盖问题运行模拟退火本质上没什么区别？

Markov 链的收敛速度

讨论了这么久 MCMC 方法，一个重要的问题是，到底要迭代几次，我们才能得到一个足够满意的样本点呢？这里为了简化问题，只考虑有限状态链。

Mixing Time

首先，我们需要一个标准衡量 Markov 链的收敛速度。一个想法是，无论我们从何种分布出发，迭代固定次数之后，当前的前缀平均分布 \(\mathbf a(t)\) 与目标分布 \(\pi\) 的足够接近。

关键就是，如何衡量两个分布的是否相近呢？

我们这里使用 \(L_1\) 距离或者总变差距离。

总变差距离的定义是对于分布 \(f,g\)，\(\Vert f-g\Vert_{TV}=\frac{1}{2}\sup_A |f(A)-g(A)|\)，即对于同一个事件，概率相差的最大值。

另一个本质完全相同的定义 \(L_1\) 距离，\(\Vert f-g\Vert_1=\sum_{x\in \Omega} |f(x)-g(x)|\)，即把样本空间中每个原子事件的概率看作一个向量之后，两个向量之间的 \(L_1\) 距离。可以发现 \(L_1\) 距离就是总变差距离的两倍。

对于不可约正常返的 Markov 链，定义 \(\epsilon\)-Mixing Time 为一个最小的非负整数 \(t\)，无论起始分布如何，迭代 \(t\) 次之后，当前的前缀平均分布 \(\mathbf a(t)\) 与稳定分布 \(\pi\) 一定满足 \(\Vert \mathbf a(t)-\pi \Vert_1\le \epsilon\)。

导流率

什么东西在限制一个 Markov 链迟迟不收敛，也就是所谓的 "Mixing" 呢？

比如说带权无向图上的随机游走，可以发现，阻止概率扩散到整个图上的结构是“瓶颈”，即这张图可以分割成两个部分，使得这两个部分之间的连接的边权总和很少。

考虑用数值描述一个 Markov 链的“瓶颈”程度。对于 Markov 链，将其状态集合划分成两个集合 \(S\cup \bar S\)。定义 \(\pi(T)=\sum_{x\in T} \pi_x\)，则该划分的导流率 (conductance) 为：

\[\Phi(S)=\frac{\sum_{x\in S,y\in \bar S} \pi_x P_{x,y}}{\min\left\{\pi(S),\pi(\bar S)\right\}}=\frac{\sum_{y\in \bar S,x\in S} \pi_y P_{y,x}}{\min\left\{\pi(S),\pi(\bar S)\right\}}=\Phi(\bar S) \]

（中间的等号是因为在平稳分布下，一个集合流出的概率一定等于流入的概率。 \(\Phi(S)=\Phi(\bar S)\) 说明了导流率是对于一个划分良定义的。）

定义整个 Markov 链的导流率 \(\Phi\) 为所有非平凡（两边非空）划分的导流率最小值。

Mixing Time 的简单下界

注意到对于一个状态集合划分 \(S\cup \bar S\)，\(\pi(S)\le \frac{1}{2}\)，如果我们令初始分布为所有的在 \(S\) 中的状态 \(x\) 概率为 \(\frac{\pi_x}{\pi(S)}\)，否则为 \(0\)。

这样的话容易发现对于服从该初始分布的一个样本点，每一次迭代进入 \(\bar S\) 的概率为 \(\Phi(S)\)，如果没有成功进入 \(\bar S\) 那么该样本点的服从概率分布仍然是初始分布。因此，该样本点首次进入 \(\bar S\) 的期望步数为 \(\frac{1}{\Phi(S)}\)。

Mixing Time 的意义就是为了保证 MCMC 的采样的迭代次数。因此我们知道在这种情况下 MCMC 至少需要迭代 \(\Omega(\frac{1}{\Phi})\) 次。所以，对于常数 \(\epsilon\)，Mixing Time 具有下界 \(\Omega(\frac{1}{\Phi})\)。

（总感觉这个书上的论证有一点点说明不充分，不过还是能让我们很直观的感受到这个结论的正确性）

对称链的 Mixing Time 上界

重头戏来喽！

由于 MCMC 算法构造的 Markov 链都满足细致平衡条件，所以如果我们仅仅考虑对称链的 Mixing Time，就可以得到一个形式非常好的上界。

具体地，对于不可约对称链，我们可以证明 \(\epsilon\)-Mixing Time 拥有上界 \(O(-\frac{\ln(\pi_{\min})}{\Phi^2\epsilon^3})\)，其中 \(\pi_{\min}\) 是稳定向量的所有分量中的最小值。

设 \(t=-\frac{c\ln(\pi_{\min})}{\Phi^2\epsilon^3}\)，其中 \(c\) 是某个合适的常数。我们需要证明 \(\Vert \mathbf a(t)-\pi\Vert_1\le \epsilon\)。

设 \(v_i=\frac{\mathbf a(t)_i}{\pi_i}\)，对 \(v_i\) 从大到小排序，按照排序的顺序对于所有状态重标号，使得状态的标号满足 \(v_1\ge v_2\ge \cdots\ge v_{i_0}> 1\ge v_{i_0+1}\ge v_{i_0+2}\ge \cdots\)。称 \(v\) 值比一大，即下标在 \(i_0\) 前的一部分状态叫重状态，否则叫做轻状态。

根据 \(L_1\) 距离的定义，我们有：

\[\Vert \mathbf a(t)-\pi\Vert_1=\sum_{i} |v_i-1|\pi_i=2\sum_{i\le i_0} (v_i-1)\pi_i=2\sum_{i>i_0} (1-v_i)\pi_i \]

不妨考虑反证法，假设 \(\sum_{i>i_0} (1-v_i)\pi_i> \frac{\epsilon}{2}\)，利用这个条件去证明 \(\sum_{i\le i_0} (v_i-1)\pi_i\le \frac{\epsilon}{2}\) 导出矛盾。

考虑如何估计 \(\sum_{i\le i_0} (v_i-1)\pi_i\)。这个和式可以看作求阶梯状网格的面积，我们换一种考虑方式，将这个网格转置，得到这个求和的另一个形式：

\[\sum_{i\le i_0} (v_i-1)\pi_i=(v_{i_0}-1)\sum_{j\le i_0}\pi_j+\sum_{i< i_0} (v_{i}-v_{i+1}) \sum_{j\le i} \pi_j \]

这里如果记 \(\gamma_i=\sum_{j\le i} \pi_j\)，我们考虑将这个和式的 \(i_0\) 项分成连续的 \(k\) 组，第 \(i\) 组的左端点为 \(l_i\)，并且假设 \(l_{k+1}=i_0+1\), 那么有：

\[\begin{aligned} \sum_{i\le i_0}(v_i-1)\pi_i&=(v_{i_0}-1)\gamma_{i_0}+\sum_{i<i_0}(v_i-v_{i+1})\gamma_i\\ &\le (v_{l_k}-1)\gamma_{i_0}+\sum_{i<k}(v_{l_i}-v_{l_{i+1}})\gamma_{l_{i+1}-1}\\ &\le \sum_{i\le k}(v_{l_i}-v_{l_{i+1}})\gamma_{l_{i+1}-1} \end{aligned} \]

我们具体如何进行分组呢？从左到右确定所有 \(l_i\)，首先 \(l_1=1\)，假设 \(l_1,l_2,\dots,l_{x}\) 都确定了，那么 \(l_{x+1}\) 在满足如下条件的前提下尽可能大：

\[l_{x+1}\le i_0+1,\gamma_{l_{x+1}-1}-\gamma_{l_x}\le \frac{\epsilon \Phi\gamma_{l_x}}{4} \]

不断地取最大的边界，直到合法的 \(l_{x+1}=i_0+1\)。

这样的分组方式，一共会分出多少组呢？除了最后一个组，所有组 \(x\) 都满足 \(\gamma_{l_{x+1}}-\gamma_{l_x}>\frac{\epsilon \Phi \gamma_{l_x}}{4}\)，即 \(\gamma_{l_{x+1}}>(1+\frac{\epsilon\Phi}{4})\gamma_{l_x}\)。由于 \(\forall i,\gamma_i\le 1\)，且 \(\gamma_1\ge \pi_{\min}\)，因此组的个数是 \(O(\log_{(1+\frac{\epsilon\Phi}{4})} \frac{1}{\pi_{\min}})=O(-\frac{\ln \pi_{\min}}{\epsilon\Phi})\) 级别的（这里再次用到 \(\ln(1+x)\sim x)\)。

那么在同一组中，\((v_{l_i}-v_{l_{i+1}})\gamma_{l_{i+1}-1}\) 的上界如何确定呢？

设 \(l,r\) 分别为当前组的左右边界，即对于当前考虑的组 \(i\)，\(l=l_{i},r=l_{i+1}-1\)。我们有：

\[\sum_{i=1}^l \mathbf a(t)_i-[\mathbf a(t)P]_i=\frac{1}{t}\sum_{i=1}^l \mathbf p(0)_i -\mathbf p(t)_i\le \frac{1}{t} \]

（其中第二步我们已经在 Markov 链的收敛性一节证明过了）

考虑这个式子的实际意义，即将 \(\mathbf a(t)\) 当做 Markov 链的输入分布，经过一次迭代之后，考虑前 \(l\) 个数的分布总概率流出了多少。由流出量的另一个计算方式，我们得出流出量的下界：

\[\begin{aligned} &\sum_{i\le l,j>l} \mathbf a(t)_i P_{i,j}-\mathbf a(t)_j P_{j,i}\\ =&\sum_{i\le l,j>l} \mathbf \pi_i P_{i,j} v_i-\mathbf \pi_j P_{j,i} v_j\\ =&\sum_{i\le l,j>l} \mathbf \pi_j P_{j,i} (v_i-v_j)\\ \ge&\sum_{i\le l,j>r} \mathbf \pi_j P_{j,i} (v_i-v_j)\\ \ge&(v_l-v_{r+1})\sum_{i\le l,j>r} \mathbf \pi_j P_{j,i} \end{aligned} \]

（其中第二个等号运用了对称链的细致平衡条件。）

而由于反证法的假设 \(\sum_{i>i_0} (1-v_i)\pi_i> \frac{\epsilon}{2}\)，我们有 \(\sum_{i>i_0} \pi_i>\frac{\epsilon}{2}\)，如果我们将前 \(l\) 个状态看成集合 \(A\)，那么 \(\min\{\pi(A),\pi(\bar A)\}\ge \frac{\epsilon\gamma_l}{2}\)（\(\epsilon\) 显然应该比 \(2\) 小，否则 Mixing Time 没有任何意义）。

根据导流率的定义，\(\sum_{i\le l,j>l} \pi_j P_{j,i}=\Phi(A)\min\{\pi(A),\pi(\bar A)\}\ge \frac{\epsilon\Phi\gamma_l}{2}\)。

而 \(\sum_{i\le l,l<j\le r} \pi_j P_{j,i}=\sum_{l<j\le r} \pi_j \sum_{i\le l} P_{j,i}\le \sum_{l<j\le r} \pi_j=\gamma_{r}-\gamma_{l}\)。注意到我们选取分组的方式保证了 \(\gamma_r-\gamma_l\le \frac{\epsilon\Phi\gamma_l}{4}\)。我们有 \(\sum_{i\le l,j>r} \pi_jP_{j,i}=\sum_{i\le l,j>l} \pi_j P_{j,i}-\sum_{i\le l,l<j\le r}\pi_j P_{j,i}\ge \frac{\epsilon\Phi\gamma_l}{4}\)。

将这个下界带入 \(A\) 的流出量下界，再结合 \(A\) 的流出量上界，我们有：

\[(v_l-v_{r+1})\cdot \frac{\epsilon\Phi\gamma_l}{4}\le \frac{1}{t} \]

因此我们知道了，\((v_l-v_{r+1})\gamma_l\) 具有上界 \(O(\frac{1}{\epsilon \Phi t})\)。

对于原来的求和 \(\sum_{i\le k}(v_{l_i}-v_{l_{i+1}})\gamma_{l_{i+1}-1}\) 每一组的求和具有上界 \(O(\frac{1}{\epsilon \Phi t})\)，一共 \(O(-\frac{\ln \pi_{\min}}{\epsilon \Phi})\) 组，代入 \(t=-\frac{c\ln \pi_{\min}}{\Phi^2\epsilon^3}\)，并选取合适的 \(c\)，我们就可以证明 \(\sum_{i\le i_0}(v_i-1)\pi_i\le \frac{\epsilon}{2}\) 与反证法假设矛盾。