第9章 · 变分蒙特卡洛#

9.2.3 Backflow wave functions#

Backflow transformations introduce a many-body effect into the orbitals themselves. The idea is to replace the coordinate of each particle by a “collective” coordinate that depends on the positions of all other particles:

The resulting determinant is built from orbitals of these backflow coordinates, and the wave function includes more correlation beyond the Jastrow factor.

9.3 Optimization techniques#

The goal in variational Monte Carlo is to minimize $E_{T}$ or, equivalently, to bring the trial state as close to the ground state as possible. We discuss three common optimization methods.

9.3.1 Stochastic reconfiguration#

The stochastic reconfiguration method uses Monte Carlo samples of the derivative of the logarithm of the trial function with respect to parameters to compute the gradient of the energy. The steepest descent or conjugate gradient method then updates the parameters.

9.3.2 Linear method#

The linear method expands the trial state $\psi_{T}$ around a reference parameter set $\boldsymbol{p}_{0}$, and solves a generalized eigenvalue equation in the space of the trial function and its derivatives:

where $| O_{0} \rangle = | \psi_{T} \rangle$ and $| O_{i} \rangle = \partial | \psi_{T} \rangle / \partial p_{i}$. The ground state solution of the generalized eigenvalue problem $H \boldsymbol{c} = E S \boldsymbol{c}$ gives an improved estimate.

9.3.3 Correlated sampling#

When comparing different parameter sets, correlated sampling uses the same set of configurations to compute energy differences. This reduces the statistical noise in the energy and the gradient.

9.4 9.4 小结#

变分蒙特卡洛方法将灵活的试探态与蒙特卡洛采样相结合来计算期望值。它给出了基态能量的上界，可用于研究其他量子蒙特卡洛方法会遭遇符号问题的大系统。关键要素是试探波函数。现代方法通常组合斯莱特行列式、贾斯特罗因子、回流变换以及其他多体关联器。参数的优化依赖于利用局域能量涨落来引导最小化的随机或确定性方法。

该方法不仅是计算基态性质的工具，也是后续章节将要讨论的扩散蒙特卡洛和格林函数蒙特卡洛等更高级技术的基础。

或者

这是一个相当强大的结果，因为任何可允许的试探函数的 Rayleigh 商都必须满足这些界限。

从Rayleigh商可以得出几个值得注意的事实。首先，如果试探态的精度为 $\mathcal{O} ( \varepsilon )$，其中 $\varepsilon$ 很小，那么变分能量 $E_{T}$ 的精度为 $\mathcal{O} ( \varepsilon^{2} )$，即精度更高。因此，使 $\left| {\psi_{T}} \right.$ 尽可能接近 $| \psi_{0} \rangle$ 可以使 $E_{T}$ 更接近 $E_{0}$。然而，硬币的另一面是，$E_{T}$ 接近 $E_{0}$ 并不一定意味着 $\vert \psi_{T} \rangle$ 能如我们所愿地准确表示 $\left| {\psi_{0}} \right.$ 的物理性质。

第二个事实，或者更准确地说是一组相关事实，与试探波函数参数的最优能量收敛有关。如果试探波函数的参数集为 $p = ( p_{1} , p_{2} , \ldots , p_{N_{\mathrm{opt}}} )$，那么我们期望增加参数数量会得到一个序列 $E_{0} \ \leq \ \cdots \ \leq \ E_{T} ( p_{1} , p_{2} , p_{3} ) \leq E_{T} ( p_{1} , p_{2} ) \leq E_{T} ( p_{1} )$。例如，如果基态在相空间的有限区域内只有一个峰，那么我们可以预期，一个由高斯的位置和宽度参数化的试探态会比仅由高斯的位置或宽度参数化的试探态得到更低的能量。如果只有一个参数，我们需要对另一个参数做出良好的猜测，才能达到使用两个参数的效果。

我们可以证明试探态随参数数量增加而系统收敛的一种情况是，我们将试探态表示为 n 个基态的线性组合。任何张成给定哈密顿量的希尔伯特空间的基都适用。如果这个空间的维度是 N，但我们构造一个试探态 $\begin{array} {r} {| \psi_{T} \rangle = \sum_{i = 1}^{n} p_{i} | \phi_{i} \rangle} \end{array}$，其中 $\{| \phi_{i} \rangle \}_{i = 1}^{N}$ 是希尔伯特空间中的某个基（不一定正交），那么对于 $n \textless N$，我们可以证明，使 Rayleigh 商最小化的参数集 $\{p_{i} \}_{i = 1}^{n}$ 满足广义特征值问题

其中，哈密顿量矩阵 H 和交叠矩阵 S 定义为

交叠矩阵是正定的。

如果基是正交归一化的，那么交叠矩阵就是单位矩阵，我们便得到标准的特征值方程。我们本可以通过在约束 $\textstyle \sum_{i} p_{i}^{2} = 1$ 下最小化二次型 $\begin{array} {r} {Q = \sum_{i j} H_{i j} p_{i} p_{j}} \end{array}$ 来得到相同的方程。设 $\lambda_{1} \ \leq \ \lambda_{2} \ \leq \ \cdot \ \cdot \ \leq \ \lambda_{n}$ 为 Q 的特征值。接下来，假设我们通过 j 个线性齐次约束将 $Q$ 简化为一个二次型 $\bar{Q}$，其特征值为 $\bar{\lambda}_{1} \leq \bar{\lambda}_{2} \leq \cdot \cdot \cdot \leq \bar{\lambda}_{n - j}$。那么我们有（Courant and Hilbert, 1965）

特别地，如果我们取 $j = 1$ 且 $p_{n} = 0$，那么二次型 $Q$ 变为其第一个主要子式的二次型，因此我们得到：第一个主要子式的第 i 个特征值至少等于原二次型的第 i 个特征值，并且至多等于原二次型的第 $( i + 1 )$ 个特征值。如果我们重复此过程，将 $( n - 1 )$ 阶子式简化为 $( n - 2 )$ 阶子式，也会得到类似的特征值交错结果。这些交错结果的重要性在于，在 Rayleigh 变分框架内，我们需要使用至少阶数为 n 的基，才能对第 n 个精确特征值进行界定。如果我们将基从 n 增加到 $n + 1$，那么有

对于一般相邻的一对精确特征值也是如此。如果我们不断增加基函数的数量，当 $n N$ 时，我们将收敛到精确结果。

在第9.2节中，我们将讨论多体尝试态选择的各种常见选项。在某些情况下，尝试态是一个包含所有参数的单一态。在其他情况下，它是基态的线性组合。还有一些情况，它可能是前两种情况的混合。在第9.3节中，我们将介绍一种针对混合情况进行最小化的数值技术。接下来，我们将讨论在变分法中使用蒙特卡洛方法的基础知识。

9.1.2 蒙特卡洛采样#

我们现在讨论如何通过蒙特卡洛方法来评估(9.4)。使用一个完备的正交归一基 $\{|C\rangle\}$，满足 $\sum_{C} | C \rangle \langle C | = I$ 和 $\langle C | C^{\prime} \rangle = \delta_{C C^{\prime}}$，我们称之为构型基，则(9.4)变为

我们将其重写为

其中

$E ( C )$ 被称为构型能量。从 $p ( C )$ 中采样构型使我们能够对变分能量应用标准的蒙特卡洛估计：

我们可以使用Metropolis-Hastings算法(2.28)或类似算法从 $p ( C )$ 中采样构型 $C_{i}$。这些算法需要用户定义的提议概率 $T ( C^{\prime} | C )$，用于将态 $| C \rangle$ 改变为 $| C^{\prime} \rangle$。我们回忆一下，Metropolis-Hastings算法然后确定所提议变化的接受概率 $A ( C^{\prime} | C )$ 为

利用相同的样本，我们可以通过下式估计其他可观测量：

其中

一般来说，这些其他估计值并非变分界限。

虽然能量变分泛函(9.1)是优化尝试态的自然起点，但其他选择也是可能的，有时会导致更方便和更有效的算法。最常见的替代方法使用与 H 和 $\left| {\psi_{T}} \right.$ 相关联的能量方差的定义（Umrigar et al., 1988），

