第8章 · Green 函数蒙特卡洛#

8.1 基本概念与形式推导#

本章我们将介绍实现所谓 Green's Function Monte Carlo (GFMC) 技术所需的基本原理和更高级的细节。在这种方法中，通过使用1.7节讨论的幂法，对给定哈密顿量的基态波函数进行随机采样。正如第5章对变分蒙特卡罗方法所做的那样，我们在希尔伯特空间中固定一个完备基组{|x-}：

其中状态被视为正交归一的。此外，我们假设对于给定的哈密顿量 $\mathcal{H} .$ ，矩阵元 $\mathcal{H}_{x^{\prime} , x} = \langle x^{\prime} | \mathcal{H} | x \rangle$ 可以对每个|x-和$| x^{\prime} \rangle$ 高效计算。对于局域哈密顿量，尽管希尔伯特空间的维数随系统大小 $L_{☉}$ 呈指数增长，但对每个 $| x \rangle$ ，非零元素 $\mathcal{H}_{\boldsymbol{x}^{\prime} , \boldsymbol{x}}$ 的数量与 $L$ 成正比；因此，给定基态 $| x \rangle$ ，所有矩阵元 $\mathcal{H}_{\boldsymbol{x}^{\prime} , \boldsymbol{x}}$ 可以在合理的计算量内求得。利用这一性质，我们可以定义一种随机算法，以统计方式执行幂法（参见1.7节）（Ceperley and Alder, 1980; Reynolds et al., 1982; Foulkes et al., 2001）。这里，基态波函数 $| \Upsilon_{0} \rangle$ 通过使用合适的投影算子从初始态 $| \Psi_{0} \rangle$ 中筛选出来：

其中  是一个对角算子，满足 $\Lambda_{x , x} ~ = ~ \lambda$ 。逼近基态 $| \Upsilon_{0} \rangle$ 的必要条件是取足够大的λ值，并选择初始态 $\left| {\Psi_{0}} \right.$ 使得 $\langle \Upsilon_{0} | \Psi_{0} \rangle \neq 0$ 。在实践中，式(8.2)可以迭代实现：

在给定基组上展开，我们得到：

其中 $\Psi_{n} ( x ) = \langle x | \Psi_{n} \rangle$ 而

就是所谓的格林函数（尽管它并不对应于真正的格林函数）。在格点上，GFMC 技术特别简单，并且首次被用于研究方格上海森堡模型的基态性质（Trivedi and Ceperley, 1989, 1990）。

在幂法的统计实现中，我们可能会倾向于将式(8.4)解释为随机变量的主方程，其中 $\Psi_{n} ( x )$ 表示第n次迭代的概率分布，而 $\mathcal{G}_{x^{\prime} , x}$ 是转移概率。然而，为了最终获得这种量子演化的统计解释，必须阐明一些重要问题。首先，在这种方法中，格林函数 $\mathcal{G}_{x^{\prime} , x}$ 的所有矩阵元必须非负，这是定义转移概率所要求的。就对角元 ${\mathcal G}_{x , x}$ 而言，这没有问题，因为我们总可以定义足够大且为正的λ值使得 $\mathcal{G}_{x , x} ~ \geq ~ 0$ 。相反，要求非对角元非负是非常不平凡的，实际上只有某些模型满足这一条件；如果对于某些对 $( x^{\prime} , x )$ 有 $\mathcal{G}_{x^{\prime} , x} \ < \ 0$ ，我们就说存在符号问题，这将在第10章中讨论。接下来，我们假设哈密顿量使得对所有对 $( x , x^{\prime} )$ 都有 $\mathcal{G}_{x^{\prime} , x} \geq 0$ 。此外，即使格林函数非负，通常我们仍有：

因此，$\mathcal{G}_{x^{\prime} , x}$ 不能被解释为转移概率，因为它没有归一化。注意，方程(8.6)表明在方程(8.4)的投影技术中，波函数的归一化并不守恒。尽管如此，我们总可以将格林函数分解为两个因子的乘积：一个代表真实转移概率的随机矩阵 $p_{x^{\prime} , x}$，以及一个实数 $b_{x}$，它是格林函数的归一化因子：

其中 $b_{x}$ 由下式给出：

然后矩阵 $p_{x^{\prime} , x}$ 可以通过下式得到：

我们现在想要设计一个简单的马尔可夫过程，从中可以得到方程(8.4)中波函数的演化。由于格林函数不能直接解释为转移概率，仅用构型 $x_{n}$（标识迭代步骤 n 时的量子态 |x-）的描述是不充分的。考虑尺度因子 $b_{x}$ 存在的最简单方法是在投影方法的统计描述中增加一个权重 $w_{n}$。因此，马尔可夫过程由一个二元组 $( x_{n} , w_{n} )$ 定义；根据方程(8.8)和(8.9)的分解，其演化描述如下：

1 ) \mathrm{以概率} p_{x^{\prime} , x_{n}} \mathrm{生成} x_{n + 1} = x^{\prime} ,\tag{8.10}

2 ) \mathrm{更新权重为} w_{n + 1} = w_{n} b_{x} .\tag{8.11}

注意，这里新的构型 $x_{n + 1}$ 可以与 $x_{n}$ 相同，因为它是从所有满足 $p_{x^{\prime} , x_{n}} \neq 0$ 的构型中选出的，包括当 $p_{x_{n} , x_{n}} \neq 0$ 时的 $x_{n}$ 本身（这种情况发生在 ${\mathcal G}_{x , x} ~ = ~ \lambda - {\mathcal H}_{x , x} ~ \neq ~ 0$ 时）。马尔可夫过程描述了构型 $x_{n}$ 及其权重 $w_{n}$ 的扩散，其转移概率如下：

给定旧的二元组 $(x, w)$，该定义给出了新的二元组 $( x^{\prime} , w^{\prime} )$。因此，对应于概率密度 $\mathcal{P}_{n} ( x , w )$ 的主方程为：

在下文中，假设变量 w 的积分限从 −∞ 到 +∞，而概率密度 $\mathcal{P}_{n} ( x , w )$ 在那些沿马尔可夫链无法达到的 w 值（例如，$w < 0$）处为零。在此形式下，波函数 $\Psi_{n} ( x )$ 可以通过对权重积分概率密度 $\mathcal{P}_{n} ( x , w )$ 乘以 w 得到。实际上，从主方程(8.13)和方程(8.12)的转移概率定义，我们有：

将两边乘以 $w^{\prime}$ 并对 $w^{\prime}$ 积分，我们得到：

这正是方程(8.4)的原始迭代过程，只要进行以下等同：

因此，经过大量迭代 n 后，$\mathcal{P}_{n} ( x , w )$ 收敛到一个平衡分布 $\mathcal{P}_{\mathrm{eq}} ( x , w )$，该分布决定了基态波函数 $\Upsilon_{0} ( x )$

下面，我们将描述马尔可夫过程的一种实际实现方法，它允许我们在根据式 (8.4) 的迭代过程滤除 $| \Psi_{0} \rangle$ 中存在的高能成分后，计算基态 $| \Upsilon_{0} \rangle$ 上的可观测量。

8.2 单行走者技术#

在这里，我们解释如何用最简单的方法实现式 (8.10) 和 (8.11) 的马尔可夫过程。为此，我们定义该随机过程的基本元素——所谓的行走者。行走者由一个标记组态|x>的指标 $x$ 和一个权重 $w$ 决定，权重关联于波函数的振幅，如式 (8.16) 所示。行走者通过具有离散迭代时间 $n$ 的马尔可夫过程改变其组态和权重：二元组 $( x_{n} , w_{n} )$（表示 $n$ 时刻的行走者）按照 $\mathcal{P}_{n} ( x , w )$ 分布。最重要的是，行走者在统计意义上决定了量子态 $\Psi_{n} ( x )$

