来源论文: https://arxiv.org/abs/2606.20522v1 生成时间: Jun 19, 2026 06:59
0. 执行摘要
在一维量子多体物理学中,代数衰减(或幂律衰减)的长程相互作用广泛存在于诸如等效自旋链、RKKY 耦合以及里德堡原子和极性分子等量子模拟器体系中。然而,传统的无限矩阵乘积态(iMPS)和矩阵乘积算符(MPO)算法在处理长程相互作用时面临本质困难:MPO 结构天然适用于表达有限程相互作用或指数衰减相互作用。传统的权宜之计是采用**指数之和(Sum of Exponentials, SOE)**来拟合代数尾部,将其转化为有限极点的有理近似。尽管这种方法在一定空间窗口内具有高精度,但其本质上依然是在红外区(长程极限下)引入了人为的指数截断(有限极点代用品),从而混淆了系统真实的物理临界性质,并在变分优化中引入了不可控的 Hamiltonian 律定偏差。
由阿姆斯特丹大学的 Qi Yang 提出的最新研究,展示了一种创新的无截断变分方法。该方法完全不修改目标长程 Hamiltonian,而是通过在转移矩阵通道中直接对代数相互作用的无限序列进行求和。这一操作将长程两点关联函数的变分计算转化为作用在连通转移矩阵 $\widetilde{T}_A$ 上的非有理矩阵函数 $F_{\alpha,Q}(\widetilde{T}_A)$ 的求值。通过结合 Krylov 子空间方法(Arnoldi 过程)进行低成本、高精度的矩阵函数作用(Matrix-Function Action)评估,并利用 Fréchet 共轭技术与隐式定点微分(基于 Daleckii-Krein 公式的分裂差分)实现稳定、精确的逆向模式自动微分,该研究彻底打通了变分 iMPS 直接求解长程幂律体系的计算瓶颈。这一方法在长程自由费米子、Haldane-Shastry 链以及长程横场 Ising 模型等标杆体系上得到了完美验证,展示出其在捕捉真实物理临界行为方面超越传统 SOE-MPO 的绝对优势。
1. 核心科学问题、理论基础、技术难点与方法细节
1.1 核心科学问题与传统 SOE 方案的缺陷
考虑一维代数衰减相互作用的两点哈密顿量项:
$$ J(r) = \frac{\cos(Qr)}{r^\alpha} $$其中,$\alpha$ 为衰减指数,$Q$ 控制振荡频率($Q=0$ 对应非振荡长程作用,$Q\neq 0$ 对应反铁磁或不共度调制长程作用)。在无限大系统(iMPS)下,由于系统平移对称性,能量密度的变分评估需要计算如下无限求和:
$$ E_{\text{long}} = \sum_{r \ge 1} J(r) \langle X_i Y_{i+r} \rangle $$传统的做法是将幂律衰减 $r^{-\alpha}$ 拟合为有限个($K$ 个)指数因子的叠加(即 SOE-MPO 方案):
$$ \frac{1}{r^\alpha} \approx \sum_{k=1}^K c_k \lambda_k^r $$这一方法成功地将无限长程求和转化为了一个具有辅助键维数 $\chi_{\text{MPO}} \propto K$ 的标准 MPO。然而,在物理临界点附近,这一拟合带来了严峻的技术缺陷:
- 红外截断与人为交叉(Crossover):任何有限极点 $K$ 的 SOE 在渐近极限 $r \to \infty$ 下都不可避免地退化为最慢的指数衰减 $\lambda_{\max}^r$。这意味着 SOE 近似本质上是在长程极限下将幂律哈密顿量替换为了短程哈密顿量。这会人为引入一个特征关联长度 $\xi_{\text{SOE}} \propto -1/\ln \lambda_{\max}$,在接近物理临界区时,它会与系统本身的物理关联长度混淆,从而污染临界指数的提取。
- Hamiltonian 律定偏差:在进行 iMPS 的变分键维数 $D$ 的外推时,SOE-MPO 优化的是一个“被修改”的哈密顿量,因而增加 $D$ 只会改善对该“有理代用品”的描述,而无法消除代用品与原真实代数哈密顿量之间的本征偏差。为了得到干净的临界行为,必须进行 $K \to \infty$ 和 $D \to \infty$ 的双重外推,计算成本极高且极不稳定。
1.2 理论基础:转移矩阵函数(Transfer-Matrix Functions)
本工作提出了一种完全不同的范式:保持目标代数哈密顿量精确不变,利用 iMPS 的连通两点关联函数天生由转移矩阵的幂次生成的性质,直接在转移通道中完成对代数尾部的解析求和。
对于平移对称的单位置元 iMPS(其三阶张量为 $A^s$),定义标准转移矩阵 $\mathcal{T}_A$ 为一个作用在 $D^2$ 维虚拟双线空间(Virtual Double Space)上的线性映射:
$$ \mathcal{T}_A(W) = \sum_{s} A^s W (A^s)^\dagger, \quad \mathcal{T}_A \equiv T^I_A $$设其左右最大特征向量为 $\langle L|$ 和 $|R\rangle$,满足 $\langle L|\mathcal{T}_A = \langle L|$,$\mathcal{T}_A|R\rangle = |R\rangle$,且归一化 $\langle L|R\rangle = 1$。为了分离关联函数中的常数贡献,定义投影算符 $\mathcal{P} = \mathbb{I} - |R\rangle\langle L|$,由此可构建连通转移矩阵(Connected Transfer Matrix):
$$ \widetilde{T}_A = \mathcal{P} \mathcal{T}_A \mathcal{P} $$定义算符 $O = X, Y$ 的一位置信期望值 $\langle O \rangle = \langle L|\mathcal{T}^O_A|R\rangle$,以及中心化算符 $O - \langle O \rangle \mathbb{I}$ 的连通读出(Readout)左向量 $\langle L_X|$ 与源(Source)右向量 $|R_Y\rangle$:
$$ \langle L_X| = \langle L|\mathcal{T}_A^X \mathcal{P}, \quad |R_Y\rangle = \mathcal{P}\mathcal{T}_A^Y |R\rangle $$根据转移矩阵的代数性质,距离为 $r \ge 1$ 的两点关联函数可以精确表示为:
$$ \langle X_i Y_{i+r} \rangle = \langle X \rangle \langle Y \rangle + \langle L_X| \widetilde{T}_A^{r-1} |R_Y\rangle $$将此代入长程求和式,并定义特征标量核函数 $F_{\alpha,Q}(z)$:
$$ F_{\alpha,Q}(z) = \sum_{n=0}^\infty \frac{e^{iQ(n+1)}}{(n+1)^\alpha} z^n = \frac{\text{Li}_\alpha(e^{iQ} z)}{z} $$其中 $\text{Li}_\alpha(x)$ 是多重对数函数(Polylogarithm)。由线性叠加原理,真实的幂律相互作用能量密度可精确写为:
$$ E_{\alpha,Q}[X, Y] = J \text{Re} \left[ \text{Li}_\alpha(e^{iQ}) \langle X \rangle \langle Y \rangle + \langle L_X| F_{\alpha,Q}(\widetilde{T}_A) |R_Y\rangle \right] $$这一等式的物理意义非凡:长程相互作用尾部的无限求和被完美浓缩为非有理矩阵函数 $F_{\alpha,Q}(\widetilde{T}_A)$ 在源向量 $|R_Y\rangle$ 上的单次矩阵动作! 没有引入任何截断,没有拟合极点,在固定的变分键维数 $D$ 下,该表达式是恒等式。
1.3 技术难点与破解方案
难点一:高维非正规矩阵的函数计算
由于转移矩阵 $\widetilde{T}_A$ 具有 $D^2 \times D^2$ 的巨大维度(当 $D=256$ 时,维度高达 $65536$),对其进行密集的完全对角化计算需要 $O(D^6)$ 的时间和 $O(D^4)$ 的空间复杂度,这在计算上是不可承受的。更为关键的是,由于 iMPS 转移矩阵一般是非正规矩阵(Non-normal Matrix),其特征向量之间可能高度接近线性相关,这导致直接光谱对角化在数值上极度不稳定,容易导致伪谱(Pseudospectra)效应和数值发散。
破解方案:Krylov 子空间投影(Arnoldi 算法) 我们不显式求解完全的 $\widetilde{T}_A$,而是通过 Arnoldi 算法,从源向量 $|R_Y\rangle$ 出发,构建一个低维(维度 $m \ll D^2$)的克罗夫正交基矩阵 $V_m$,将高维转移矩阵投影压缩为一个 $m \times m$ 维的上 Hessenberg 矩阵 $H_m$:
$$ H_m = V_m^\dagger \widetilde{T}_A V_m $$通过对低维矩阵 $H_m$ 进行双正交对角化,求得 Ritz 值 $\theta_i$ 以及相应的左右特征向量 $u_i$ 和 $\tilde{u}_i$。转移矩阵函数的作用即可通过在克罗夫子空间内进行光谱读出来高精度逼近:
$$ \langle L_X| F_{\alpha,Q}(\widetilde{T}_A) |R_Y\rangle \approx \sum_{i=1}^m F_{\alpha,Q}(\theta_i) \langle L_X| V_m | u_i \rangle \langle \tilde{u}_i | V_m^\dagger |R_Y\rangle $$单次 Arnoldi 迭代只需进行矩阵-向量乘法 $\widetilde{T}_A v$,其时间复杂度仅为 $O(D^3)$,总成本压缩至 $O(m D^3)$,与常规 MPO 算法处于同一数量级。
难点二:通过 Krylov 投影和非有理函数的变分梯度回传
变分能量最小化依赖于计算关于 iMPS 张量 $A$ 的精确梯度。由于计算链条中包含了 Krylov 子空间投影、非有理多重对数函数求值、以及可能的非对称对角化,传统的通过计算图自动微分(AD)会面临极大的不稳定性,特别是在 Arnoldi 递归中回传梯度会导致内存暴增以及严重的累积误差。
破解方案:Fréchet 共轭与隐式定点微分 为了实现数学上严谨、数值上高度稳定的梯度回传,本研究摒弃了穿透克罗夫迭代图的传统微分法,而是将矩阵动作 $y = F(M)b$ 视作一个独立的微分原语,通过求解其 Fréchet 导数算符 的共轭作用(Adjoint Action)来实现反向传播。
根据 Daleckii-Krein 定理,矩阵函数 $F(M)$ 的一阶变分可以由以下分块矩阵的函数求值给出:
$$ F\begin{pmatrix} M & \delta M \\ 0 & M \end{pmatrix} = \begin{pmatrix} F(M) & \mathcal{L}_F(M; \delta M) \\ 0 & F(M) \end{pmatrix} $$其中 $\mathcal{L}_F(M; \delta M)$ 即为 Fréchet 导数。对于输入的伴随算符(梯度共轭)$\bar{y}$,其向源向量 $b$ 和矩阵 $M$ 的拉回(Pullback)公式分别为:
$$ \bar{b} = F(M)^\dagger \bar{y} $$$$ \langle \bar{M}, \delta M \rangle = \langle \bar{y}, \mathcal{L}_F(M; \delta M) b \rangle $$在数值实现中,我们构建独立的左、右克罗夫子空间,将高维 Fréchet 映射投影为小维度的标量分裂差分矩阵(Daleckii-Krein Table):
$$ F^{[1]}_{ij} = \begin{cases} \frac{F(\theta^L_i) - F(\theta^R_j)}{\theta^L_i - \theta^R_j}, & \theta^L_i \neq \theta^R_j \\ F'(\theta^L_i), & \theta^L_i = \theta^R_j \end{cases} $$利用该分裂差分矩阵在 Ritz 基下进行逐元素矩阵乘法,可以极其稳定、高精度地重构出关于转移矩阵和张量的全部变分梯度。这彻底避免了 Arnoldi 展开过程中的病态不稳定性。
2. 关键 Benchmark 体系、计算数据与性能展示
为了全面验证该转移矩阵函数(TMF)算法的优越性,研究者在三个极具代表性的长程多体物理体系中进行了严苛的数值标定。
2.1 长程自由费米子(带有 Jordan-Wigner 弦)
长程无自旋费米子的哈密顿量为:
$$ H_c = \sum_i \sum_{r \ge 1} \frac{\cos(Qr)}{r^\alpha} (c^\dagger_i c_{i+r} + c^\dagger_{i+r} c_i) $$在经过 Jordan-Wigner 变换映射到自旋链表示时,费米子算符的非局部性会引入一个复杂的宇称算符弦(Parity String):
$$ c^\dagger_i c_{i+r} = \sigma^-_i \left( \prod_{m=i+1}^{i+r-1} \sigma^z_m \right) \sigma^+_{i+r} $$这对应于在转移通道中,普通的转移矩阵被替换为带有宇称算符的弦转移矩阵(String Transfer Matrix) $\mathcal{T}_A^\Pi$。该模型在热力学极限下具有精确的解析色散关系 $\varepsilon(k)$,因此是验证算法精度的绝佳平台。
标定数据(图 3 展示):
- 能量收敛性:对于调制频率 $Q=\sqrt{3}$ 和 $Q=\sqrt{5}$,在不同的幂律指数 $\alpha = 1, 1.5, 2$ 下,变分能量误差 $E(D) - E_{\text{exact}}$ 随着变分键维数 $D$(自 16 外推至 128)的增加而呈现完美的幂律下降,在 $D=128$ 时能量精度已轻松达到 $10^{-5}$ 至 $10^{-7}$。
- 纠缠熵标度分析:临界体系下的纠缠熵应满足共形场论(CFT)的有限纠缠标度律 $S_D \approx \frac{c_{\text{eff}}}{6} \ln \xi_D$。数值拟合结果显示:
- 对于 $Q=\sqrt{3}$,拟合得到的有效中心电荷 $c_{\text{eff}} \approx 2$(具体为 $1.99(5)$),这与物理上该参数下费米海具有 2 个独立费米口袋(Fermi Pockets)的理论预期完全契合。
- 对于 $Q=\sqrt{5}$,拟合得到 $c_{\text{eff}} \approx 1$(具体为 $0.986(4)$),对应系统具有 1 个费米口袋。这极为精确地印证了算法在处理带有非局部弦算符的长程非有理尾部时的完美保真度。
2.2 逆平方 Heisenberg 模型族(含 Haldane-Shastry 点)
研究自旋-1/2 逆平方 Heisenberg 长程模型:
$$ H_H(Q) = \frac{1}{4} \sum_i \sum_{r \ge 1} \frac{\cos(Qr)}{r^2} \vec{\sigma}_i \cdot \vec{\sigma}_{i+r} $$当 $Q=0$ 时,该模型即为著名的 Haldane-Shastry (HS) 链,具有解析的可积基态。对于 $Q \neq 0$,该族模型不存在解析解。
标定数据(图 4、图 5 展示):
HS点($Q=0$)的能量极高精度:利用 TMF 方法,在变分键维数 $D=512$ 下进行优化,算得的能量与精确热力学能量 $E_{HS} = -\pi^2/24$ 相比,残差仅为:
$$ E_\infty - E_{HS} = +2.2 \times 10^{-10} $$这一精度惊人地达到了双精度浮点数的数值极限,直接确立了转移矩阵函数方法的无误差律定属性。
两点关联函数点对点吻合:变分求得的自旋-自旋连通关联函数 $r |C^S(r)|$ 与 Gebhard-Vollhardt 解析解 $3 \text{Si}(\pi r)/(4 \pi r)$ 展现了完美的一致性(如图 4(c) 所示),无任何肉眼可见偏差。
非共度 $Q \neq 0$ 的非平凡物理探索:对于 $Q=0.5, 1, 1.5$,有限纠缠标度分析表明它们均具有 $c_{\text{eff}} \approx 1$ 的共形临界行为,且关联长度 $\xi_D$ 随 $D$ 呈无饱和清洁增长。这一结果首次排除了该族模型被额外调制尺度驱动离共形临界区的可能性,为探索新型长程物理开辟了道路。
2.3 长程铁磁横场 Ising 链(与 SOE-MPO 的深度对决)
为了彻底揭示避免 Hamiltonian 律定偏差的实际物理意义,研究者对比了 TMF 方法与传统 SOE-MPO 拟合在长程铁磁横场 Ising 模型临界点下的行为:
$$ H = -\sum_i \sum_{r \ge 1} \frac{\sigma^x_i \sigma^x_{i+r}}{r^\alpha} - \Gamma \sum_i \sigma^z_i $$在指数 $\alpha = 1.5$ 处,该系统具有长程平均场行为。量子蒙特卡洛(QMC)给出的高精度临界点为 $\Gamma_c = 4.75999$。在这个固定的、已知物理临界点处,对比两种算法计算得到的关联长度 $\xi_D$ 和补偿关联函数 $r^\eta |C(r)|$(理论预期 $\eta = 0.75$ 处应出现长程平台)。
关键性能对比(图 6 核心展示):
| 物理指标 | 传统 SOE-MPO 方案 (极点数 $K=5 \sim 20$, 拟合截止 $R_{\max}=10^4$) | 本文转移矩阵函数 (TMF) 方案 (无拟合无截断) |
|---|---|---|
| 关联长度 $\xi_D$ 随键维 $D$ 的标度 | 随着 $D$ 的增加,$\xi_D$ 发生严重的人为饱和(Artificial Saturation)。例如对于 $K=5$ 或 $K=10$,$\xi_D$ 在很小的 $D$ 处就停止增长。 | $\xi_D$ 随着 $D$ 增加展现出完美、无截断的指数幂律上升,直达计算瓶颈。 |
| 补偿关联函数 $r^\eta \|C(r)\|$ 的行为 | 关联函数仅在极短的距离内遵循幂律,随后在 $r \sim K$ 相关的尺度处发生提前的指数级快速衰减,彻底丢失了长程物理性质。 | 关联函数在极宽的空间尺度(从 $r=1$ 直达 $r=1000$ 级别)内,完美地贴合在 $0.14$ 的物理临界平台上。 |
| 临界性质诊断的科学偏置 | 无法在固定 $\Gamma_c$ 下判断系统是否真正临界。为了掩盖这一由于有理极点拟合导致的偏差,必须对每一个 $K$ 重新调节临界点 $\Gamma_c(K)$,从而引入严重的系统误差。 | 能够纯净、无偏地在真实 QMC 临界点处观察到极其完美的共形场论和标度行为,彻底解耦了张量变分截断与哈密顿量截断。 |
3. 代码实现细节、复现指南与开源生态
本算法已在 Julia 语言生态下完全开源。其核心设计思想是利用 Julia 强大的多态性与张量网络生态,无缝接入现有的变分流形优化框架。
3.1 核心开源库与技术栈
论文关联的开源复现仓库及依赖链如下:
- 核心复现仓库:transfer-matrix-imps (GitHub)
- 张量代数基础设施:TensorKit.jl — 用于处理带有任意 Abel(如 $U(1)$)和非 Abel(如 $SU(2)$)对称性的稠密和对称张量体系。
- iMPS 变分流形算法库:MPSKit.jl — 提供了尖端的、基于 Riemannian 梯度流的无限尺寸基态变分流形算法(VUMPS)。
- 流形优化器:OptimKit.jl 及 TensorKitManifolds.jl — 实现了基于等距等价类的 Riemannian L-BFGS 变分能量最小化控制。
3.2 算法核心复现逻辑与伪代码剖析
在变分更新步骤中,最核心的模块是矩阵函数动作的计算。以下展示了如何在 Julia 中利用 KrylovKit.jl 对转移矩阵函数及其 Fréchet 导数进行正向评估与伴随反向传播:
using KrylovKit
using GSL # 用于高精度评估标量多重对数函数 Li_α
# 1. 定义标量核函数 F_α,Q(z) 及其导数
function polylog_F(α, Q, z)
if abs(z) < 1e-12
return exp(im * Q)
end
# 使用 GSL 库或多倍精度库计算 Li_α(e^{iQ} * z)
return gsl_sf_polylog(α, exp(im * Q) * z) / z
end
function d_polylog_F(α, Q, z)
# 利用复现附录中的解析导数公式:
# F'_{α,Q}(z) = [ Li_{α-1}(e^{iQ} z) - Li_α(e^{iQ} z) ] / z^2
term1 = gsl_sf_polylog(α - 1, exp(im * Q) * z)
term2 = gsl_sf_polylog(α, exp(im * Q) * z)
return (term1 - term2) / (z^2)
end
# 2. 正向克罗夫子空间评估
function evaluate_matrix_function_action(T_A, b, m, α, Q)
# T_A 是一个高维的连通转移矩阵线性映射
# b 是源右向量 |R_Y>
# m 是克罗夫子空间最大维数
# 使用 Arnoldi 过程构建克罗夫子空间
V_m, H_m, β, L_m = arnoldi(T_A, b; maxiter = m)
# 对低维上 Hessenberg 矩阵 H_m 进行对角化
θ, U_right = eigen(H_m)
U_left = inv(U_right)
# 在 Ritz 基下计算标量函数值
F_θ = [polylog_F(α, Q, t) for t in θ]
# 重构变分作用结果 y = F(T_A)b
# y ≈ V_m * U_right * diag(F_θ) * U_left * (V_m^† * b)
coeff = U_right * (F_θ .* (U_left[:, 1] * norm(b)))
y = V_m * coeff
return y, θ, V_m, H_m
end
# 3. 基于 Daleckii-Krein 的反向 Fréchet 共轭计算
function compute_frechet_adjoint(θ_L, θ_R, y_bar, V_L, V_R, α, Q)
# 计算分裂差分矩阵 (Daleckii-Krein Table)
n_L, n_R = length(θ_L), length(θ_R)
F1 = zeros(ComplexF64, n_L, n_R)
for i in 1:n_L, j in 1:n_R
diff = θ_L[i] - θ_R[j]
if abs(diff) < 1e-10 # 避免退化特征值的 Lorentz 展宽发散
F1[i, j] = d_polylog_F(α, Q, θ_L[i])
else
F1[i, j] = (polylog_F(α, Q, θ_L[i]) - polylog_F(α, Q, θ_R[j])) / diff
end
end
# 得到低维子空间伴随算符 H_bar
# 进而将其拉回并提升至高维转移矩阵空间 M_bar
# M_bar_full ≈ V_L * H_bar * V_R^†
# ... (细节参见附录 C 的两步提升逻辑)
end
3.3 复现指南
要想完美重构论文中的全部插图数据,用户可遵循以下流程:
- 克隆官方 GitHub 仓库并使用 Julia 激活环境:
git clone https://github.com/qi-yang-physics/transfer-matrix-imps.git cd transfer-matrix-imps julia --project=. - 实例化项目以自动下载、配置包括
TensorKit和MPSKit在内的全部依赖项:using Pkg; Pkg.instantiate() - 运行对应的基准测试脚本。例如,要复现图 2 中关于 Haldane-Shastry 模型的 Krylov 维度收敛性曲线,直接执行:
julia --project=. scripts/haldane_shastry_krylov_convergence.jl - 优化长程横场 Ising 模型并比较 SOE 与 TMF 的临界行为(图 6):该脚本会自动调用 python 环境中的
julia --project=. scripts/ising_tf_comparison.jlTeNPy模块进行 SOE 多极点拟合作为基线,并在 Julia 端运行无截断转移矩阵优化,最终将数据存为 HDF5 并使用内置脚本绘制成高品质矢量图。
4. 关键引用文献与局限性批判评论
4.1 关键引用文献
在学术谱系上,这项突破性工作深深扎根于并整合了以下几项基石性研究:
- 长程相互作用 MPO 的经典 SOE 拟合:
- Crosswhite et al., Phys. Rev. B 78, 035116 (2008) [16];
- Pirvu et al., New J. Phys. 12, 025012 (2010) [17];
- Zaletel et al., Phys. Rev. B 91, 165112 (2015) [18]。
- 评价:这些文献奠定了长程相互作用张量网络计算的标准范式(SOE-MPO),是本工作致力于超越的经典基准。
- 矩阵函数的数学理论与 Krylov 近似:
- N. J. Higham, “Functions of Matrices: Theory and Computation” (2008) [34]。
- 评价:该权威著作提供了本工作计算转移矩阵函数的数学根基(包括 Fréchet 导数、非正规矩阵稳定对角化等)。
- 张量网络伴随微分与 dominant eigensolver 自动求导:
- Xie, Liu, Wang, Phys. Rev. B 101, 245139 (2020) [39];
- Daleckii & Krein, Am. Math. Soc. Transl. 47, 1 (1965) [71]。
- 评价:Xie 等人的工作打通了 iMPS 主特征向量的自动微分瓶颈,而 Daleckii-Krein 公式则是本工作解决 Ritz 值简并、确保微分不发散的关键工具。
4.2 局限性批判性评论
虽然转移矩阵函数(TMF)方法在物理精确性上树立了新的里程碑,但在实际科学计算中,它仍具有不容忽视的局限性,技术作者对此进行了以下深入的批判性反思:
极弱衰减指数下的红外发散控制($\alpha \le 1$ 的困境) 当衰减指数 $\alpha \le 1$ 时(例如极长程库仑力相互作用或极小幂律系统),在非振荡情况 $Q=0$ 下,标量常数项 $\text{Li}_\alpha(1)$ 会发生本质红外物理发散。尽管这一发散在物理上可以通过引入中性背景电荷(背景抵消法)或适当的能量重整化来消除,但在转移矩阵函数的数值处理中,这会导致常数和读出算符在投影前的系数极其巨大,造成严重的数值精度截断和溢出风险。如何将中性化背景电荷的处理天然嵌入到转移矩阵内核中,是该算法亟待突破的一大方向。
临界 gap 关闭导致的 Krylov 维度 $m$ 的规模膨胀 在强临界区域内,iMPS 的纠缠急剧增加,其物理关联长度 $\xi_D$ 趋于无限。在转移矩阵的表现上,这意味着其最大本征值(主极点)与其子最大本征值之间的 Gap 极度收窄(转移通道能隙关闭)。这在数值数学上会导致 Arnoldi 投影收敛速度显著下降。为了保持给定的计算精度,必须将克罗夫子空间维度 $m$ 扩展到 $128$ 甚至 $256$ 以上。这会显著增加变分计算的内存消耗(需要存储大量高维 Krylov 基向量 $V_m$),使得算法在大键维数 $D$ 下的成本红利有所打折。
多位置元与非均匀长程体系的几何拓展难度 本工作目前仅完整推导并实现了单位置元(Single-site)平移对称 iMPS 体系。虽然论文指出将其推广至大胞元(Multi-site unit cell)在理论上是平铺直叙的,但在实际中,当大胞元具有内部非均匀局部相互作用时,连通转移矩阵将变成一个空间几何异质的算符乘积。此时,转移矩阵将失去其统一的标量核 $F_{\alpha,Q}(z)$ 形式,每一个局部键都需要构建专属的、包含多种空间路径交叉的矩阵函数,这会导致代码实现的几何复杂度呈指数级增加。
5. 补充深度物理与数理分析:Fréchet 微分定理的代数证明
为了给变分量子化学与张量网络方向的科研人员提供最扎实的数理支持,本节补充给出本方法得以成立的核心数理基石——Fréchet 导数的分块矩阵表示公式(Eq. B3)的严格数学推导与物理内幕。
5.1 物理内幕:为什么传统的求导在张量网络中失效?
在变分张量网络(如 VUMPS 或 iDMRG)中,我们要对如下泛函进行极小化:
$$ E[A] = \langle L_X[A]| F(\widetilde{T}_A) |R_Y[A]\rangle $$如果我们直接对基于 Krylov 迭代生成的张量求导(即沿计算图追溯三对角化或 Arnoldi 的每一步):
$$ \frac{\partial E}{\partial A} = \sum_{k=1}^m \frac{\partial E}{\partial v_k} \frac{\partial v_k}{\partial A} $$由于 Arnoldi 过程中包含大量的正交化操作(Gram-Schmidt 过程),其导数链条会包含无数的分母项 $1/\|v_k\|$。当基矢接近正交极限或 Ritz 值发生微小简并时,这些分母会引入极大的不稳定性。这就是为什么流传的张量自动微分在遇到特征值求解器(Eigensolver)和 Krylov 投影时频频崩溃的本征原因。
5.2 严格证明:分块矩阵函数的幂级数展开
Fréchet 导数的巧妙之处在于,它通过代数重构,将不稳定的动力学求导转化为了一个静态的分块矩阵函数求值问题。我们要证明如下关系式(即论文中的 Eq. B3):
$$ F \left( \begin{pmatrix} M & \delta M \\ 0 & M \end{pmatrix} \right) = \begin{pmatrix} F(M) & \mathcal{L}_F(M; \delta M) \\ 0 & F(M) \end{pmatrix} $$证明: 假设函数 $F(z)$ 具有解析的泰勒展开式(对于非有理多重对数函数,在其收敛域内完全满足):
$$ F(z) = \sum_{n=0}^\infty a_n z^n $$定义分块矩阵 $X$ 为:
$$ X = \begin{pmatrix} M & \delta M \\ 0 & M \end{pmatrix} $$我们来计算其 $n$ 次幂 $X^n$ 的代数形式。对于 $n=1$:
$$ X^1 = \begin{pmatrix} M & \delta M \\ 0 & M \end{pmatrix} $$对于 $n=2$:
$$ X^2 = \begin{pmatrix} M & \delta M \\ 0 & M \end{pmatrix} \begin{pmatrix} M & \delta M \\ 0 & M \end{pmatrix} = \begin{pmatrix} M^2 & M \delta M + \delta M M \\ 0 & M^2 \end{pmatrix} $$对于 $n=3$:
$$ X^3 = \begin{pmatrix} M^2 & M \delta M + \delta M M \\ 0 & M^2 \end{pmatrix} \begin{pmatrix} M & \delta M \\ 0 & M \end{pmatrix} = \begin{pmatrix} M^3 & M^2 \delta M + M \delta M M + \delta M M^2 \\ 0 & M^3 \end{pmatrix} $$通过数学归纳法,易证对于任意正整数 $n$,其 $n$ 次幂的右上方分块具有如下严格代数形式:
$$ X^n = \begin{pmatrix} M^n & \sum_{k=0}^{n-1} M^k \delta M M^{n-1-k} \\ 0 & M^n \end{pmatrix} $$现在将此代入泰勒展开式中,计算矩阵函数 $F(X)$:
$$ F(X) = \sum_{n=0}^\infty a_n X^n = \begin{pmatrix} \sum_{n=0}^\infty a_n M^n & \sum_{n=1}^\infty a_n \left( \sum_{k=0}^{n-1} M^k \delta M M^{n-1-k} \right) \\ 0 & \sum_{n=0}^\infty a_n M^n \end{pmatrix} $$显然,左上角和右下角分块精确对应 $F(M)$。而东北角(右上角)分块为:
$$ [F(X)]_{12} = \sum_{n=1}^\infty a_n \sum_{k=0}^{n-1} M^k \delta M M^{n-1-k} $$根据定义,这恰好是算符一阶导数的定义。设 $\epsilon$ 为无穷小量,计算 $F(M + \epsilon \delta M)$ 的一阶展开:
$$ F(M + \epsilon \delta M) = \sum_{n=0}^\infty a_n (M + \epsilon \delta M)^n = \sum_{n=0}^\infty a_n \left( M^n + \epsilon \sum_{k=0}^{n-1} M^k \delta M M^{n-1-k} + O(\epsilon^2) \right) $$$$ = F(M) + \epsilon \left[ \sum_{n=1}^\infty a_n \sum_{k=0}^{n-1} M^k \delta M M^{n-1-k} \right] + O(\epsilon^2) $$由此立得:
$$ \mathcal{L}_F(M; \delta M) \equiv \lim_{\epsilon \to 0} \frac{F(M + \epsilon \delta M) - F(M)}{\epsilon} = \sum_{n=1}^\infty a_n \sum_{k=0}^{n-1} M^k \delta M M^{n-1-k} $$这与分块矩阵东北角元素 $[F(X)]_{12}$ 完全恒等。证毕。
5.3 这一代数结果对计算的颠覆性贡献
该数学证明为数值梯度计算提供了一条无比优美的通道:为了求得连通转移矩阵函数的变分梯度 $\bar{M}$,我们完全不需要对转移矩阵进行物理上的微扰和求导。只需把当前变分得到的张量 $A$、梯度伴随 $\bar{y}$ 以及通过 Krylov 正交化得到的左右投影子空间基矢带入,直接计算这个双线性的分块矩阵函数,并通过计算其低维 Ritz 基下的 Daleckii-Krein Table(分层求导),即可一揽子、无误差、无物理发散地导出张量 $A$ 的最新变分下降方向。这堪称应用数理方法在强关联张量网络领域的一次绝美演出!