该泛函实际上与我们用于推导能量变分泛函(9.2)的泛函相同。它明确指出了这样一个事实：随着 $| \psi_{T} \rangle$ 收敛到基态，与瑞利商相关的能量方差收敛到零。这一事实被称为零方差原理。¹ 零方差原理也有一个蒙特卡洛含义：随着 $\mid \psi_{T} \rangle$ 收敛到基态，构型能量 $E ( C )$ 收敛到基态能量 $E_{0}$，从而变得与构型无关。如果它与构型无关，那么它的方差为零。

零是(9.8)的一个下界，并且它具有已知的下界，这是基于方差泛函的变分计算与基于能量泛函(9.1)的变分计算之间的一个区别特征。虽然能量泛函存在下界，但它是未知的。

我们可以很容易地证明

其中

这意味着方差计算需要与基于（9.4）的方法相同的能量估计子和采样，这并不令人意外，因为瑞利商（Rayleigh quotient）正是使能量方差最小化的能量估计子。形式上，当 $| \psi_{T} \rangle \to | \psi_{0} \rangle$ 时，两个泛函给出相同的 $E_{T}$ 和 $\left| {\psi_{T}} \right.$，即它们都得到 $E_{0}$ 和 $\left| {\psi_{0}} \right.$。但在实践中，两者的计算效率可能相差很大。

理想情况下，变分蒙特卡洛方法会使用一个现成的 $\left| {\psi_{T}} \right.$，该态应尽可能好地逼近 $| \psi_{0} \rangle$。运行模拟后，任务就完成了。通常，一个好的猜测态会包含一组可调参数。如果参数数量较少，可以通过图形化方法轻松找到最小值。如果参数数量较多，则需要结合蒙特卡洛方法和确定性方法。在这种结合中，能量的估计与波函数的优化同时进行；也就是说，我们不是先优化波函数再估计能量。

在下一节中，我们将讨论在晶格模型的变分蒙特卡洛模拟中已被证明有用的不同类别的尝试态（trial states）。然后，在第9.3节中，我们将讨论优化尝试态的过程。随着尝试态参数化变得越来越复杂，这一过程也会变得更加有用。

9.2 尝试态（Trial states）#

尝试态的使用赋予了变分蒙特卡洛方法（就此而言，任何变分计算）一种“种瓜得瓜，种豆得豆”的特性。这种特性表现在两个方面。首先，如果我们知道基态的对称性，并构造出与之匹配的尝试态，那么通过增加尝试态的参数化程度，我们可以使变分能量尽可能接近基态能量。能有多接近取决于我们在构造尝试态上投入了多少精力。对于量子晶格问题，我们通常只知道哈密顿量的对称性，

并关心这些对称性中是否有一个会发生破缺。然而，变分计算本身不会破缺对称性。

变分计算的一个令人困扰的问题是，即使基态的不同试探波函数（Ansätze）可能是正交的，它们的变分能量也可能几乎相同。那么哪一个才是正确的态呢？不幸的是，即使能量相差很大，我们也不能断定能量较低的态一定代表基态。为基态选择合适的试探波函数需要额外的信息。通常，我们会构造尝试态来包含我们对基态的一切确切了解，然后进行试验。这些确切信息可能仅仅是基态是实函数且平方可积。

从事连续介质中原子、分子和固体电子结构计算的研究人员，确信其基态是自旋单态（spin singlet），他们通过开发强大的程序来优化包含大量参数的尝试态，从而精细地打磨了变分蒙特卡洛方法。虽然晶格上的变分蒙特卡洛模拟与连续介质中的模拟面临许多共同的挑战，但前者还需要计算比基态能量更多的东西，这增加了额外的复杂性。探测物质的新奇相态是一种需求；识别可能的量子相变是另一种需求。满足这两种需求都需要非常精确地计算关联函数。此外，除了电子性质模型外，还需要考虑量子自旋和玻色子的晶格模型。因此，晶格模拟需要更多类别的尝试态。在电子连续介质中，直到最近，一直使用的是斯莱特-贾斯特罗（Slater-Jastrow）尝试态，但正如我们将要看到的，这种情况正在发生变化。我们以斯莱特-贾斯特罗类作为比较和对比的参照点来开始讨论尝试态，然后讨论几类适用于晶格模型的尝试态。其中一些函数的形式类似于斯莱特-贾斯特罗构造，但另一些则引出了相当不同的概念。

9.2.1 斯莱特-贾斯特罗态（Slater-Jastrow states）#

对原子、分子和固体电子结构进行的变分蒙特卡洛计算通常使用斯莱特-贾斯特罗尝试态。在第一量子化表示中，该状态具有如下形式：

其中 $R$ 代表 $N$ 个电子的空间坐标 $( r_{1} , r_{2} , \ldots , r_{N} )$，$J ( R )$ 是贾斯特罗函数（Jastrow function），而 $D ( R )$ 是一个斯莱特行列式（Slater determinant）。一种常见的贾斯特罗函数形式为：

其中 $r_{i j} = | r_{i} - r_{j} |$，尽管目前也在使用更复杂的形式。斯莱特行列式是自旋向上行列式和自旋向下行列式的乘积：

其中 $\psi_{n}^{\uparrow}$ 和 $\psi_{m}^{\downarrow}$ 分别是自旋向上和自旋向下电子的单粒子轨道。行列式参数中的 $\psi_{n}^{\sigma} ( r_{i \sigma} )$ 表示由这些轨道构建的一个 $N_{\sigma} \times N_{\sigma}$ 矩阵。更一般地，

在特定的应用中，我们固定斯莱特-贾斯特罗参数，使得尝试态满足平移和旋转对称性。一个基本的对称性是泡利不相容原理（Pauli exclusion principle）所要求的，关于同自旋电子交换的坐标反对称化。斯莱特行列式显然遵守这一点，但为了使完整的尝试态也如此，贾斯特罗因子必须在此类交换下是对称的（玻色子型的）。

行列式具有多电子问题单电子近似的形式。我们尝试使用手头物理系统的最佳单粒子轨道来构建它们。原则上，使用相对任意的函数形式是可能的。实践中，我们使用某个平均场理论的解。平均场理论，例如哈特里-福克（Hartree-Fock）理论，无法充分捕捉电子关联。贾斯特罗因子的目的正是弥补这一缺陷。例如，贾斯特罗因子的一个重要用途是捕捉由电子与原子核相互作用引起的波函数尖点（cusp）处的关联能。

通常，假设的函数形式（如式 (9.9)）的参数决定了贾斯特罗函数，而斯莱特行列式包含了对单粒子轨道进行建模的参数。$E_{T}$ 和 $\psi_{T}$ 都继承了对贾斯特罗因子和斯莱特行列式参数的非线性依赖关系。如果我们使用多态行列式试探 (9.10)，唯一表面上线性的依赖关系是对 $d_{k}$ 系数的依赖。

9.2.2 古兹维勒投影态#

哈伯德模型#

对于电子晶格问题，古兹维勒尝试态（Gutzwiller trial state）的形式与斯莱特-贾斯特罗尝试态类似。古兹维勒为了配合他提出的、我们现在称为哈伯德模型（Hubbard model）的铁磁性多电子模型而提出了这个状态（Gutzwiller, 1965），

古兹维勒尝试态的形式为：

其中投影算符 $P_{G}$ 定义为：

这里，$g = \exp{( - \eta )}$，而$\left| {\psi_{0}} \right.$是基态（ground state）的某种估计。通过调整 $g$，我们可以调整双占据对$| \psi_{0} \rangle$中基态能量估计的贡献。根据$P_{G}$的定义，我们看到当$g = 1$时，它不投影掉任何双占据，但当$g = 0$时，它投影掉全部双占据。原则上，我们可以进行替换$g g_{i}$。在实践中，这很少做。

我们经常使用某种哈特里-福克（Hartree-Fock）近似来指导选择定义$\left| {\psi_{0}} \right.$的轨道。在哈特里-福克中，我们将势能项（这些项在产生和湮灭算符中是四阶的）替换为二阶的项。例如，我们可以将（9.11）近似为

一个博戈留波夫变换（Bogoliubov transformation，见 Blaizot and Ripka, 1986）通过定义旧算符的线性组合作为新的产生和湮灭算符，形式上将新的二次型对角化。然而，新的二次型包含未知的标量量，比如上述公式中的$\langle n_{i \sigma} \rangle$。在这种近似的确定性用法中，这些标量构成我们自洽确定的参数集。而在变分蒙特卡洛计算中，我们不是自洽地选择它们的值，而是通过调整它们来最小化能量。