其中 $\langle \langle \dots \rangle \rangle$ 表示对以组态 $( x_{n} , w_{n} )$ 进行的（无限）多次独立马尔可夫链实现的统计平均。我们想指出，在实际操作中很少会真正确定波函数 $\Psi_{n} ( x )$，因为在多体系统中希尔伯特空间呈指数增长，波函数的信息甚至无法存储在计算机内存中。此外，为了对组态获得足够的误差限，每个组态必须被多次访问，这意味着计算时间呈指数增长。值得注意的是，量子蒙特卡洛方法基于这样一个事实：为了可靠估计能量或关联函数，并不需要精确的波函数统计信息。

在实际中，我们可以使用一个行走者，从初始条件 $x_{0}$（对应给定的多体组态 $\left| x_{0} \right.$）和权重 $w_{0} = 1$ 开始，从而有：

然后，我们利用式 (8.10) 和 (8.11) 演化行走者。对于格点上的通用哈密顿量，马尔可夫过程可以非常容易地实现，因为对于给定的 $x_{n}$，随机矩阵 $p_{x^{\prime} , x_{n}}$ 中非零元素的数量很少，并且随着格点数量 $L$ 的增长而适度增加。因此，为了定义 $x_{n + 1}$，只需将区间 [0, 1) 划分为与 $x_{n}$ 以非零概率 $p_{x^{\prime} , x_{n}}$ 相连的所有可能 $\{x^{\prime} \}$ 对应的更小区间，并生成一个介于 0 和 1 之间的随机数 $\xi$。如第 3.5 节所述，$\xi$ 将位于上述某个区间内，命中对应特定 $x^{\prime}$ 区间的概率恰好等于 $p_{x^{\prime} , x_{n}}$，见图 8.1。权重 $w_{n + 1}$ 只需将 $w_{n}$ 乘以 $b_{x_{n}}$ 即可得到。

最后，我们讨论仅对组态 $x$ 的边际概率的演化：

由于在单行走者实现中，组态 $x_{n}$ 的演化与其权重无关（如式 (8.10) 和 (8.11) 所示），因此可以很容易写出 $\Pi_{n} ( x )$ 的主方程：

图8.1 将区间 $[ 0 , 1 )$ 划分为若干个长度为概率 $p_{x_{i}^{\prime} , x_{n}}$ 的子区间，其中 $i$ 标记所有可能的组态 $\breve{x^{\prime}}$ 且对应非零 $p_{x_{i}^{\prime} , x_{n}}$ （这里显示有 5 个条目的情况）。然后，生成一个随机数 $0 \leq \xi < 1$，并选择新的组态 $x_{n + 1}$ （此处 $i = 3$），因为 $\xi$（用箭头表示）落在相应的区间内，详见第 3.5 节。

因此，在平衡态下，边缘概率由随机矩阵 $p_{x^{\prime} , x}$ 特征值 $\lambda_{0} = 1$ 的右特征向量给出，即：

确实，我们有：

其中我们只利用了矩阵 $\mathcal{G}_{x^{\prime} , x}$ 是对称矩阵这一事实。那么，对于遍历的哈密顿量，这意味着遍历的矩阵 $p_{x^{\prime} , x}$，概率 $\Pi_{n} ( x )$ 会收敛（以指数速率）到 $\Pi_{\mathrm{eq}} ( x )$，无论初始条件如何，详见第 3.8 节。

8.2.1 基态能量#

现在让我们讨论如何计算物理可观测量，并首先考虑基态能量 $E_{0}$。由于对于足够大的 $n$，有 $\mathcal{HG}^{n} | \Psi_{0} \rangle \propto E_{0} | \Upsilon_{0} \rangle$，我们得到：

利用式 (8.18) 中 $\Psi_{n} ( x ) = \langle x | {\mathcal{G}}^{n} | \Psi_{0} \rangle$ 的形式，得到：

其中我们将局域能量定义为：

因此，基态能量可以通过进行若干独立的 $n$ 步计算（$n$ 足够大以使初始状态 $\left| \Psi_{0} \right\rangle$ 中的基态波函数被过滤出来），然后将它们平均得到。

一种更高效、更直接的计算方法是进行一次长度为 $N_{\mathrm{sim}} \gg n$ 的模拟；经过热化时间后，马尔可夫过程中的组态 $x$ 将根据式 (8.22) 的 $\Pi_{\mathrm{eq}} ( x )$ 达到平衡；然后，我们可以设想从每一步 $n-p$ 开始一个长度为 $p$ 的投影技术。在这个意义上，初始概率不再由式 (8.19) 给出，而是由：

经过 $p$ 步后，每个投影过程的累积权重为：

然后，通过在累积权重中取足够大的 $p$ 值，可以估计基态能量：

在实际操作中，为避免数值上溢/下溢，将每个马尔可夫步的因子 $b_{x_{n}}$（甚至更好的 $b_{x_{n}} / \bar{b}$，其中 $\bar{b}$ 是 $b_{x_{n}}$ 平均值的估计）存储起来（例如存储在硬盘上）而不进行大量乘法运算是很有用的。在后处理分析中，我们可以计算能给出收敛能量的最小可能 $p$ 值所对应的 $G_{n}^{p}$，并通过方程(8.29) 评估 $E_{0}$。

这一简单过程给出了基态能量的精确估计。然而，仍需注意几点。首先，在计算误差棒时，必须谨记：将每次迭代作为投影技术的起点时，我们得到的是相关测量值，这需要如第3.11节所述的分块方法。最重要的是，根据方程(8.26)，能量的计算不满足零方差性质，因为随机量 $e_{L} ( x )$ 不依赖于任何变分猜测。从这个意义上说，引入重要性抽样对于减少统计波动至关重要。

8.3 重要性抽样#

现在让我们描述如何在GFMC技术中实现重要性抽样。这可以通过引入所谓的引导函数 $\Psi_{G} ( x )$ 来实现，该函数必须是实数（见下文），并且对所有构型 $| x \rangle$ 均非零。通常，它被选为最佳变分波函数或与其非常接近的波函数。将方程(8.4) 两边乘以 $\Psi_{G} ( x^{\prime} )$，我们得到：

然后，通过定义带有重要性抽样的格林函数：

我们得到波函数

满足与方程(8.4) 相同的演化，只需将 $\mathcal{G}$ 替换为 $\tilde{\mathcal{G}}$：

我们想指出一个重要事实：方程(8.31) 表示一个相似变换（在第3.8节讨论马尔可夫链时也曾使用）。这种变换不会改变矩阵的谱，即 $\mathcal{G}$ 和 $\tilde{\mathcal{G}}$ 具有相同的特征值，但会平凡地改变特征向量：如果 $\upsilon(x)$ 是 $\mathcal{G}$ 的特征向量，那么 $\tilde{\upsilon} ( x ) = \Psi_{G} ( x ) \upsilon ( x )$ 是 $\tilde{\mathcal{G}}$ 的（右）特征向量。注意，通常 $\tilde{\mathcal{G}}_{x^{\prime} , x}$ 不再是对称的。

当 $\tilde{\mathcal{G}}_{x^{\prime} , x}$ 对所有对 $( x^{\prime} , x )$ 都是实数且非负时，可以应用与方程(8.7) 相同的分解，用条件概率和实数系数表示（具有实格林函数这一事实意味着引导函数也必须是实数）：

其中：

