来源论文: https://arxiv.org/abs/2606.07825v1 生成时间: Jun 09, 2026 18:52

投影逆迭代(PII):神经网络量子态基态计算的特征值方法深度解析

0. 执行摘要

自 Carleo 和 Troyer 于 2017 年首次将神经网络量子态(Neural Quantum States, NQS)引入量子多体物理以来,利用深度学习方法求解薛定谔方程(Schrödinger Equation)已成为凝聚态物理和量子化学领域最具前景的路径之一。然而,高维非线性变分流形上的参数优化仍然是一个严峻的技术瓶颈。传统的优化算法(如 Adam、一阶梯度下降)通常表现不佳,而作为行业黄金标准的随机重构(Stochastic Reconfiguration, SR)方法(在机器学习界等价于自然梯度下降,Natural Gradient Descent),其收敛速度极度依赖于体系的光谱间隙(Spectral Gap, $\Delta = E_1 - E_0$)。在面对强挫折磁体、临界点附近的自旋系统以及强关联电子材料时,光谱间隙趋近于零(即小能隙甚至无能隙体系),SR 会遭遇灾难性的“减速”现象,甚至陷入难以收敛的轨道振荡。

为了解决这一核心瓶颈,来自苏黎世联邦理工学院(ETH Zurich)和柏林工业大学(TU Berlin)的研究团队(Hang Zhang, Jannes Nys, Marius Zeinhofer 等)在最新论文中提出了**投影逆迭代(Projected Inverse Iteration, PII)算法。PII 颠覆了将基态求解视为单纯“能量最小化”的传统思维,转而直接将其作为特征值问题(Eigenvalue Problem)进行求解。通过将经典数值线性代数中的逆迭代法(Inverse Iteration)通过 Galerkin 投影 映射到变分流形的切空间,PII 实现了具有能隙不敏感性(Gap-Insensitivity)**的高效二阶收敛。在理论上,PII 消除了收敛率对体系谱展宽(Spectral Spread)的依赖;在计算复杂度上,PII 保持了与 SR 完全相同的多项式级扩展标度。在一系列具有挑战性的二维自旋体系(包括高度受挫的二维 $J_1-J_2$ 模型)的 Benchmark 测试中,PII 展现出了远超 SR 的收敛速度与数值稳定性。这一方法不仅是量子多体计算优化的重大突破,也为深度学习在一般非线性变分特征值问题中的应用树立了新的里程碑。


1. 核心科学问题、理论基础、技术难点与方法细节

1.1 核心科学问题:为什么传统优化方法会在小能隙体系中失效?

在变分蒙特卡洛(Variational Monte Carlo, VMC)框架下,我们将波函数参数化为 $\psi_\theta$,其中 $\theta \in \mathbb{R}^P$ 是神经网络的权重与偏置。寻找基态的经典表述是最小化 Rayleigh 商(Rayleigh Quotient):

$$\min_{\theta} E(\psi_\theta) = \frac{\langle \psi_\theta | \hat{H} | \psi_\theta \rangle}{\langle \psi_\theta | \psi_\theta \rangle} \quad (1.1)$$

在优化实践中,SR 通过引入量子几何张量(Quantum Geometric Tensor, QGT)作为预条件子(Preconditioner),来校正非线性流形上的梯度方向:

$$S(\theta)_{ij} = \langle \partial_{ heta_i} \hat{\psi}_\theta | \partial_{ heta_j} \hat{\psi}_\theta \rangle \quad (1.2)$$$$\theta_{k+1} = \theta_k - \eta S^{-1}( heta_k) \nabla_\theta E(\psi_{ heta_k}) \quad (1.3)$$

然而,从函数空间的几何视角来看,SR 本质上是单位球面 $\mathbb{S}$ 上的黎曼梯度下降(Riemannian Gradient Descent)。根据黎曼梯度下降的收敛性分析(详见论文 Theorem 2),其泛函收敛率 $\rho_{\text{SR}}$ 受以下公式控制:

$$\rho_{\text{SR}} = \left| 1 - \frac{1}{2} \frac{E_1 - E_0}{E_{\max} - E_0} \right| = \left| 1 - \frac{\Delta}{2\Gamma} \right| \quad (1.4)$$