我们经常做出各种假设来减少参数数量。例如，如果我们假设（9.14）中的$\langle n_{i \sigma} \rangle$与 $i$ 和 $\sigma$ 无关，我们就得到了限制性哈特里-福克近似，它只有单个未知量 $\langle n \rangle$。在这种情况下，轨道是平面波，对于$| \psi_{0} \rangle$我们选择费米海（Fermi sea）$| F S \rangle$。为了构造它，我们定义产生算符$\alpha_{k \sigma}^{\dagger}$，用于将自旋为σ的电子放入波数为 $k$ 的平面波态。这些算符只是原始算符的傅里叶变换：

费米海是

在限制性哈特里-福克中，轨道与 $\langle n \rangle$ 无关。为了调整 $\langle n \rangle$，我们只需要截断（9.15）中的乘积，使得被占据的最低能量态的数量等于电子数。那么我们的变分态就只有 $g$ 作为可调参数。在这种情况下，变分蒙特卡洛的能力仅限于无需进一步近似地估计能量，例如避免了常用的古兹维勒近似（Gutzwiller approximation，见 Gutzwiller, 1965; Fazekas, 1999）。

使 $\langle n_{i \sigma} \rangle$ 与 $i$ 无关，留下 $\langle n_{\sigma} \rangle$ 未知，就得到了非限制性哈特里-福克近似。非限制性哈特里-福克近似对于研究磁性基态很有用。这种态具有明确的波数 $K$，表示具有长程磁序的铁磁、反铁磁或非公度态。在许多方面，自旋密度波（spin-density wave）图像更有用。例如，对于反铁磁自旋密度波，我们取 $\left| {\psi_{0}} \right.$（Yokoyama and Shiba, 1987b）为

其中

并且

且 $E_{k} = \sqrt{\varepsilon_{k}^{2} + \Delta^{2}}$（$\Delta$ 是一个变分参数）。对于二维正方晶格，$K \stackrel{\cdot}{=} ( \pi , \pi )$

当我们期待超导配对关联时，我们从另一种类型的哈特里-福克近似中获得指导。这里，我们采用超导的BCS理论，并取 $\left| {\psi_{0}} \right.$

其中

且 $E_{k} = \sqrt{\varepsilon_{k}^{2} + \Delta_{k}^{2}}$ 。通常对 $\Delta_{k}$ 的参数化做得很简单，其 $k$ 依赖性反映特定的分波。例如，对于二维中的 $d$ 波对称性，$\Delta_{k} = \Delta ( \cos k_{x} - \cos k_{y} )$

我们强调，在(9.18)和(9.20)中，$\Delta$ 都是一个变分参数，而非像哈特里-福克理论中那样的序参量。在后两个试控态中，$\Delta$ 是一个可通过调整来有望捕获基态物理并使模拟更高效的参数。与哈特里-福克不同，$\Delta$ 的非零值并不代表长程有序。长程有序由代表序参量的算符的非零期望值来标志。

除了序参量之外，变分蒙特卡洛方法还可用于估计一系列其他物理量。例如，可以估计准粒子权重和色散（Paramekanti 等，2001；Nave 等，2006；Yang 等，2007）、费米面形变（Himeda and Ogata，2000）、条纹相（Himeda 等，2002）以及谱函数（Paramekanti 等，2001；Yunoki 等，2005）。Edegger 等人（2007）提供了相关工作的简短综述。

除了包含必要的物理之外，试控波函数还应易于采样，并易于用于估算构型能量和其他可观测量。采样通常使用 Metropolis 算法（Shiba，1986；Yokoyama and Shiba，1987a）。一个构型态 $|C\rangle$ 是状态 $| R \rangle =$ $| r_{1 \uparrow} , r_{1 \downarrow} , r_{2 \uparrow} , . . . , r_{N \downarrow} \rangle ~ = ~ c_{r_{1 \uparrow}}^{\dag} c_{r_{1 \downarrow}}^{\dag} \cdot \cdot \cdot c_{r_{N \downarrow}}^{\dag} | 0 \rangle$；也就是说，它由电子位置和自旋参数化。我们接受或拒绝一个提议的（位置或自旋）变化，取决于是否满足

其中 $\zeta$ 是 [0, 1] 上的均匀随机数。对于哈伯德模型，我们几乎总是使用

这暗示了我们提议的构型变化类型。例如，我们可以随机选择一个电子位置，然后等概率地选择是移动该电子还是翻转其自旋。我们以等概率将该电子移动到近邻格点上。如果我们的试控态是 (9.12)，那么（Yokoyama and Shiba，1987a）

其中 $\Delta N_{D}$ 是双占据的变化量，并且如果我们将自旋为 $\sigma$ 的电子从格点 $r_{j}$ 移动到 $r_{j} + a_{\mathrm{\ell}}$，我们有

我们按类似于(9.7)的方式评估其他期望值。

计算量大的步骤是计算行列式的比值。如果没有一个技巧能将成本降低至少等于定义行列式的矩阵阶数倍，这将是计算上不可接受的。这个技巧与我们在第7章讨论行列式方法的更新过程中使用的技巧相同。（另见 Hammond et al. (1994)，附录B。）简而言之，计算并存储矩阵 $D_{\sigma}$ 的逆。然后，

因为我们知道 $D_{\sigma}^{- 1}$ 且蒙特卡罗移动是局域的，$\delta D_{\sigma}$ 只有一行和一列非零。因此上述行列式计算毫不费力。然后，如果移动被接受，我们使用附录G中给出的 Sherman-Morrison-Woodbury 公式（Press et al., 2007，第2章）来更新 $D_{\sigma}^{- 1}$。

如(9.19)所述，$|BCS\rangle$ 超出了刚才描述的算法，因为它缺乏固定的电子数，因此无法表示为斯莱特行列式。形式上，我们通过将试波函数投影到固定电子数来修正这种情况：

更仔细地检查这个问题，我们注意到

其中 $\phi_{k} = \nu_{k} / u_{k}$。展开指数并只保留 $P= N/2$ 次幂（即配对的数量，Bouchoud et al., 1988），我们得到

其中

如果 $\phi_{k} = \exp ( i k \ell )$，那么

这被解释为在方向 $\ell$ 上的一系列价键的和。

两个电子可以配对成单重态或三重态。取 $\phi ( r )$ 为 $\phi_{k}$ 的傅里叶变换，对于单重态和 $S_{z} = 0$ 的三重态，我们可以证明^{2}

其中 $\phi ( r_{i \uparrow} - r_{j \downarrow} )$ 是一个 $P \times P$ 矩阵的元素（Bouchoud et al., 1988; Gros, 1988）。

态

其中 $\begin{array} {r} {P_{G}^{0} = \prod_{i} ( 1 - n_{i \uparrow} n_{i \downarrow} )} \end{array}$ 就是Anderson（1987）提出的著名共振价键态。Anderson提出，当海森堡模型母体绝缘基态中预先存在的单重态配对在系统进行空穴掺杂后变为可移动的带电超导配对时，共振价键态（RVB state）就会导致超导。对电子（和量子自旋）系统此类基态的研究，在许多方面已不仅是对有用尝试态的建议（Edegger等人，2007），更成为了一种物质理论。

电子系统的尝试态日益复杂，伴随而来的是变分参数数量的增加（Sorella, 2005; Tahara and Imada, 2008; Neuscamman等人，2012）。绝大部分注意力都集中在揭示高温超导体的基态上。在高温超导体中，反铁磁性的同时存在导致越来越多地使用能够容纳两种长程有序的尝试态。例如（Giamarchi and Lhuillier, 1990, 1991），

其中产生算符由（9.17）定义，系数由（9.20）定义。实际上，各种附加投影算符可以加在 $\left| {\psi_{T}} \right.$ 之前。通常，这些投影算符会引入额外的变分参数。相关例子见 Sorella 等人 (2002); Sorella (2005); Watanabe 等人 (2006); Weber 等人 (2006); Tahara and Imada (2008)。