然后，可以类比式 (8.10) 和式 (8.11) 来设计相应的马尔可夫过程。

通过使用重要性采样变换，可以实现两个重要事实。首先，当引导波函数 $\Psi_{G} ( x )$ 是正确基态波函数时，无需统计波动即可得到基态能量 $E_{0}$ 的精确值（见下文）。此外，引导函数可以用作解决符号问题的手段：在某些情况下，原始格林函数 $\mathcal{G}_{x^{\prime} , x}$ 并不满足所有非对角矩阵元的非负性条件，而这一事实可以通过包含一个适当选取的 $\Psi_{G} ( x )$ 来调整。正方晶格上的反铁磁最近邻海森堡模型是最简单的例子，其中格林函数可能为负，即对于具有两个最近邻自旋方向相反（上下取向）的组态 $| x \rangle$ 和 $| x^{\prime} \rangle$，有 $\mathcal{G}_{x^{\prime} , x} = - J / 2$。一旦考虑一个具有所谓马歇尔符号 (Marshall, 1955) 的引导函数：

其中 $N_{A , \uparrow} ( x )$ 是组态 |x- 中一个子晶格上的自旋向上数量，则带有重要性采样的格林函数对所有组态对都是非负的。

在存在重要性采样的情况下，第 8.2 节的所有推导均可以重复进行。例如，式 (8.14) 被替换为：

此外，式 (8.18) 变为：

最后，我们可以证明随机矩阵 $\tilde{p}_{x^{\prime} , x}$ 的右本征向量（特征值为 1）由下式给出：

8.3.1 基态能量与关联函数#

我们现在可以计算存在重要性采样时的基态能量。再次强调，对于足够大的 n，我们有：

利用以下事实：

结合式 (8.39) 并考虑引导函数是实函数，我们得到：

其中我们定义了带重要性采样的局域能量：

注意，如果 $| \Psi_{G} \rangle$ 是哈密顿量具有本征值 E 的本征矢，那么所有组态 $| x \rangle$ 的 $\tilde{e}_{L} ( x ) = E$，这意味着不存在统计涨落。当然，对于一般的关联哈密顿量来说，这是一种非常不现实的情况；然而，考虑重要性抽样的优势在于，每当引导函数接近某个本征态（例如精确的基态）时，统计涨落会相对于没有重要性抽样的情况大幅降低。

在实际操作中，和之前一样，我们可以执行一次长度为 $N_{\mathrm{sim}} \gg n$ 的模拟，并在热化时间之后，设想从根据方程(8.40)的 $\tilde{\Pi}_{\mathrm{eq}} ( x )$ 分布的组态开始执行 p 步的投影技术。在这种情况下，从第 $n - p$ 次迭代累积的 p 步权重为：

基态能量由下式给出：

通过使用所谓的正向行走技术，GFMC 方法还可以高效地计算在基组 {|x-} 中对角的局域算符 $\mathcal{O}$ 的期望值：

其中 $\mathcal O ( x )$ 是对应于组态 $| x \rangle$ 的本征值。实际上，对于大的 n 和 m，基态上的真实期望值可以写为：

对于这些算符，马尔可夫链可以轻松修改，以考虑在选定的第 $n - m$ 步施加算符 $\mathcal{O}$。事实上，在统计意义上，这完全等价于修改权重：

使得：

和之前一样，可以通过使用方程(8.45)的累积权重在单次运行中计算右端的量：

从而给出当 $m = p / 2$ 且 $p$ 足够大时的基态期望值。这里的关键点是，式(8.48)的分子和分母都在单次模拟中计算，即使用同一条马尔可夫链。

计算一般算符平均值的另一种方法是通过对式(8.48)的分子和分母分别进行两次独立模拟。这种方法的优势在于不仅可以测量对角算符，还可以测量对于 $\vert x^{\prime} \rangle \neq \vert x \rangle$ 满足 $\mathcal{O}_{x^{\prime} , x} = \langle x^{\prime} | \mathcal{O} | x \rangle \neq 0$ 的一般算符。在这种情况下，我们首先进行一次模拟来评估分母。然后，我们取平衡后的行走者配置 $( x_{n} , \tilde{w}_{n} )$，并在给定迭代步应用算符 $\mathcal{O}$。当算符具有非对角元时，这并不简单地等同于像式(8.49)那样缩放权重。然而，我们可以使用一种随机方法，从通过非零矩阵元 $\mathcal{O}_{x^{\prime} , x}$ 与 x 相连的配置中仅选择一个配置。最后，从所选配置开始，再进行 m 步传播。

只要 $\mathcal{O}_{x^{\prime} , x} > 0$，这种随机方法可以通过以下几个步骤实现：

1. 对于式(8.51)的分母，使用原始的有 n 步的马尔可夫链来计算累积权重 ${\tilde{G}}_{n}^{p}$。

2. 对于式(8.51)的分子：

考虑原始马尔可夫链在迭代步 $n - m$ 时的相应配置 $x_{n - m} \equiv x$ 和权重 ${\tilde{G}}_{n - m}^{p}$。

_• 计算 $\mathcal{O} ( x_{n - m} ) \equiv \mathcal{O} ( x )$

_• 以正比于 $\mathcal{O}_{{x^{\prime}} , x} \Psi_{G} ( x^{\prime} ) / \Psi_{G} ( x )$ 的概率选择一个随机的新配置 $x^{\prime}$。

_• 最后，从 $x^{\prime}$ 启动一条独立的马尔可夫链进行 m 步，并持续将权重从 ${\tilde{G}}_{n - m}^{p}$ 累积到 $\left( {\tilde{G}}_{n}^{p} \right)^{\prime}$。注意，在这种情况下，由于前向传播进行了 m 步，分子和分母上的权重是不同的。

我们要强调一个重要事实：矩阵元中符号任意的一般算符 ${\mathcal{O}}_{:}$，总可以写成两个具有正定矩阵元的算符之差 $\mathcal{O} = \mathcal{O}^{+} - \mathcal{O}^{-}$。因此，上述方法可以分别应用于 $\mathcal{O}^{+}$ 和 $\mathcal{O}^{-}$。

最后，我们提到，使用式(8.48)且 $m = 0$ 时，可以得到所谓的混合平均值 $\mathcal{O}_{\mathrm{MA}}$，它是精确量子平均值的有偏估计：

这种混合平均值的计算对于任何类型的算符都是可能的，而不仅仅是局域算符。对于作用在基态上的算符（即 $\mathcal{O} | \Upsilon_{0} \rangle = \mathcal{O}_{0} | \Upsilon_{0} \rangle$，其中 $\mathcal{O}_{0}$ 是特征值），混合平均估计是精确的（正如我们在基态能量中看到的）。对于所有其他算符，可以使用一种近似方案来评估基态期望值，该方案在状态 $| \Psi_{G} \rangle$ 接近 $| \Upsilon_{0} \rangle$ 时有效。事实上，如果 $| \Psi_{G} \rangle = | \Upsilon_{0} \rangle + \epsilon | \Upsilon^{\prime} \rangle$（其中 $| \Upsilon^{\prime} \rangle$ 是归一化的且与 $\left| \Upsilon_{0} \right\rangle$ 正交），则在 $O ( \epsilon^{2} )$ 阶精度下有以下关系成立：

8.4 连续时间极限#

