基于 Metropolis 的朗之万算法

1. 未经调整的朗之万算法

未调整的朗之万算法 (ULA) 是一个离散时间马尔可夫链 ${\mathbf{U}}_n$，它是对普通朗之万扩散 ${\mathbf{L}}_t$的自然离散化。

任何使用$$\Vert P_{\mathrm{M}}^n(\mathbf{x},\cdot )-\pi \Vert \to 0$$的简单算法都可以通过这种方式构造，例如 Parisi (1981) 或 Grenander 和 Miller (1994) 所描述的。

我们将看到，该算法可能具有一些不理想的收敛性质，尽管由于其实现可能比某些更稳健的替代方案需要较少的计算开销，因此它仍可能具有实际价值。

为了形成这个链，给定${\mathbf{U}}_{n-1}$，我们只需根据以下公式构造 ${\mathbf{U}}_n$:
$$N\left({\mathbf{U}}_{n-1}+\frac{1}{2}\nabla \log \pi \left({\mathbf{U}}_{n-1}\right),h{I}_k\right)$$正如 Besag (1994) 所指出的，这个链仅能近似维持 $\pi$ 的不变性：例如，如果 $\pi$ 本身在$\mathbb{R}$上是 $N(0,1)$ ，那么当 $h=2$ 时，我们有 $U_n\sim N(0,2)$ ，这显然表明如果离散化步长$h$ 如此粗䊁，那么我们会得到立即"收敛"，但却是到一个完全不期望的分布。

ULA 链实际上可能表现得相当糟糕：例如，即使原始扩散是指数遍历的，它可能会收敛但并非几何快速收敛，或者更为惊人的是，它实际上可能是一个瞬态链，尽管 ${\mathbf{L}}_t$ 具有非常良好的不变分布。

2. 基于 Metropolis 的朗之万算法 (MALA)

2.1. MH算法

这些算法首先考虑一个候选转移核，其密度为 $q(\mathbf{x},\mathbf{y})$，其中 $\mathbf{x},\mathbf{y}\in \mathrm{X}$，用于生成在 $X$ 上演化的离散时间马尔可夫链的潜在转移。在此，我们通常将 $X$ 视为 ${\mathbb{R}}^k$ 的子集，并配备了 Borel $\sigma$-代数 $\mathcal{B}$，同时 $\pi (\mathbf{y})$ 和 $q(\mathbf{x},\mathbf{y})$ 都是相对于 Lebesgue 测度 ${\mu }^{\text{Leb }}$ 的密度，尽管更一般的形式化也是可能的。