尽管许多变分蒙特卡洛计算都是在单带哈伯德模型上进行的，但也有一些研究了其他模型，如周期安德森模型（periodic Anderson model）（Shiba, 1986; Shiba and Fazekas, 1990; McQueen and Wang, 1991; Oguri and Asahata, 1992; Yanagisawa, 1999; Watanabe and Ogata, 2009）和三带哈伯德模型（Oguri and Asahata, 1992; Oguri 等人，1994）。正如 Edegger 等人（2007）的综述所证明的，紧接在哈伯德模型之后的受关注对象是 $t J$ 模型，我们接下来将讨论它。该模型连接了 Gutzwiller 投影态和 RVB 态的主题，这是下一节的主题。实际上，Anderson（1987）正是针对 tJ 模型提出了（9.22）作为其基态。

9.2.2 tJ 模型#

当排斥性在位库仑相互作用很大时，除了使用 $P_{G}^{0}$ 从 $\left| {\psi_{0}} \right.$ 中投影掉所有双占据外，我们还可以将（哈伯德）哈密顿量替换为在 $U \gg t$ 极限下与之等价的哈密顿量。“等价”指的是，该哈密顿量与原始哈密顿量具有相同的本征值和相同的本征态，最多差一个整体的幺正变换。然后，我们使用适合等价哈密顿量物理的尝试态进行分析。

根据（9.11）和（9.21），我们很容易确信，对于 $t = 0$，双占据格点数 $N_{D}$ 是一个好量子数，并且对于给定的电子密度，具有最小 $N_{D}$ 的状态是能量最低的状态。这些状态高度简并，在半满情况下简并度最高，此时每个格点有一个电子。将讨论限制在半满或更低的电子密度，我们注意到 $N_{D} ~ = ~ 0$ 标志着能量最低的状态。在构造等价哈密顿量时，我们将动能视为打破 $N_{D} = 0$ 状态间简并的微扰，并通过执行二阶简并微扰论来推导等价哈密顿量。

我们可以通过进行幺正变换最轻松地执行此微扰论

其中 $S$ 是某个厄米算符，并且

将 $H = H_{0} + H_{1}$ 写出，我们取 $H_{0}$ 为库仑项，$H_{1}$ 为动能项。现在我们将此变换展开为 $S$ 的幂级数，

并通过选择 $S$ 使得

我们得到一个有效哈密顿量

更一般地，如果我们对由某量子数 $n$ 标记的态进行微扰，那么在矩阵表示中，我们取

其中 $m \neq n .$ 。对于手头的问题，态 $n$ 是没有双占据格点的那些态，而态 $m$ 是具有一个或多个双占据格点的态。只有一个双占据格点的态占主导。在二次量子化表示中，

经过一些代数运算，我们得到

其中

并且

这里

其中 $\sigma^{\alpha}$ 是泡利自旋矩阵。交换相互作用为 $J = 4 t^{2} / U^{2}$

有效哈密顿量 (9.23) 就是 tJ 模型。它只能在一个没有任何格点被双占据的空间中运作。如前所述，省略的一项使其在大 $U$ 极限下只是近似等价于哈伯德模型。在半满填充时，晶格上向上和向下的电子数目相等，每个格点一个电子，此时 $t J$ 模型约化为自旋 $\frac{1}{2}$ 反铁磁海森堡模型

在半满状态下引入空穴掺杂后，每个格点有零个或一个电子（上旋或下旋），其动能具有物理上的解释：空穴在反铁磁涨落的背景中跳跃。

9.2.3 价键态#

随着tJ模型的引入，价键态变得重要起来，它既代表了一种新的绝缘基态（自旋液体），也作为海森堡模型量子蒙特卡罗研究的一种便捷变分基。安德森在高温超导出现之前就提出，这类态是新型绝缘基态的表述（Anderson, 1973）。仍然在高温超导出现之前，Fazekas和Anderson提出，它们是受挫晶格上反铁磁海森堡模型基态的表述（Fazekas and Anderson, 1974）。在某种意义上，这些态已成为强关联电子系统中打破自旋旋转对称性和平移对称性的新颖基态的范式。

安德森提出的新颖基态是自旋液体态，一种具有能隙且仅在相反取向电子自旋间存在短程有序的态。一维中的一些量子自旋模型具有RVB自旋液态基态（Majumdar and Ghosh, 1969a,b; Affleck et al., 1987, 1988），即价键态的线性组合。无能隙自旋液体和RVB晶体只是关于新颖RVB基态的额外提议中的一小部分。这些态通常缺乏序参量。它们似乎特别适用于像三角晶格这样的受挫晶格上的模型。这里，我们仅将这类态视为变分蒙特卡罗模拟的便捷基。

在价键基中，尝试态是价键态的叠加，

价键态是所有可能配对的电子束缚成单重态的直积态。这些态是超完备的，即

它们共有 $N ! / ( [ ( N / 2 ) ! ]^{2} ( N / 2 + 1 ) )$ 个。我们的主要目标是详述关于该基的几个事实，这些事实在变分蒙特卡罗模拟中是有用的。

Marshall (1955)的一个定理指出，在二部晶格上的哈伯德模型和海森堡模型的基态必须具有总自旋 S 为零的数值。对于海森堡模型，标准基 $\{| C \rangle \}$ 是 $\{| s_{1} , s_{2} , \ldots , s_{n} \rangle \}$ 这个集合，包含 $2^{N}$ 个构型，每个构型是每个晶格点上自旋 $-\frac{1}{2}$ 算符的 z-分量的可能本征态。通过选择一半格点自旋向上、另一半自旋向下的态，很容易选出总 $S^{z} = 0$ 的态，但使用这些态的基集可能会将更高总自旋的贡献包含到尝试态展开中。价键基是一个具有给定总自旋值的基。我们只专注于构建那些 $S = 0$ 的态。

我们处理如下形式的直积态

其中 $| ( i , j ) \rangle$ 是晶格点 i 和 j 之间的一个自旋单重态键。在二部晶格上，其约定是

这里 i 是 A 子格点，j 是 B 子格点。这个约定固定了单重态的相位。如果没有这个相位约定，我们就必须在格点上自旋交换时跟踪基态的符号。该约定还将基中态的数量减少到 $( N / 2 ) !$，但基仍然是超完备的。这种超完备性源于一般关系式

对于单重态乘积态。该关系表明，如果单重态对键合同一亚晶格上的位点，我们可以将其重新表示为单重态键合不同亚晶格上位点的态之和。不难证明：

由第一个方程可知，一对的总自旋大小为零，即$( \vec{S}_{i} + \vec{S}_{j} )^{2} | ( i , j ) \rangle = 0$

变分蒙特卡洛（Variational Monte Carlo）

图9.1 两个价键态(a) |V 和(b) $\vert V^{\prime} \rangle$ 以及它们的交叠(c) $\langle V \vert V^{\prime} \rangle$ 。A亚晶格位点为黑点。正文中左侧边缘的位点从下到上编号为1、5、9和13。

图9.1a和b展示了$4 \times 4$方格上的两个可能价键态。A亚晶格位点为黑点。为方便讨论，我们将这些态表示为$^{3}$：

和

这些数组的顶行是A亚晶格位置。底行是与之相连的B位点。这种连接就是价键。该键不一定是物理键。给定基矢中的一个，比如|V，我们可以通过置换其底行的$N / 2$个数字来生成所有其他基矢。根据我们对每个单重态对归一化的选择，

更一般地，

后一个表达式表明，对于$\langle V | V \rangle$，$N_{\mathrm{loops}} ~ = ~ N / 2$。该表达式也需要一些解释。

图9.1c是对应于$\langle V \vert V^{\prime} \rangle$的图表。它是图9.1a和b图表的叠合。我们看到这些键形成了$N_{\mathrm{loops}}$个环。在当前情况下$N_{\mathrm{loops}} = 3$且$N = 16$，所以$\langle V | V^{\prime} \rangle = 1 / 32$。如果我们通过下式定义|V的伴随态：

其中我们将B位点放在顶行，则可以通过交错排列(9.26)和(9.29)来表示$\langle V|V^{\prime} \rangle$：

在这种表示中，A位点的列对应右矢的键，而B位点的列对应左矢的键。选取顶行的一个位点，我们在底行找到其成键位点。将该位点带回顶行，我们可以找到与其成键的位点。重复此过程直到我们回到起始位点，便描画出重叠图中的一个环。很容易证明，在这种特定情况下只有三个环。

除了计算重叠以外，识别和计数环之所以重要，是因为