到目前为止，我们看到了幂法的一种统计实现方式，它基于方程(8.10)和(8.11)的马尔可夫过程。正如我们在第8.1节中讨论的，鉴于GFMC方法的概率性质，格林函数必须是非负的。因此，方程(8.2)中的常数λ必须取足够大的值，以使所有对角元素 $\mathcal{G}_{x , x} \geq 0$。这里，为了简化符号，我们不对各种量加波浪号，假设我们正在使用重要性抽样。对角元素非负的要求通常需要一个非常大的移位λ，该值也随着团簇尺寸的增加而增大；此外，在某些情况下，这个移位无法先验给出（例如，当考虑固定节点近似时，参见第10章）。如果对于选定的 $\lambda$，发现负的对角元素 ${\mathcal G}_{x , x}$，我们需要重新开始一次新的模拟，并采用更大的 $\lambda$ 值，这会浪费大量计算时间。为了避免这个问题，我们可以使用一个非常大的 $\lambda$ 值；然而，这种选择会降低算法的效率，因为停留在同一构型的概率

变得非常接近1，从而意味着非常大的相关时间。

下面，我们将展示，在不损失效率的情况下，可以解决使用大λ的问题。实际上，给定 $p_{d}$，游走者连续k步停留在同一构型 |x⟩（即进行k次对角移动）然后转变到 $\vert x^{\prime} \rangle \neq \vert x \rangle$ 的概率 (k) 简单由下式给出：

正如我们在第3.3节中所述，我们可以通过抽取一个在[0, 1)上均匀分布的随机数 $\xi$ 来对这个离散分布进行抽样；那么对角移动次数 $k$ 由以下不等式确定：

这导致

这些关系等价于取：

其中 $[ \dots ]$ 表示整数部分。这里我们假设 $0 \leq \xi < 1$，这样k总是有限的。

这个想法是将M次基本迭代的集合视为单个马尔可夫步，对于这个步长，$\mathcal{Q} ( k )$ 由方程(8.56)对于 $k = 0 , \ldots , M - 1$ 给出，且 $\mathcal{Q} ( M ) = p_{d}^{M}$（表示在M次尝试中没有非对角移动的可能性）。然后，在M次迭代过程中，我们可以通过下式对 (k) 进行抽样：

其中 $M_{\mathrm{left}}$ 是完成该马尔可夫步所需剩余迭代次数。开始时，$M_{\mathrm{left}} = M$，然后我们可以迭代地应用方程(8.60)并记录 $M_{\mathrm{left}}$。游走者的权重 w 根据k次对角移动进行更新：

如果 $k < M_{\mathrm{left}}$，则根据由下式定义的转移概率 $t_{x^{\prime} , x}$ 随机抽取一个新构型：

最后，如果$k < M_{\mathrm{left}}$，则将$M_{\mathrm{left}}$更新为$M_{\mathrm{left}} - ( k + 1 )$（k次对角移动和一次非对角移动），如此我们可以继续使用式(8.60)，直到$M_{\mathrm{left}} = 0$，此时马尔可夫步骤结束，见图8.2。我们强调，马尔可夫步骤由M个基本步骤组成；因此，在这种情况下，计算式(8.45)中$G_{n}^{p}$所使用的相应因子被替换为沿这M步累积的权重。

有趣的是，这种方法可以轻松推广到$\lambda \to \infty$的情况，只需通过在λ中增加M来实现，即$M = \lambda \beta$，其中$\beta$表示精确传播子$\exp(-\beta \mathcal{H})$的虚时间演化，并以统计方式应用：

图8.2 基于对式(8.60)中对角移动次数进行采样的马尔可夫过程示意图。此处，基本步骤数为 $M = 18$：构型未改变的步骤用叉号标记，而（两个）构型发生改变的步骤用实心点标记。为了计算可观测量，只有马尔可夫步骤结束时（即18个基本过程后）的权重被写入文件，以构建式(8.45)中的累积权重 $G_{n}^{p}$。当 $\lambda \to \infty$ 时，类似情况成立，权重仅存储在由虚时间 $\tau$ 演化描述的马尔可夫步骤结束时。

为此，我们只需记录完成长度为$\beta$的虚时间传播所剩余的对应时间$\beta_{\mathrm{left}}$。当$\lambda \to \infty$时，$p_{d}$变为：

于是，式(8.60)变为：

其中我们利用了$\ln p_{d} \approx [ e_{L} ( x ) - \mathcal{H}_{x , x} ] / \lambda$。在虚时间$\beta$上，我们迭代应用这个公式，并记录完成长度为$\beta$的马尔可夫步骤所剩余的剩余时间$\beta_{\mathrm{left}}$。初始时，$\beta_{\mathrm{left}} = \beta$，然后使用式(8.65)抽取$\tau$；权重根据下式更新：

如果$\tau < \beta_{\mathrm{left}}$，则根据式(8.62)采用新构型，并将剩余时间更新为$\beta_{\mathrm{left}} - \tau$。当$\beta_{\mathrm{left}} = 0$时，马尔可夫步骤完成。

这种方法的一个重要方面是，使用虚时间传播子$\exp ( - \beta \mathcal{H} )$进行的演化不会引入任何Trotter误差（Trotter, 1959；Suzuki, 1976a,b），这在格点上是可以实现的。

8.5 多游走子表述#

在实际应用中，当式(8.45)中的投影次数$p$过大时，前述章节描述的单游走子表述会变得不稳定。实际上，累积权重$G_{n}^{p}$是$p$个随机变量的乘积，因此具有巨大的涨落，如第2.7节所示：

其中 $\sigma^{2}$ 是 ln $b_{x}$ 的方差。（这里我们同样倾向于使用简化记号，不在不同量上添加波浪线，但依然假设我们是在使用重要性抽样。）

为了克服大幅度涨落的问题，我们讨论一种同时传播一组 $N_{\mathrm{w}}$ 个行走子的方法，这些行走子由其权重 $w_{\alpha , n}$ 和构型 $x_{\alpha , n}$ 定义，对应于每个马尔可夫迭代步 n 且 $\alpha = 1 , \ldots , N_{\mathrm{w}}$。给定 $N_{\mathrm{w}}$ 个行走子，我们用一对向量 $( \mathbf{x}_{n} , \mathbf{w}_{n} )$ 表示相应的构型和权重。如果让它们独立演化，对于前述的大涨落没有改进（例如，我们只是将累积权重的方差减小了 $1 / N_{\mathrm{w}}$ 倍）：

相反，通过在权重方差变得过大之前对行走子进行迭代重构，可以在模拟过程中使所有行走子的权重大致相等，从而大幅减少评估 $G_{n}^{p}$ 时每个单独因子出现的统计误差，这意味着：

这使得公式 (8.67) 中的涨落呈指数级减小。因此，通过取 $N_{\mathrm{w}} \approx p$，该过程允许我们对 p 步的投影进行稳定传播。

现在我们来讨论这种被称为分支的重构方案。主方程 (8.38) 可以很容易地推广到多个独立行走子的情形。实际上，如果概率分布 $\mathcal{P}_{n} ( {\mathbf{x}} , {\mathbf{w}} )$ 的演化没有任何限制，每个行走子与其他行走子都不相关，那么我们有：

类似于之前公式 (8.39) 的情况，我们可以将迭代步 n 演化后的状态定义为：

表 8.1. 重构方案。这里 $\overline{{w}}_{n}$ 由式 (8.72) 定义，新漫步者根据式 (8.73) 的概率 $p_{\alpha , n}$ 选取，见图 8.3。