这里,$\Delta = E_1 - E_0$ 是第一激发态与基态之间的光谱间隙,而 $\Gamma = E_{\max} - E_0$ 是算符 $\hat{H}$ 的总光谱展宽。由此可见:

  1. 能隙依赖性:当体系接近量子相变临界点、或者存在竞争有序相时,能隙 $\Delta \to 0$,导致 $\rho_{\text{SR}} \to 1$,收敛极其缓慢。
  2. 光谱展宽依赖性:随着系统尺寸 $N$ 的增大,最大特征值 $E_{\max}$ 呈线性或超线性增长,导致 $\Gamma$ 极其庞大。此时,即使能隙 $\Delta$ 适中,比值 $\Delta / \Gamma$ 也会变得极小,强行压低了收敛速率。这一机制完美解释了为什么在大尺寸 NQS 模拟中,SR 的学习率 $\eta$ 必须微调到极小($\eta \sim 10^{-3}$),且极其容易因不稳定性而崩溃。

1.2 理论基础:非参数化逆迭代(Inverse Iteration)

为了绕过上述瓶颈,数值线性代数中公认的最强基态(最小特征值对应特征态)求解算法是逆迭代法(Inverse Iteration)。给定一个接近基态能量的估算值(位移参数 $\tau \approx E_0$),逆迭代通过算符的逆来极大地放大基态分量:

$$|\hat{\psi}_{k+1}\rangle = \frac{(\hat{H} - \tau \hat{I})^{-1}|\hat{\psi}_k\rangle}{\|(\hat{H} - \tau \hat{I})^{-1}|\hat{\psi}_k\rangle\|} \quad (1.5)$$

为了将其与近代流形优化结合,我们可以将式 (1.5) 改写为球面 $\mathbb{S}$ 上的黎曼更新格式

$$|\hat{\psi}_{k+1}\rangle = \frac{|\hat{\psi}_k\rangle - |d_k angle}{\||\hat{\psi}_k\rangle - |d_k angle\|} \quad (1.6)$$

其中,更新方向 $|d_k angle$ 位于当前状态 $\hat{\psi}_k$ 的切空间 $T_{\hat{\psi}_k}\mathbb{S}$ 内(即满足 $\langle d_k | \hat{\psi}_k \rangle = 0$):

$$|d_k angle = |\hat{\psi}_k angle - \frac{(\hat{H} - \tau \hat{I})^{-1}|\hat{\psi}_k angle}{\langle \hat{\psi}_k | (\hat{H} - \tau \hat{I})^{-1} | \hat{\psi}_k \rangle} \quad (1.7)$$

逆迭代法最吸引人的物理与数学特性在于其收敛率 $\rho_{\text{PII}}$(详见 Theorem 1):

$$\rho_{\text{PII}} = \left| \frac{E_0 - au}{E_1 - au} \right| \quad (1.8)$$

关键结论:$\rho_{\text{PII}}$ 完全独立于最大能量特征值 $E_{\max}$ (即不受总谱展宽 $\Gamma$ 的影响)!并且,当位移 $\tau \to E_0$ 时,$\rho_{\text{PII}} \to 0$,实现超线性甚至二次收敛。即使在能隙 $\Delta \to 0$ 的极限情况下,只要我们能合理设定位移 $\tau$,逆迭代依然能保持极高质量的稳健收敛。

1.3 技术难点:希尔伯特空间的指数墙与物理算符逆的不可解性

尽管逆迭代在理论上具有无与伦比的优越性,但在多体物理的实际应用中存在两个绝望的技术难点:

  1. 指数墙问题:在 $N$ 个自旋-$1/2$ 组成的系统中,希尔伯特空间维度为 $2^N$(例如 $10 \times 10$ 晶格,维度为 $2^{100} \approx 1.2 \times 10^{30}$),任何显式存储或直接计算自旋状态的操作都是绝对不可能实现的。
  2. 算符逆的不可积性:算符 $(\hat{H} - \tau \hat{I})^{-1}$ 是高度非局域的,即使借助蒙特卡洛采样,也无法直接评估一个“逆算符”作用在波函数上的结果。

1.4 方法细节:通过 Galerkin 投影攻克难关

本工作的核心突破在于,作者指出:通过将黎曼更新方向 (1.7) 投影到变分流形的切空间上,可以将无法直接求解的特征值逆算符,转化为在参数向量空间中完全 tractable(可处理)的线性预条件系统

设变分流形在当前参数 $\theta$ 处的局部线性切空间为:

$$\mathcal{V}_\theta = \text{span} \{ \partial_{ heta_i} \hat{\psi}_\theta \}_{i=1}^P \subset T_{\hat{\psi}_\theta} \mathbb{S} \quad (1.9)$$

由于直接计算切空间中的正交投影极度困难,作者巧妙地定义了一个特殊的、依赖于哈密顿量和位移 $\tau$ 的双线性形式(即新型内积):

$$b(v, w) = \langle v | \hat{H} - \tau \hat{I} | w \rangle \quad (1.10)$$