根据密度 $q(\mathbf{x},\cdot )$ 生成的“候选转移”到 $\mathbf{y}$ 被接受的概率为 $\alpha (\mathbf{x},\mathbf{y})$，其表达式为
$$(1)\alpha (\mathbf{x},\mathbf{y})=\{\begin{array}{ll}\min \left\{\frac{\pi (\mathbf{y})}{\pi (\mathbf{x})}\frac{q(\mathbf{y},\mathbf{x})}{q(\mathbf{x},\mathbf{y})},1\right\}& \pi (\mathbf{x})q(\mathbf{x},\mathbf{y})>0\\ 1& \pi (\mathbf{x})q(\mathbf{x},\mathbf{y})=0\end{array}$$因此，Hastings 链的实际转移（我们记作 ${\Phi }_n$）根据转移概率密度的规律 $P$ 进行，其转移概率密度为
$$(2)p(\mathbf{x},\mathbf{y})=q(\mathbf{x},\mathbf{y})\alpha (\mathbf{x},\mathbf{y}),\mathbf{y}\ne \mathbf{x}$$且保持在同一点的概率为
( 3 ) r ( x ) = P ( x , { x } ) = ∫ q ( x , y ) [ 1 − α ( x , y ) ] d y (3)\quad r(\mathbf{x})=P(\mathbf{x},\{\mathbf{x}\})=\int q(\mathbf{x}, \mathbf{y})[1-\alpha(\mathbf{x}, \mathbf{y})] \mathrm{d} \mathbf{y} (3)r(x)=P(x,{x})=∫q(x,y)[1−α(x,y)]dy通过选择这样的 $\alpha$，我们有 $\pi$ 是不变测度：即满足 $\pi (A)=\int \pi (\mathbf{x})P(\mathbf{x},A)\mathrm{d}\mathbf{x},\mathbf{x}\in \mathrm{X},A\in \mathcal{B}$。

只要链具有适当的不可约性和非周期性，那么标准结果表明，定义为 $n$ 步转移概率 $P^n(\mathbf{x},A)=P\left({\Phi }_n\in A\mid {\Phi }_0=\mathbf{x}\right)$ 对于每个 $n\ge 1$ 而言，$\mathbf{x}\in \mathrm{X},A\in \mathcal{B}$，在全变差范数下收敛于 $\pi$：即对于几乎所有 $\pi$ 上的 $\mathbf{x}$，
$$(4)\Vert P^n(\mathbf{x},\cdot )-\pi \Vert :=\frac{1}{2}\underset{A\in \mathcal{B}}{\sup }\left|P^n(\mathbf{x},A)-\pi (A)\right|\to 0$$

2.2. 基于 Metropolis 的朗之万算法 (MALA)

根据 Besag (1994) 的建议，我们引入了进一步的修改，并遵循 (1) 和 (2) 式的结构，构造了基于 Metropolis 的朗之万算法 (MALA)。

这是一个 Hastings-Metropolis 链 ${\mathbf{M}}_n$，它使用 ULA 来构造候选链。因此，在给定 ${\mathbf{M}}_{n-1}$ 的情况下，${\mathbf{U}}_n$ 首先被设为如下分布的变量：
$$N\left({\mathbf{M}}_{n-1}+\frac{1}{2}h\nabla \log \pi \left({\mathbf{M}}_{n-1}\right),h{I}_k\right)$$将此提议密度记为 $q\left({\mathbf{M}}_{n-1},{\mathbf{U}}_n\right)$。接下来执行接受/拒绝步骤，接受 ${\mathbf{U}}_n$ 的概率为
$$(5)\alpha \left({\mathbf{M}}_{n-1},{\mathbf{U}}_n\right)=1\wedge \frac{\pi \left({\mathbf{U}}_n\right)q\left({\mathbf{U}}_n,{\mathbf{M}}_{n-1}\right)}{\pi \left({\mathbf{M}}_{n-1}\right)q\left({\mathbf{M}}_{n-1},{\mathbf{U}}_n\right)}$$如果 ${\mathbf{U}}_n$ 被接受，则设 ${\mathbf{M}}_n={\mathbf{U}}_n$，否则令 ${\mathbf{M}}_n={\mathbf{M}}_{n-1}$。通过 Hastings 构造，如 (2) 和 (3) 式，MALA 链收敛于 $\pi$，其意义是
$$\Vert P_{\mathrm{M}}^n(\mathbf{x},\cdot )-\pi \Vert \to 0$$对于几乎所有 $\pi$ 上的 $\mathbf{x}$，其中我们写作 $P_{\mathrm{M}}^n(\mathbf{x},A)=P\left({\mathbf{M}}_n\in A\mid {\mathbf{M}}_0=\mathbf{x}\right)$：这遵循于链在 Roberts 和 Tweedie (1996) 中明确为 ${\mu }^{\text{Leb }}$-不可约且非周期性的结果。作为我们结果的一个次要但有用的副产品，我们展示了在几何遍历的情况下，收敛性也适用于所有起始点。

寻找几何速率收敛且适用于每个起始点的条件。

1. 当 ULA 是瞬态时，MALA 不是指数遍历的；
2. 在 ULA 不是瞬态的情况下，MALA 通常是几何遍历的，意味着它可以较快地收敛到目标分布。如果目标分布的尾部比指数分布更重（即，目标分布在远离中心的区域衰减得比指数分布更慢），那么这种快速的几何收敛性可能会受到影响。

3. Metropolis 调整的朗之万截断算法（MALTA）

最后，我们简要提到对算法进行的一个简单调整，旨在尝试结合随机游走 Metropolis 算法和“目标”朗之万候选 ULA 的最佳特性。我们称此算法为 MALTA（Metropolis 调整的朗之万截断算法）。这个修订算法涉及用截断候选分布替换第一个 ULA 近似：
$${\mathbf{T}}_n\sim N\left({\mathbf{M}}_{n-1}+R\left({\mathbf{M}}_{n-1}\right),h{I}_k\right)$$其中漂移项现在为：
$$R\left({\mathbf{M}}_n\right)=\frac{D\nabla \log \pi (\mathbf{x})}{2(D\vee |\nabla \log \pi (\mathbf{x})|)}$$其中 $D>0$ 是某个常数。

然后，调整候选跳跃 ${\mathbf{T}}_n$ 以确保正确的平稳分布成立，如在 (5) 中所示。

使用 MALTA，链具有更加稳健的几何遍历性。我们不对 MALTA 进行详细分析，仅指出本文以及 Roberts 和 Tweedie (1996) 中使用的方法可以很容易地应用于该算法的分析。