关键点是物理信息蕴含于分布 $\mathcal{P}_{n} ( {\mathbf{x}} , {\mathbf{w}} )$ 的一阶矩（关于权重）中，该矩给出了波函数（例如，含重要性抽样的 $\Psi_{G} ( x ) \Psi_{n} ( x )$）。因此，我们可以定义一个重构过程，在不改变式 (8.71) 统计平均值的前提下改变概率分布。这可以通过一个马尔可夫步骤实现，该步骤将构型 $( \mathbf{x}_{n} , \mathbf{w}_{n} )$ 转变为 $( \mathbf{x}_{n}^{\prime} , \mathbf{w}_{n}^{\prime} )$。特别地，以下重构能够很好地完成这一任务 (Calandra Buonaura and Sorella, 1998)：

1. 取新权重等于所有旧权重的平均值：

2. 以正比于旧构型权重的概率 $p_{\alpha , n}$ 从旧构型中选择新构型：

分支过程如示意图表8.1和图8.3所示。原则上，新构型的选择可以通过抽取 $N_{\mathrm{w}}$ 个独立随机数 $z_{\alpha}$（均匀分布在 [0, 1) 内）并应用第3.5节的策略来完成。然而，一种更好的重构方案是通过抽取 $N_{\mathrm{w}}$ 个相关随机数 $z_{\alpha} = ( \alpha + \xi - 1 ) / N_{\mathrm{w}}$ 来实现，其中 $\xi$ 是单一随机数，均匀分布在 [0, 1) 内（且 $\alpha = 1 , \ldots , {{N_{\mathrm{w}}}}$）。然后，这组数被用来选择新构型，尽管它们相关，但能给出无偏结果。注意，当所有权重 $p_{\alpha , n}$ 相等时，这种改进方案不会改变游走者群体；而在一般情况下，存活游走者的数量是最优的。这次重构之后，新的 $N_{\mathrm{w}}$ 个游走者集合具有相同的权重 $\overline{{w}}_{n}$，并且大多数权重很小且无关紧要的游走者已被剔除，因为它们不会根据 $p_{\alpha , n}$ 被选中。

图8.3 对于 $N_{\mathrm{w}} = 6$ 个游走者的分支过程示意图。这里，每个构型 $x_{\alpha , n}$（用实心点表示）可能被杀死（用叉号表示）或继续其演化（用箭头表示）。在后一种情况下，构型可能根据方程(8.73)的概率产生它的一个或多个副本。重构每隔 $n_{\mathrm{bra}}$ 个马尔可夫步进行一次（此处中间步骤未显示）。前向行走技术也在此展示：在最左边步，应用一个（对角线）算子，其值为 $\alpha = 1 , \ldots , 6$ 的 ${\mathcal{O}}_{\alpha}$（马尔可夫步的索引省略）；在最右边步，前向行走结束（在 $m = 4 \ \mathrm{steps}$ 之后），并且算子的值根据分支过程而改变。

这种重构起到了与传统分支方案相同的稳定效果（Ceperley and Alder, 1980; Reynolds et al., 1982; Foulkes et al., 2001），但优点是游走者数量在模拟过程中保持固定。事实上，传统分支方案的一个问题是控制游走者群体和权重的波动，这可能会在模拟中引入一些偏差。

在方程(8.72)和(8.73)的重构中，游走者之间引入了某种相关性，然而，我们可以严格证明：在概率 $\mathcal{P}_{n} ( {\mathbf{x}} , {\mathbf{w}} )$ 上对方程(8.71)的统计平均等于对新概率 $\mathcal{P}_{n}^{\prime} ( \mathbf{x}^{\prime} , \mathbf{w}^{\prime} )$ 进行的统计平均，即我们有：

事实上，分支过程由下式描述：

其中转移概率 $T_{\mathrm{B}} ( \mathbf{x}^{\prime} , \mathbf{w}^{\prime} | \mathbf{x} , \mathbf{w} )$ 由下式给出：

因此，分支过程后的新波函数由下式给出：

这里，为了方便，我们单独挑出求和中的每一项，即 $w_{\gamma}^{\prime} \delta_{x , x_{\gamma}^{\prime}}$，并对所有其他变量 $x_{\alpha}^{\prime} \left( w_{\alpha}^{\prime} \right)$（其中 $\alpha \neq \gamma$）求和（积分）。此项的贡献（记为 $\left[ \Psi_{G} ( x ) \Psi_{n}^{\prime} ( x ) \right]_{\gamma} $）由下式给出：

然后，通过施加条件 $\begin{array} {r} {w_{\gamma}^{\prime} = 1 / N_{\mathrm{w}} \sum_{\alpha} w_{\alpha} \colon} \end{array}$ 的 delta 函数，之前方程中出现的因子 $\textstyle \sum_{\alpha} w_{\alpha}$ 在对 $w_{\gamma}^{\prime}$ 积分和对 $x_{\gamma}^{\prime}$ 求和后消去：

即 $\left[ \Psi_{G} ( x ) \Psi_{n}^{\prime} ( x ) \right]_{\gamma} = \Psi_{G} ( x ) \Psi_{n} ( x ) / N_{\mathrm{w}}$ 。对所有 $\gamma$ 求和后，我们证明了式 (8.74) 的结论。

实际上，式 (8.75) 的重配置方案每 $n_{\mathrm{bra}}$ 次幂迭代执行一次。我们想强调的是，采用较大的 $n_{\mathrm{bra}}$ 意味着一个相当低效的分支过程，因为游走者的权重差异会很大，只有少数能在重配置中幸存下来（这里的低效性在于，我们花费时间演化那些最终对平均结果无贡献的游走者）。相反，$n_{\mathrm{bra}}$ 过小也不是最优的，因为在这种情况下所有游走者都将不受影响，我们执行的重配置并不会改变分布 $\mathcal{P}_{n} ( {\mathbf{x}} , {\mathbf{w}} )$。

在每次重配置后，存储平均权重 $\overline{{w}}_{n}$ 并将所有游走者 $\alpha = 1 , \ldots , N_{\mathrm{w}}$ 的权重设为 $w_{\alpha , n}^{\prime} = 1$（而不是取 $w_{\alpha , n}^{\prime} = \overline{{w}} )$ 是很有用的。累积权重 $G_{n}^{p}$ 的值可以通过追踪 $p$ 次重配置过程中权重的演化来恢复，其形式为：

这对应于幂法迭代 $p \times n_{\mathrm{bra}}$ 次的应用。

第 8.4 节的连续时间方法也可以应用；在这种情况下，分支方案可以在长度为 $\tau_{\mathrm{bra}}$ 的演化之后进行，即通过 $\exp ( - \tau_{\mathrm{bra}} \mathcal{H} )$ 传播之后进行。

一旦分支过程完成，就无法再获得边际概率的平衡分布的显式形式：

事实上，通过重复与之前相同的步骤，我们可以证明，重配置后新的边际概率由下式给出：

该式依赖于行走子数量 $N_{\mathrm{w}}$，并且当 $N_{\mathrm{w}}$ 增加时（即当第一个括号中的项趋近于常数时，根据中心极限定理，见第 2.7 节），可以恢复方程 (8.77) 的波函数。

8.5.1 基态能量与关联函数#

与单行走子技术的情况类似，基态能量可以通过单次模拟计算得到，其中构型根据 $\Pi_{\mathrm{eq}} ( x )$（未知显式形式）平衡，并考虑方程 (8.80) 的累积权重：

其中 $n$ 表示第 $n$ 次重配置过程；局域能量可以在重配置过程之后立即评估，此时所有行走子具有相同权重：

或者，为了获得更好的统计误差，可以在重配置之前，适当考虑每个行走子的权重来评估：

图8：4 报告了 $N_{\mathrm{w}}$ 取两个值时收敛性的示例，使用的是第 8.4 节的连续时间算法。这里展示了 $6 \times 6$ 团簇（可通过 Lanczos 对角化得到精确解）上方格子 Heisenberg 模型的基态能量情况。注意，$p = 0$ 的起始点对于 $N_{\mathrm{w}} = 10$ 和 50 是不同的，后一种情况更接近精确结果。