我们寻求一个投影更新矢量 $|d_k^{\text{proj}}\rangle = \sum_j \xi_j |\partial_{ heta_j}\hat{\psi}_\theta\rangle \in \mathcal{V}_\theta$,使得它在双线性形式 $b$ 下最逼近无限维的逆迭代更新方向 $|d_k\rangle$(即进行所谓的 Galerkin 投影):

$$\sum_j b(\partial_{ heta_i}\hat{\psi}_\theta, \partial_{ heta_j}\hat{\psi}_\theta) \xi_j = b(\partial_{ heta_i}\hat{\psi}_\theta, |d_k\rangle) \quad \forall i=1,\dots,P \quad (1.11)$$

现在,让我们见证奇迹的发生。将双线性形式 (1.10) 和逆迭代更新方向 (1.7) 代入上式左右两端:

  • 左端项: $$b(\partial_{ heta_i}\hat{\psi}_\theta, \partial_{ heta_j}\hat{\psi}_\theta) = \langle \partial_{ heta_i}\hat{\psi}_\theta | \hat{H} - \tau \hat{I} | \partial_{ heta_j}\hat{\psi}_\theta \rangle = H(\theta)_{ij} - \tau S(\theta)_{ij} =: Q(\theta)_{ij} \quad (1.12)$$ 这里我们定义了全新的 PII 预条件矩阵 $Q(\theta) = H(\theta) - \tau S(\theta)$。其中 $S(\theta)$ 为标准的重叠矩阵,而 $H(\theta)_{ij} = \langle \partial_{ heta_i}\hat{\psi}_\theta | \hat{H} | \partial_{ heta_j}\hat{\psi}_\theta \rangle$ 为哈密顿算符在切空间基底下的表示。
  • 右端项: $$b(\partial_{ heta_i}\hat{\psi}_\theta, |d_k\rangle) = \langle \partial_{ heta_i}\hat{\psi}_\theta | \hat{H} - \tau \hat{I} | \left[ \hat{\psi}_k - \frac{(\hat{H} - \tau \hat{I})^{-1}\hat{\psi}_k}{\langle \hat{\psi}_k | (\hat{H} - \tau \hat{I})^{-1} \hat{\psi}_k \rangle} \right] \rangle$$ 由于哈密顿算符的厄米性,且 $\hat{\psi}_k$ 的切向导数满足正交性 $\langle \partial_{ heta_i}\hat{\psi}_\theta | \hat{\psi}_k \rangle = 0$,上式中的“不可求逆”项 $(\hat{H} - \tau \hat{I})^{-1}$ 被精确消除(消去)!最终只剩下: $$b(\partial_{ heta_i}\hat{\psi}_\theta, |d_k\rangle) = \langle \partial_{ heta_i}\hat{\psi}_\theta | \hat{H} | \hat{\psi}_k \rangle - 0 - \langle \partial_{ heta_i}\hat{\psi}_\theta | \hat{\psi}_k \rangle \times (\dots) = \frac{1}{2} \partial_{ heta_i} E(\psi_\theta) \quad (1.13)$$

因此,在参数空间中,这一投影过程最终极简地收敛为一个可解的线性代数方程组:

$$Q(\theta) \xi = \frac{1}{2} \nabla_\theta E(\psi_\theta) \quad (1.14)$$

而参数的更新方程为:

$$\theta_{k+1} = \theta_k - \eta_k Q^{-1}( heta_k) \left[ \frac{1}{2} \nabla_\theta E(\psi_{ heta_k}) \right] \quad (1.15)$$

通过这一令人叹为观止的 Galerkin 投影技术,作者成功将一个需要对哈密顿算符求逆的无限维物理问题,转化为了只需在参数向量空间中解线性系统 $Q\xi = g$ 的代数问题


2. 关键 Benchmark 体系、计算所得数据与性能分析

为了全面验证 PII 的优越性,论文在不同尺度、不同物理性质的体系上设计了极其严苛的 Benchmark 测试。

2.1 体系 1:$4 \times 4$ 二维横场伊辛模型(TFIM)——能隙依赖性深度探测

横场伊辛模型(Transverse-Field Ising Model, TFIM)是检验优化算法性质的经典试金石:

$$\hat{H}_{\text{TFIM}} = -J \sum_{\langle i,j \rangle} \hat{\sigma}_i^z \hat{\sigma}_j^z - h \sum_i \hat{\sigma}_i^x \quad (2.1)$$

作者选取了两种截然不同的物理耦合区间进行对比,对应不同大小的能隙 $\Delta = E_1 - E_0$:

  1. 小能隙(强关联临界区):$h/J = 2$,此时精确对角化(ED)算出的能隙仅为 $\Delta = 0.015$。由于系统处于高度纠缠与近简并态,优化极其困难。
  2. 大能隙(顺磁相区):$h/J = 4$,对应的能隙 $\Delta = 3.25$,优化相对简单。