利用这一结果，我们可以计算能量、自旋-自旋关联函数和交错磁化强度，这些是反铁磁自旋-1/2海森堡模型的主要关注对象。

在价键基中，尝试态的形式为

如Liang等人（1988）所建议，我们假设

并进一步引用Marshall（1955）的另一个定理，即对于任意二分格点上的海森堡模型，我们可以取 $f_{i_{1} , j_{1} , i_{2} , \ldots , i_{N / 2}}^{} > 0$ ，从而可取 $h ( i , j ) > 0$ 。我们还假设 $h ( i , j ) ~ = ~ h ( | i - j | )$ 。对于h的选择有多种方案。例如，Liang等人（1988）通过调整几个最近邻距离上的 $h ( i , j )$ 来优化 $| \psi_{T} \rangle$ ，然后通过假设函数 $h ( i , j ) \sim | i - j | ^{\alpha}$ 的幂律形式并调整幂指数α来描述剩余部分。重要的事实在于，该尝试态同时支持短程有序和长程有序，具体取决于衰变是短程还是长程的。这种灵活性是价键基的一个优势。

图9.2 提出的蒙特卡洛移动。

由于价键基是非正交且过完备的，我们不能直接使用第9.1.2节中描述的变分蒙特卡洛方法。相反，对于那些适用Marshall正性定理的模型，我们利用(9.1)和(9.25)式，并写出

其中

这表明我们可以用(9.30)来评估构型能量$H_{W^{\prime}}$，并利用$f_{V} f_{V^{\prime}} \langle V | V^{\prime} \rangle$作为采样函数。遵循Liang等人的做法，我们采用Metropolis方法进行采样。他们的建议是选择两个A子格点，比如$a$和$^{b ,}$，然后更新与这些格点相关联的$| V \rangle$或$\vert V^{\prime} \rangle$。这里我们选择$| V \rangle$。如果其在$a$和$^ b$格点处的键末端是格点$c$和$d$，那么提议的移动是通过交换末端来创建一个新状态$| V^{\prime \prime} \rangle$。我们在图9.2中展示了这一移动。利用$P_{V V^{\prime}}$，我们根据下式接受该提议：

我们可以将重叠比简化为

这显著提高了效率。提议的移动要么将所有格点最初位于同一环上时使环数增加1，要么在它们不在同一环上时使环数减少1。判断属于哪种情况只需查表即可。在接受移动后，我们需要更新这一结构。Sandvik 和 Evertz (2010) 提出了另一种方法，有望获得更高的计算效率。

9.2.4 张量网络态#

我们现在简要讨论矩阵乘积态（matrix product state, MPS）和投影纠缠对态（projected entangled pair state, PEPS）。我们使用“张量网络态”作为统称术语，并且目前正在积极开发一个不断增长的关联态列表，主要用于扩展密度矩阵重整化群（density matrix renormalization group, DMRG）方法（White, 1992）的适用范围。这些态与前几小节讨论的态的一个显著区别在于，它们强调在构建基态$\left| {\psi_{0}} \right.$的估计量$\left| {\psi_{T}} \right.$过程中增加纠缠度，而不是像Gutzwiller态（9.12）那样事后通过作用于$| \psi_{0} \rangle$的投影算子来增加纠缠。

到目前为止，我们使用的多体试探波函数的基展开形式为

其中基态$\vert \phi_{i_{1}} , \phi_{i_{2}} , \ldots , \phi_{i_{N}} \rangle$是一个张量积态

采用这一ansatz后，我们系统的空间由与不同格点相关联的N个子空间复合而成。每个子空间由基$\{| \phi_{i_{j}} \rangle \}$张成。例如，对于量子自旋系统，$| \phi_{i_{j}} \rangle$是$S_{i}^{z}$自旋-$\frac{1}{2}$算符的两个本征态|↑⟩和|↓⟩之一。通常，每个子空间的维数相同。对于自旋-$\frac{1}{2}$的例子，维数为2。以下我们将假设它为d。

物理态在这种复合基下的展开很常见，部分原因是系统动力学往往是局域的。对于格点模型，以这种方式构建的基具有随格点数量增加而自然扩展的额外便利。然而，随着格点数量的增加，基态的数量以及展开式(9.31)中的项数会呈指数增长。张量网络态在展开式(9.31)中也使用了张量积基。新的要素在于这些基态振幅的构造方式，进而也是试验态参数化的构造方式。

让我们先做一个对比。从广义Hartree-Fock视角出发构建试验态，会将多体问题替换为一个有效的单体问题，其在波数空间中的基态非常紧凑：

这个态没有纠缠。各种Jastrow因子和其他投影仪可能会引入纠缠，但只有在构建$| \psi_{0} \rangle$之后才能实现。新的思路是，通过用不同于平均场构造$c_{i_{1} , i_{2} , \ldots , i_{n}} \approx c_{i_{1}} c_{i_{2}} \cdot \cdot \cdot c_{i_{N}} ,$的方式来近似$c_{i_{1} , i_{2} , \ldots , i_{n}}$，从而在先验上给基态引入更多的纠缠。

我们从一维的MPS开始。如果我们的系统具有开边界条件，我们定义MPS为

更紧凑地，

对于每个格点上的每个态，我们关联一个不同的矩阵，但边界上的格点除外，我们为其关联行向量和列向量。由于开边界条件下系统性质不均匀，每个格点上的矩阵各不相同。如果系统具有周期性边界条件，我们在每个格点放置相同的矩阵，并定义

通过利用每个格点上与A相关的额外自由度，我们可以在相邻格点之间纠缠这些状态，并在多个格点之间协调这种纠缠。向量和矩阵的阶通常用D表示，称为键维数。这个参数，以及施加物理或假定的对称性，控制了尝试状态参数化的程度。

对于MPS，如同其他张量网络态一样，我们构造展开系数 $c_{i_{1} , i_{2} , \dots , i_{N}}$ 不仅是为了反映晶格模型的对称性，也是为了最大化量子纠缠。我们期望展开系数 $c_{i_{1} , i_{2} , \dots , i_{N}}$ 的结构也能反映动力学的局域性，例如，当任意两个格点之间的距离超过几个晶格间距时，这些系数会很小。如果是这样，小的展开系数会降低相关基态的重要性。这一期望是在RVB基下对尝试状态的展开系数进行成对分解的合理性依据。在该基组中，子系统是格点对上的自旋态，而不是单个格点上的自旋态；并且不涉及所有可能的自旋态，而只存在自旋单态。RVB态是一类特定的张量网络态。自旋单态是纠缠态。

矩阵乘积态自然地在一维系统中作为DMRG方法（Takasaki et al., 1999）的定点基出现。这种确定性的基态方法对于在一维晶格或准一维系统（如少腿梯子和Cayley树）上定义的、具有开边界条件和短程相互作用的模型效果最佳。关于该方法为何对这些几何构型和边界条件效果如此出色，存在一些物理上的见解。这些见解围绕面积律（area law）（Verstraete et al., 2008; Cirac and Verstraete, 2009）展开。

面积律适用于基态和第一激发态之间存在能隙且具有短程相互作用的系统。它指出，如果我们考虑一个由最近邻粒子（比如自旋）组成的块A，那么与该块的基态约化密度矩阵 $S = - \mathrm{Tr} \rho_{A}$ $\rho_{A}$ 相关的信息论熵尺度与块边界上的粒子数成正比，而非块内部的粒子数。PEPS可以约化为MPS。在这两种情况下，我们在每个格点上放置一个对象，其下标数量等于其最近邻格点的配位数。在二维中，格点的配位数高于一维。同样重要的是，所需基态的数量随系统大小的多项式规模增长。这种缩放规律表明，原则上我们可以从这些态中构造出一个高效的基组。使其比目前更高效是一个活跃的研究课题。在许多情况下，尽管是多项式规模，但缩放仍然很差。

这些基态也揭示了DMRG的变分性质（Ostlund and Rommer, 1995; Dukelsky et al., 1998; Takasaki et al., 1999）。因此，它们可以作为变分尝试波函数（variational Ansätze）用于变分蒙特卡洛计算。DMRG的细节在Schollwock (2005)中有综述，而关于张量积态及其与DMRG联系的技术细节在McCulloch (2007)、Verstraete et al. (2008)以及Cirac and Verstraete (2009)中有所论述。这里，我们将更为谨慎，首先介绍矩阵乘积基形式背后的动机，然后非常简要地讨论张量网络态在变分蒙特卡洛模拟中的应用。在变分蒙特卡洛模拟中使用MPS和PEPS是一个快速发展的课题。