图 8.4 $6 \times 6$ 团簇（方格子）上 Heisenberg 模型的每站点能量随投影指数 $p$ 的变化，见方程 (8.83)。显示了 $N_{\mathrm{w}} = 10$（圆圈）和 $N_{\mathrm{w}} = 50$（方块）两种情况。使用 $\tau = 0.1$ 的连续时间方法。虚线表示通过 Lanczos 对角化得到的精确结果。

让我们通过讨论如何计算关联函数来结束这一部分。对于满足方程 (8.47) 条件的局域算符，我们可以使用单条马尔可夫链同时计算方程 (8.48) 的分子和分母，如第 8.3 节中单行走子技术所述。然而，为了控制前向行走技术，在第 $n$ 次迭代中，对于 $\alpha = 1 , \ldots , N_{\mathrm{w}}$ 的一组测量值 $O_{\alpha , n}$，必须在每次前向方向发生的重配置过程之后进行修改。事实上，在具有方程 (8.76) 重配置的马尔可夫过程中，分母可以通过构建累积权重 $G_{n}^{p}$ 轻松计算；而对于分子，分支方案必须进行修正，因为第 $n$ 次迭代中算符的应用会改变后续迭代的重配置，从而修改权重。在实践中，每次分支之后，只需记录对应于重配置中存活的行走子的可观测量值。这可以通过在每次分支后存储每个行走子的算符值 ${\mathcal{O}}_{\alpha , n}$ 以及一个表 $\beta = j_{n} ( \alpha )$ 来实现，该表给出了从中生成新行走子 $\alpha$ 的旧行走子 $\beta$。然后，经过 $m$ 次重配置步骤后的贡献由函数 $j_{n} ( \alpha )$ 的递归应用给出：

必须强调，前行走技术（forward walking technique）在采用少量行走子且考虑较大传播步长时可能不稳定：实际上，只有少数可观测量 ${\mathcal{O}}_{\alpha , n}$ 的值能在前行走中幸存，最终当 $m \to \infty$ 时，仅剩一个值（见图8.3）。因此，该方法仅能在不稳定性出现之前、传播步长足够小的条件下给出可靠结果。

关于随重配置步数收敛性的算例展示在图8.5中，采用第8.4节的连续时间算法，行走子数 $N_{\mathrm{w}} = 50$。此处展示的是正方晶格海森堡模型（仍为 $6 \times 6$ 团簇，可精确求解）实空间中的自旋-自旋关联函数：

其中 $L$ 为格点数，$S_{i}^{z}$ 为格点 i 处自旋算符的 z 分量。

图8.5 $6 \times 6$ 团簇（正方晶格）海森堡模型中公式(8.87)的自旋-自旋关联 $S ( r )$ 随重配置步数 m 的函数关系（见公式(8.86)）。考虑了沿团簇对角线的三个独立距离（从上至下距离递增）。此处使用 $N_{\mathrm{w}} = 50$ 和 $\tau = 0.1$ 的连续时间方法。虚线表示通过兰佐斯对角化获得的精确结果。

8.6 实际实现#

最后，我们简要总结格林函数蒙特卡罗算法在实际实现中的重要步骤。

1. 计算起始时的初始化。#

生成 $N_{\mathrm{w}}$ 个随机配置 $\vert x_{\alpha} \rangle$ 并初始化权重 $w_{\alpha} = 1$。与变分蒙特卡罗方法类似，配置可存储在一组向量 $\mathtt{iconf}(L, N_{\mathrm{w}})$ 中，其元素给出每个行走子 $\alpha = 1, \ldots, N_{\mathrm{w}}$ 在格点 $i = 1, \ldots, L$ 上的局域态；以及 $\mathtt{kel}(2L, N_{\mathrm{w}})$ 中，其非零元素给出每个格点上定义采样配置 |x-（公式(5.44)）的字符串中产生算符的位置。

_• 验证初始配置非奇异，即 $\langle x_{\alpha} | \Psi_{\mathrm{J}} \rangle \neq 0$，与变分蒙特卡罗方法中的讨论类似。

_• 与变分蒙特卡罗方法一样，计算公式(5.33)的表格以及所有需要用重要性采样计算 $\tilde{\mathcal{G}}$ 并进行时间演化的格林函数。对于行列式情况，只需公式(5.80)的静态格林函数；对于普发阶情况，则需要标准格林函数(5.129)以及反常格林函数(5.130)和(5.131)。

2. GFMC投影（采用连续时间方法）。#

_• 对每个行走子 α 独立操作，在与其 $\left| x_{\alpha} \right\rangle$ 通过随机矩阵 $t_{x^{\prime}, x}$（公式(8.62)）连接的所有配置中，选择新配置 $| x_{\alpha}^{\prime} \rangle$。

_• 对每个行走子 α，使用公式(8.66)更新权重因子 $w_{\alpha}$。

_• 更新贾斯特罗因子 (5.34) 的表格以及所有格林函数，即行列式的公式(5.79)或普发阶的公式(5.125)和(5.128)。

_• 经过长度为 $\tau_{\mathrm{bra}}$ 的时间传播后，根据公式(8.72)和(8.73)对行走子进行重配置（亦见表8.1）。

3. 观测量计算#

每次进行分支重配置时（之前或之后），都可以计算观测量。对于基态能量的计算，建议将公式(8.72)的平均权重 $w$ 和公式(8.84)（分支前）或公式(8.85)（分支后）的局域能量写入硬盘，然后利用公式(8.45)计算累积权重，取 $p$ 的最小值。对于局域观测量的计算，建议保存所有组态 iconf $( L , N_{\mathrm{w}} )$ 以及表格 $\beta = j_{n} ( \alpha )$（在每次 $n$ 分支过程中获得），以便执行前向行走子技术所需的公式(8.86)重排操作。

爬行量子蒙特卡洛#

9.1 一种简单的路径积分技术#

在第8章中，我们描述了格林函数蒙特卡洛（GFMC）方法，该方法利用幂法过滤出基态波函数。其中，通过引入由电子组态 |x- 和权重 $w$ 标记的行走子，以统计方式实现投影技术。然后，我们看到可以通过考虑多个行走子并迭代地以近似相同的权重传播它们，来解决权重指数增长的问题。尽管这类蒙特卡洛方法被广泛使用，但当模拟所需的行走子数量过大时，该方法可能变得低效。这通常发生在用于重要性采样的波函数不够精确的情况下。此外，对于大量前向投影，可能需要极其大量的行走子。基于这些原因，引入爬行量子蒙特卡洛（RQMC）算法，作为一种在路径积分表示中简单高效地执行投影方法的方式（Baroni and Moroni, 1998; Moroni and Baroni, 1999）。RQMC的主要优点在于，与变分蒙特卡洛方法类似，它只需要一个不带任何权重的行走子，从而允许我们计算相关函数，即使对于较大的投影时间也是如此。然而，在这种情况下，自相关时间可能会变得非常大。

RQMC方法的目标是定义一个马尔可夫过程，使得我们可以采样：

其中 $| \Psi_{\mathrm{var}} \rangle$ 是给定的初始态（例如最佳的变分拟设），随后对其执行 $N$ 步投影（如同通常情况，$\Lambda$ 是对角算子，满足 $\Lambda_{x , x} = \lambda$，其中 $\lambda$ 是足够大的数，以保证收敛到基态）。需要注意的是，与由公式(8.2)定义的 GFMC 方法不同，这里的投影针对固定的 $N$ 步数进行（然后计算与初始态的重叠）。在 RQMC 中，$\mathcal{Z}$ 扮演着伪配分函数的角色，用于定义经典平均。事实上，$\mathcal{Z}$ 的一个自然表示是利用离散路径积分表示（其中对不同路径的积分被代之以对组态的求和）；然后，物理量将使用这个伪配分函数写成经典平均的形式。

