第1章 · 引言#

量子蒙特卡罗方法（Quantum Monte Carlo method）简单来说就是应用于量子问题的蒙特卡罗方法。量子蒙特卡罗方法与经典蒙特卡罗方法的区别在于，需要先付出努力将量子问题转化为适合蒙特卡罗模拟的形式。正是在进行这种转换时，问题的量子性质才会显现出来——不仅通过物理变量的非对易性（noncommutativity）和波函数对称化或反对称化等显而易见的问题，还通过符号问题（sign problem）等不那么明显的问题表现出来。几乎在所有情况下，这种转换都会用量子自由度替代经典自由度，而蒙特卡罗方法实际上就是应用于这些经典自由度。后续章节将介绍并解释许多已成功应用于各类量子问题的量子蒙特卡罗方法。在第1章和第2章中，我们将重点讨论什么是蒙特卡罗方法以及它为何有用。

1.1 蒙特卡罗方法#

蒙特卡罗方法并非一种具体技术，而是一种解决过于复杂无法解析求解、或数值计算强度太大无法确定性求解的问题的通用策略。通常，一种具体策略会融合多种不同的蒙特卡罗技术。在可能是最早使用"蒙特卡罗"一词的期刊文章中，Metropolis和Ulam（1949）讨论了这一策略。概括地说，

蒙特卡罗方法是一种迭代随机过程（iterative stochastic procedure），该过程与某个函数的定义关系一致，使得可以在不完全确定该函数的情况下对其进行估计。

这与口语化说法"一种使用随机数的方法"大相径庭。让我们逐部分剖析这一定义，一个关键点将浮现出来。

Ulam和Metropolis当时是在介绍研究微分方程和积分微分方程的一种统计方法的动机和总体描述。这些方程就是他们所说的"某个函数的定义关系"，而"函数"则是这些方程的解。这个函数当然先验未知。几年后，Metropolis、Rosenbluth、Rosenbluth、Teller和Teller（1953）提出了一种研究平衡统计力学的统计方法。那里的定义关系是物理量在玻尔兹曼分布上的热力学平均（thermodynamic average）。函数就是物理量，未知量则是它的平均值。Ulam和Metropolis提出并由我们转述的蒙特卡罗方法的一般描述，涵盖了统计力学应用以及比最初的微分和积分微分方程多得多的应用。

该方法是一种"迭代随机过程"，意味着该过程被反复应用。表面上，其中一个目的是产生大量测量值，对这些测量值进行统计分析，并借助大数定律（law of large numbers）加以验证。当然，这种应用正是Ulam和Metropolis所设想的蒙特卡罗方法的用途之一。对于许多其他应用，这个短语的含义更为微妙。从一个步骤到下一个步骤的移动依赖于随机过程，从而创建一个事件链，并在此基础上进行统计研究。通常，这个事件链是一个马尔可夫链（Markov chain）。

随机数通过"随机过程"这一表述而间接引入。除了对第一批计算机的运算速度感到惊讶外，Ulam和Metropolis还惊讶于这种运算可以生成在零到一区间上均匀分布的"随机数"，从而使各种从几乎任何分布中抽样的方案成为可能。

然而，蒙特卡罗策略的决定性特征来自陈述的其余部分：蒙特卡罗方法估计某个函数，而"不完全确定它"。对于Ulam和Metropolis来说，这意味着不必逐点获得满足微分和积分微分方程的函数值。更广泛地说，预见到在统计力学问题中的应用，我们转述说蒙特卡罗方法可以在不需要完全确定函数的情况下进行估计。由于不完全确定函数，蒙特卡罗方法因此只能提供一个估计值。正如我们将要讨论的，这个估计值可以具有高度的确定性。

蒙特卡洛方法并非唯一不需要完全确定函数的方法。确定性方法通常也无需如此。例如，一个标准的数值问题是估计一个多元函数的积分。确定性算法使用离散化方案，通过在每个维度上仅评估 N 个点的被积函数来估计积分。因此，这些方法的计算量随 $N^{D}$ 增长，其中 D 是维度（变量）数量。因此，高维问题的计算量随维度呈指数增长。正如我们将看到的，许多蒙特卡洛方法通过仅需在有限点上评估函数来打破“维数灾难”，而这些点的数量对维度的依赖性很弱。这是关键点。