为了引出MPS的形式，我们首先设想系统是由两个子系统组合而成，它们分别具有标准正交基 $\{| a_{i} \rangle \}$（张成一个子空间）和标准正交基 $\{| b_{j} \rangle \}$（张成另一个子空间）。然后，创建一个标准的张量积基，我们可以写出

然而，施密特定理（Nielsen and Chaung，2000，第2.5节）指出，存在一种更紧凑的表示形式：如果态$| \psi \rangle$位于两个子系统的复合空间中（这两个子系统的空间维数不一定相同），那么对于这些子系统，存在一组正交归一基$\{| u_{i} \rangle \}$和$\{| \nu_{j} \rangle \}$，使得

其中$\textstyle \sum_{i} \lambda_{i}^{2} = 1$。$\lambda_{i}$是非负实数。非零$\lambda_{i}$的数目称为施密特数。它等于两个子空间中较小空间的维数。如果施密特数等于1，则$| \psi \rangle = | u \rangle \otimes | \nu \rangle$是非纠缠的（纯态）。如果施密特数大于1，则该态是纠缠的。因此，施密特数是$| \psi \rangle$中纠缠程度的一种度量。施密特分解一次仅适用于两个子系统，并且其唯一性仅在任意子系统基的酉变换意义下成立。

现在我们将施密特分解应用于张量积展开式(9.31)，首先将由态$\{| \phi_{i_{1}} \rangle \}$张成的空间作为一个子系统，由态$\{| \phi_{i_{2}} , \phi_{i_{3}} , \ldots , \phi_{i_{N}} \rangle \}$张成的空间作为另一个子系统。根据施密特定理，我们得到

由于基集$\{| \phi_{i_{1}} \rangle \}$和$\{| \psi_{\alpha_{1}} \rangle \}$张成相同的空间，我们可以变换回原始基$| \phi_{i_{1}} \rangle$，从而写成

接下来，我们对$| \psi_{\alpha_{1}}^{[ 2 , N ]} \rangle$应用施密特分解，并再次变换主导态，得到

重复这些步骤直到遍历每个位置，会得到一个看起来杂乱的表达式，但其形式与(9.32)相同。这些步骤产生了用矩阵乘积展开中的系数表示的$c_{i_{1} , i_{2} , \ldots , i_{n}}$的参数化形式。显然，这种参数化可以非常丰富。

当我们在变分蒙特卡洛计算中使用张量网络态时，基本算法与第9.1.2节中先前描述的相同，其中

它们在使用上的一个区别在于计算概率比率 $P ( C^{\prime} ) / P ( C )$ 的计算成本，因为计算矩阵乘积会产生开销。利用迹运算的循环特性、存储部分矩阵乘积并对其进行更新，而不是整个乘积字符串，这一点非常重要。关于矩阵元素、其乘积以及算符期望值的高效计算的详细信息，可参考 Sandvik and Vidal (2007)， Schuch et al. (2008)， Sfondrini et al. (2010)，以及 Wang et al. (2011a) 等文章。这个领域仍在发展，效率可能会进一步提高。

通常，通过利用张量网络概念的灵活性来实现效率。没有必要将纠缠限制在最近邻对之间。我们可以将自旋纠缠在一个最近邻的自旋plaquette中。在这方面，我们还可以将展开系数分解为对乘积（Changlani et al., 2009; Marti et al., 2010; Boguslawski et al., 2011; Neuscamman et al., 2011; Neuscamman et al., 2012），沿袭 Liang et al. (1988) 的方式，但单重态对的放置具有规律性；或者分解为plaquette乘积（Sandvik, 2008; Mezzacapo et al., 2009; Mezzacapo, 2011; Wang et al., 2011b）；或者分解为字符串乘积（Schuch et al., 2008; Sfondrini et al., 2010）。我们既可以包含单重态，也可以包含三重态（Neuscamman et al., 2012）。我们还可以将投影算子重新引入到计算中。目前已有用于玻色子、电子（Corboz et al., 2010a,b; Neuscamman et al., 2012）、量子自旋系统，以及最近用于连续介质中电子系统（Marti et al., 2010; Boguslawski et al., 2011; Neuscamman et al., 2012）的变分蒙特卡洛模拟的张量网络方案。虽然我们的介绍主要从一维角度出发，但大部分计算工作都聚焦于二维。

随着波函数参数化的增强，对试探能量或其能量方差的极小化成为模拟中非常重要的一部分（Schuch et al., 2008; Sfondrini et al., 2010; Wang et al., 2011b）。在某些情况下，特定的优化技术能很好地适用于特定的张量网络ansätze。在其他情况下，优化过程遵循与我们接下来将要讨论的类似步骤。本讨论的重点是一种近期提出的通用试探态方法，适用于参数集规模较小和较大、以及有无Jastrow型投影算子的各种问题。Neuscamman et al. (2012) 的工作就使用了这种方法。

9.3 试探态优化#

试探态优化最明显的方法是在一组蒙特卡洛样本上极小化变分能。实际上，直到大约1988年，这仍是常用的方法。大致来说，该方法是先从一个依赖于一组参数 $p = ( p_{1} , p_{2} , . . . , p_{N_{\mathrm{opt}}} )$ 的好的试探态 $| \psi_{T} ( p ) \rangle$ 开始，然后使用该态的模平方作为分布函数，对一组足够大的构型 $\{C_{1} , C_{2} , \ldots , C_{N_{\mathrm{config}}} \}$ 进行变分蒙特卡洛采样，使得和 $\sum_{i} | \psi_{T} ( C_{i} , p ) \rangle / N_{\mathrm{config}}$ 能够较好地近似基态。接下来，在固定构型系综的情况下，连续生成多组新参数，每组相对于前一组只有微小变化。使用一种“贪婪”蒙特卡洛方法，根据变分能是升高还是降低来接受或拒绝一个新状态。

这种方法的缺点在于，构型系综中的蒙特卡洛样本数量必须比待优化的参数数量高出多个数量级。因此，只能优化少数几个参数。之所以需要这种不平衡，是因为多参数试探波函数具有拟合数据的灵活性，可能导致参数集上的系综变分能低于真实能量。然而，当这些参数用于评估新独立抽取的构型系综上的变分能时，可能会得到高得多的能量。人们发现，最小化构型能量的方差（Bartlett 等人，1935；Coldwell，1977；Umrigar 等人，1988）而非最小化构型能量可以克服这个问题。在1988年至2001年期间，进行这种方差最小化是标准做法。

现在已有更复杂的方法，可以有效优化能量，或者更佳地，优化能量与局部能量方差的任意线性组合（Lin 等人，2000；Nightingale 和 Melik-Alaverdian，2001；Schautz 和 Filippi，2004；Sorella，2005；Umrigar 和 Filippi，2005；Toulouse 和 Umrigar，2007，2008；Umrigar 等人，2007）。我们在这里只讨论这些方法中最有效的一种，最初由 Nightingale 和 Melik-Alaverdian (2001) 提出，用于同时优化基态和激发态的试探波函数。他们的试探波函数对参数具有线性和非线性依赖性。他们通过最小化能量来优化线性参数，通过最小化方差来优化非线性参数。Nightingale 方法后来被扩展，允许通过能量最小化（Toulouse 和 Umrigar，2007；Umrigar 等人，2007）来同时优化线性和非线性参数，然后进一步扩展以优化能量和方差的任意线性组合（Toulouse 和 Umrigar，2008）。我们主要关注这最后一种优化方法，并将其称为“线性方法”（linear method），尽管它也优化非线性参数。刚才提到的各种方法都是在连续介质中的多电子问题背景下发展起来的，因此针对的是Slater-Jastrow形式的试探态。所以在文献中，一些算法细节可能会侧重于最优地处理这种形式，但基本方法通常具有普适性，实际上也已被应用于晶格问题（Changlani 等人，2009）。

9.3.1 线性方法#

至此，采用开创性论文中简洁的符号表示会更为方便。特别地，我们将构型基 $\{| C \rangle \}$ 视为N个粒子标记位置 $R = ( r_{1} , r_{2} , \ldots , r_{N} )$ 的基集 $\{| R \rangle \}$