与变分法和 GFMC 方法（见第5章和第8章）类似，我们在多体希尔伯特空间中选定一个完备基 {|x-}：

其中，这些状态被取为正交归一化的。此外，我们假设哈密顿量在该基下的矩阵元 $\mathcal{H}_{x^{\prime} , x} = \langle x^{\prime} | \mathcal{H} | x \rangle$ 对于每个 |x- 和 $| x^{\prime} \rangle$ 都可以高效计算。对于局域哈密顿量，非零矩阵元 $\mathcal{H}_{\boldsymbol{x}^{\prime} , \boldsymbol{x}}$ 的数量相比其总数极小；事实上，对于给定的 |x-，非零矩阵元的数量约为系统大小 L 的量级，而非整个希尔伯特空间的大小。因此，对于给定的状态 |x-，所有矩阵元 $\mathcal{H}_{\boldsymbol{x}^{\prime} , \boldsymbol{x}}$ 都可以在合理的计算量下完成。

通过在所选正交归一基中插入 $N + 1$ 个完备性关系，伪配分函数 $\mathcal{Z}$ 可以写成沿 N 个时间切片的路径积分形式：

这里，权重 $\mathcal{W} ( x_{N} , \ldots , x_{0} )$ 定义为：

其中 $\Psi_{\mathrm{var}} ( x ) ~ = ~ \langle x | \Psi_{\mathrm{var}} \rangle$，并且与 GFMC 方法类似，我们定义了带有重要性抽样的格林函数 $\mathcal{G}_{x^{\prime} , x}$：

此后，为简化记号，我们不再对带有重要性抽样的各量加上波浪号。

与 GFMC 方法类似，为了使得权重 $\mathcal{W} ( x_{N} , \ldots , x_{0} )$ 非负，必须确保格林函数的所有矩阵元 $\mathcal{G}_{x^{\prime} , x}$ 非负。当某对 $( x^{\prime} , x )$ 使得 $\mathcal{G}_{x^{\prime} , x} < 0$ 时，我们就面临符号问题，这将在第10章中讨论。

9.1.1 基态能量与关联函数#

局域算符的基态能量和关联函数可以用伪配分函数 $\mathcal{Z}$ 来表示。例如，基态能量可以（对于足够大的 N）通过下式估计：

其中，与 GFMC 方法一样，局域能量定义为：

通过注意到式 (9.6) 中 $\mathcal{H}$ 与 $\mathcal{G}$ 对易，我们可以得到对基态能量的略微更好估计：

类似地，对于在该基下对角的局域算符，即 $\mathcal{O} | x \rangle = \mathcal{O} ( x ) | x \rangle$，其中 $\mathcal O ( x )$ 是对应于组态 $| x \rangle$ 的本征值，我们有：

因此，量子期望值被写成一个经典概率 $\mathcal{W} ( x_{N} , \ldots , x_{0} )$ 的形式，该概率依赖于由投影过程中从 $x_{0}$ 到 x_N 的所有组态定义的爬虫（或聚合物）：

图9.1展示了爬虫的示意图。

9.2 一种简单的组态采样方法#

为了计算基态能量和关联函数，我们可以考虑一个基于Metropolis算法的马尔可夫过程，用于对经典权重 $\mathcal{W} ( \mathcal{R} )$ 进行采样。我们定义两个基本移动 $\mathcal{R} \mathcal{R}^{\prime}$，将爬虫向右或向左移动，并约定右移（左移）记作 $d = + 1$ $( d = - 1 )$，见图9.1。在最简单的RQMC采样中，每一步随机选择方向 $d$，左右移动的概率相等。我们记 $T ( \mathcal{R}^{\prime} , d^{\prime} | \mathcal{R} , d )$ 为定义Metropolis算法中提议移动的转移概率。由于新方向 $d^{\prime}$ 是随机选择的，左右移动概率相等，与 $d$ 的实际值无关，因此有：

其中 $t^{( d^{\prime} )} ( \mathcal{R}^{\prime} | \mathcal{R} )$ 定义了在给定方向 $d^{\prime}$ 下爬虫 $\mathcal{R} \mathcal{R}^{\prime}$ 的变化。接受概率 $a^{( d^{\prime} )} ( \mathcal{R}^{\prime} | \mathcal{R} )$ 可通过利用细致平衡条件得到：

在该方程中，$d = - d^{\prime}$，因为反向移动只有通过改变传播方向 $d$ 才能实现。然后，在Metropolis算法中，接受概率由下式给出：

其中

在这种方法中，爬虫的概率分布沿马尔可夫迭代 $n$ 按以下方式演化：

其中全转移概率由下式给出：

注意，全转移概率满足细致平衡：

这可以通过使用式(9.11)和(9.12)轻松验证。经过热化时间后，例如对于较大的 $n$，$\mathcal{P}_{n} ( \mathcal{R} , d )$ 收敛到所需的平衡概率 ${\mathcal{P}}_{\mathrm{eq}} ( {\mathcal{R}} ) \equiv{\mathcal{W} ( \mathcal{R} ) / ( 2 \mathcal{Z} )}$，参见第3.8节。

接下来，我们像在GFMC方法中那样分解格林函数 $\mathcal{G}$，即 $\mathcal{G}_{x^{\prime} , x} = p_{x^{\prime} , x} b_{x}$，使得经典权重变为：

图9.1 由从 $x_{0}$ 到 $x_{N}$（此处 $N = 6$）的所有组态定义的爬虫 $\mathcal{R}$ 的示意图。新爬虫 $\mathcal{R}^{\prime}$ 是通过执行左/右移动从旧爬虫生成的，在此过程中添加一个新的头部/尾部，并丢弃旧的尾部/头部。

接下来我们详细讨论这两种右移和左移。首先考虑 $d = + 1$ (右移)。在这种情况下：

其中 $x_{T}$ 是在爬虫最右侧的试探移动。该过程的转移概率为：

新爬虫 $\mathcal{R}^{\prime}$ 对应的权重为：

因此，在实施 Metropolis 算法所需的比值 $\mathcal{W} ( \mathcal{R}^{\prime} ) / \mathcal{W} ( \mathcal{R} )$ 中，几乎所有因子都相互抵消，得到：

其中我们利用了 $\mathcal{G}_{x^{\prime} , x} = p_{x^{\prime} , x} b_{x}$ 是带有重要性抽样的格林函数，且满足 ${\mathcal G}_{y , x} = {\mathcal G}_{x , y} {\Psi}_{\mathrm{var}}^{2} ( y ) / {\Psi}_{\mathrm{var}}^{2} ( x )$ 。相比之下，对于 $d = - 1$ (左移)，我们有：

现在 $x_{T}$ 是在爬虫最左侧的试探移动。该过程的转移概率为：

新爬虫 $\mathcal{R}^{\prime}$ 对应的权重为：

因此，进入 Metropolis 算法的比值 $\mathcal{W} ( \mathcal{R}^{\prime} ) / \mathcal{W} ( \mathcal{R} )$ 具有简单形式：