变分波函数采用非线性 Vision Transformer(ViT)架构(4层,隐藏层维度 $d=20$,2个注意力头)。

核心计算数据与收敛对比(图2):

  • 在小能隙体系中($h/J=2, \Delta = 0.015$)
    • PII(参数设置:$\eta = 1.0, \lambda = 0.01, \tau = E_0$):仅需 10次迭代,能量的相对误差 $\epsilon = |(E_{\text{ViT}} - E_{\text{ED}}) / E_{\text{ED}}|$ 就急剧下降至 $10^{-4}$ 以下,在 30 次迭代内彻底收敛至极高精度的基态($\epsilon \sim 10^{-6}$)。
    • SR(参数设置:$\eta = 0.01, \lambda = 10^{-4}$,已做最精细超参微调):收敛极其缓慢,呈现典型的渐进性减速,在 80 次迭代后仍有明显的相对误差,且在小能隙下出现严重的阻尼震荡。
    • Adam(一阶梯度基线,$\eta = 0.01$):几乎完全失效,停留在相对误差 $\epsilon \sim 10^{-2}$ 处无法前行。
  • 在大能隙体系中($h/J=4, \Delta = 3.25$)
    • PII 与小能隙表现完全一致,表现出完美的能隙不敏感性(Gap-Insensitivity)
    • SR 此时由于 $\Delta$ 变大,收敛有所加快,但其速度依然显著慢于 PII。

这组数据强有力地支持了论文的理论预测:PII 能够彻底抹平不同物理能隙之间的收敛速度差异。

2.2 体系 2:位移参数 $\tau$ 的“欠冲(Undershooting)”效应与采样噪声抑制(图3)

理论上,当 $\tau = E_0$(精确基态能量)时,PII 收敛最快。然而,在 VMC 的实际运行中,由于蒙特卡洛采样的存在,哈密顿矩阵元 $H(\theta)_{ij}$ 包含统计噪声。当 $\tau \to E_0$ 且采样数 $M$ 较小时,预条件矩阵 $Q(\theta) = H(\theta) - \tau S(\theta)$ 会由于噪声扰动而失去正定性,导致优化崩塌。为了提升数值鲁棒性,作者引入了“位移欠冲(Undershooting)”技术,即选择 $\tau < E_0$(注意:在负能量系统如自旋模型中,这意味着令 $\tau$ 大于 $E_0$ 的绝对值,例如 $\tau = 1.2 E_0$)。

在 $4\times 4$ TFIM ($h/J=2$) 体系上进行的不同采样数 $M$ 的消融实验显示:

  • 无欠冲($\tau = E_0$):当采用大采样量 $M=16384$ 时,PII 收敛完美;但当采样量缩减至 $M=4096$ 或 $M=512$ 时,由于采样噪声,优化轨道在后期发生剧烈崩塌,无法稳定收敛。
  • 有欠冲($\tau = 1.2 E_0$):即使在极小样本量 $M=512$ 下,优化也展现出极强的稳定性,相对误差极其平滑地降至 $\sim 10^{-5}$。这证明了增加位移偏置能有效平移 $Q$ 矩阵的谱分布,使其在随机噪声下保持良态(Positive Definite)。
  • 无噪声确定性极限(Full-Sum):在不采用采样、直接进行全希尔伯特空间求和的消融实验中,调节 $\tau = \alpha E_0$($\alpha = 1, 1.2, 2, 4$)的收敛曲线严格遵循理论公式 (1.8) 的预测。这完美厘清了物理机制(确定性收敛率)与随机涨落(统计稳定性)之间的权衡关系。

2.3 体系 3:$10 \times 10$ 大尺寸二维自旋体系的卓越表现(图4)

为了证明 PII 在科学计算前沿(不可精确求解的指数级维度)的实用价值,作者在 $10\times 10$ 晶格(希尔伯特空间维度 $2^{100}$)上挑战了三个凝聚态物理学标志性的 Hamiltonian:

  1. 横场伊辛模型(TFIM, $h/J=2$):stoquastic,使用实自旋 ViT($P \approx 12,000$ 变分参数)。
  2. 海森堡反铁磁模型(Heisenberg AFM, $J_2=0$):stoquastic,使用实自旋 ViT。
  3. 受挫反铁磁 $J_1-J_2$ 模型($J_2/J_1 = 0.5$):存在极其严重的符号问题(Sign Problem),地面态具有高度复杂的非平凡自旋相位结构。波函数参数化必须使用复数形式的 ViT 架构(4层,隐藏维数 $d=60$,10个头,参数量巨大)。