蒙特卡洛方法是一个强大但相对简单的想法。这一想法的种子源于乌拉姆对纸牌游戏“通关”的挫败感。在游戏中的任何中间时刻，他都因难以计算自己是否会赢或输而感到沮丧，但他意识到，如果他玩得足够多，就能可靠地估计获胜概率。第一台电子计算机的高速计算能力促使他将这个想法发展成一种解决数学问题的新方法。他与冯·诺依曼分享想法后，他们共同开发了一种随机方法来模拟中子通过裂变材料的输运过程。这种输运过程本质上是随机的：由平均自由程定义，中子在特定距离内不发生碰撞的概率是确定的，而在碰撞时，散射截面则规定了碰撞是弹性、非弹性、裂变还是吸收的概率。第一个蒙特卡洛方法是一个随机过程，可在计算机上实现，它模仿了自然界中发生的事件。它通过迭代和随机的方式求解类似玻尔兹曼输运的问题，而从未直接诉诸玻尔兹曼方程。

梅特罗波利斯-乌拉姆论文是在该方法首次应用后几年才发表的。它指出随机程序在解决问题（即使是那些没有随机类比的问题）时的普遍实用性。梅特罗波利斯因创造了“蒙特卡洛”这个名字而受到赞誉（Metropolis, 1985, 1987）。选择这个名字并非是为了将这种方法与摩纳哥这个著名赌场的赌博游戏联系起来，而是梅特罗波利斯对乌拉姆的调侃，因为乌拉姆有个叔叔偶尔会“需要”去赌场。蒙特卡洛方法被视为理论家们偶尔不得不求助的一种程序。但可能未曾预料到的是，我们求助于它的频率竟然如此之高。

1.2 量子蒙特卡洛方法#

梅特罗波利斯和乌拉姆（1949）将预见蒙特卡洛方法可用于解决量子问题的功劳归于费米。据他们所述，费米指出，这种新的计算方式使我们能够通过引入虚时依赖关系，利用变换得到定态薛定谔方程（Schrödinger equation）的基态：

\nabla^{2} \psi \left( x , y , z \right) = \left[ E - V ( x , y , z ) \right] \psi \left( x , y , z \right) $$

通过变换引入虚时依赖：

$$ u \left( x , y , z , \tau \right) = \psi \left( x , y , z \right) e^{- E \tau} , $$

因此 $u ( x , y , z , \tau )$ 满足类扩散方程：

$$ {\frac{\partial u} {\partial \tau}} = \nabla^{2} u - V u . $$

已知该方程可以用一组带权重的粒子进行蒙特卡洛表达，其中每个粒子独立执行随机游走，同时其权重会受到点 $( x , y , z )$ 处 V 值的乘法作用。在这类蒙特卡洛模拟中，粒子的总权重以由特征值 E 控制的速率在（虚）时间上呈指数衰减，而粒子的空间分布则提供了 $\psi ( x , y , z )$（定态薛定谔方程的解）的估计。我们稍后将证明，这个特征对（eigenpair）对应于基态解。事实上，费米的变换是所有基态量子蒙特卡洛方法的起点。然而，必须等计算机发展之后，蒙特卡洛方法才能应用于非平凡的量子问题。

费米的观察等价于对含时薛定谔方程进行变换：

$$ i \frac{\partial \psi} {\partial t} = - \nabla^{2} \psi + V \psi\tag{1.1} $$

通过变量替换 $i t = \tau$，从实时间变换到虚时间：

$$ {\frac{\partial \psi} {\partial \tau}} = \nabla^{2} \psi - V \psi .\tag{1.2} $$

正如费曼后来指出的 (Feynman and Hibbs, 1965)，相同的变量替换（在虚时间中采用周期性边界条件）将量子力学的实时路径积分形式变换为有限温度量子统计力学。有限温度量子蒙特卡罗模拟通常在虚时间中进行。我们将在后续详细阐述这一观点。眼下我们有一个更具体、更直接的任务——阐述蒙特卡罗方法及其依赖它的必要性。

1.3 经典蒙特卡罗#

为了便于解释为什么有时需要蒙特卡罗模拟，我们从外磁场中的一维伊辛模型（Ising model）入手讨论。该模型并非量子模型，但很可能是具有许多相互作用自由度的最简单模型。