与之前一样，我们利用了 $\mathcal{G}_{x^{\prime} , x} = p_{x^{\prime} , x} b_{x}$ 是带有重要性抽样的格林函数，且满足 $\mathcal{G}_{y , x} = {\mathcal G}_{x , y} \Psi_{\mathrm{var}}^{2} ( y ) / \Psi_{\mathrm{var}}^{2} ( x)$ 。

现在我们可以简化出现在公式(9.14)的 Metropolis 算法中的项。对于右移，将其逆向操作 $\mathcal{R}^{\prime} \to \mathcal{R}$ 的是一个左移，且 $x_{T} = x_{N}$ ：

这里我们利用了在第一次右移后，爬虫 $\mathcal{R}^{\prime}$ 的最左侧组态是 $x_{N - 1}$ 这一事实。因此，我们得到：

类似地，对于左移，将其逆向操作 $\mathcal{R}^{\prime} \to \mathcal{R}$ 的是一个右移，且 $x_{T} = x_{0}$ ：

这里我们利用了爬虫 $\mathcal{R}^{\prime}$ 的最右侧组态是 $x_{1}$ 这一事实。然后，我们得到：

总之，Eq. (9.28) 和 (9.30) 完整定义了在标准 Metropolis 算法 (9.13) 中接受或拒绝新提议的爬虫 $\mathcal{R}^{\prime}$ 的规则。

9.3 弹跳算法#

前述更新过程的主要问题在于存在长相关时间。实际上，每一步都选择不同的方向会导致爬虫只晃动其头部和尾部（即最左和最右位置的少数构型），而身体部分（即所有内部构型）却保持不变。当考虑大量时间切片 N 时，这个问题尤为严重。弹跳算法通过在一个方向（右或左）上执行多个步骤后再改变方向，实现了较小的相关时间（Pierleoni and Ceperley, 2005）。在实践中，变量 d 不再随机采样，而仅仅在 Eq. (9.13) 中的移动被拒绝时才改变符号。这样一来，细致平衡条件不再满足。然而，马尔可夫链的一般理论（Meyer, 2000）保证了只要满足遍历性，就会存在唯一的平衡分布。接下来，我们将证明弹跳算法以 $\mathcal{W} ( \mathcal{R} )$ 为平衡分布，但不会讨论收敛到该分布的问题。

弹跳算法的完整转移概率为：

其中 $B^{( d )} ( {\mathcal{R}} )$ 可由归一化条件确定：

由此可得：

转移概率 $K_{B} ( \mathcal{R}^{\prime} , d^{\prime} | \mathcal{R} , d )$（不满足细致平衡条件）同时决定了新的爬虫 $\mathcal{R}^{\prime}$ 和方向 $d^{\prime}$ 的新值。例如，如果移动 $\mathcal{R} \to \mathcal{R}^{\prime}$ 被接受（第一项），则 $d^{\prime} = d$；如果移动被拒绝（第二项），则 $d^{\prime} = - d$。

于是，主方程为：

它以 $\mathcal{P}_{\mathrm{eq}} ( \mathcal{R} , d ) \equiv \mathcal{W} ( \mathcal{R} ) / ( 2 \mathcal{Z} )$ 为平稳解。事实上，利用 Eq. (9.31) 中转移概率的显式形式，并令 $\mathcal{P}_{n} ( \mathcal{R} , d ) = \mathcal{W} ( \mathcal{R} ) / ( 2 \mathcal{Z} )$，我们得到：

然后，使用 Eq. (9.12) 并令 $d = - d^{\prime}$，我们有：

最后，利用 Eq. (9.33)，我们得到：

证明了 $\mathcal{W} ( \mathcal{R} ) / ( 2 \mathcal{Z} )$ 是由 $K_{B} ( \mathcal{R}^{\prime} , d^{\prime} | \mathcal{R} , d )$ 确定的随机过程的一个平稳解这一事实。

9.4 连续时间极限#

与 GFMC 方法类似，我们可以考虑格林函数 ${\mathcal{G}} = ( \Lambda - {\mathcal{H}} )$ 中 $\lambda \to \infty$ 的极限，例如，通过考虑精确传播子 $\exp ( - \beta \mathcal{H} )$ 在虚时下的统计演化（见第 8.4 节）。这可以在没有任何 Trotter 误差的情况下完成（Trotter, 1959; Suzuki, 1976a,b）。在 RQMC 方法中，我们可以令 $\beta = N \tau$，并考虑有 N 个时间片的路径积分表示，每个时间片对应一个使用 $\exp ( - \tau \mathcal{H} )$ 的投影操作。然后，之前的更新方案可以很容易地通过传播最左侧（或最右侧）的构型 $\left| x_{N} \right. \left( \left| x_{0} \right. \right)$ 来实现，即使用公式 (8.62) 的跃迁概率进行时间 τ 的传播，直到最终达到 $\vert x_{T} \rangle$。最后，使用 Metropolis 算法来决定接受或拒绝该移动（即整个 τ 时间的传播）；在这种情况下，$r^{( d )} ( \mathcal{R}^{\prime} | \mathcal{R} )$ 仍然由公式 (9.28) 或 (9.30) 给出，但分子和分母中的权重是在相应时间片中累积的，这与公式 (8.66) 类似。

图9：2 展示了一个关于时间片数 N 收敛性的例子。此处，我们展示了在 $6 \times 6$ 团簇（方晶格）上海森堡模型的基态能量（该系统的精确解可以通过 Lanczos 对角化获得）。

图 9.2 $6 \times 6$ 团簇（方晶格）海森堡模型每格点能量随时间片数 N 的变化关系（考虑 $\tau = 0.1$ 的连续时间方法）。虚线表示通过 Lanczos 对角化获得的精确结果。

9.5 实际实现#

最后，我们简要总结一下在实际实现爬行量子蒙特卡洛算法中的重要步骤。

1. 计算开始时的初始化#

_• 为所有时间片生成相同的随机构型 $\left| x_{i} \right.$。这些构型可以存储在一组向量中：iconf$( L , N + 1 )$，其元素给出每个格点 $i = 1 , \ldots , L$ 和时间片上的局域态；以及 kel$( 2 L , N + 1 )$，其非零元素给出每个格点上定义采样构型 $| x \rangle$（见公式 (5.44)）的字符串中产生算符的位置。

_• 验证初始构型非奇异，即 $\langle x_{i} | \Psi_{\mathrm{J}} \rangle \neq 0$，这与变分蒙特卡洛方法中讨论的类似。

_• 与变分蒙特卡洛方法一样，计算用于快速更新 Jastrow 因子和所有必要的格林函数的表格（公式 (5.33)）。对于行列式情况，只需要公式 (5.80) 中的静态格林函数；而对于 Pfaffian 情况，则需要标准格林函数 (5.129) 以及反常格林函数 (5.130) 和 (5.131)。

2. RQMC 投影（采用连续时间方法）#

在给定方向 $d = \pm 1$ 内，使用公式 (8.62) 的跃迁概率，将最左侧（最右侧）构型 $\left| x_{N} \right. \left( \left| x_{0} \right. \right)$ 在虚时 τ 内传播，直到达到试探构型 $| x_{T}^{\prime} \rangle$。类似于公式 (8.66)，累积用于 Metropolis 算法的对应权重，以计算 $r^{( d )}$。接受或拒绝试探移动。如果试探移动被接受，则保持相同的 d 用于下一步；否则，改变方向 $d \to - d$。

_• 每当试探移动被接受时，更新构型、Jastrow 因子表格 (5.34) 以及所有格林函数，即对于行列式情况为公式 (5.79)，或对于 Pfaffian 情况为公式 (5.125) 和 (5.128)。

3. 计算观测量。#

观测量（例如，基态能量）可以每隔 $O ( L )$ 步计算一次。