关键性能数据(图4):

  • TFIM ($10 \times 10$):PII 仅需 20次迭代 能量便彻底饱和并收敛到极佳基态;而 SR 需要近 200次迭代 才能勉强接近相同能量水平。
  • 海森堡模型 ($10 \times 10$):PII 在 30次迭代 内完全收敛;SR 则面临长达 150次迭代 的漫长平台期,且收敛后的最终基态能量依然略逊于 PII。
  • $J_1 - J_2$ 模型 ($10 \times 10$, 最具挑战性的强关联受挫磁体):PII(参数 $\tau = E_0, \eta = 1.0, \lambda = 10^{-1}$)在 40次迭代 内将单格点能量推进至 $\sim -0.495 J_1$ 的极低水平,完全展现出强有力的二阶跳跃收敛特征;而 SR 方法($\eta = 0.004$)在前 100 次迭代中收敛速率极慢,极易陷入局部亚稳态。

这一系列令人振奋的前沿数据表明,PII 在计算复杂度不增加的前提下,将高维强关联体系的收敛效率提升了整整一个数量级!


3. 代码实现细节、复现指南与开源生态

3.1 预条件矩阵 $Q(\theta)$ 的高效 Monte Carlo 估计数学物理表述

要在 VMC 中实用化 PII,最核心的一步是如何通过对 Born 概率分布 $p_\theta(x) = |\psi_\theta(x)|^2 / \sum_{x'} |\psi_\theta(x')|^2$ 进行采样,来高效估计复杂的 $H(\theta)_{ij}$ 和 $S(\theta)_{ij}$ 矩阵。以下为核心算符的 Monte Carlo 估计子公式:

定义对数导数(Log-Derivative)和局部能量(Local Energy):

$$J_i(x) = \partial_{ heta_i} \log \psi_\theta(x) \quad (3.1)$$$$E_L(x) = \frac{[\hat{H}\psi_\theta](x)}{\psi_\theta(x)} \quad (3.2)$$

对于 $M$ 个独立采样点 $\{x^{(m)}\}_{m=1}^M$,我们定义两个极其关键的辅助矩阵 $O \in \mathbb{R}^{2M \times P}$ 和 $A \in \mathbb{R}^{2M \times P}$(利用实部与虚部分割技术将复数系统实数化,以统一处理复数波函数):

$$O_{nj} = \frac{1}{\sqrt{M}} \text{ReIm}_n \left( J_j(x^{(m)}) - \mathbb{E}_p[J_j] \right) \quad (3.3)$$$$A_{nj} = \frac{1}{\sqrt{M}} \text{ReIm}_n \left( E_{L,\theta_j}(x^{(m)}) + E_L(x^{(m)})(J_j(x^{(m)}) - \mathbb{E}_p[J_j]) \right) \quad (3.4)$$

这里 $\text{ReIm}_n(\cdot)$ 的物理含义是:对于样本序号 $n = 1, \dots, M$ 取其实部,对于 $n = M+1, \dots, 2M$ 取其虚部。从而将复线性方程变换为 $2M$ 维度的实矩阵运算。

借助这些高度优化的辅助矩阵,能量梯度 $\nabla_\theta E$ 以及 PII 算符估计子写为极其简洁的矩阵外积形式:

$$\nabla_\theta E = 2 (O^T e) \quad (3.5)$$$$S(\theta) = O^T O \quad (3.6)$$$$H(\theta) = O^T A \quad (3.7)$$

因此,PII 的预条件子 $Q(\theta)$ 表示为:

$$Q(\theta) = O^T A - \tau O^T O \quad (3.8)$$

3.2 极小空间隐式求解技术:MinPII (Push-Through Identity)

对于具有大型神经网络架构的 NQS,其参数量 $P$(通常在 $10^4 \sim 10^5$ 量级)通常远大于单步蒙特卡洛采样数 $2M$(通常在 $10^3 \sim 10^4$ 之间)。这意味着矩阵 $Q(\theta) \in \mathbb{R}^{P \times P}$ 极其庞大,显式构建并求逆的计算成本高达 $\mathcal{O}(P^3)$,是完全不可接受的。

作者在此处无缝移植了近代优化算法中的推通恒等式(Push-Through Identity),将其命名为 MinPII

$$(O^T A - \tau O^T O + \lambda I_P)^{-1} O^T e = O^T \left( A O^T - \tau O O^T + \lambda I_{2M} \right)^{-1} e \quad (3.9)$$

重大物理意义:MinPII 将需要求逆的矩阵维度从变分参数空间 $P \times P$ 暴降到了样本空间 $2M \times 2M$。通过先在低维样本空间中解线性系统 $\left( A O^T - \tau O O^T + \lambda I_{2M} \right) y = e$,再将其投影回切空间 $\xi = O^T y$,计算开销被直接压缩至 $\mathcal{O}(\min(M^2P, P^2M))$。这与目前业界最高效的 MinSR 方法完全属于同一复杂度标度,保证了在大系统、深层网络下的高能物理扩展性。

