第6章 · 变分波函数的优化#
6.1 引言#
在第一个将变分蒙特卡洛方法应用于多体系统的研究中，所考虑的波函数相对简单，仅包含少数几个参数；因此，只需要对少数几个能量值进行简单的拟合过程，就可以获得这些参数的最优值（McMillan, 1965）。在蒙特卡洛框架中包含多个变分参数的可能性是相对较新的成就。事实上，直到九十年代初，人们仍然将优化限制在非常少（例如，最多两到三个）的变分参数上，这些参数通过试错法计算能量或通过拟合方法来确定。如今，人们对实现快速而准确的方案以优化大量参数来描述关联多体波函数的兴趣日益增长。
在下文中，我们用向量 ${\pmb{\alpha}} = ( \alpha_{1} , \ldots , \alpha_{p} )$ 表示描述波函数的所有 $p$ 个变分参数：
$$ \Psi_{\alpha} ( x ) = \langle x | \Psi_{\alpha} \rangle ,\tag{6.1} $$其中 $\{\vert x \rangle \}$ 表示希尔伯特空间中的一个完备基集（假定为正交归一化的）。优化的第一个突破在于认识到，一旦获得了 $\Psi_{\alpha} ( x )$ 的抽样，就可以利用重新加权技术来计算具有不同参数的波函数（即 $\Psi_{\alpha + \delta \alpha} ( x )$）的能量或方差（Umrigar et al., 1988）。这种方法被证明对于少量粒子（例如，几个原子和分子）是非常高效和稳定的，特别是在采用方差最小化时，而对于大量粒子则存在相当大的缺陷。在重新加权技术被应用的几年之后，在多个参数的随机优化技术方面取得了进一步进展。基本思想是针对参数的任意微小变化计算能量差，即评估能量导数。结果表明，这些量对应于定义明确的估计量，其方差仅随粒子数多项式增长（例如，最多二次增长）。最简单的方法基于最速下降最小化（Harju et al., 1997）。更精细的算法包括所谓的随机重构（Sorella, 1998, 2001）、牛顿-拉弗森方法的近似版本（Umrigar and Filippi, 2005; Sorella, 2005）以及所谓的线性方法（Nightingale and Melik-Alaverdian, 2001; Umrigar and Filippi, 2005; Umrigar et al., 2007; Toulouse and Umrigar, 2007）。所有这些技术的关键点在于设计一种方案，即使在存在统计噪声的情况下也能保持稳定，并且能够快速收敛到全局最小值。这些迭代方案基于能量导数，在给定的蒙特卡洛模拟中必须尽可能精确地计算这些导数。
我们强调，一种稳健且高效的波函数优化方法也将改善投影技术的收敛性，因为投影技术关键依赖于所使用的引导波函数的质量（参见第8章和第9章）。特别地，固定节点方法受益于优化的波函数，因为通过系统改进尝试态，这一近似变得精确（参见第10章）。
6.2 波函数优化的重新加权技术#
在此，我们简要讨论如何使用第3.2节中描述的重新加权技术，通过单次模拟获得多组变分参数 ${\pmb{\alpha}}$ 的能量或方差。特别地，我们考虑基于方差优化的方法，该方法由Umrigar等人（1988）引入后，被广泛使用了数年。假设一组构型 $\{x_{i} \}$（$i = 1 , \ldots , N$）是从变分态 $\Psi_{\alpha} ( x )$ 的平方中抽样得到的。当然，对应于该具有α的"初始"猜测的能量和方差由下式给出：
$$ E_{\alpha} = \frac{\langle \Psi_{\boldsymbol{\alpha}} | \mathcal{H} | \Psi_{\boldsymbol{\alpha}} \rangle} {\langle \Psi_{\boldsymbol{\alpha}} | \Psi_{\boldsymbol{\alpha}} \rangle} \approx \frac{1} {N} \sum_{i = 1}^{N} e_{L , \boldsymbol{\alpha}} ( x_{i} ) ,\tag{6.2} $$$$ \sigma_{\Psi_{\alpha}}^{2} = \frac{\langle \Psi_{\alpha} | ( \mathcal{H} - E_{\alpha} )^{2} | \Psi_{\alpha} \rangle} {\langle \Psi_{\alpha} | \Psi_{\alpha} \rangle} \approx \frac{1} {N} \sum_{i = 1}^{N} \left[ e_{L , \alpha} ( x_{i} ) - E_{\alpha} \right]^{2} ,\tag{6.3} $$其中局部能量由式 (5.7) 给出：
$$ e_{L , \alpha} ( x_{i} ) = \frac{\langle x_{i} | \mathcal{H} | \Psi_{\alpha} \rangle} {\langle x_{i} | \Psi_{\alpha} \rangle} .\tag{6.4} $$然后，利用式 (3.10)，结合下式：
$$ \mathcal{R}_{\alpha + \delta \alpha_{k}} ( x ) = \left| \frac{\Psi_{\alpha + \delta \alpha_{k}} ( x )} {\Psi_{\alpha} ( x )} \right| ^{2} ,\tag{6.5} $$我们可以（在同一组构型 $\{x_{i} \}$ 下）估算对应于（略微）不同变分参数 $\pmb{\alpha} + \delta \pmb{\alpha}$ 的波函数的能量和方差：
$$ E_{\alpha + \delta \alpha} = \frac{\sum_{i = 1}^{N} e_{L , \alpha + \delta \alpha} ( x_{i} ) \mathcal{R}_{\alpha + \delta \alpha} ( x_{i} )} {\sum_{i = 1}^{N} \mathcal{R}_{\alpha + \delta \alpha} ( x_{i} )} ,\tag{6.6} $$$$ \sigma_{\Psi_{\alpha + \delta \alpha}}^{2} = \frac{\sum_{i = 1}^{N} \left[ e_{L , \alpha + \delta \alpha} ( x_{i} ) - E_{\alpha + \delta \alpha} \right]^{2} \mathcal{R}_{\alpha + \delta \alpha} ( x_{i} )} {\sum_{i = 1}^{N} \mathcal{R}_{\alpha + \delta \alpha} ( x_{i} )} .\tag{6.7} $$因此，通过一次使用变分参数 α 的模拟，就可以“重构” ${\pmb{\alpha}}$ 邻域内不同取值对应的能量和方差。为了找到最优参数集，我们可以应用最速下降法（在给定的构型集 $\{x_{i} \}$ 下，这些是变分参数的非线性函数），针对 δα 最小化式 (6.6) 的能量或式 (6.7) 的方差。在重加权方法中，方差最小化的优势在于它是正项之和，因此总是存在极小值（换言之，方差总是正量，即使对于有限数量的样本也是如此）。相反，执行能量最小化时，并不能保证找到有下界的结果；实际上，可能出现这样的情况：根据 $| \Psi_{\alpha} ( x ) | ^{2}$ 采样的某个构型 $x_{i}$（对应于有限的局部能量 $e_{L , \alpha} ( x_{i} )$），对于参数 $\alpha + \delta \alpha$ 可能产生非常大且负的局部能量，例如 $x_{i}$ 可能任意接近 $\Psi_{\alpha + \delta \alpha} ( x )$ 的节面，使得 $e_{L , \alpha + \delta \alpha} ( x_{i} ) \to - \infty$。
重加权技术在大粒子数情况下是不稳定的。事实上，让我们考虑向量 α 中仅有一个分量 $\alpha_{k}$ 增加了 $\delta \alpha_{k}$ 的情况；那么重加权因子由两个按指数级增大或减小的多体波函数之比决定，导致：
$$ \mathcal{R}_{\alpha + \delta \alpha_{k}} ( x ) \approx \exp ( \pm \delta \alpha_{k} N_{p} ) ,\tag{6.8} $$其中 $N_{p}$ 是粒子数。因此，对于大 $N_{p}$ 的系统，即使是变分参数的微小变化也会导致统计噪声呈指数级增长。这一事实使得在包含多个粒子的系统中找到最优参数集变得极其困难（若非不可能），因为能量或方差的变分中统计噪声非常大。因此有必要考虑替代的优化方法：一种可能是基于能量导数计算的方法。实际上，尽管基于方差最小化的方法效果令人满意，但仍有充分理由发展基于能量最小化的方法：例如，能量优化的波函数通常比方差优化的波函数能提供更好的期望值，因为前者对低能变化比后者更敏感。
6.3 能量导数#
在本节中，我们讨论计算变分能量相对于给定变分参数 $\alpha_{k}$ 的导数所需的基本要素：
$$ f_{k} = - \frac{\partial E_{\alpha}} {\partial \alpha_{k}} = - \frac{\partial} {\partial \alpha_{k}} \frac{\langle \Psi_{\alpha} | \mathcal{H} | \Psi_{\alpha} \rangle} {\langle \Psi_{\alpha} | \Psi_{\alpha} \rangle} .\tag{6.9} $$$E_{\alpha}$ 对变分参数的依赖仅仅是波函数 $| \Psi_{\alpha} \rangle$ 依赖于 α 这一事实的结果。因此，为了对 $E_{\alpha}$ 求导，方便的做法是将 $| \Psi_{\alpha} \rangle$ 针对小变化 $\alpha_{k} \ \to \ \alpha_{k} + \delta \alpha_{k}$ 进行展开。对于给定的构型 $| x \rangle$，其中 $\Psi_{\alpha} ( x )$ 是一个复数（在参数为复数的情况下，一旦我们将其实部和虚部分开考虑，就可以假设所有 $\alpha_{k}$ 都是实数），我们有：
$$ \Psi_{\alpha + \delta \alpha_{k}} ( x ) = \Psi_{\alpha} ( x ) + \delta \alpha_{k} \frac{\partial \Psi_{\alpha} ( x )} {\partial \alpha_{k}} + O ( \delta \alpha_{k}^{2} ) ,\tag{6.10} $$其中符号 $\Psi_{\alpha + \delta \alpha_{k}} ( x )$ 表示只有向量 α 中的分量 $\alpha_{k}$ 增加了 $\delta \alpha_{k}$ 。在下文中，为简单起见，我们假设对于所有构型 $\Psi_{\alpha} ( x ) ~ \neq ~ 0$ 。对于连续空间中的费米子系统，节点区域 $\Psi_{\alpha} ( x ) = 0$ 代表一个可忽略（即测度为零）的积分域。在格点上，可以通过对变分试探波函数施加微小扰动（例如添加一个小噪声项）并考虑扰动趋近于零的极限，来移除 $\Psi_{\alpha} ( x ) = 0$ 的偶然构型。于是，方程(6.10)可以形式上写为对应参数 $\alpha_{k}$ 的局部算子 $\mathcal{O}_{k}$ ，该算子由对角矩阵元 ${\mathcal{O}}_{k} ( x )$ 定义：
$$ \langle x | \mathcal{O}_{k} | x^{\prime} \rangle = \delta_{x , x^{\prime}} \mathcal{O}_{k} ( x ) ,\tag{6.11} $$$$ \mathcal{O}_{k} ( x ) = \frac{\partial \ln \Psi_{\alpha} ( x )} {\partial \alpha_{k}} = \frac{1} {\Psi_{\alpha} ( x )} \frac{\partial \Psi_{\alpha} ( x )} {\partial \alpha_{k}} ;\tag{6.12} $$这里，原则上 ${\mathcal{O}}_{k} ( x )$ 可能依赖于变分参数 ${\pmb{\alpha}}$ ；然而，为使符号简单，我们更倾向于不在局部算子中标出 α。
重要的是，${\mathcal{O}}_{k} ( x )$ 通常可以根据变分态的具体试探形式计算出来（关于 Jastrow-Slater 波函数，见第6.7节）。这样，我们可以将多体态的形式展开写为：
$$ | \Psi_{\alpha + \delta \alpha_{k}} \rangle = ( 1 + \delta \alpha_{k} \mathcal{O}_{k} ) | \Psi_{\alpha} \rangle ,\tag{6.13} $$这可以通过取上述方程两侧与$| x \rangle$的内积，并利用式(6.11)和(6.12)轻松验证。注意，对角算符$\mathcal{O}_{k}$不一定是厄米的，因为对于一般的复情况，其对角元不一定是实数。
现在让我们展示如何得到能量相对于给定变分参数的导数的显式形式。显然，变分能量$E_{\alpha}$（以及任何其他关联函数）不依赖于波函数的整体归一化（和全局相位）。换句话说，将波函数乘以任意复数常数c，即$| \Psi_{\alpha} \rangle \rightarrow c | \Psi_{\alpha} \rangle$，$E_{\alpha}$保持不变。为了利用这一性质，最好考虑显式归一化的波函数。首先我们定义：
$$ | \nu_{0 , \alpha} \rangle \equiv \frac{| \Psi_{\alpha} \rangle} {| | \Psi_{\alpha} | |} ,\tag{6.14} $$其中$| | \Psi_{\alpha} | |$表示态$| \Psi_{\alpha} \rangle$的范数。然后，我们定义一组态（对应于每个$k = 1 , \hdots , p$值）：
$$ | \nu_{k , \alpha} \rangle \equiv \left( \mathcal{O}_{k} - \overline{{\mathcal{O}}}_{k} \right) | \nu_{0 , \alpha} \rangle ,\tag{6.15} $$其中
$$ \overline{{\mathcal{O}}}_{k} = \langle \nu_{0 , \alpha} \vert \mathcal{O}_{k} \vert \nu_{0 , \alpha} \rangle = \frac{\langle \Psi_{\alpha} \vert \mathcal{O}_{k} \vert \Psi_{\alpha} \rangle} {\langle \Psi_{\alpha} \vert \Psi_{\alpha} \rangle} .\tag{6.16} $$态$| \nu_{k , \alpha} \rangle$与$| \nu_{0 , \alpha} \rangle$正交，这在使用式(6.16)时很容易验证；然而，它们既不归一化也不相互正交，即一般来说对于$\boldsymbol{k} , \boldsymbol{k}^{\prime} \neq 0$，有$\langle \nu_{k , \alpha} | \nu_{k^{\prime} , \alpha} \rangle \neq \delta_{k , k^{\prime}}$。因此，态集$| \nu_{0 , \alpha} \rangle$和$\{| \nu_{k , \alpha} \rangle \}$定义了一个半正交基。
为了在参数$\alpha_{k}$改变时计算归一化波函数，我们首先计算$| \Psi_{\alpha + \delta \alpha_{k}} \rangle$的范数：
$$ \begin{array} {r l r} & {} & {| | \Psi_{\alpha + \delta \alpha_{k}} | | ^{2} = \langle \Psi_{\alpha} | ( 1 + \delta \alpha_{k} \mathcal{O}_{k} )^{*} ( 1 + \delta \alpha_{k} \mathcal{O}_{k} ) | \Psi_{\alpha} \rangle} \\ & {} & {= | | \Psi_{\alpha} | | ^{2} \left[ 1 + 2 \Re ( \delta \alpha_{k} \overline{{\mathcal{O}}}_{k} ) + O ( \delta \alpha_{k}^{2} ) \right] .} \end{array}\tag{6.17} $$于是，我们有：
$$ \begin{array} {l} {| \nu_{0 , \alpha + \delta \alpha_{k}} \rangle = \displaystyle \frac{| \Psi_{\alpha + \delta \alpha_{k}} \rangle} {| | \Psi_{\alpha + \delta \alpha_{k}} | |} = | \nu_{0 , \alpha} \rangle + \left[ \delta \alpha_{k} \mathcal{O}_{k} - \Re ( \delta \alpha_{k} \overline{{\mathcal{O}}}_{k} ) \right] | \nu_{0 , \alpha} \rangle + O ( \delta \alpha_{k}^{2} )} \\ {\displaystyle \qquad = \left[ 1 + i \Im ( \delta \alpha_{k} \overline{{\mathcal{O}}}_{k} ) \right] | \nu_{0 , \alpha} \rangle + \delta \alpha_{k} | \nu_{k , \alpha} \rangle + O ( \delta \alpha_{k}^{2} ) , \qquad (} \end{array}\tag{6.18} $$这最终可以重新写成：
$$ \left| \nu_{0 , \alpha + \delta \alpha_{k}} \right. = \exp ( i \delta \phi ) \left[ | \nu_{0 , \alpha} \rangle + \delta \alpha_{k} | \nu_{k , \alpha} \rangle \right] + O ( \delta \alpha_{k}^{2} ) ,\tag{6.19} $$其中$\delta \phi = \Im ( \delta \alpha_{k} {\overline{{\mathcal{O}}}}_{k} )$。
利用上述表达式，可以立即得到变分能量$E_{\alpha}$相对于给定变分参数$\alpha_{k}$的导数：
$$ \begin{array} {r l} & {\frac{\partial E_{\alpha}} {\partial \alpha_{k}} = \underset{\delta \alpha_{k} \rightarrow 0} {\operatorname*{lim}} \frac{\langle \nu_{0 , \alpha + \delta \alpha_{k}} | \mathcal{H} | \nu_{0 , \alpha + \delta \alpha_{k}} \rangle - \langle \nu_{0 , \alpha} | \mathcal{H} | \nu_{0 , \alpha} \rangle} {\delta \alpha_{k}}} \\ & {\qquad = \langle \nu_{k , \alpha} | \mathcal{H} | \nu_{0 , \alpha} \rangle + \langle \nu_{0 , \alpha} | \mathcal{H} | \nu_{k , \alpha} \rangle = 2 \Re \Bigg [ \frac{\langle \Psi_{\alpha} | \mathcal{H} ( \mathcal{O}_{k} - \overline{{O}}_{k} ) | \Psi_{\alpha} \rangle} {\langle \Psi_{\alpha} | \Psi_{\alpha} \rangle} \Bigg ] .} \end{array}\tag{6.20} $$同时注意到，相位因子 $\delta \phi$ 如预期那样未出现在上述表达式中。
为了通过标准蒙特卡洛采样计算力 $f_{k}$，我们引入完备性关系得到：
$$ \begin{array} {r l r} {{f_{k} = - 2 \Re [ \frac{\sum_{x} \langle \Psi_{\alpha} | \mathcal{H} | x \rangle \langle x | ( \mathcal{O}_{k} - \overline{{\mathcal{O}}}_{k} ) | \Psi_{\alpha} \rangle} {\sum_{x} \langle \Psi_{\alpha} | x \rangle \langle x | \Psi_{\alpha} \rangle} ]}} \\ & {} & {= - 2 \Re [ \sum_{x} \frac{e_{L}^{*} ( x ) ( \mathcal{O}_{k} ( x ) - \overline{{\mathcal{O}}}_{k} ) | \Psi_{\alpha} ( x ) | ^{2}} {\sum_{x} | \Psi_{\alpha} ( x ) | ^{2}} ] ,} \end{array}\tag{6.21} $$其中 $e_{L}^{*} ( x )$ 是局部能量的复共轭（此处，为了与局部算符 $\mathcal{O}_{k}$ 采用的符号一致，我们省略了指数 $\alpha$），因为它通常是一个复值函数。那么，$f_{k}$ 可以通过下式估算：
$$ f_{k} \approx - 2 \Re \left[ \frac{1} {N} \sum_{i = 1}^{N} e_{L}^{*} ( x_{i} ) \left( \mathcal{O}_{k} ( x_{i} ) - \overline{{\mathcal{O}}}_{k} \right) \right] ,\tag{6.22} $$$$ \overline{{\mathcal{O}}}_{k} \approx \frac{1} {N} \sum_{i = 1}^{N} \mathcal{O}_{k} ( x_{i} ) .\tag{6.23} $$我们指出，在评估误差棒时，利用期望值 $\overline{{\mathcal{O}}}_{k}$ 与实际导数 $f_{k}$ 之间的相关性非常重要，这些导数额外必须使用同一组样本 $\{x_{i} \}$ 进行计算。
最后，我们强调上述导数表达式的两个主要性质，它们构成了任何有效随机最小化方法的基础：
由式 (6.22) 可知，如果波函数是哈密顿量的精确本征态，则局部能量等于相应的精确本征值 E，而与样本 $\{x_{i} \}$ 无关。这意味着 $f_{k}$ 完全为零，没有统计涨落，因此我们恢复了能量导数的零方差性质。
在一般情况下，方差始终受粒子数平方 $N_{p}^{2}$ 的限制，因为局部能量是一个广延随机变量，它与一个导数估计量（即 ${\mathcal{O}}_{k} ( x )$）相乘，而该估计量也至多是广延的。因此，两个广延量乘积的方差至多 $O ( N_{p}^{2} )$。这意味着这些表达式可用于优化具有大量粒子的波函数，因为与基于重加权技术的方法不同，估算能量导数的统计涨落随粒子数的增加至多以多项式增长。
6.4 随机重构#
式 (6.9) 中能量导数的知识使我们能够采用最速下降法 (Press et al., 2007) 来改变变分参数 $\alpha \_ =$ $( \alpha_{1} , \ldots , \alpha_{p} )$，即使 $p$ 非常大：
$$ \alpha_{k}^{\prime} = \alpha_{k} + \delta \alpha_{k} ,\tag{6.24} $$$$ \delta \alpha_{k} = \Delta f_{k} ,\tag{6.25} $$其中 $\Delta$ 是一个任意（小）常数。原则上，可以优化其值以在每次迭代中达到尽可能低的能量；然而，在大多数应用中，通常的做法是在整个最小化过程中保持 $\Delta$ 不变。然后，沿着马尔可夫链过程迭代改进变分参数。在没有噪声的情况下，最速下降法总是收敛到极小值，此时满足欧拉条件 $f_{k} = 0$。实际上，假设 $f_{k} \neq 0$，则 ${\pmb{\alpha}}^{\prime}$ 的能量由泰勒展开至 $\Delta$ 的一阶项给出：
$$ E_{\alpha^{\prime}} = E_{\alpha} + \sum_{k} {\frac{\partial E_{\alpha}} {\partial \alpha_{k}} \delta \alpha_{k}} + O ( \Delta^{2} ) = E_{\alpha} - \Delta \sum_{k} {f_{k}^{2}} + O ( \Delta^{2} ) ,\tag{6.26} $$这里我们使用了 $\partial E_{\alpha} / \partial \alpha_{k} = - f_{k}$ 和 $\delta \alpha_{k} = \Delta f_{k}$。因此，对于小的 $\Delta$，当泰勒展开的线性截断足够精确时，我们得到：
$$ \Delta E \equiv E_{\alpha^{\prime}} - E_{\alpha} = - \Delta \sum_{k} f_{k}^{2} \le 0 ;\tag{6.27} $$这里，等式成立仅当 $f_{k} = 0$。因此，该方法之所以能在大量迭代次数下收敛到最小值，正是因为能量随迭代次数单调递减。在最速下降法中，只计算能量的一阶导数，并且可以确信，平行于力 $\mathbf{f} = ( f_{1} , \ldots , f_{p} )$ 的参数微小变化 $\delta \pmb{\alpha} = ( \delta \alpha_{1} , \ldots , \delta \alpha_{p} )$ 会降低能量；唯一的问题涉及 $\Delta$ 的大小，必须取足够小的 $\Delta$ 以使方程(6.26)中的二次项可以忽略不计。
让我们讨论一个简单情况，其中所有参数都是“等价的”，并考虑一个由 $p$ 个正交归一 Slater 行列式 $| D_{k} \rangle ( \text{即} , \langle D_{k} | D_{k^{\prime}} \rangle = \delta_{k , k^{\prime}} )$ 线性组合而成的波函数：
$$ | \Psi_{\boldsymbol{\alpha}} \rangle = \sum_{k} \alpha_{k} | D_{k} \rangle ;\tag{6.28} $$那么，合理的做法是在满足以下条件下进行搜索：
$$ \delta s^{2} = \sum_{k} \delta \alpha_{k}^{2}\tag{6.29} $$在优化过程中足够小，这里假设对于小的 δs，能量对 δα 的一阶泰勒展开是足够精确的。为了强制执行方程(6.29)的约束，通常的方案是引入一个拉格朗日乘子 $\mu$ 并最小化二次型：
$$ \Delta E + \mu \delta s^{2} = \sum_{k} \left( - f_{k} \delta \alpha_{k} + \mu \delta \alpha_{k}^{2} \right) ;\tag{6.30} $$该方法给出了一个最小值条件：
$$ \delta \alpha_{k} = \frac{f_{k}} {2 \mu} ,\tag{6.31} $$这与方程(6.25)中 $\Delta = 1 / 2 \mu$ 的最速下降法兼容。虽然上述论证对于简单波函数（例如线性依赖于参数的波函数）是有意义的，但在强关联态中会出现严重困难，因为变分参数的依赖性高度非线性（例如，在 Jastrow 因子中）。在这里，给定变分参数的微小变化可能产生非常不同的波函数和物理量，而另一个参数可能对波函数影响很弱。为了克服这些困难，引入更合适的度量 $\delta s^{2}$ 来估计两个归一化（复）波函数 $| \nu_{0 , \alpha} \rangle$ 和 $| \nu_{0 , \alpha + \delta \alpha} \rangle$ 的“邻近性”非常重要：
$$ \begin{array} {r} {\delta s^{2} = \mathrm{Min}_{\delta \theta} | | \exp ( - i \delta \theta ) \nu_{0 , \alpha + \delta \alpha} - \nu_{0 , \alpha} | | ^{2} .} \end{array}\tag{6.32} $$在上述 $\delta s^{2}$ 的定义中，对相位因子 δθ 的最小化是必要的，因为我们不希望区分仅相差一个整体相位因子的两个波函数，因为它们产生相同的关联函数。换句话说，我们希望定义一个距离 $\delta s^{2}$，当波函数物理等价时它为零。然后，我们将方程(6.32)中的 $| \nu_{0 , \alpha + \delta \alpha} \rangle$ 表达式替换为通过将方程(6.19)推广到多个参数同时变化的情况而得到的表达式：
$$ | \nu_{0 , \alpha + \delta \alpha} \rangle = \exp ( i \delta \phi ) \left[ | \nu_{0 , \alpha} \rangle + \sum_{k} \delta \alpha_{k} | \nu_{k , \alpha} \rangle \right] + O ( | \delta \alpha | ^{2} ) ,\tag{6.33} $$其中 $\begin{array} {r} {\delta \phi = \sum_{k} \Im ( \delta \alpha_{k} \overline{{\mathcal{O}}}_{k} )} \end{array}$。现在，对 $\delta \theta$ 的最小化给出了 $\delta \theta \ : = \ : \delta \phi$，从而得到：
$$ \delta s^{2} = \sum_{k , k^{\prime}} \langle \nu_{k , \alpha} | \nu_{k^{\prime} , \alpha} \rangle \delta \alpha_{k} \delta \alpha_{k^{\prime}} + o \left( | \delta \alpha | ^{2} \right) .\tag{6.34} $$由于所有增量 $\delta \alpha_{k}$ 均假设为实数（如前所述，此处假设所有参数均为实数），我们对关于指标 k 和 $k^{\prime}$ 的表达式进行对称化处理，并忽略 $o ( | \delta \alpha | ^{2} )$ 项，得到：
$$ \delta s^{2} = \sum_{k , k^{\prime}} \left( \frac{\langle \nu_{k , \alpha} | \nu_{k^{\prime} , \alpha} \rangle + \langle \nu_{k^{\prime} , \alpha} | \nu_{k , \alpha} \rangle} {2} \right) \delta \alpha_{k} \delta \alpha_{k^{\prime}} .\tag{6.35} $$这样一来，我们最终可以确定一个矩阵 S，它完全决定了归一化波函数空间中的度规：
$$ S_{k , k^{\prime}} = \Re \big ( \langle \nu_{k , \alpha} | \nu_{k^{\prime} , \alpha} \rangle \big ) ,\tag{6.36} $$这意味着两个波函数之间的距离为：
$$ \delta s^{2} = \sum_{k , k^{\prime}} S_{k , k^{\prime}} \delta \alpha_{k} \delta \alpha_{k^{\prime}} .\tag{6.37} $$此时，很自然地利用 S 给出的度规来改进最速下降法。使用式 (6.37) 中给出的度规 δs² 来最小化 $\Delta E + \mu \delta s^{2}$ 能够改善向变分能量最小值的收敛性，相比简单的最速下降法，非等价参数可以以不同的尺度进行适当调整。这种方法称为随机重构 (Sorella, 1998 和 2001)。最小化 $\Delta E + \mu \delta s^{2}$ 给出：
$$ \sum_{k^{\prime}} S_{k , k^{\prime}} \delta \alpha_{k^{\prime}} = \frac{f_{k}} {2 \mu} ,\tag{6.38} $$这是关于未知向量 δα 的一组线性方程。求解这个线性系统后，我们可以更新变分参数直至收敛；如同最速下降法，我们可以将 $\Delta = 1 / ( 2 \mu )$ 设置得足够小，并在优化过程中保持不变。我们想强调一点，由于矩阵 S 是严格正定的，能量在优化过程中单调递减，具体如下：
$$ \Delta E = - \Delta \sum_{k , k^{\prime}} S_{k , k^{\prime}}^{- 1} f_{k} f_{k^{\prime}} < 0 .\tag{6.39} $$在蒙特卡洛过程中，矩阵 S 通过对 N 个构型 $\{x_{i} \}$ 的有限采样来估计，形式为：
$$ S_{k , k^{\prime}} \approx \Re \left[ \frac{1} {N} \sum_{i = 1}^{N} \left( \mathcal{O}_{k} ( x_{i} ) - \overline{{\mathcal{O}}}_{k} \right) \left( \mathcal{O}_{k^{\prime}} ( x_{i} ) - \overline{{\mathcal{O}}}_{k^{\prime}} \right) \right] ;\tag{6.40} $$力也以类似方式计算，参见式 (6.22)。因此，线性系统 (6.38) 的解受到统计误差的影响，即使在达到收敛后，也会导致最终变分参数 $\{\alpha_{k} \}$ 的统计波动。在这种情况下，进行多次迭代是方便的，以便通过平衡后的平均来获得变分参数的准确值。
通过考虑 $\Delta 0$ 的极限，优化过程简化为一组近似由朗之万动力学描述的随机微分方程（见第4章）：
$$ \frac{d \pmb{\alpha} ( t )} {dt} = \mathbf{S}^{- 1} \mathbf{f} [ \pmb{\alpha} ( t ) ] + \pmb{\eta} ( t ) ,\tag{6.41} $$其中我们明确指出力依赖于变分参数。此处，我们忽略统计噪声中力分量之间的相关性，因为这不会改变定性讨论。噪声由以下公式描述：
$$ \langle \eta_{k} ( t ) \rangle = 0 ,\tag{6.42} $$$$ \langle \eta_{k} ( t ) \eta_{k^{\prime}} ( t^{\prime} ) \rangle = 2 T_{\mathrm{noise}} \delta_{k , k^{\prime}} \delta ( t - t^{\prime} ) .\tag{6.43} $$在蒙特卡罗优化中，参数 $T_{\mathrm{noise}}$ 取决于 f 和 S 均通过使用 N 个构型的有限采样进行评估这一事实。我们想强调，矩阵 S 的求逆会引入可控偏差，该偏差在 $N \to \infty$ 时消失。这是因为由方程 (6.40) 得到的矩阵元 $S_{k , k^{\prime}}$ 是围绕相应精确值分布的随机变量，其涨落为 $O ( 1 / \sqrt{N} )$。因此，根据中心极限定理，通过增加 N（即 $T_{\mathrm{noise}} \propto 1 / N$）可以将噪声调至零。达到收敛所需的迭代次数对“温度” $T_{\mathrm{noise}}$ 的依赖性较弱（参见第4章）；但它严重依赖于系统的能量景观。N 的最优值是能够使参数涨落在所需精度内的最小值。需要注意的是，在朗之万模拟中平均的变分参数接近真实能量最小值，但可能受随 $T_{\mathrm{noise}}$ 趋于零的偏差影响，这是由于能量景观中存在非二次项。因此，获得精确结果的最佳策略是通过系统性地增加 N 来平均变分参数。
图6.1 半填充（即 $N_{e} = L$）下 $L = 162$ 格点，$U / t = - 1$ 的吸引哈伯德模型的结果。波函数是通过对具有位点s波配对 $\Delta_{k} = \Delta_{s}$ 的BCS波函数 (1.91) 施加密度-密度Jastrow因子 (1.65) 得到的；Jastrow赝势的数量为29。上图中：利用随机重配置，$\Delta_{s}$ 作为蒙特卡罗步数的函数进行优化，此处 $N = 5000$。下图中：固定所有变分参数计算的变分能量，对应最优值 $\Delta_{s} = 0.04$（左图）和 $\Delta_{s} = 0$（右图）。
作为随机重配置收敛性的一个例子，我们在图6.1中报告了半填充下 $U / t = - 1$ 的吸引哈伯德模型的情况，其中 $N_{e} = L = 162$。变分波函数由具有位点s波配对 $\Delta_{k} = \Delta_{s}$ 的BCS态 (1.91) 给出，并辅以密度-密度Jastrow因子 (1.65)。对于该电子-电子相互作用值，稳定了一个小的 $\Delta_{s}$ 值，即 $\Delta_{s} \approx 0.04$。我们想强调，此项带来的能量增益为 $\Delta E = 0.00007 ( 2 )$，远小于模拟中存在的能量涨落；尽管如此，利用随机优化仍可能精确确定该变分参数，这展示了此类蒙特卡罗方法的能力。
6.4.1 协方差性质#
随机重配置优化的一个重要性质是，当多体态的参量化方式改变但形式不变时，它具有相同的效率。这一事实代表了相对于简单最速下降法的重要改进。例如，让我们考虑波函数包含Gutzwiller因子 $g$ 和逸度项 $f$ 的情况（当模拟中费米子数量不固定时适用）；对于费米子情况，其中 $n_{i} = n_{i , \uparrow} + n_{i , \downarrow}$ 且 $n_{i \sigma}^{2} = n_{i \sigma}$，关联项可以用两种等价方式写出：
$$ \exp \left( - \sum_{i} \frac{g} {2} n_{i}^{2} - \sum_{i} f n_{i} \right) = \exp \left( - \sum_{i} f^{\prime} n_{i} \right) \prod_{i} \left[ 1 + ( g^{\prime} - 1 ) n_{i \uparrow} n_{i \downarrow} \right] ,\tag{6.44} $$其中 $g^{\prime} = \exp ( - g )$ 和 $f^{\prime} = f + g / 2$。选择一种参数化（用 $g$ 和 $f$）或另一种（用 $g^{\prime}$ 和 $f^{\prime}$）不一定在优化过程中给出相同的效率（尽管在平衡时，除了与局部极小值相关的问题外，它们必须给出相同的能量）。在此，我们证明，在随机重构技术中，当改变量子态的参数化方式时，最小化过程具有相同的有效性。
一般情况下，考虑变分参数的任意变换：
$$ \alpha_{k}^{\prime} = \alpha_{k}^{\prime} ( \alpha_{1} , \ldots , \alpha_{p} ) ,\tag{6.45} $$该变换给出相同的量子态，即 $\Psi_{\alpha^{\prime}} ( x ) = \Psi_{\alpha} ( x )$。然后，我们希望施加条件：在优化一步之后，使用 $\alpha$ 或 $\pmb{\alpha}^{\prime}$ 的参数化方式应获得相同的改进。条件 $\Psi_{\alpha^{\prime} + \delta \alpha^{\prime}} ( x ) = \Psi_{\alpha + \delta \alpha} ( x )$ 要求：
$$ \delta \alpha_{k}^{\prime} = \sum_{k^{\prime}} \frac{\partial \alpha_{k}^{\prime}} {\partial \alpha_{k^{\prime}}} \delta \alpha_{k^{\prime}} = ( \mathbf{J} \delta \pmb{\alpha} )_{k} ,\tag{6.46} $$其中我们定义了雅可比矩阵 $\begin{array} {r} {J_{k , k^{\prime}} = \frac{\partial \alpha_{k}^{\prime}} {\partial \alpha_{k^{\prime}}}} \end{array}$。在随机重构技术中，上述关系确实成立。事实上，参与能量导数计算的局部算符 $\mathcal{O}_{k}^{\prime} ( x )$ 由下式给出：
$$ \mathcal{O}_{k}^{\prime} ( x ) = \frac{1} {\Psi_{\alpha^{\prime}} ( x )} \frac{\partial \Psi_{\alpha^{\prime}} ( x )} {\partial \alpha_{k}^{\prime}} = \sum_{k^{\prime}} \frac{\partial \alpha_{k^{\prime}}} {\partial \alpha_{k}^{\prime}} \left( \frac{1} {\Psi_{\alpha} ( x )} \frac{\partial \Psi_{\alpha} ( x )} {\partial \alpha_{k^{\prime}}} \right) = \sum_{k^{\prime}} J_{k , k^{\prime}}^{- 1} O_{k^{\prime}} ( x ) ,\tag{6.47} $$这里我们利用了 $\Psi_{\alpha^{\prime}} ( x ) = \Psi_{\alpha} ( x )$ 以及逆变换的雅可比矩阵恰好是直接变换的雅可比矩阵的逆矩阵这一性质。这样一来，新的力 $\mathbf{f^{\prime}}$ 和新的矩阵 $\mathbf{S}^{\prime}$ 为：
$$ \mathbf{f^{\prime}} = \mathbf{J}^{- 1} \mathbf{f} ,\tag{6.48} $$$$ \mathbf{S}^{\prime} = \mathbf{J}^{- 1} \mathbf{S} \mathbf{J}^{- 1} .\tag{6.49} $$因此，使用新参数化 ${\pmb{\alpha}}^{\prime}$ 的迭代步骤为：
$$ \begin{array} {r} {\delta \pmb{\alpha}^{\prime} = \Delta ( \mathbf{S}^{\prime} )^{- 1} \mathbf{f}^{\prime} = \Delta \left( \mathbf{J}^{- 1} \mathbf{S} \mathbf{J}^{- 1} \right)^{- 1} \mathbf{J}^{- 1} \mathbf{f} = \mathbf{J} \delta \pmb{\alpha} ,} \end{array}\tag{6.50} $$这正好是我们一开始所要求的条件。这一事实可称为方法的协变性，因为（正如广义相对论中一样）矩阵 S 代表我们非线性空间（波函数空间）的度量，优化的方向是协变的，不依赖于波函数表示的方式。这一性质非常重要，因为它意味着通过参数间的任何非线性变换都无法改进该方法。
6.4.2 尺度不变的规范化方法#
在随机优化过程中，矩阵 S 的求逆是一个基本要素，即使考虑有限采样 N 个构型，它也是一个非负矩阵。然而，其特征值可能相差几个数量级，因此，进行预处理以稳定求逆过程总是值得的，即考虑如下尺度变换：
$$ S_{k , k^{\prime}}^{\mathrm{pc}} = \frac{S_{k , k^{\prime}}} {\sqrt{S_{k , k} S_{k^{\prime} , k^{\prime}}}} .\tag{6.51} $$通过这种方式，我们首先求解线性方程：
$$ \sum_{k^{\prime}} S_{k , k^{\prime}}^{\mathrm{pc}} \delta \alpha_{k^{\prime}}^{\mathrm{pc}} = \frac{f_{k}^{\mathrm{pc}}} {2 \mu} ,\tag{6.52} $$其中缩放后的未知量为 $\delta \alpha_{k^{\prime}}^{\mathrm{pc}} = \sqrt{S_{k^{\prime} , k^{\prime}}} \delta \alpha_{k^{\prime}}$，缩放后的力为 $f_{k}^{\mathrm{pc}} = f_{k} / \sqrt{S_{k , k}}$；然后，通过简单的重新缩放 $\delta \alpha_{k} = \delta \alpha_{k}^{\mathrm{pc}} / \sqrt{S_{k , k}}$ 即可得到原始的 δα。
当矩阵 $\mathbf{S}^{\mathrm{pc}}$ 获得非常小甚至可忽略的特征值时，就会出现一个非常严重的问题，这意味着将 $( \mathbf{S}^{\mathrm{pc}} )^{-1}$ 作用于力 $\mathbf{f}^{\mathrm{pc}}$ 可能是不稳定的，特别是考虑到这两个量都是在随机过程中评估的。例如，如果波函数的参数化存在冗余，则力的某个特定线性组合将为零，而矩阵 $\mathbf{S}^{\mathrm{pc}}$ 也将是奇异的。例如，通过对 $\mathbf{S}^{\mathrm{pc}}$ 进行对角化，我们得到：
$$ \lambda_{i} ~ ( \delta \alpha^{\mathrm{pc}} )_{i}^{U} = ( f^{\mathrm{pc}} )_{i}^{U} ,\tag{6.53} $$其中 $\lambda_{i}$ 是 $\mathbf{S}^{\mathrm{pc}} = \mathbf{U}^{\dagger} \boldsymbol{\Lambda} \mathbf{U}$ 的特征值，$( \pmb{\delta \alpha}^{\mathrm{pc}} )^{U} = \mathbf{U} \pmb{\delta \alpha}^{\mathrm{pc}}$，以及 $( \mathbf{f}^{\mathrm{pc}} )^{U} = \mathbf{U} \mathbf{f}^{\mathrm{pc}}$ 对于冗余参数，$( f^{\mathrm{pc}} )_{k}^{U} = \lambda_{k} = 0$，这使得问题的解不确定。因此，显然需要某种形式的正则化。一个简单而有效的可能性是仅修改对角元：
$$ S_{k , k^{\prime}}^{\mathrm{pc} , \epsilon} = S_{k , k^{\prime}}^{\mathrm{pc}} + \epsilon \delta_{k , k^{\prime}} .\tag{6.54} $$这定义了一个正则矩阵，对于任意 $\epsilon > 0$ 都可以安全地求逆；这样，我们就消除了变分参数之间可能的依赖关系。实际上，式(6.53)被修改为：
$$ \begin{array} {r} {( \lambda_{i} + \epsilon ) ( \delta \alpha^{\mathrm{pc}} )_{i}^{U} = ( f^{\mathrm{pc}} )_{i}^{U} ,} \end{array}\tag{6.55} $$这意味着，当 $( f^{\mathrm{pc}} )_{i}^{U} \approx 0$ 且 $\lambda_{i} \ll \epsilon$ 时，$( \delta \alpha^{\mathrm{pc}} )_{i}^{U} \approx 0$。因此，引入有限 $\epsilon$ 的效果是保持 $( f^{\mathrm{pc}} )_{i}^{U} \approx 0$ 的参数方向固定。由于 $( f^{\mathrm{pc}} )_{i}^{U} \approx \sqrt{\lambda_{i}}$，忽略这个方向会导致优化误差约为 $\sqrt{\epsilon}$ 的量级。
尽管这种正则化不能保持上述的协方差性质，但它保留了在通过任意常数 $\xi_{k}$ 缩放参数时获得的变换子集：
$$ \alpha_{k}^{\prime} = \xi_{k} \alpha_{k} ;\tag{6.56} $$事实上，我们可以立即验证这一性质，因为在这种情况下 J 是对角矩阵，因此式(6.51)中的 $\mathbf{S}^{\mathrm{pc}}$ 在对矩阵 S 进行变换(6.49)时保持不变。因此，这种正则化是尺度不变的。 $\epsilon$ 的值越小，趋近最低变分能量的速度就越快，但出现不稳定的概率也越高（如果矩阵 S 是奇异的）；通常的做法是使用 $\epsilon \approx 10^{-3}$，以获得合理的效率和精度。
6.4.3 利用信噪比的优化#
至此，我们想讨论一个非常简单的论点，以解释上述方法为何是一种高效的随机优化过程。广义力分量 $f_{k}$（定义于式(6.22)）是一个随机变量，具有相应的噪声，其值可以与其他分量的相关性一起通过统计方法估算。实际上，协方差矩阵可以通过第3.10节中描述的刀切法或自助法来评估：
$$ \sigma_{k , k^{\prime}}^{2} \approx \langle \langle f_{k} f_{k^{\prime}} \rangle \rangle - \langle \langle f_{k} \rangle \rangle \langle \langle f_{k^{\prime}} \rangle \rangle .\tag{6.57} $$现在让我们考虑参数空间中的一个特定方向，该方向由 $\pmb{\alpha} ( \Delta ) = \pmb{\alpha}_{0} + \Delta \pmb{\tau}$ 描述，其中 $\Delta$ 参数化了参数沿向量 $\boldsymbol{\tau}$ 方向从 ${\pmb{\alpha}}_{0}$ 出发的线性变化。这意味着沿方向 $\boldsymbol{\tau}$ 的广义力 $f_{\tau}$ 由下式给出：
$$ f_{\boldsymbol{\tau}} = - \frac{\partial E_{\alpha}} {\partial \Delta} = - \sum_{\boldsymbol{k}} \frac{\partial E_{\alpha}} {\partial \alpha_{k}} \frac{\partial \alpha_{k}} {\partial \Delta} = \mathbf{f} \cdot \boldsymbol{\tau} .\tag{6.58} $$其标准差为：
$$ \sigma_{\tau} = \sqrt{\sum_{k , k^{\prime}} \tau_{k} \sigma_{k , k^{\prime}}^{2} \tau_{k^{\prime}}} .\tag{6.59} $$现在，我们可以寻找这样一个方向：在该方向上，信号（即导数的实际值）与其标准差相比尽可能大。如果信噪比的平方
$$ \Sigma^{2} ( \tau ) = \frac{\sum_{k , k^{\prime}} \tau_{k} f_{k} f_{k^{\prime}} \tau_{k^{\prime}}} {\sum_{k , k^{\prime}} \tau_{k} \sigma_{k , k^{\prime}}^{2} \tau_{k^{\prime}}}\tag{6.60} $$远大于1（例如，如果沿给定方向的导数非零且超过3个标准差），则沿此方向改变参数可以保证能量降低。在这种情况下，通过将 $\pmb{\alpha}_{0}$ 改变为 $\pmb{\alpha}_{0} + \Delta \pmb{\tau}$（其中 $\Delta$ 足够小）将降低能量。从这个意义上说，一个简单的标准最大化过程得出“最佳方向”由下式给出：
$$ \tau_{k} = \Delta_{\tau} \sum_{k^{\prime}} \left( \sigma^{2} \right)_{k , k^{\prime}}^{- 1} f_{k^{\prime}} ,\tag{6.61} $$其中 $\Delta_{\tau}$ 是一个任意常数，因为 $\Sigma^{2} ( \pmb{\tau} )$ 的最大值在按任意常数缩放解时不会改变。在实际应用中，对于精确的变分波函数，期望协方差矩阵接近于矩阵 $\mathbf{S}$，因为力分量 $f_{k}$ 的涨落主要由 $( \mathcal{O}_{k} - \overline{{\mathcal{O}}}_{k} )$ 的涨落给出，而局域能量的涨落可以忽略，这导致：
$$ \sigma_{k , k^{\prime}}^{2} \simeq S_{k , k^{\prime}} .\tag{6.62} $$因此，方程(6.61)和(6.62)解释了为什么随机重构方法是一种特别有效的执行蒙特卡洛优化的方式。事实上，根据方程(6.38)进行的极小化（近似地）遵循了最大化信噪比的方向，这对于稳定的随机优化至关重要。
6.5 作为投影技术的随机重构#
在这里，我们证明，对于实波函数的情况，随机重构方法等价于一种虚时投影技术，该技术被限制在由 $| \nu_{0 , \alpha} \rangle$ 和 $\{| \nu_{k , \alpha} \rangle \}$ 定义的子空间中，如方程(6.14)和(6.15)所定义。为了证明这一论断，让我们定义一个波函数，它由变分波函数 $| \Psi_{\alpha} \rangle$ 通过应用单步幂法（见第1.7节）得到：
$$ | \Psi_{\Delta} \rangle = ( 1 - \Delta \mathcal{H} ) | \Psi_{\alpha} \rangle ,\tag{6.63} $$其中 $\Delta$ 是一个小常数，用于相对于 $| \Psi_{\alpha} \rangle$ 改善能量。这里我们考虑虚时演化（即 $\Delta$ 为实数），但所有结果也适用于实时传播（即 $\Delta \to i \Delta$），这为在变分蒙特卡洛框架内模拟量子动力学奠定了基础(Carleo等人，2012)。然而，在这种情况下，我们不能假设变分参数是实数，因为经过动力学的前几步后，它们就会变为复数。因此，我们需考虑一般的复数量子态，它们是变分参数的解析函数。然后很容易证明，即使考虑复变分 $\{\delta \alpha_{k}\}$，方程(6.14)、(6.15)和(6.19)仍然保持不变。一般来说，对于 $\Delta$ 的任意有限值，态 $| \Psi_{\Delta} \rangle$ 将不同于任何可用原始变分形式参数化的态，即对于任何参数选择 ${\pmb{\alpha}}^{\prime}$，都有 $| \Psi_{\Delta} \rangle \neq | \Psi_{\alpha^{\prime}} \rangle$。然而，当 $\Delta$ 足够小时，$| \Psi_{\Delta} \rangle$ 非常接近 $| \Psi_{\alpha} \rangle$，我们可以探讨是否可以通过参数化 ${\pmb{\alpha}}^{\prime} = {\pmb{\alpha}} + \delta{\pmb{\alpha}}$ 来精确表示 $\left| \Psi_{\Delta} \right.$。特别地，我们可以最小化 $\left| \Psi_{\Delta} \right.$ 与 $e^{i \delta \theta} | \Psi_{\alpha + \delta \alpha} \rangle$（其中 $\delta \theta$ 是一个通用角度，需考虑实时传播）之间的“距离”。使用方程(6.14)和(6.15)的记号，我们得到：
$$ | \Psi_{\Delta} \rangle = | | \Psi_{\alpha} | | ( 1 - \Delta \mathcal{H} ) | \nu_{0 , \alpha} \rangle ,\tag{6.64} $$$$ | \Psi_{\alpha + \delta \alpha} \rangle = e^{i \delta \phi} | | \Psi_{\alpha + \delta \alpha} | | \left[ | \nu_{0 , \alpha} \rangle + \sum_{k} \delta \alpha_{k} | | \nu_{k , \alpha} \rangle \right] ,\tag{6.65} $$其中 $| | \Psi_{\alpha + \delta \alpha} | | = C | | \Psi_{\alpha} | |$ 且 $C \approx 1$，见方程(6.17)。距离的最小化可通过将 $| \Psi_{\Delta} \rangle$ 投影到半正交基张成的子空间来实现：
$$ \langle \nu_{0 , \alpha} | \Psi_{\Delta} \rangle = e^{i \delta \theta} \langle \nu_{0 , \alpha} | \Psi_{\alpha + \delta \alpha} \rangle ,\tag{6.66} $$$$ \langle \nu_{k , \alpha} | \Psi_{\Delta} \rangle = e^{i \delta \theta} \langle \nu_{k , \alpha} | \Psi_{\alpha + \delta \alpha} \rangle ,\tag{6.67} $$这导出了关于 $\delta \alpha_{k}^{\mathbf{\Upsilon}} \mathbf{\dot{s}}$ 的以下线性方程组：
$$ 1 - \Delta \langle \nu_{0 , \alpha} | \mathcal{H} | \nu_{0 , \alpha} \rangle = C e^{i ( \delta \theta + \delta \phi )} ,\tag{6.68} $$$$ - \Delta \langle \nu_{k , \alpha} | \mathcal{H} | \nu_{0 , \alpha} \rangle = C e^{i ( \delta \theta + \delta \phi )} \sum_{k^{\prime}} \langle \nu_{k , \alpha} | \nu_{k^{\prime} , \alpha} \rangle \delta \alpha_{k^{\prime}} .\tag{6.69} $$这里，第一个方程定义了归一化常数 $C = 1 + O ( \Delta )$，这是无关紧要的（对于实数 $\Delta$ 有 $\delta \theta = - \delta \phi$），而第二个方程可以改写为：
$$ \sum_{k^{\prime}} \tilde{S}_{k , k^{\prime}} \delta \alpha_{k^{\prime}} = \frac{\tilde{f}_{k}} {2 \mu} ,\tag{6.70} $$其中 $\mu ~ = ~ C e^{i \theta} / \Delta$，且力分量 $\tilde{f}_{k}$ 和矩阵 $\tilde{S}_{k , k^{\prime}}$ 通常为复数，由下式给出：
$$ \tilde{f}_{k} = - 2 \langle \nu_{k , \alpha} | \mathcal{H} | \nu_{0 , \alpha} \rangle = - 2 \frac{\langle \Psi_{\alpha} | \left( \mathcal{O}_{k} - \overline{{\mathcal{O}}}_{k} \right)^{*} \mathcal{H} | \Psi_{\alpha} \rangle} {\langle \Psi_{\alpha} | \Psi_{\alpha} \rangle} ,\tag{6.71} $$$$ \tilde{S}_{k , k^{\prime}} = \left. \psi_{k , \alpha} \middle \vert \nu_{k^{\prime} , \alpha} \right. = \frac{\left. \Psi_{\alpha} \middle \vert \left( \mathcal{O}_{k} - \overline{{\mathcal{O}}}_{k} \right)^{*} \left( \mathcal{O}_{k^{\prime}} - \overline{{\mathcal{O}}}_{k^{\prime}} \right) \middle \vert \Psi_{\alpha} \right.} {\left. \Psi_{\alpha} \middle \vert \Psi_{\alpha} \right.} .\tag{6.72} $$我们想强调，式(6.70)的投影方法与我们在式(6.38)中推导的最小化算法一致，仅当变分波函数 $\Psi_{\alpha} ( x )$ 为实数时成立。
最后我们指出，即使在投影技术中，矩阵 $\tilde{\bf S}$ 也需要一些正则化处理，以避免其求逆时出现过小的特征值。如果我们希望以 $\Delta^{2}$ 的精度积分短时动力学，那么正则化截断需按 $\epsilon \propto \Delta^{4}$ 进行缩放。这在原则上不是问题，但使得高精度的时间积分变得不可行，因为我们必须将样本数增加至 $N \propto 1 / \epsilon^{2}$ ，才能获得 $O ( \epsilon )$ 的统计精度；否则，误差会大于目标精度。幸运的是，如果我们只关注变分波函数的优化，而不需要像实时量子动力学那样精确描述真实时间传播，那么这个问题并不构成障碍。
6.6 线性方法#
当优化问题中不存在能量尺度差异过大的情况时（例如晶格模型或通过赝势消除了核心电子高能标的从头算电子计算），随机重构方法表现出色。然而，在某些特殊情况下，即使在随机重构方案中，收敛到最小值的速度仍可能极其缓慢。加速任何优化方案的一种成熟方法是考虑能量相对于变分参数的精确二阶导数；这可以通过牛顿-拉夫森方案实现，该方案极其高效，只需很少的步数就能收敛到能量最小值。遗憾的是，黑塞矩阵的精确求值：
$$ H_{k , k^{\prime}} = \frac{\partial^{2} E_{\alpha}} {\partial \alpha_{k} \partial \alpha_{k^{\prime}}}\tag{6.73} $$在实践中相当繁琐。实际上，这需要计算变分波函数在参数α变化下的泰勒展开至二阶：
$$ \Psi_{\alpha + \delta \alpha} ( x ) \approx \Psi_{\alpha} ( x ) + \sum_{k} \delta \alpha_{k} \frac{\partial \Psi_{\alpha} ( x )} {\partial \alpha_{k}} + \frac{1} {2} \sum_{k , k^{\prime}} \delta \alpha_{k} \delta \alpha_{k^{\prime}} \frac{\partial^{2} \Psi_{\alpha} ( x )} {\partial \alpha_{k} \partial \alpha_{k^{\prime}}} ;\tag{6.74} $$如前所述，参见式(6.11)和(6.12)，我们可以用局部算符 $\mathcal{O}_{k , k^{\prime}}$ 来表达二阶导数。
$$ \langle x | \mathcal{O}_{k , k^{\prime}} | x^{\prime} \rangle = \delta_{x , x^{\prime}} \mathcal{O}_{k , k^{\prime}} ( x ) ,\tag{6.75} $$$$ \mathcal{O}_{k , k^{\prime}} ( x ) = \frac{1} {\Psi_{\alpha} ( x )} \frac{\partial^{2} \Psi_{\alpha} ( x )} {\partial \alpha_{k} \partial \alpha_{k^{\prime}}} ;\tag{6.76} $$然而，这些算符的实际求值通常并不容易（尽管对于Jastrow参数而言是可行的）。
与其采用完整的黑塞方法，一个极大的简化出现在我们使用波函数展开至线性阶时（即仅考虑一阶导数），然后计算哈密顿量的期望值，即使超出线性范围：
$$ E_{\alpha + \delta \alpha} = \frac{\sum_{{k , k^{\prime} = 0}}^{p} z_{k}^{*} z_{k^{\prime}} \langle \nu_{k , \alpha} | \mathcal{H} | \nu_{k^{\prime} , \alpha} \rangle} {\sum_{{k , k^{\prime} = 0}}^{p} z_{k}^{*} z_{k^{\prime}} \langle \nu_{k , \alpha} | \nu_{k^{\prime} , \alpha} \rangle} ,\tag{6.77} $$其中我们明确指出求和包含 $k = 0$ 和 $k^{\prime} = 0$ 的项；那么对于 $k > 0$ 有 $z_{k} ~ = ~ \delta \alpha_{k}$ ，而 $z_{0} = 1$ 。这种方法被称为线性方法 (Umrigar et al., 2007; Toulouse and Umrigar, 2007)。在这种简化方法中，我们需要计算矩阵元：
$$ H_{k , k^{\prime}} = \langle \nu_{k , \alpha} | \mathcal{H} | \nu_{k^{\prime} , \alpha} \rangle ,\tag{6.78} $$$$ \overline{{S}}_{k , k^{\prime}} = \langle \nu_{k , \alpha} | \nu_{k^{\prime} , \alpha} \rangle ,\tag{6.79} $$其中 $\overline{{S}}_{0 , 0} = 1$，对于 $k > 0$ 有 $\overline{{S}}_{k , 0} = \overline{{S}}_{0 , k} = 0$，而对于 $k$ 和 $k^{\prime} > 0$ 有 $\overline{{S}}_{k , k^{\prime}} = \langle \nu_{k , \alpha} | \nu_{k^{\prime} , \alpha} \rangle$。然后，可以通过对每个 $z_{k}$ 最小化能量来得到变化量 $\delta \alpha_{k}$。这一过程导出了一个广义本征值方程（因为态 $\{| \nu_{k} \rangle \}$ 之间并不正交）：
$$ \sum_{k^{\prime} = 0}^{p} H_{k , k^{\prime}} z_{k^{\prime}} = E \sum_{k^{\prime} = 0}^{p} \overline{{S}}_{k , k^{\prime}} z_{k^{\prime}} .\tag{6.80} $$在上述广义本征值方程的所有可能右本征向量中，取对应的（实）本征值 E 最低的那个并不总是方便的。实际上，当矩阵特别大或模拟噪声特别大时，最好取与 $| \nu_{0 , \alpha} \rangle$ 具有最大重叠 $| z_{0} | ^{2}$ 的本征向量；这样，修正量 $\delta \alpha_{k} = z_{k} / z_{0}$ 较小，且对 $| \Psi_{\alpha + \delta \alpha} \rangle$ 的展开仍然可控。
矩阵元 $\overline{{S}}_{k , k^{\prime}}$ 和 $H_{k , k^{\prime}}$ 可以通过蒙特卡罗采样来评估：
$$ \overline{{S}}_{k , k^{\prime}} \approx \frac{1} {N} \sum_{i = 1}^{N} \left( \mathcal{O}_{k} ( x_{i} ) - \overline{{\mathcal{O}}}_{k} \right)^{*} \left( \mathcal{O}_{k^{\prime}} ( x_{i} ) - \overline{{\mathcal{O}}}_{k^{\prime}} \right) ,\tag{6.81} $$$$ H_{k , k^{\prime}} \approx \frac{1} {N} \sum_{i = 1}^{N} \left( \mathcal{O}_{k} ( x_{i} ) - \overline{{\mathcal{O}}}_{k} \right)^{\ast} \frac{\langle x_{i} | \mathcal{H} \left( \mathcal{O}_{k^{\prime}} - \overline{{\mathcal{O}}}_{k^{\prime}} \right) | \Psi_{\alpha} \rangle} {\langle x_{i} | \Psi_{\alpha} \rangle} ,\tag{6.82} $$其中我们定义了 $\mathcal{O}_{0} ~ = ~ \mathbb{I}$ 和 $\overline{{\mathcal{O}}}_{0} ~ = ~ 0$，以便包含 $k = 0 \mathrm{or} k^{\prime} = 0$ 的情况。方程(6.82)右侧的因子可以用局域能量的导数来表示，这些导数可以在变分蒙特卡罗技术中高效计算。注意，哈密顿矩阵 H 仅在平均意义上是厄米的，因为其厄米特性被有限配置数样本中存在的统计涨落所破坏。因此，很自然地会考虑对矩阵进行对称化以恢复其厄米性，例如通过考虑 $( H_{k , k^{\prime}} + H_{k^{\prime} , k}^{*} ) / 2$。然而，我们的经验表明，不进行这种对称化而直接求解上述方程更为方便，以维持所谓的“强”零方差性质。实际上，如果半正交基中态的一个特定线性组合
$$ \vert \varphi \rangle = \sum_{k = 0}^{p} z_{k} \left( \mathcal{O}_{k} - \overline{{\mathcal{O}}}_{k} \right) \vert \Psi_{\alpha} \rangle\tag{6.83} $$是系统(6.80)的一个具有本征值 $E_{\varphi}$ 的右本征态，那么我们有：
$$ \begin{array} {r l} & {\displaystyle \sum_{k^{\prime} = 0}^{p} H_{k , k^{\prime}} z_{k^{\prime}} \approx \frac{1} {N} \sum_{i = 1}^{N} ( \mathcal{O}_{k} ( x_{i} ) - \overline{{\mathcal{O}}}_{k} )^{*} \frac{\langle x_{i} | \mathcal{H} \sum_{k^{\prime} = 0}^{p} z_{k^{\prime}} \left( \mathcal{O}_{k^{\prime}} - \overline{{\mathcal{O}}}_{k^{\prime}} \right) | \Psi_{\alpha} \rangle} {\langle x_{i} | \Psi_{\alpha} \rangle}} \\ & {\quad \quad \quad \quad \quad \quad = E_{\varphi} \frac{1} {N} \sum_{i = 1}^{N} ( \mathcal{O}_{k} ( x_{i} ) - \overline{{\mathcal{O}}}_{k} )^{*} \frac{\langle x_{i} | \sum_{k^{\prime} = 0}^{p} z_{k^{\prime}} ( \mathcal{O}_{k^{\prime}} - \overline{{\mathcal{O}}}_{k^{\prime}} ) | \Psi_{\alpha} \rangle} {\langle x_{i} | \Psi_{\alpha} \rangle}} \\ & {\quad \quad \quad \approx E_{\varphi} \sum_{k^{\prime} = 0}^{p} \overline{{S}}_{k , k^{\prime} z_{k^{\prime}}} ,} \end{array}\tag{6.84} $$这意味着式 (6.80) 以零方差得到满足：通过用相同的构型估计 $H_{k , k^{\prime}}$ 和 $\overline{{S}}_{k , k^{\prime}}$，我们得到了线性系统的一个右本征向量，该向量对应于由式 (6.83) 假设给出的精确本征态。这是一个“强”零方差原理，因为即使当 $| \Psi_{\alpha} \rangle$ 不是哈密顿量的精确本征态时，它也能被满足。
我们希望对线性方法做一些最后的思考，该方法在许多情况下可能相对于随机重构方法带来实质性的改进。尽管它的实现比随机重构技术更困难，并且每次迭代需要额外的计算时间，但达到收敛结果所需的总迭代次数通常比随机重构技术大幅减少。遗憾的是，无法量化使用线性方法时预期能节省的实际计算时间。事实上，在量子蒙特卡洛模拟中，即使使用收敛较慢的方法，也能获得稳定、准确且快速的优化，因为每次迭代可能只需要少量样本；因此，用随机重构技术执行数千次迭代并不构成严重问题，特别是当变分参数数量足够少时。
在所有这些优化技术中，为达到优化所涉及矩阵的有意义精度所需的样本数量，应远大于矩阵的主要维度（即参数数量），否则采样将产生秩亏缺的奇异（即“脏”）矩阵。一个实用的经验法则是取样本数 $N$，使得：
$$ N \geq 10 \times p ,\tag{6.85} $$其中 $p$ 是必须优化的参数数量。
作为线性方法的一个示例，我们展示了在 $L = 98$ 格点团簇上半填充 $(N_{e} ~ = ~ L)$ 的排斥性哈伯德模型，其中 $U / t = 4$ 的情况。此时，变分态由 Jastrow-Slater 波函数给出，其中 Jastrow 项包含密度-密度项（见式 (1.65)），而 Slater 行列式由自由电子轨道构成。待优化的变分参数是各独立距离对应的 Jastrow 赝势（包括在位 Gutzwiller 因子）。在图 6.2 中，我们报告了前十次迭代（每次迭代对应 $N = 80000$ 步）的能量和少数 Jastrow 参数（Gutzwiller 项和最小四个距离对应的项）。初始参数设为零。我们想强调，使用线性方法，即使像 $N \approx 4$ 这样非常小的步数也足以收敛到平衡值，对于对应长距离的 Jastrow 项也是如此。
图 6.2 排斥性哈伯德模型在 $U / t = 4$、$L = 98$ 格点半填充（即 $N_{e} = L$ ）下的结果。波函数是通过将密度-密度 Jastrow 因子 (1.65) 应用于由自由电子轨道生成的 Slater 行列式得到的；Jastrow 赝势总数为 19（包括在位 Gutzwiller 项）。使用线性方法，报告了 Jastrow 赝势（Gutzwiller 项和最小 3 距离对应的项）的优化过程（上图）。同时展示了优化过程中的变分能量（下图）。每个优化步的独立蒙特卡罗尝试次数为 $N = 80000$
最后，我们列出一些仍在研究中的问题，这些问题在线性方法中尚未得到解决：
• 当某些参数被限制为实数时（例如，通常我们希望使用实数的 Jastrow 因子），“强”零方差性质不成立。相反，如果计算限制为实数的 Ansatz 和实数哈密顿量，则此方法不存在这一缺陷，因为在这种情况下，$H_{k , k^{\prime}}$ 和 $\overline{{S}}_{k , k^{\prime}}$ 都是实数，相应的特征向量也是实数。
• 在远离最小值时，由 $\delta \alpha_{k} = z_{k} / z_{0}$ 给出的变分参数变化可能相当大，有时甚至会增加能量（因为变化可能超出线性展开有效的范围）。在这种情况下，最好先用随机重构技术开始优化，然后仅在接近最小值时应用线性方法。请注意，随着系统规模的增大，这个问题可能变得越来越突出。
6.7 Jastrow-Slater 情况下的导数计算#
在此，我们展示在 Jastrow-Slater Ansatz 情况下如何计算波函数的导数（参见式 (6.12)）：
$$ | \Psi_{\mathrm{J}} \rangle = \mathcal{I} | \Phi_{0} \rangle ,\tag{6.86} $$其中
$$ \mathcal{T} = \exp \left( - \frac{1} {2} \sum_{i , j} \nu_{i , j} n_{i} n_{j} \right) ,\tag{6.87} $$而 $\left| {\Phi_{0}} \right.$ 是非相互作用 (BCS) 哈密顿量的基态：
$$ {\mathcal{H}}_{0} = \sum_{I , J} t_{I , J} d_{I}^{\dagger} d_{J} .\tag{6.88} $$在下文中，我们考虑平移不变的情况，其中 Jastrow 赝势和 $\mathcal{H}_{0}$ 中的参数仅依赖于格点间的距离，即 $\nu_{i , j} ~ ~ \nu_{k}$ 和 $t_{I , J} \ \to \ t_{k}$ 。推广到平移对称性被破坏的情况是直接的。如果 $\left| {\Phi_{0}} \right.$ 不是通过辅助哈密顿量定义，而是由给定且经过适当参数化的单粒子轨道定义，则可以为 Slater 部分设计不同的方法。
让我们从 Jastrow 因子的情况开始，这非常简单。考虑到平移不变性，我们有：
$$ \mathcal{T} = \exp \left( - \frac{1} {2} \sum_{k} \nu_{k} \sum_{i} n_{i} \sum_{s_{i} ( k )} n_{s_{i} ( k )} \right) ,\tag{6.89} $$其中 $s_{i} ( k )$ 表示与格点 i 距离为 k 的所有格点，例如，在一维空间中，$s_{i} ( k )$ 包括左右两个格点（对于偶数个格点，最大距离对应单个格点）；在二维空间中，连接格点的数量取决于晶格的点群，也取决于距离（例如，在正方晶格上，可以是 8、4、2 和 1）。由式 (6.12) 可得：
$$ \mathcal O_{k} ( x ) = - \frac{1} {2} \frac{\langle x | \left( \sum_{i} n_{i} \sum_{s_{i} ( k )} n_{s_{i} ( k )} \right) | \Psi_{\mathfrak I} \rangle} {\langle x | \Psi_{\mathfrak I} \rangle} = - \frac{1} {2} \langle x | \left( \sum_{i} n_{i} \sum_{s_{i} ( k )} n_{s_{i} ( k )} \right) | x \rangle ,\tag{6.90} $$其中最后一个等式是因为 $n_{i}$ 在基 $\{\vert x \rangle \}$ 中是对角的，该基中电子在晶格上具有确定的位置（和自旋），参见第 5.4 节。因此，导数的计算对应于一组简单的对角关联函数的求值。
接下来，我们考虑对于行列式形式（参见第 5.6 节）的变分 Ansatz，其关于无相互作用态 $\left| {\Phi_{0}} \right.$ 中包含参数的导数。那么，一个微小的变化 $t_{k} \ \ t_{k} + \delta t_{k}$ 对应于辅助哈密顿量的修改，即 $\mathcal{H}_{0} \mathcal{H}_{0} + \delta t_{k} \mathcal{V}^{( k )}$，其中 $\mathcal{V}^{( k )}$ 是一个两体算符，包含所有与 $t_{k}$ 相关的跳跃项。
$$ {\mathcal V}^{( k )} = \sum_{I , J} {V}_{I , J}^{( k )} d_{I}^{\dagger} d_{J} ,\tag{6.91} $$其中矩阵 $V_{I , J}^{( k )}$ 当跳跃 $t_{k}$ 连接（不连接）格点 I 和 J 时等于 1（0）。在下文中，我们假设 $\mathcal{H}_{0}$ 的基态是非简并的，这意味着所有对应本征值 $\varepsilon_{\alpha} \leq \varepsilon_{F}$ 的轨道均被占据，而其他轨道均为空。那么，无相互作用波函数的变化可以通过扰动理论计算：
$$ | \Phi_{0} ( t_{k} + \delta t_{k} ) \rangle = | \Phi_{0} ( t_{k} ) \rangle - \frac{\delta t_{k}} {\mathcal{H}_{0} - E_{0}} \left( \mathcal{V}^{( k )} - \langle \mathcal{V}^{( k )} \rangle_{0} \right) | \Phi_{0} ( t_{k} ) \rangle + O ( \delta t_{k}^{2} ) ,\tag{6.92} $$其中 $\langle \mathcal{V}^{( k )} \rangle_{0} ~ = ~ \langle \Phi_{0} ( t_{k} ) | \mathcal{V}^{( k )} | \Phi_{0} ( t_{k} ) \rangle$，$E_{0}$ 是无相互作用基态能量。通过将扰动 $\mathcal{V}^{( k )}$ 在对角化哈密顿量 $\mathcal{H}_{0}$ 的基中表示，可以轻松求出导数，参见式 (5.42)：
$$ \mathcal{V}^{( k )} = \sum_{\alpha , \beta} \overline{{V}}_{\alpha , \beta}^{( k )} \phi_{\alpha}^{\dagger} \phi_{\beta} ,\tag{6.93} $$其中
$$ \overline{{V}}_{\alpha , \beta}^{( k )} = \sum_{I , J} U_{I , \alpha}^{*} V_{I , J}^{( k )} U_{J , \beta} ,\tag{6.94} $$U 是用来对角化 $\mathcal{H}_{0}$ 的幺正矩阵，参见式 (5.42)。因此，我们得到：
$$ | \Phi_{0} ( t_{k} + \delta t_{k} ) \rangle = | \Phi_{0} ( t_{k} ) \rangle - \delta t_{k} \sum_{\alpha \neq \beta} \frac{\overline{{V}}_{\alpha , \beta}^{( k )}} {\varepsilon_{\alpha} - \varepsilon_{\beta}} \phi_{\alpha}^{\dagger} \phi_{\beta} | \Phi_{0} ( t_{k} ) \rangle + O ( \delta t_{k}^{2} ) ;\tag{6.95} $$这里，我们使用了态 $\phi_{\alpha}^{\dagger} \phi_{\beta} | \Phi_{0} ( t_{k} ) \rangle$ 是 $\mathcal{H}_{0}$ 的本征态，本征值为 $E_{\alpha , \beta} = E_{0} + \varepsilon_{\alpha} - \varepsilon_{\beta}$。因此，在扰动理论内：
$$ \frac{| \Phi_{0} ( t_{k} + \delta t_{k} ) \rangle - | \Phi_{0} ( t_{k} ) \rangle} {\delta t_{k}} = - \sum_{\alpha \neq \beta} \frac{\overline{{V}}_{\alpha , \beta}^{( k )}} {\varepsilon_{\alpha} - \varepsilon_{\beta}} \phi_{\alpha}^{\dagger} \phi_{\beta} | \Phi_{0} ( t_{k} ) \rangle + O ( \delta t_{k} ) ;\tag{6.96} $$局部算符 ${\mathcal{O}}_{k} ( x )$ 由取 $\delta t_{k} \to 0$ 的极限得到:
$$ \mathcal{O}_{k} ( x ) = - \sum_{\alpha \neq \beta} \frac{\overline{{V}}_{\alpha , \beta}^{( k )}} {\varepsilon_{\alpha} - \varepsilon_{\beta}} \frac{\langle x | \phi_{\alpha}^{\dagger} \phi_{\beta} | \Phi_{0} ( t_{k} ) \rangle} {\langle x | \Phi_{0} ( t_{k} ) \rangle} ;\tag{6.97} $$最后一项可通过将算符 $\phi_{\alpha}^{\dagger}$ (和 $\boldsymbol{\phi}_{\beta} )$ 用 $d_{I}^{\dagger}$ (和 $d_{J} )$ 重新表示来计算:
$$ \mathcal{O}_{k} ( x ) = - \sum_{I , J} W_{I , J}^{( k )} \frac{\langle x | d_{I}^{\dagger} d_{J} | \Phi_{0} ( t_{k} ) \rangle} {\langle x | \Phi_{0} ( t_{k} ) \rangle} ,\tag{6.98} $$其中
$$ W_{I , J}^{( k )} = \sum_{\alpha , \beta}^{\prime} U_{I , \alpha} \ \frac{\overline{{{V}}}_{\alpha , \beta}^{( k )}} {\varepsilon_{\alpha} - \varepsilon_{\beta}} \ U_{J , \beta}^{*} ;\tag{6.99} $$这里，带撇的求和遍历占据轨道 $( \varepsilon_{\beta} \leq \varepsilon_{F} )$ 和未占据轨道 $( \varepsilon_{\alpha} > \varepsilon_{F} )$。
对于 Pfaffian 情形，参见式 (5.106)，可通过前述步骤的直接推广得到。此处，我们还需包含含有两个产生或湮灭算符的微扰:
$$ {\mathcal V}^{( k )} = \sum_{I , J} W_{I , J}^{( k )} d_{I}^{\dagger} d_{J}^{\dagger} + h . c . ,\tag{6.100} $$这对应于配对项的一个微小变化（对于平移不变情形，$\Delta_{j , i}^{\tau , \sigma} \to \Delta_{k}^{\tau , \sigma} )$。然后，通过逆推式 (5.107) 的 Bogoliubov 变换:
$$ d_{I} = \sum_{\alpha} \left( u_{I , \alpha}^{*} \Phi_{\alpha} + \nu_{I , \alpha} \Phi_{\alpha}^{\dagger} \right) ,\tag{6.101} $$微扰将同时包含标准项和反常项，即 $\Phi_{\alpha}^{\dagger} \Phi_{\beta}$ (类似于式 (6.93) 中的项)，以及 $\Phi_{\alpha}^{\dagger} \Phi_{\beta}^{\dagger}$ 和 $\Phi_{\alpha} \Phi_{\beta}$。其余推导过程与行列式情形所示类似。
最后，我们展示如何计算变分波函数对定义反流关联的参数 $\eta_{k}$ 的导数。在我们的表述中，“关联”轨道线性地依赖于这些参数，例如参见式 (5.98) 的最简单情形。对数导数可通过矩阵 M 的通用公式计算:
$$ \frac{\partial} {\partial \eta_{k}} \ln \operatorname*{det} \mathbf{M} = \frac{\partial} {\partial \eta_{k}} \mathrm{Tr} \ln \mathbf{M} = \mathrm{Tr} \left( \mathbf{M}^{- 1} \frac{\partial \mathbf{M}} {\partial \eta_{k}} \right) ,\tag{6.102} $$其中 $\partial \mathbf{M} / \partial \eta_{k}$ 可通过将矩阵 M 的每个元素对参数 $\eta_{k}$ 求导轻松得到。然后，局部算符 ${\mathcal{O}}_{k} ( x )$ 通过将式 (5.99) 中定义具有反流关联的态的矩阵 $\tilde{\mathbf{U}}^{b}$ 作为式 (6.102) 中的 M 来获得。注意逆矩阵不会增加额外计算时间，因为它在 Metropolis 算法中已被计算。