该模型的能量由下式给出：

$$ E = - \sum_{i} \left( J s_{i} s_{i + 1} + H s_{i} \right) ,\tag{1.3} $$

其中 $J > 0$ 是（铁磁）交换积分，i 标识晶格位点，局域伊辛变量 $s_{i}$ 取值为 $\pm 1$，而 H 是外磁场。晶格共有 N 个位点，并假设采用周期性边界条件；即，假设 $s_{i} = s_{i + N}$。

对于一个相互作用多体问题，通常会提出两个问题：随着 N 趋近热力学极限，其零温性质是什么？其有限温度性质又是什么？对于外磁场中的伊辛模型，零温态，即基态（ground state），其能量为 $- N ( J \pm H )$，所有伊辛变量取 $+ 1 \ \mathrm{o r -} 1$，具体取决于 H 的符号。在零场下，基态是二重简并的，自旋全正或全负时能量 -NJ 最低。知道了能量和状态就足以确定所有基态性质。因此，我们重点关注有限温度性质。

对于有限温度，一维伊辛模型的配分函数（partition function）$Z$ 是已知的 (Thompson, 1972)：

$$ Z = \lambda_{+}^{N} + \lambda_{-}^{N} ,\tag{1.4} $$

其中

$$ \lambda_{\pm} = e^{J / k T} \cosh ( H / k T ) \pm \left[ e^{2 J / k T} \sinh^{2} ( H / k T ) + e^{- 2 J / k T} \right]^{1 / 2} ,\tag{1.5} $$

因此自由能 $F = - k T$ ln Z 也是已知的。通过微分自由能，可以得到模型的所有平衡态热力学性质。

通常，感兴趣的有限温度性质是状态方程，重点关注系统在相变（phase transitions）附近的行为。对于相互作用自旋模型，出现磁性长程有序（long-range order）是可能的相变。统计力学告诉我们，要研究相变，需要将系统外推到热力学极限，并在自由能或导出的热力学函数中寻找奇异行为。

本模型的热力学极限很容易得到。由(1.4)和(1.5)式，我们发现每个位点的自由能为：

$$ F / N = - k T \ln \lambda_{+} , $$

每个位点的磁化强度（magnetization）为：

$$ M / N = - {\frac{\partial} {\partial H}} \left( F / N \right) = \sinh ( H / k T ) \bigl [ \sinh^{2} ( H / k T ) + e^{- 4 J / k T} \bigr ]^{- 1 / 2} . $$

在零场下，我们看到所有非零温度下的磁化强度均为零。这排除了模型自发磁化的可能性，这一特性将一维伊辛模型与更高维度的伊辛模型区分开来。

能够写出配分函数的精确表达式就相当于该模型已被解决。方程(1.4)和(1.5)是通过转移矩阵方法（transfer matrix technique）得到的（Thompson, 1972）。虽然该方法也能得到其他模型的精确解，例如零场二维伊辛模型，但统计力学通常需要提供良好的近似或优秀的数值模拟。现在我们来讨论高效蒙特卡罗模拟伊辛模型面临的基本挑战。

对于蒙特卡罗方法，热力学最便捷的处理方式是根据微观变量可能构型的玻尔兹曼概率（Boltzmann probability）计算平均值。如果 (X) 代表与伊辛变量相关的某个物理量，我们可以通过下式计算温度依赖的平均值（热期望值，thermal expectation value）：

$$ \langle X ( T ) \rangle = {\frac{\displaystyle{\sum_{C} X ( C ) e^{- E ( C ) / k T}}} {\displaystyle{\sum_{C} e^{- E ( C ) / k T}}}} = {\sum_{C} X ( C ) {\frac{e^{- E ( C ) / k T}} {Z}}} .\tag{1.6} $$

这里，(C) 代表伊辛变量的一个构型；即 (C = ( s_{1} , s_{2} , \ldots , s_{N} ))，(X(C)) 和 (E(C)) 分别为该构型下 (X) 和能量 (E) 的值，求和范围是所有可能构型的集合。指数项 (e^{- E ( C ) / k T}) 是构型 (C) 的玻尔兹曼因子，配分函数 (\begin{array} {r} {Z = \sum_{C} e^{- E ( C ) / k T}} \end{array}) 是其归一化常数。因此我们可以将(1.6)改写为：