3.3 核心复现伪代码与算法设计

以下是基于论文核心逻辑,在 Python 框架中进行 PII 单步参数更新的典型代码架构(基于 JAX 及自动微分实现):

import jax
import jax.numpy as jnp

def pii_step(params, samples, model_apply, hamiltonian_local_energy, tau, lmbda, eta):
    """
    执行单步投影逆迭代(PII)参数更新
    params: 神经网络变分参数 (P,)
    samples: 蒙特卡洛采样的构型 (M, N_sites)
    model_apply: 神经网络计算对数波函数值 log_psi(params, x)
    hamiltonian_local_energy: 计算每个样本的局部能量 E_L(x)
    tau: 能量位移参数 (float)
    lmbda: Tikhonov 正则化常数 (float)
    eta: 学习率 (float)
    """
    M = samples.shape[0]
    
    # 1. 计算局部能量和对数波函数对参数的雅可比矩阵 (Jacobian)
    # 使用 jax.vmap 和 jax.jacobian 实现高度并行计算
    log_psi_fn = lambda p, x: model_apply(p, x)
    jac_fn = jax.vmap(jax.jacobian(log_psi_fn), in_axes=(None, 0))
    
    # J_raw 形状: (M, P) - 复数矩阵
    J_raw = jac_fn(params, samples)
    # E_L_val 形状: (M,) - 复数矩阵
    E_L_val = hamiltonian_local_energy(params, samples)
    
    # 2. 估计 E_L 对参数的导数 (E_L,j)
    # 论文指出:E_L,j = d(E_L)/d_theta_j,可通过 JAX 的双重微分或 VJP 高效获得
    local_energy_fn = lambda p, x: hamiltonian_local_energy(p, x)
    jac_EL_fn = jax.vmap(jax.jacobian(local_energy_fn), in_axes=(None, 0))
    J_EL = jac_EL_fn(params, samples) # (M, P)
    
    # 计算均值以进行中心化
    E_mean = jnp.mean(E_L_val)
    J_mean = jnp.mean(J_raw, axis=0)
    
    # 3. 构建实数化的 O 矩阵和 A 矩阵
    J_centered = J_raw - J_mean
    E_diff = E_L_val - E_mean
    
    # 分离实部与虚部,拼装成 2M 行的实数辅助矩阵
    O_real = jnp.real(J_centered) / jnp.sqrt(M)
    O_imag = jnp.imag(J_centered) / jnp.sqrt(M)
    O = jnp.concatenate([O_real, O_imag], axis=0) # (2M, P)
    
    E_diff_real = jnp.real(E_diff) * 2.0 / jnp.sqrt(M)
    E_diff_imag = jnp.imag(E_diff) * 2.0 / jnp.sqrt(M)
    e = jnp.concatenate([E_diff_real, E_diff_imag], axis=0) # (2M,)
    
    # A 矩阵项: E_L,j + E_L * (J_j - E[J_j])
    A_term = J_EL + E_L_val[:, None] * J_centered
    A_real = jnp.real(A_term) / jnp.sqrt(M)
    A_imag = jnp.imag(A_term) / jnp.sqrt(M)
    A = jnp.concatenate([A_real, A_imag], axis=0) # (2M, P)
    
    # 4. 采用 Push-Through Identity (MinPII) 解低维线性系统 (假设 2M < P)
    # 构建核算子:K = A_O^T - tau * O_O^T + lambda * I_(2M)
    # 注意:A @ O.T 的尺寸为 (2M, 2M)
    K = (A @ O.T) - tau * (O @ O.T) + lmbda * jnp.eye(2 * M)
    
    # 解线性方程组 K * y = e
    y = jnp.linalg.solve(K, e)
    
    # 将解投影回切空间,获得参数更新方向 xi
    xi = O.T @ y # (P,)
    
    # 5. 执行参数更新
    new_params = params - eta * 0.5 * xi
    return new_params

3.4 关键开源工具与复现仓库链接

  • 官方代码实现开源仓库projected_inverse_iteration
  • 核心底层计算框架:基于 NetKet 3(量子多体深度学习物理社区的顶梁柱,集成丰富的自旋与费米子哈密顿量算符以及 MCMC 采样技术)与 JAX(提供极为强悍的自动微分与 GPU 高性能并行加速)。

4. 关键引用文献与局限性批判性讨论