接下来，我们对投影到这个基集的量使用函数表示法，例如，$\psi_{T} ( R ) =$ $\langle R | \psi_{T} \rangle$。最后，我们去掉下标 T，使得 $\psi ( R ) = \psi_{T} ( R )$。我们注意到，与构型能量的情况类似，当 $\psi ( R )$ 成为精确的本征态时，$E ( R )$ 是常数且等于该态的能量。因此，可以通过最小化构型能量的方差来优化试探态。在连续介质中，构型能量被称为局部能量（local energy）。

在每个优化步骤中，我们围绕当前参数集 $p_{0}$ 将试探波函数展开至 $\Delta p = p - p^{0}$ 的线性阶

其中 $\psi_{0} ( R ) = \psi ( R , p_{0} )$ 是当前的波函数（还不是基态），并且

是波函数相对于参数的 $N_{\mathrm{opt}}$ 个导数。这个线性化波函数成为变分原理中用于能量的尝试态。在无限多个蒙特卡洛样本下，使用(9.33)的线性化波函数计算能量时，使能量最小化的最优参数变化量 $\Delta p$ 是所得广义特征值方程的最低特征值解：

其中 $H$ 和 $S$ 是由当前波函数及其导数 $\left\{\psi_{0} , \psi_{1} , \psi_{2} , \ldots , \psi_{N_{\mathrm{opt}}} \right\}$ 构成的 $( N_{\mathrm{opt}} + 1 )$ 维基下的哈密顿矩阵和重叠矩阵。这里 $\Delta p_{0} = 1$ 。抽取有限数量的蒙特卡洛样本，我们将这些矩阵估计为：

其中

且

在无限样本下，$\langle j \rangle = 0$ 且

因此$H_{i j}$是对称的，但在有限样本上并非如此。由于真实哈密顿量矩阵是对称的，并且非对称广义特征值问题可能产生复特征值，我们或许会认为对矩阵进行对称化处理是可取的。然而，这种非对称性实际上是有利的，因为它产生了强零方差性质（Nightingale和Melik-Alaverdian，2001），即$\Delta p$每个分量的方差不仅在$\psi_{T} ( R , p )$变得精确时消失，而且在仅$\psi_{\mathrm{lin}} ( R , p )$变得精确时也消失，也就是说，在基态张成哈密顿量不变子空间的极限下消失。在实践中，非对称哈密顿量将参数波动降低一到两个数量级。

9.3.2 更新参数#

优化算法的一部分现已就位。在迭代的每一步，我们围绕$p_{0}$对$\psi$进行线性化，并通过蒙特卡洛采样估计$H_{i j}$和$S_{i j}$。然后，利用现成的软件，我们求解广义特征值问题，得到具有最低特征值的特征向量$\Delta p$。通过求解式(9.34)获得参数变化$\Delta p$后，下一个问题是：我们如何利用它们来更新试探波函数中的参数？

将当前参数增加$\Delta p$，即$p_{0} \gets p_{0} + \Delta p$的最简单方法对于线性参数有效，但如果式(9.33)的线性近似不够好，则对于非线性参数无效。调整ψ中的参数使其适配优化的线性形式是一种可能的策略。Umrigar及其合作者提出了一种更简单的替代方案，该方案利用了试探态归一化的自由度。ψ归一化选择的任意性导致了参数解的任意性。对于线性参数，得到的是相同物理波函数但归一化不同，而对于非线性参数，它在每个优化步骤都会改变预测的波函数。

现在我们考虑归一化不同的波函数$\overline{{\psi}} ( R , p ) = N ( p ) \psi ( R , p )$，其中要求$\overline{{\psi}} ( R , p_{0} ) = \psi ( R , p_{0} ) \equiv \psi_{0} ( R )$，并且$N ( p )$仅依赖于非线性参数。那么，在$p = p_{0}$处$\overline{{\psi}} ( p )$的导数为

$\overline{{\psi}} ( p )$的线性近似为

由于我们通过在同一变分空间中的优化得到$\overline{{\psi}}_{\mathrm{lin}}$和$\psi_{\mathrm{lin}}$，它们彼此成正比，

$\begin{array} {r} {\mathrm{so} c = 1 / ( 1 - \sum_{i = 1}^{N_{\mathrm{opt}}} \Delta p_{i} N_{i} )} \end{array}$，且$\Delta \bar{p}$与$\Delta p$通过统一尺度变换关联

此外，由于重缩放因子可以在 $- \infty$ 到 $\infty$ 之间任意取值，归一化选择不仅会影响参数变化的幅度，还会影响其符号。

我们还需讨论归一化因子 $N_{i}$ 的参数梯度选择。Sorella (2005) 在另一种优化方法中采用了一个好的选择：选择 $N_{i}$ 使得试探态 ${\overline{{\psi}}}_{i}$ 的参数导数与 $\psi_{0}$ 正交。Umrigar 及其合作者采用了一个更通用的选择：他们选择 $N_{i}$ 使得

其中 $\xi$ 是介于 0 和 1 之间的常数，从而得到

其中 $\begin{array} {r} {D = \sqrt{1 + 2 \sum_{j}^{\mathrm{nonlin}} S_{0 j} \Delta p_{j} + \sum_{j k}^{\mathrm{nonlin}} S_{j k} \Delta p_{j} \Delta p_{k}}} \end{array}$，求和仅对非线性参数进行。$\xi = 0$ 对应 Sorella 的选择，而 $\xi = 1 / 2$ 则强制要求 $\overline{{\psi}}_{\mathrm{lin}}$ 的归一化与 $\psi_{0}$ 相同。

在上述描述中，我们遵循了 Umrigar 等人 (2007) 的公式化方法：先通过求解原始基矢下的哈密顿矩阵和重叠矩阵的广义特征值问题来计算 $\Delta p$，然后利用 (9.36) 计算 $\Delta \bar{p}$。另一种方法是，我们可以遵循 Toulouse 和 Umrigar (2007) 的公式化方法，直接从半正交基矢下的哈密顿矩阵和重叠矩阵构造 $\Delta \bar{p}$。这两种方法导致相同的参数变化。然而，如果使用了稳定化技术（将在后面讨论），这两种方法会导致略有不同的参数变化。

9.3.2 牛顿法#

在讨论使用线性方法最小化能量方差之前，我们做一个重要的延伸，揭示线性方法本质上是一种稳定化的牛顿法，在实践中通常比标准牛顿法收敛更快。

牛顿法是一种常用的多变量极小化方法（Press 等人，2007，第 9 章）。它首先考虑函数 $E ( p )$ 关于某点 $p_{0}$ 的截断泰勒级数展开

其中 g 是 $E$ 的梯度向量：

h 是 $E$ 的海森矩阵：

最小化这个二次表达式，我们得到最优参数由下式给出

由于 $E ( p )$ 一般不是 $p$ 的二次函数，因此需要迭代⁴

1. 求解 $h \cdot \Delta p = - g$ 得到 $\Delta p$

2. $p_{0} \leftarrow p_{0} + \Delta p ,$

直到 $\Delta p$ 的范数满足某个容差条件（算法30）。

关于试态优化，Umrigar 和 Filippi (2005) 利用在某些情况下两个随机变量 $X$ 和 $Y$ 的协方差 $\operatorname{cov}(X,Y) = \langle (X - \langle X \rangle)(Y - \langle Y \rangle) \rangle = \langle XY \rangle - \langle X \rangle \langle Y \rangle$ 的波动远小于乘积 $\langle XY \rangle$ 的波动这一观察，给出了梯度 $g$ 的零方差表达式以及海森矩阵 $h$ 的减方差表达式。这些零方差表达式通常能将效率提高两到三个数量级。

显然，我们可以使用牛顿法来最小化不仅是能量，还有局部能量的方差以及两者的任意线性组合。我们只需将 (9.38) 中的表达式替换为所需线性组合的二次近似即可。

算法30 牛顿法。 输入：$p_{0}$ 和 $g_{\mathrm{tol}}$，以及计算 $\begin{array} {r} {g_{k} = \frac{\partial E} {\partial p} | _{p_{k}}} \end{array}$ 和 $\begin{array} {r} {h_{k} = \frac{\partial^{2} E} {\partial p^{2}} | _{p_{k}} .} \end{array}$ 的过程 $k = 0$ ; 重复 求解 $h_{k} \delta p_{k} = - g_{k}$ 得到 $\delta p_{k}$ $p_{k + 1} = p_{k} + \delta p_{k}$ $k \gets k + 1 ;$ 直到 $\| g_{k + 1} \| < g_{\mathrm{tol}} .$ 返回 $p_{k + 1} .$