$$ \langle X ( T ) \rangle = \sum_{C} X ( C ) p ( C ) ,\tag{1.7} $$

其中

$$ p ( C ) = {\frac{e^{- E ( C ) / k T}} {Z}}\tag{1.8} $$

是构型 (C) 的概率，且 (\begin{array} {r} {\sum_{C} p ( C ) = 1} \end{array})

在实际应用中，(X) 可以代表能量：

$$ \langle E ( T ) \rangle = \sum_{C} E ( C ) p ( C ) , $$

磁化强度：

$$ \langle M ( T ) \rangle = \sum_{C} \bigg ( \sum_{i} s_{i} ( C ) \bigg ) p ( C ) , $$

或者自旋-自旋关联函数（spin-spin correlation function）：

$$ \left. s_{i} s_{j} ( T ) \right. = \sum_{C} s_{i} ( C ) s_{j} ( C ) p ( C ) . $$

以此方式计算平均值的核心任务在于生成所有构型 (C) 并计算每个构型的概率 (p ( C ))。为了完成后者，我们必须计算配分函数，这就需要计算所有的玻尔兹曼因子。对于我们的 N 位点模型，可能构型的数量为 (2^{N})，因此这项任务的复杂性呈指数增长，很快就使配分函数的计算变得不可行。我们需要一种方法，能在不完全求解问题的情况下，提供可靠的解决方案。我们需要蒙特卡罗方法。

针对伊辛模型这类问题，蒙特卡罗方法并不生成所有构型，而是仅以正确的概率选取 (M) 个代表性构型。它将(1.7)替换为：

$$ \langle X \rangle_{\mathrm{MC}} = \frac{1} {M} \sum_{i = 1}^{M} X ( C_{i} ) .\tag{1.9} $$

如果所有 (X ( C_{i} )) 的测量值在统计上独立，那么蒙特卡罗方法还可以提供与该平均值相关的统计误差 (\sigma_{X}) 的估计值：

$$ \sigma_{X} \approx{\sqrt{{\frac{1} {M - 1}} \left( {\frac{1} {M}} \sum_{i = 1}^{M} X^{2} ( C_{i} ) - \left( {\frac{1} {M}} \sum_{i = 1}^{M} X ( C_{i} ) \right)^{2} \right)}} .\tag{1.10} $$

理想情况下，$M \ll 2^{N}$，即我们在不完全求解问题的情况下估计平均值。

我们现在开始定义具体的蒙特卡洛方法，这些方法允许我们生成（1.9）式中出现的位形 $C_{i}$。我们需要回答的一个基本问题是，如果我们无法在合理的计算机时间内生成计算配分函数所需的所有位形，蒙特卡洛方法如何以正确的概率选择位形？¹ 正如我们将看到的，存在一些方法可以从非归一化的概率密度中选出具有正确概率的位形。最著名的方法称为Metropolis算法，这是我们在前面提到的由Metropolis、Rosenbluth、Rosenbluth、Teller和Teller完成的成果。我们在第2.5.1节讨论该算法。在讨论了这个算法和其他多个蒙特卡洛主题后，我们将在第4章和第5章回到伊辛模型。在后一章中，我们将其量子化，然后开发第一个用于模拟量子自旋模型的量子蒙特卡洛算法。

总体而言，我们关注格点上的量子多体模型的马尔可夫链蒙特卡洛算法。并非所有蒙特卡洛算法都是马尔可夫链算法。有些算法只是盲目地掷飞镖。马尔可夫链的概念在第2章中定义和讨论。它产生随机的事件序列，我们对其执行统计分析。格点模型是指诸如用于相互作用量子自旋的海森堡模型（Heisenberg model）和用于相互作用电子的哈伯德模型（Hubbard model）等模型，其中自旋和电子仅存在于格点位置。用于格点上相互作用电子和玻色子的零温算法与用于连续介质中相互作用电子和玻色子的零温算法非常相似。对于格点模型，存在更多种类的有限温算法。

第1章 · 引言#

1.1 蒙特卡罗方法#

1.2 量子蒙特卡洛方法#

1.3 经典蒙特卡罗#

推荐阅读#