4.1 本工作立足的重要学术基石

  1. 神经网络量子态的开山之作:Carleo, G. & Troyer, M. Solving the quantum many-body problem with artificial neural networks. Science 355, 602–606 (2017). (定义了 VMC + NQS 的基本框架)。
  2. 随机重构(SR)方法的经典引入:Sorella, S. Green function monte carlo with stochastic reconfiguration. Phys. Rev. Lett. 80, 4558 (1998). (提出了量子物理中的自然梯度法)。
  3. 无限维几何视角的二阶优化框架:Müller, J., Zeinhofer, M. et al. Position: Optimization in SciML should employ the function space geometry. ICML 2024. (为 PII 的 Galerkin 投影提供了核心数学思想)。
  4. Kaczmarz 加速算法的基础:Goldshlager, G. et al. A kaczmarz-inspired approach to accelerate the optimization of neural network wavefunctions. J. Comput. Phys. 516, 113351 (2024). (启发了 PII-SPRING 算法的演进)。

4.2 局限性与前沿挑战性评论

尽管 PII 表现出无与伦比的收敛优势,但作为一项先驱性工作,在实际推广和工业/科研落地上面临以下极具启发性的技术挑战:

  1. 位移参数 $\tau$ 的超参选择困境

    • 物理硬伤:由公式 (1.8) 可知,PII 的收敛率取决于 $\tau$。若不小心将 $\tau$ 设得高于真实基态能 $E_0$(且越过了第一激发能),PII 会受到数值吸引力的拉扯,直接收敛到第一激发态甚至更高阶的激发态,而无法求出真正的基态。这使得用户需要对基态能量有一个相对准确的下界(Lower Bound)预估。
    • 自适应的尴尬:虽然可以采用自适应位移 $\tau_k = E(\psi_{\theta_k})$(演变为瑞利商迭代,Rayleigh Quotient Iteration),但在 VMC 优化前期,变分能量起伏极大,这种自适应很容易在前期直接将系统锁死在亚稳态,数值极不稳定。
  2. 估计子包含一阶导数,对采样噪声更敏感

    • 在构建 $A_{nj}$ 矩阵时,我们需要显式求取局部能量对参数的导数 $E_{L,\theta_j}(x)$。这引入了哈密顿算符 $\hat{H}$ 的二次相互作用项。在复杂的强关联自旋晶格或量子化学分子中,这个导数的方差通常比普通的对数导数 $J_i(x)$ 的方差大得多。这导致 PII 的采样开销可能在某些极度无序的系统中显著上升,需要更精妙的统计消噪和正则化处理。
  3. 非厄米物理算符与复数流形的正则化阵痛

    • 虽然理论上哈密顿量 $\hat{H}$ 是厄米的,但在实际 Monte Carlo 抽样估计下,$H(\theta) = O^T A$ 仅具有渐进厄米性(Asymptotic Hermiticity)。在有限样本数下,非厄米扰动可能导致 $Q$ 矩阵产生复杂的虚部特征值,导致参数更新轨道失稳,因此对正定正则化项(如 Tikhonov $Q + \lambda I$)的微调仍然具有一定的技术门槛。
  4. 费米子体系符号问题的阴霾

    • 在量子化学和强关联电子材料中,费米子波函数的反对称性会引发灾难性的自旋/电荷符号问题。PII 只是一个代数层面的“解方程”优化器,无法在物理层面上消灭符号问题所带来的 Monte Carlo 采样权重的指数级坍塌。这需要后续研究将 PII 与神经网络(如 FermiNet 或 Pfap神经网络)进行更底层的架构融合。

5. 补充拓展:深度解读 PII 优化方法的数学美学与衍生算法

为了帮助研究人员更透彻地理解 PII,我们将从流形微分几何的视角对现有主流 NQS 优化算法进行一次大一统的高屋建瓴式梳理,并详细推导 PII 极其重要的物理加速变体:PII-SPRING

5.1 黎曼优化视角的“大一统”图景

在现代科学机器学习(SciML)中,最具有美学冲击力的发现之一是:所有的经典优化方法,都只是在特定空间、选择特定内积、对特定泛函优化方向进行 Galerkin 投影的产物。我们可以用下表进行完美的总结(对应论文 Table 1 的深化表述):