9.3.3 线性方法与牛顿法之间的联系#

我们现在建立线性方法与牛顿法之间的联系（Toulouse 和 Umrigar, 2008）。我们首先将能量 (9.38) 重新表达为一种更具启发性的分块矩阵形式

其中 $\bar{S}$ 是重叠矩阵，并且

是在 $\xi = 0$ 的半正交基下的 $N_{\mathrm{opt}} \times N_{\mathrm{opt}}$ 哈密顿矩阵。对 $p$ 最小化该能量得到一个广义特征值方程

这与线性方法中的情况类似，只不过我们得到了一个对称的哈密顿矩阵。

9.3.4 能量方差优化#

在建立了线性方法与牛顿法之间的联系之后，我们现在修改线性方法来优化局部能量的能量和方差的线性组合。类似于 (9.38)，所需的一切是对方差进行二次近似

最小化导出以下广义特征值方程

Toulouse 和 Umrigar (2008)给出了$g_{V}$和$h_{V}$的显式表达式。

9.3.5 稳定化#

在线性方法和牛顿法中，如果能量的二次近似太差，或者用于评估H和S矩阵的蒙特卡洛样本数量太少，更新后的参数可能不稳定。我们可以通过以下方式稳定这两种方法：牛顿法通过在Hessian矩阵的对角元素上添加一个正常数 $a_{\mathrm{diag}} \geq 0$ 来稳定，这样提议的移动会变小，并向最速下降方向旋转。我们在线性方法中通过将 $a_{\mathrm{diag}}$ 添加到哈密顿矩阵的对角线上来达到同样的效果：$H_{i j} \to H_{i j} + a_{\mathrm{diag}} \delta_{i j} ( 1 - \delta_{i 0} )$ 。随着 $a_{\mathrm{diag}}$ 增大，参数变化 $\Delta p$ 会变小并旋转。我们可以在每个优化步骤中自动调整 $a_{\mathrm{diag}}$ 的值，具体如下：一旦我们计算出矩阵H和S，我们使用三个相差10倍的 $a_{\mathrm{diag}}$ 值来预测三个新的波函数。然后，我们使用关联采样（correlated sampling）进行一次短暂的蒙特卡洛运行，以比能量本身更高的精度计算这些波函数的能量差。利用这些关联能量，我们计算出接近最优的 $a_{\mathrm{diag}}$ 值。所谓关联采样，是指我们不通过独立的蒙特卡洛运行来估计相关量，而是使用来自单次蒙特卡洛运行的 $N_{\mathrm{config}}$ 个蒙特卡洛构型来估计它们。其优点是这些量的相对误差远小于单个量的绝对误差。在我们的例子中，相关量是三个波函数的能量。我们在附录N中讨论关联采样。

9.3.6 线性优化方法和牛顿优化方法总结#

现在我们以一种强调其共同特征的方式来总结线性方法（算法31）和牛顿方法（算法32）。两种方法都是迭代的，并在两次蒙特卡洛运行之间交替进行。具体来说，在线性方法和牛顿法中，我们首先使用一种蒙特卡洛运行来估计重要的矩阵。在线性方法中，我们使用蒙特卡洛采样来计算重叠矩阵和哈密顿矩阵。然后，利用哈密顿矩阵和重叠矩阵，我们求解一个广义特征值方程以获得提议的参数变化。在牛顿法中，为了获得提议的参数变化，我们使用蒙特卡洛采样来计算梯度向量和Hessian矩阵，然后求解一个线性方程组。我们为三个不同的 $a_{\mathrm{diag}}$ 值求解特征值问题或线性系统，以获得三组提议的新参数。然后，我们进行第二种蒙特卡洛运行，该运行采用关联采样（附录N）来获得目标函数（能量和能量方差的线性组合）的差值，其精度高于单个值。利用这些计算出的差值，我们获得接近最优的 $a_{\mathrm{diag}}$ 值，然后使用该值获得下一个迭代循环的参数。

算法31 线性化优化方法。 输入：参数 $p$ 和容差 $\delta .$ 重复 进行一次蒙特卡洛运行以估计哈密顿矩阵 $H_{i j}$ 和重叠矩阵 $S_{i j}$ 用三个 $a_{\mathrm{diag}}$ 值求解广义特征值问题，以计算三组重新缩放的参数增量； 进行一次包含关联采样的短期蒙特卡洛运行，以计算这三个波函数的相对能量，从而找到接近最优的 $a_{\mathrm{diag}}$ 使用接近最优的 $a_{\mathrm{diag}}$，设置 $p = p_{0} + \Delta \bar{p}$ 直到 $| \Delta \bar{p} | < \delta .$ 返回优化后的 $p .$

算法 32 牛顿优化方法。 输入：参数 $p$ 和容差 $\delta .$ 重复 进行蒙特卡洛运行以估计Hessian矩阵 $h_{i j}$ 和梯度向量 $g_{i}$； 用三个 $a_{\mathrm{diag}}$ 值解线性方程，计算三组重新缩放参数增量； 用关联采样进行短蒙特卡洛运行，以计算这三个波函数的相对能量，从而找到接近最优的 $a_{\mathrm{diag}}$； 利用接近最优的 $a_{\mathrm{diag}}$，设置 $p = p_{0} + \Delta \bar{p}$； 直到 $| \Delta \bar{p} | < \delta$ 返回优化后的 $p .$

延伸阅读#

C. Gros, "Physics of projected wavefunctions," Ann. Phys. 189, 53 (1989).

B. L. Hammond, W. A. Lester, Jr., and P. J. Reynolds, Monte Carlo Methods for ab initio Quantum Chemistry (Singapore: World-Scientific, 1994), 第2章、第5章和附录B。

J. Toulouse and C. J. Umrigar, "Full optimization of Slater-Jastrow wave functions with application to the first-row atoms and homonuclear diatomic atoms," J. Chem. Phys. 128, 174101 (2008).

J. I. Cirac and F. Verstraete, "Renormalization and tensor product states in spin chains and lattices," J. Phys. A: Math. Theor. 42, 504004 (2009).

S. Sorella, "Variational Monte Carlo and Markov chains for computational physics," in Strongly Correlated Systems, Springer Series in Solid-State Physics 176, ed. A. Avella and F. Mancini (Heidelberg: Springer-Verlag, 2013).

习题#

9.1 如果尝试态精确到 $\mathcal{O} ( \varepsilon )$ ，证明Rayleigh商给出的能量估计精确到 $\mathcal{O} ( \varepsilon^{2} )$ 。

9.2 如果A是N阶Hermitian矩阵，B是由如下分块形式定义的Hermitian矩阵：

其中a是N维向量，α是标量。证明B的特征值 $\lambda_{1} \le \lambda_{2} \le \dots \le \lambda_{N + 1}$ 与A的特征值交错分布。如果A正定，则B也正定的充要条件是什么？

9.3 考虑一个无周期性边界条件的双位点电子系统（$\mathrm{H}_{2}$ 分子），哈密顿量分别为海森堡反铁磁模型、tJ模型和正U哈伯德模型。

1. 针对每种情况，解析地求出模型的特征值和特征向量，对于每个位点一个电子的填充以及根据模型相应更少的填充情况。

2. 在大U极限下，将哈伯德模型的能量和态映射到tJ模型的能量和态。

3. 在半满填充下，将tJ模型的能量和态与海森堡模型的能量和态联系起来。

9.4 验证Gutzwiller投影关系式(9.13)。

9.5 Yokoyama和Shiba (1988)不是将BCS变分波函数投影到固定N，而是通过变量替换实现相同结果：

适当定义c和d粒子的真空态 $\mathinner{| {\bar{0}} \rangle}$，使得 $c_{k} | \bar{0} \rangle = d_{k} | \bar{0} \rangle = 0$ ，则可得：

$\mathinner{| {\bar{0}} \rangle}$ 与原真空态 $\mathinner{| {0} \rangle}$ 有何关系？

9.6 推导两个价键态重叠的表达式(9.28)，以及自旋-自旋关联函数期望值的表达式(9.30)。