参数空间更新算法函数空间(希尔伯特空间)优化器模型采用的双线性投影内积 $b(v, w)$对应在数值线性代数中的特征值求解算子
Stochastic Reconfiguration (SR)球面 $\mathbb{S}$ 上的黎曼 $L^2$ 梯度下降$b(v, w) = \langle v \vert w \rangle$ ($L^2$ 内积)幂迭代法(Power Iteration)(参见附录 F,SR 可自然重写为哈密顿量的移位幂迭代过程)
Projected Inverse Iteration (PII) (本工作)球面 $\mathbb{S}$ 上的全局化黎曼牛顿法$b(v, w) = \langle v \vert \hat{H} - \tau \hat{I} \vert w \rangle$移位逆迭代(Shifted Inverse Iteration)
Rayleigh Gauss-Newton (RGN)球面 $\mathbb{S}$ 上的经典黎曼牛顿法$b(v, w) = \langle v \vert \hat{H} - E(\psi) \hat{I} \vert w \rangle$ (黎曼海森矩阵)瑞利商迭代(Rayleigh Quotient Iteration)
Wasserstein Quantum Monte Carlo (WQMC)概率分布空间上的 Wasserstein 梯度流$\dot{\rho} = - \nabla_W E$ (不受球面流形约束)

这一宏大的数学图景揭示了为什么 SR 会在能隙变小时表现糟糕:因为它在投影时使用了完全不包含哈密顿谱信息的简单 $L^2$ 内积,从而对哈密顿量的能量地貌“一无所知”;而 PII 巧妙地将投影内积本身与哈密顿能量算符 $\hat{H} - \tau\hat{I}$ 绑定,从而使得每一次参数迭代都预载了最本质的量子能域物理结构信息。

5.2 极致样本效率提升:Kaczmarz 加速版 PII-SPRING 算法

在蒙特卡洛采样预算严重受限(例如分子轨道庞大,只能使用极小 batch size $M$ 采样)的场景下,PII 面对的线性方程组 (1.14) 通常处于高度欠定(Underdetermined)状态,此时求解极易因不确定性噪声而发散。受优化界 SPRING 算法的启发,作者推导了 PII-SPRING 算法,其代数核心是从最小二乘惩罚项锚定的视角重新审视优化轨迹。

我们将 Galerkin 投影在最小二乘意义下写为损失函数最小化问题:

$$\xi^* = \arg\min_{\xi \in \mathbb{R}^P} \mathcal{L}(\xi) \quad \text{with} \quad \mathcal{L}(\xi) = \frac{1}{2} \left\| \sum_j |\partial_j \hat{\psi}_k\rangle \xi_j - |d_k\rangle \right\|_b^2 \quad (5.1)$$

将其展开,即得到标准二次型:

$$\mathcal{L}(\xi) = \frac{1}{2} \xi^T Q(\theta_k) \xi - \frac{1}{2} \xi^T \nabla_\theta E(\psi_k) + \text{const} \quad (5.2)$$

为了通过引入历史更新步长 $\xi_{k-1}$ 的相关性来平滑随机采样涨落(即利用 Kaczmarz 代数投影法),我们不再惩罚参数更新矢量的绝对范数 $\|\xi\|^2$,而是改用锚定惩罚项

$$\xi^* = \arg\min_{\xi \in \mathbb{R}^P} \left[ \mathcal{L}(\xi) + \frac{\lambda}{2} \|\xi - \mu \xi_{k-1}\|_2^2 \right] \quad (5.3)$$

这里,$\mu \in [0, 1]$ 是动量衰减常数,其充当了历史梯度的物理惯性。通过极值条件 $\nabla_\xi [\dots] = 0$,我们能够立刻导出具有极强抗噪能力的 PII-SPRING 更新方程

$$(Q(\theta_k) + \lambda I_P)(\xi - \mu \xi_{k-1}) = \frac{1}{2} \nabla_\theta E(\psi_k) - \mu Q(\theta_k) \xi_{k-1} \quad (5.4)$$

同样地,我们可以利用推通恒等式(MinPII-SPRING)进行高速隐式求逆:

$$\xi^* = \mu \xi_{k-1} + O^T \left( A O^T - \tau O O^T + \lambda I_{2M} \right)^{-1} O^T \left[ \frac{1}{2} e - \mu (A - \tau O) \xi_{k-1} \right] \quad (5.5)$$

物理性能增益分析(图 G3):

在 $10\times 10$ TFIM ($h/J=2$) 体系、极度嘈杂的小样本($M=4000$)严苛测试中:

  • 标准 MinPII:在优化后期由于极小样本引入的重构偏置,相对误差横盘在 $10^{-4}$ 处,并夹杂大量的刺骨噪声(随机尖峰)。
  • MinPII-SPRING(参数 $\mu = 0.8$):通过复用上一代参数空间的特征几何信息,以极度平滑的轨迹直冲 $\epsilon \sim 10^{-6}$ 的物理收敛极限,展现了难以置信的样本节省能力。

这一拓展加速算法彻底打消了物理学界关于“二阶特征值优化算法由于导数项较多必然极度耗费样本”的固有偏见,为超大变分参数 NQS 迈向复杂真实分子及三维关联电子晶格铺平了坚实的道路。