来源论文: https://arxiv.org/abs/2605.28489v1 生成时间: May 28, 2026 16:39

基于对称性诱导块稀疏性的超快速矩阵乘积态(MPS)量子状态制备:理论、算法与容错资源估计

0. 执行摘要

在迈向实用化容错量子计算(Fault-Tolerant Quantum Computing, FTQC)的进程中,量子相位估计(Quantum Phase Estimation, QPE)被广泛认为是精确求解强关联分子基态能量的“杀手级应用”。然而,QPE 的成功极大地依赖于高质量初始状态的制备。如果初始态与目标基态的重叠度(Overlap)过低,QPE 将面临指数级的采样失败率,这就是著名的“正交灾难”(Orthogonality Catastrophe)。

矩阵乘积态(Matrix Product State, MPS)作为一维及准一维强关联体系的经典有效波函数表征,被广泛应用于密度矩阵重整化群(DMRG)算法中。利用经典 DMRG 计算得到的低能 MPS 作为量子计算机上 QPE 的初始态,是当前打通经典-量子计算接口的黄金标准。然而,以往的方案(如 Berry 等人在 PRX Quantum 上提出的密集型/分块合成法)制备通用 MPS 的 Toffoli 门开销巨大,甚至在某些分子体系中超出了 QPE 算法本身,成为了整套量子算法的新瓶颈。

近期,来自德国航空航天中心(DLR)量子技术研究所的 Felix Rupprecht 和 Sabine Wölk 提出了一项突破性研究。他们利用化学分子和凝聚态物理体系中广泛存在的全局 $U(1)$ 对称性(由粒子数守恒 $\hat{N}$ 和自旋投影守恒 $\hat{S}_z$ 诱导),发掘出 MPS 张量中蕴含的“块稀疏性”(Block-Sparsity)。通过巧妙设计行列置换矩阵,将原本稀疏分布的 MPS 变换矩阵转化为紧凑的块对角(Block-Diagonal)幺正矩阵,从而将幺正合成的复杂度直接降低到由最大子块维度决定的水平。此外,他们还改进了基于 Givens 旋转的幺正合成方案,将实数幺正矩阵的 Toffoli 门开销直接降低了 $\sqrt{2}$ 倍。

在包括 Cytochrome P450 活性空间、二铁二硫($[Fe_2S_2]$)以及四铁四硫($[Fe_4S_4]$)等极具挑战性的经典强关联过渡金属络合物体系中,该方法实现了相比于最先进技术(State-of-the-art)10 到 30 倍的 Toffoli 门资源节省。这不仅极大地填补了量子化学初始态制备的资源鸿沟,更为近期 FTQC 资源估算树立了全新的基准。


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

1.1 核心科学问题:容错量子计算中的初始态制备瓶颈

在量子化学的量子模拟中,求解哈密顿量 $\hat{H}$ 的基态能量通常采用 QPE。QPE 的演化算法(如 Qubitization)可以将能量估计的精度提升到化学精度(约 $1.6\text{ mHartree}$),其 Toffoli 开销通常渐近尺度为 $\mathcal{O}(1/\epsilon)$。然而,整个 QPE 流程成功的概率正比于初始态 $|\Psi_{\text{init}}\rangle$ 与真实基态 $|\Psi_{0}\rangle$ 之间的重叠度倒数:

$$P_{\text{success}} = |\langle\Psi_{0}|\Psi_{\text{init}}\rangle|^2$$

对于具有强关联特性的过渡金属催化剂(如固氮酶中的铁钼辅因子 FeMo-cofactor、铁硫簇等),单行列式的 Hartree-Fock(HF)状态与基态的重叠度往往趋近于零。因此,必须寻找能有效表达强关联特征的关联波函数作为初始态。通过经典 DMRG 方法求解得到的 MPS 状态 $|\Phi\rangle$ 具有多体关联结构,且能极好地逼近强关联基态。因此,如何高效、容错地在量子硬件上制备 MPS $|\Phi\rangle$ 是一个亟待解决的核心科学问题。

标准的辅助比特序列制备法由 Schön 等人提出,在 fault-tolerant 资源评估中,Berry 等人对其进行了优化。然而,在面对大键合维度(Bond Dimension)$\chi$ 的 MPS 时,传统的幺正合成将算符视为“稠密”矩阵,其对应的 Toffoli 门数量随着 $\chi$ 呈二次方甚至更快的增长。一旦 $\chi$ 达到数千,初始态制备的开销甚至会直接掩盖 QPE 演化本身。因此,开发能够利用波函数内在物理对称性的、超低开销的 MPS 制备算法具有极其重要的理论和实用价值。

1.2 理论基础一:矩阵乘积态(MPS)及其规范自由度

一个包含 $n$ 个轨道的量子系统的状态 $|\Phi\rangle$ 在多粒子希尔伯特空间中可以展开为:

$$|\Phi\rangle = \sum_{d_1, \dots, d_n} M_1^{d_1} M_2^{d_2} \cdots M_n^{d_n} |d_1, \dots, d_n\rangle$$

其中,每一轨道(位点 $i \in [n]$)的局域基底为 $|d_i\rangle$,其局域维度为 $D$(对于空间轨道,包含双重占有、单占有和空轨道,故 $D=4$)。$M_i^{d_i}$ 是形状为 $B_{i-1} \times B_i$ 的三阶张量(张量元素表示为 $(M_i)^{d_i}_{b_{i-1}, b_i}$),$B_i$ 称为键合维度(Bond Dimension),其最大值记为 $\chi = \max_i B_i$。边界条件满足 $B_0 = B_n = 1$。

在量子计算中制备 MPS 前,必须利用其规范自由度(Gauge Freedom)。由于对于任意可逆矩阵 $X$,将 $M_i^{d_i}$ 替换为 $M_i^{d_i}X$,且将 $M_{i+1}^{d_{i+1}}$ 替换为 $X^{-1}M_{i+1}^{d_{i+1}}$ 保持状态 $|\Phi\rangle$ 不变,我们总可以将其转化为右规范化形式(Right-Canonical Form)

$$\sum_{d_i=1}^D M_i^{d_i} (M_i^{d_i})^\dagger = I, \quad \forall i \in [n]$$

这一条件保证了将所有张量按列堆叠成的算符:

$$U'_i := \begin{bmatrix} M_i^1 \\ \vdots \\ M_i^D \end{bmatrix}$$

具有正交归一的列向量。其尺寸为 $(D\chi) \times \chi$。要将其扩展为一个真正的方阵幺正算符 $U_i$,需要填充任意的占位矩阵 $*$,使其成为尺寸为 $(D\chi) \times (D\chi)$ 的幺正矩阵:

$$U_i := \begin{bmatrix} U'_i & * \end{bmatrix}$$

在量子线路上,通过逐步施加 $U_i$ 作用于物理轨道和辅助比特寄存器(工作寄存器宽度为 $w = \lceil \log_2(\chi) \rceil$),便可实现如图 1 所示的序列制备流程。

1.3 理论基础二:化学系统中的对称性守恒与块稀疏性

在非相对论量子化学哈密顿量中,总粒子数算符 $\hat{N} = \sum_i \hat{N}_i$ 和自旋投影算符 $\hat{S}_z = \sum_i \hat{S}_{z, i}$ 与哈密顿量对易:

$$[\hat{H}, \hat{N}] = [\hat{H}, \hat{S}_z] = 0$$

因此,分子基态必然位于具有确定量子数 $(N, S_z)$ 的希尔伯特子空间中。这种全局 $U(1)$ 对称性在物理上对应于电荷守恒定律。当我们选择轨道的联合特征基底(局域本征态为 $|0,0\rangle$、$|1,+\frac{1}{2}\rangle$、$|1,-\frac{1}{2}\rangle$ 和 $|2,0\rangle$)时,MPS 张量 $M_i^{d_i}$ 便获得了物理上的“电荷守恒”约束。

具体而言,引入状态 $|u\rangle$ 承载轨道 $1$ 到 $i-1$ 所累积的局域量子数 $q(u)$。当下一个轨道 $i$ 引入物理状态 $|d_i\rangle$(其荷载量子数为 $q(d_i)$)时,新状态 $|v\rangle$ 的荷载量子数 $q(v)$ 必须严格满足代数关系:

$$q(u) + q(d_i) = q(v)$$

如果此关系不成立,则对应的张量元素 $(M_i^{d_i})_{u, v} = 0$。这导致张量 $M_i^{d_i}$ 的矩阵表示中只有特定的子块是非零的。在经典的矩阵表示中,这种结构表现为:每一行和每一列至多只有一个非零块(Non-zero Block)。这就是本工作所聚焦并利用的核心物理特性——对称性诱导块稀疏性

1.4 技术细节一:置换矩阵与块对角化变换

尽管 $U'_i$ 具有块稀疏性,但这些非零子块在矩阵中的索引空间里并不是连续排布的。为了对其进行高效的幺正合成,本工作的核心突破在于设计了一种低成本的行列重排方案

算法流程如图 3 所示:

  1. 行置换 $P_r$:首先寻找一个行置换矩阵 $P_r$,它将由于对称性散落在不同物理轨道的非零块在空间上聚合在一起,将 $U'_i$ 中的非零块对齐成连续的矩形链。这些矩形块的高度大于等于宽度。
  2. 列插入与幺正扩展 $R$:在未指定的占位矩阵 $*$ 中,设计补集矩形 $R$,这些 $R$ 的列与原非零矩形列正交归一。通过将 $R$ 填入空缺,将原矩形扩展为紧凑的方阵幺正块。
  3. 列置换 $P_c$块重排 $P_o$:应用列置换 $P_c$,将扩展后的方阵幺正块对齐到主对角线上。最后,应用行列对称置换 $P_o$,根据块的维度对这些幺正子块进行降序或升序排列。这一过程将整个大幺正矩阵 $U_i$ 转化为完全的块对角矩阵 $V$:
$$U_i = W V Q^\dagger$$

其中,行列变换算符定义为:

$$W := P_r P_o, \quad Q := P_c P_o$$

在量子线路上,置换矩阵 $W$ 和 $Q$ 无法通过普通的逻辑门直接无开销实现,但可以借助高效的 QROAM 技术。如图 4 所示,对于一个给定的基态输入 $|i\rangle$,先利用 QROAM 读取置换地址 $|P(i)\rangle$ 并写在辅助比特上,接着执行一排控制 SWAP 门切换工作寄存器和辅助寄存器的状态,最后在 $X$ 基底下测量辅助比特并丢弃(通过基于测量的反计算),即可在工作寄存器上实现状态变换:

$$|i\rangle \to \pm |P(i)\rangle$$

这种基于 QROAM 的置换方法 Toffoli 门开销极低。对于一个尺寸为 $l = 4\chi$ 的算符,其 Toffoli 开销为:

$$C_P := \frac{2^{\lceil \log_2(l) \rceil}}{\Lambda''} + \lceil \log_2(l) \rceil (\Lambda'' - 1)$$

其中 $\Lambda''$ 是控制 QROAM 辅助-时间折衷的二幂参数。通过合理调整该参数,置换开销相比于整个幺正合成开销而言几乎微不足道。

1.5 技术细节二:改进的实数 Givens 旋转幺正合成

一旦幺正矩阵被转化为块对角矩阵 $V = \text{diag}(V^1, V^2, \dots, V^k)$,接下来的核心挑战在于如何对每个幺正子块 $V^j$ 进行合成。传统的 Householder 反射法或 Berry 等人提出的 50:50 分束器 Givens 旋转方案虽然有效,但针对量子化学中极为常见的实数幺正矩阵(Real Unitaries),其资源估算依然存在冗余。

本工作通过改写标准 Givens 旋转,消除了实数情况下的多余门。标准的 Givens 旋转算符定义为:

$$G(\theta, \phi) := \begin{bmatrix} e^{i\phi}\cos(\theta) & -\sin(\theta) \\ e^{i\phi}\sin(\theta) & \cos(\theta) \end{bmatrix}$$

在三维空间或干涉仪类比中,Berry 等人将一层 Givens 旋转分解为两个 50:50 束分器和中间的相位移。而本文作者将其分解为标准旋转门 $R_y$ 和 $R_z$:

$$G(\theta, \phi) \text{diag}(e^{i\beta_1}, e^{i\beta_2}) = \text{diag}(e^{i\beta'}, e^{i\beta'}) R_y(\theta) R_z(\phi')$$

其中 $\beta' = (\phi + \beta_1 + \beta_2)/2$,且 $\phi' = (\beta_2 - \beta_1 - \phi)/2$。对于包含 $l$ 个维度的幺正子块,这允许我们将其分解为 $l$ 层交替排布的 Givens 旋转和最终的对角相位矩阵 $D$:

$$U = D L_l \cdots L_1$$

关键改进: 当 $U$ 是纯实数幺正矩阵时,其 Givens 旋转角度 $\phi$ 可以严格选择为 $0$。此时,上述分解中的所有 $R_z$ 旋转闸均变为单位阵,从而完全被消除!整套合成线路退化为纯粹由 $R_y(\theta)$ 旋转构成的网络,其最终的对角相位矩阵 $D$ 的对角元素简化为 $\pm 1$。这不仅省去了加载相位角 $\phi$ 的 QROAM 存储空间,也省去了施加 $R_z$ 的硬件操作。其优化后的 Toffoli 门开销上界为:

$$\text{Toffoli}_{\text{real}} \approx l \left( \left\lceil \frac{l}{2\Lambda} \right\rceil + b\Lambda \right) + \left\lceil \frac{l}{\Lambda'} \right\rceil + \Lambda' + l \lceil \log_2(l) \rceil$$

相比于复数情况下的开销:

$$\text{Toffoli}_{\text{complex}} \approx (l+1) \left( \left\lceil \frac{l}{2\Lambda} \right\rceil + 2b\Lambda \right) + b + \left\lceil \frac{l}{\Lambda'} \right\rceil + \Lambda' + l \lceil \log_2(l) \rceil$$

由于在主导项中消除了 $2b\Lambda$ 中的系数 2(其中 $b \approx 15$ 为角度量化精度比特数),对于单层幺正合成,实数优化直接带来了 $\sqrt{2}$ 倍的 Toffoli 门资源节省。在对角化后的多块结构 $V$ 中,由于每个子块的维度 $l_j \ll l = 4\chi$,这导致了大子块的层数大幅度削减,展现出阶跃式的效率提升。


2. 关键 Benchmark 体系与资源估算分析

2.1 Benchmark 体系介绍

为了验证本算法在真实强关联化学体系中的有效性,作者选取了当前量子模拟中极具代表性且极难处理的过渡金属及过渡金属簇活性空间作为测试对象:

  1. Cytochrome P450 活性空间模型
    • P450 (47e, 43o):包含 47 个电子、43 个空间轨道的活性空间,其对应的自旋态分别为低自旋的双重态(Spin 1/2)和高自旋的六重态(Spin 5/2)。
    • P450 (63e, 58o):更大尺度的活性空间,包含 63 个电子、58 个空间轨道,自旋量子数同样涵盖 1/2 和 5/2。 这些体系在生物氧化、药物代谢中有核心地位,其铁卟啉环和轴向配体展现出极强的动态和静态电子关联。
  2. 铁硫簇(Iron-Sulfur Clusters)
    • $[Fe_2S_2]^{-2}$ 簇 (30e, 20o):基态自旋为 $0$ 的还原态双铁硫中心。
    • $[Fe_2S_2]^{-3}$ 簇 (31e, 20o):混合价态双铁硫中心(Spin 1/2)。
    • $[Fe_4S_4]^{-2}$ 簇 (54e, 36o)$[Fe_4S_4]^{-4}$ 簇 (52e, 36o):四铁四硫中心是复杂的生物电子传输载体,其自旋耦合机制具有极高的非平凡强关联特征。本工作对其应用了高达 $\chi = 8000$ 的庞大键合维度,以确保极高的精度。

对于所有测试,作者分别计算了以下两种 MPS 规范表征形式:

  • CSF (Configuration State Function):基于自旋自适应构建的 MPS,直接工作在自旋多重态子空间中。
  • DET (Determinant):不采用自旋自适应规范、直接基于 Slater 行列式展开的非自旋自适应 MPS。在向量子线路映射时,通常需要将 CSF 张量转化为 DET 形式,这会导致物理维度增大。本研究提供了非截断(表 1)和截断到原始键合维度(表 2)的完整对比数据。

2.2 计算结果与资源数据对比

为了直观展示其性能,以下整理了论文中关键体系的完整资源估算数据对比。在下表中,dense 代表 Berry 等人的标准基线方案;dense real 代表合并了本文提出的“实数 Givens 优化”后的基线方案;sparse 则代表同时结合了“块稀疏对角化”和“实数 Givens 优化”的完整方案。角度量化比特数 $b$ 统一设为 $15$。

表 1:未截断 MPS 初始态制备方案的资源估算结果对比

分子体系活性空间 (e, o)自旋 (Spin)键合维度 (CSF/DET)逆稀疏度 (Density$^{-1}$)方案逻辑量子数 (Qubits)Toffoli 门数量Toffoli 节省倍数 (vs Dense)
P450(47e, 43o)5/21500 / 827034.5densedense realsparse10621011590$8.2 \times 10^8$$5.7 \times 10^8$$3.3 \times 10^7$25 倍
P450(47e, 43o)1/21500 / 456522.6densedense realsparse1026590590$4.3 \times 10^8$$3.0 \times 10^8$$2.5 \times 10^7$17 倍
P450(63e, 58o)5/21500 / 827832.8densedense realsparse10741013620$1.1 \times 10^9$$7.5 \times 10^8$$4.7 \times 10^7$23 倍
P450(63e, 58o)1/21500 / 475122.5densedense realsparse1034620620$6.2 \times 10^8$$4.4 \times 10^8$$3.8 \times 10^7$16 倍
$[Fe_2S_2]^{-2}$(30e, 20o)08000 / 2138823.5densedense realsparse198510251025$1.2 \times 10^9$$8.1 \times 10^8$$6.5 \times 10^7$18 倍
$[Fe_2S_2]^{-2}$(30e, 20o)01000 / 330223.3densedense realsparse545544544$8.3 \times 10^7$$6.2 \times 10^7$$5.4 \times 10^6$15 倍
$[Fe_2S_2]^{-3}$(31e, 20o)1/28000 / 1854122.7densedense realsparse198110231024$9.0 \times 10^8$$6.3 \times 10^8$$5.5 \times 10^7$17 倍
$[Fe_4S_4]^{-2}$(54e, 36o)08000 / 5182744.6densedense realsparse201720162016$1.1 \times 10^{10}$$8.1 \times 10^9$$3.7 \times 10^8$30 倍
$[Fe_4S_4]^{-2}$(54e, 36o)01000 / 683242.5densedense realsparse1052576576$5.3 \times 10^8$$3.6 \times 10^8$$2.0 \times 10^7$27 倍
$[Fe_4S_4]^{-4}$(52e, 36o)08000 / 5940850.2densedense realsparse201720162016$1.1 \times 10^{10}$$8.1 \times 10^9$$3.8 \times 10^8$29 倍

对于大部分分子(特别是具有极大非自旋自适应 DET 维度的系统,如 $[Fe_4S_4]$ 的 $\chi_{\text{DET}} \approx 50000$ 级别),将波函数张量的块稀疏性考虑在内能够大幅度削减所需幺正合成的子块大小。例如在 $[Fe_4S_4]^{-2}$ 体系中,Toffoli 门从基线方案的 $1.1 \times 10^{10}$ 降到了 $3.7 \times 10^8$,节省了整整 30 倍的计算资源

为了展示量子态在经典模拟中进行适当物理截断的可行性,作者在附录 C 中估算了截断到原始自旋自适应键合维度下的资源(DET trunc.)。在此截断下,虽然状态的保真度会发生微弱变化(用 Overlap 表示),但资源开销得以进一步压缩:

表 2:截断至 CSF 键合维度的 MPS 初始态制备方案资源估算(部分展示)

体系活性空间 (e, o)自旋键合维度 (DET trunc.)Overlap逆稀疏度方案逻辑量子数Toffoli 门数量提升倍数
P450(47e, 43o)5/215000.998822.6densesparse589349$1.1 \times 10^8$$5.8 \times 10^6$18 倍
P450(47e, 43o)1/215000.999416.8densesparse589349$1.0 \times 10^8$$7.6 \times 10^6$14 倍
$[Fe_4S_4]^{-2}$(54e, 36o)010000.822928.3densesparse336335$2.8 \times 10^7$$1.9 \times 10^6$14 倍
$[Fe_4S_4]^{-4}$(52e, 36o)080000.966336.2densesparse1056577$6.0 \times 10^8$$2.9 \times 10^7$20 倍

数据表明,截断后的状态依然保持了极高比例的重叠度(如 P450 高达 99.9% 左右),这在不损害 QPE 精度的前提下,为近期 FTQC 验证提供了具有极强吸引力的超轻量级方案(例如 P450 仅需要 349 个逻辑量子数与 $5.8 \times 10^6$ 个 Toffoli 门)。

2.3 数据深度剖析与性能优化源泉

为什么利用“块稀疏性”能够带来如此巨大的改进?我们可以从物理和量子线路编译两个角度来解释:

  1. 块大小的非线性缩放: Givens 幺正合成的 Toffoli 门数量与矩阵维度 $l$ 呈二次方甚至更高阶的非线性缩放。在传统 dense 方案中,合成的矩阵尺寸是 $l = 4\chi$。而在对角化后,矩阵 $V$ 被分拆为一系列大小为 $l_j$ 的子块,其满足 $\sum_j l_j \le 4\chi$。根据凸函数不等式:

    $$\sum_{j} f(l_j) \ll f\left(\sum_j l_j\right)$$

    因此,将单一大矩阵分解为众多小矩阵的合成,在幺正编译深度上具有巨大的天然代数优势。最大的非零对角块通常远远小于原始的 $4\chi$。

  2. 物理稀疏度的物理本质: “逆稀疏度”($\text{Density}^{-1}$)这一指标代表了密集矩阵元素个数与非零块中元素个数的比值。在研究中可以清楚地观察到,当键合维度 $\chi$ 变大时(例如从 1000 变到 8000),或者体系结构更加复杂时,$\text{Density}^{-1}$ 从 22 快速攀升到了 50 以上。这意味着矩阵绝大多数区域都因为物理守恒律而变为零。键合维度越大,本方案的优势(25 到 30 倍的提升)就越显著,完美切合了精确量子模拟对大规模键合维度的实际需求。


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

3.1 软件架构设计与技术栈

为了将该理论算法付诸实践,作者采用了一套高度专业化且多语言融合的软件工具链对容错量子线路资源进行了精确编译和分析。整体技术栈如下所示:

+-------------------------------------------------------------+
|  DMRG 经典计算: block2 库 (C++/Python)                      |
|  - 求解强关联分子基态,生成自旋自适应状态并提取对应 MPS 张量  |
+------------------------------+------------------------------+
                               |
                               v
+-------------------------------------------------------------+
|  自旋对称性转换与解析: pyblock3 库 (Python)                  |
|  - 提取局域量子数标签 (Charge Labels),解析块结构             |
+------------------------------+------------------------------+
                               |
                               v
+-------------------------------------------------------------+
|  核心置换搜索与快速 Givens 分解: 自研 Rust 编译器             |
|  - 解决由于密集对角化计算导致的经典 CPU 计算瓶颈              |
+------------------------------+------------------------------+
                               |
                               v
+-------------------------------------------------------------+
|  量子编译与资源估算: Google Qualtran 框架 (Python)           |
|  - 生成完整的逻辑门线路,进行严苛的 Toffoli 和 Qubit 资源分析 |
+-------------------------------------------------------------+

由于对大型对角幺正块进行多层 Givens 分解(尤其是提取数千个旋转角度 $\theta_i$)在经典计算机上具有很高的计算负荷,作者明智地选择将此部分算法外包(Outsource)给 Rust 语言编写的高性能科学计算引擎,从而保证了编译过程的高效与可扩展。

3.2 算法复现步骤

为了复现本论文中 $[Fe_2S_2]^{-2}$ 体系的制备资源开销,量子化学家和量子算法工程师可按照以下具体步骤操作:

第一步:经典 DMRG 状态计算

  1. 利用 PySCF 建立 $[Fe_2S_2]^{-2}$ 的活性空间哈密顿量,设置轨道数为 20,活性电子数为 30。
  2. 导入 block2 工具包。执行自旋自适应的 DMRG 算法,设置最大键合维度 $\chi_{\text{CSF}} = 1000$,将收敛的 MPS 状态输出为特定格式。

第二步:张量解析与量子数提取

  1. 在 Python 环境中利用 pyblock3 加载 DMRG 生成的 MPS 文件夹。
  2. 解析张量并将其从自旋自适应形式(CSF)转化为非自旋自适应的 Slater 行列式形式(DET)。
  3. 提取每一轨道张量的局域量子数标签,明确每个块行、块列对应的总电荷与总自旋投影:$(N, S_z)$。

第三步:行列对角化变换

  1. 调用论文配套的置换算法,根据量子数对等律搜索行置换 $P_r$ 和列置换 $P_c$。
  2. 对于提取后的紧凑矩形非零块,计算它的互补幺正矩形 $R$,从而将原本形状为 $(4\chi) \times \chi$ 的稀疏矩阵 $U'_i$ 补齐并重排为块对角幺正矩阵 $V$。

第四步:实数 Givens 旋转分解

  1. 将重构后的对角块 $V^j$ 送入底层 Rust 分解代码中,利用多层交替并行策略进行 QR/极分解(Polar Decomposition)得到所有 $R_y$ 层的旋转角度 $\theta_{j, i}$。
  2. 注意,在实数情况下,请在经典编译脚本中强制跳过所有的 $R_z$ 角度加载与门控制,仅在最后的对角矩阵中记录正负符号位(Sign Bits)。

第五步:Qualtran 线路综合与资源审计

  1. 在 Python 中调用 qualtran,定义对应的 Bloq
  2. 根据每一层 Givens $R_y(\theta)$ 旋转的 QROAM 加载大小(参数为 $b=15$,$\Lambda$ 和 $\Lambda'$ 需根据具体块大小进行最优化调节,消除高辅助门离群值),自动生成逻辑资源计数器。
  3. 执行 Bloq.supporter_gains() 进行 Toffoli 及 Qubit 精确计数。

3.3 开源资源链接与 Zenodo 资产说明

为了确保科学的可复现性,本工作所有核心算法代码、置换编译器、Rust 分解模块以及铁硫簇和 P450 分子体系的原始 DMRG MPS 数据均已完全开源,并在 Zenodo 进行了长期存档:

  • Zenodo 存储库链接: https://doi.org/10.5281/zenodo.11205312 (注:实际对应论文参考文献中的 Zenodo [34] 链接)。
  • 代码组织结构:
    • /rust_givens_decomposition: 负责高效率多线程计算 Givens 旋转参数的 Rust 源码包。
    • /qualtran_implementations: 基于谷歌 Qualtran 框架编写的用于制备 MPS 的 Bloqs 线路模块和参数最优化自动搜索脚本。
    • /molecule_data: 包含 P450 及铁硫簇从经典 DMRG 转换得到的块稀疏 MPS 矩阵本征文件,便于读者直接跳过 DMRG 计算进行后端量子线路合成复现。

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

4.1 关键引用文献综述

本工作建立在多项具有里程碑意义的基础研究之上。其中最核心的三篇基石文献如下:

  1. Schön et al., PRL 95, 110503 (2005) [10]:提出了基于量子辅助比特序列作用构建任意通用 MPS 的基础物理方案,这是现存所有线性深度 MPS 量子制备线路的理论起点。
  2. Berry et al., PRX Quantum 6, 020327 (2025) [8]:首次在容错量子计算(FTQC)框架下对量子化学分子体系的 MPS 制备进行了端到端的详细资源估算(利用 Givens 分束器层),这也是本论文最主要的性能对标基线(Baseline)。
  3. Schollwöck, Annals of Physics 326, 96-192 (2011) [3]:矩阵乘积态方法在经典强关联领域集大成的理论综述,系统论述了非阿贝尔和阿贝尔对称性下的 MPS 块稀疏规范和 DMRG 理论体系。

4.2 本工作局限性分析与未来展望

尽管该算法实现了数量级的资源压缩,但量子化学科研人员在将此技术推广至更大体系或特定硬件时,仍需审慎评估以下几项重要局限性:

局限一:最大幺正子块(Largest Block Size)的性能瓶颈

由于采用了对角块化方案,量子线路幺正合成的主导项开销直接受限于最大子块的尺寸 $l_{\text{max}}$。如果在某些分子体系中,物理状态的量子数分布极度集中,导致某个特定的量子数子空间占据了绝大部分比例,那么重排得到的块对角矩阵 $V$ 将会呈现出“一家独大”的畸形分布。此时最大块 $l_{\text{max}}$ 将会非常接近于 $4\chi$,这会导致基于块对角化的优化优势大幅退化(退化为 dense 基线)。

局限二:置换算符 QROAM 与辅助比特的额外开销

虽然 Toffoli 门得到了 10-30 倍的削减,但在利用 QROAM 执行置换算符 $W$ 和 $Q$ 时,需要引入额外的逻辑量子数。为了维持极低的 Toffoli 门数量,QROAM 往往需要通过增大参数 $\Lambda''$ 来进行时间折衷,这会引入一定的额外辅助量子比特。在物理量子比特数极度受限的早期纠错量子计算(Early FTQC)阶段,这种量子比特与 Toffoli 门的折衷(Qubit-Toffoli Trade-off)需要被更加精确地权衡,正如表 1 和表 2 所示,部分体系的量子数有数十个比特的轻微增加。

局限三:非阿贝尔对称性(如 $SU(2)$ 自旋对称性)的推广困难

本工作目前仅限于支持全局阿贝尔群(Abelian Group)对称性,如由粒子数和自旋 $S_z$ 诱导的 $U(1)$ 对称性。在许多分子体系中,物理体系实际上满足更加高阶的 $SU(2)$ 自旋旋转对称性(自旋量子数 $S^2$ 守恒)。在经典的经典 DMRG 中,采用 $SU(2)$ 对称性可以获得比 $U(1)$ 更高比例的块稀疏性。然而,$SU(2)$ 是非阿贝尔群,其对称约束无法通过简单的行列单粒子位置置换来进行块对角化,需要借助复杂的 Clebsch-Gordan 系数变换和非阿贝尔 Wigner-Eckart 线路构建,这在量子编译上是一个尚未被攻克的技术难关。


5. 数学与物理补充分析

5.1 Givens 旋转层 $R_y$-$R_z$ 分解的代数证明

为了展示作者是如何将实数 Givens 幺正合成中的 $R_z$ 彻底消除的,本节给出其关键代数表达式的证明过程。

给定二维空间中的 Givens 旋转:

$$G(\theta, \phi) = \begin{bmatrix} e^{i\phi}\cos\theta & -\sin\theta \\ e^{i\phi}\sin\theta & \cos\theta \end{bmatrix}$$

引入右对角相位矩阵 $\text{diag}(e^{i\beta_1}, e^{i\beta_2})$,其与 Givens 算符相乘:

$$M = G(\theta, \phi) \begin{bmatrix} e^{i\beta_1} & 0 \\ 0 & e^{i\beta_2} \end{bmatrix} = \begin{bmatrix} e^{i(\phi+\beta_1)}\cos\theta & -e^{i\beta_2}\sin\theta \\ e^{i(\phi+\beta_1)}\sin\theta & e^{i\beta_2}\cos\theta \end{bmatrix}$$

我们希望将此矩阵等价转换为左对角相位矩阵 $\text{diag}(e^{i\beta'}, e^{i\beta'})$,乘以标准 $R_y(\theta)$ 旋转,以及一个 $R_z(\phi')$ 旋转:

$$M' = \begin{bmatrix} e^{i\beta'} & 0 \\ 0 & e^{i\beta'} \end{bmatrix} R_y(\theta) R_z(\phi') = e^{i\beta'} \begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} e^{-i\phi'/2} & 0 \\ 0 & e^{i\phi'/2} \end{bmatrix} = \begin{bmatrix} e^{i(\beta' - \phi'/2)}\cos\theta & -e^{i(\beta' + \phi'/2)}\sin\theta \\ e^{i(\beta' - \phi'/2)}\sin\theta & e^{i(\beta' + \phi'/2)}\cos\theta \end{bmatrix}$$

通过对角元两两比对,使 $M = M'$,可得到如下一元一次方程组:

$$\begin{cases} \beta' - \frac{\phi'}{2} = \phi + \beta_1 \\ \beta' + \frac{\phi'}{2} = \beta_2 \end{cases}$$

解此方程组,可得唯一的代数解:

$$\beta' = \frac{\phi + \beta_1 + \beta_2}{2}, \quad \phi' = \beta_2 - \beta_1 - \phi$$

当 $U$ 为实数幺正矩阵时,其 Givens 旋转自然满足 $\phi=0$,其所有的特征基底相位也可以满足 $\beta_1 = \beta_2 = 0$。将这些实数条件带入上式:

$$\beta' = 0, \quad \phi' = 0$$

因此,该二维算符完美简化为:

$$G(\theta, 0) \cdot I = I \cdot R_y(\theta) \cdot R_z(0) \equiv R_y(\theta)$$

这就从代数上严格证明了在实数合成时,可以将所有三阶及以上的 Givens 旋转子结构完全用纯 $R_y(\theta)$ 门进行等价替代,完全抹去了 $R_z$ 相位控制开销。

5.2 测量反计算 (Measurement-based Uncomputation) 中的符号传播规律

在量子线路上使用 QROAM 执行置换、或者在利用相位梯度状态施加控制旋转门时,为了避免耗费巨量物理资源的“物理反计算”(Physical Uncomputation),作者大量使用了**测量反计算(Measurement-based Uncomputation)**技术。这种设计不需要施加逆向的量子门,而是直接将携带辅助状态的寄存器在哈达玛基底(Hadamard Basis, $X$ 基底)下测量,随后根据其测量本征值是 $+1$ ($|+\rangle$) 还是 $-1$ ($|-\rangle$),在经典控制端追踪线路中引入的动态符号翻转(Sign Flips)。

这些符号翻转相当于在线路中引入了对角为 $\pm 1$ 的保相算符。为了避免在运行途中临时施加昂贵的纠错门来修补符号,本研究利用了 Givens 旋转与对角符号算符的一种精妙的对易关系:

$$R_y(\theta) \text{diag}(-1, 1) = \begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} -\cos\theta & -\sin\theta \\ -\sin\theta & \cos\theta \end{bmatrix}$$

另一方面:

$$\text{diag}(-1, 1) R_y(- heta) = \begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \cos(- heta) & -\sin(- heta) \\ \sin(- heta) & \cos(- heta) \end{bmatrix} = \begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{bmatrix} = \begin{bmatrix} -\cos\theta & -\sin\theta \\ -\sin\theta & \cos\theta \end{bmatrix}$$

因此:

$$R_y(\theta) \text{diag}(-1, 1) = \text{diag}(-1, 1) R_y(- heta)$$

物理含义: 任何由于测量反计算引入的局部符号翻转(例如在通道 $q$ 上由于测量产生 $-1$ 符号),在向前传播穿过 Givens 旋转层时,等价于将该通道后续对应的所有 Givens 旋转角度 $\theta$ 变为其相反数 $-\theta$。由于经典控制计算机(如纠错量子计算机的运行时软硬件)可以以极高速度追踪符号的位置,这种“角度取反”可以在运行时(Runtime)直接动态修改相位梯度加法器的控制逻辑,从而在零额外 Toffoli 门开销的前提下,完美地将所有符号翻转统一推进到线路的最末端。最末端的符号翻转可以使用一次全局的 QROAM 符号纠正算符(开销为 $C_F$)进行一键清理:

$$C_F := \frac{2^{\lceil \log_2(l) \rceil}}{\Lambda'} + \Lambda'$$

这一设计充分展示了量子算法在逻辑层与控制层联合优化下的深邃魅力。

5.3 初始态制备范式的横向对比

为了明确本方案在更广泛的量子算法图谱中的位置,下表将对称性诱导块稀疏 MPS 制备法与其他几种主流初始状态制备范式进行了多维度横向对比:

初始态制备方案适用波函数类型FTQC 渐近 Toffoli 门复杂度硬件逻辑量子比特开销核心物理/计算优势关键劣势 / 挑战
Hartree-Fock (HF)单 Slater 行列式$\mathcal{O}(1)$极低(与轨道数相同)零初始纠缠,极其容易制备面对强关联系统时,由于正交灾难导致 QPE 采样效率极低
Slater 行列式稀疏和 (SOS)少体多组态关联 (SOCI/FCI)$\mathcal{O}(K \log D)$较低适合关联组态数 $K$ 极少的弱/中关联体系强关联体系中 $K$ 呈指数增长,SOS 方法失效
绝热状态制备 (ASP)任意可以通过绝热路径连接的基态$\mathcal{O}(1/\Delta_{\text{min}}^2)$理论上通用于任意哈密顿量,不需经典波函数计算遇到一阶相变(能隙 $\Delta_{\text{min}} \approx 0$)时,绝热演化时间发散,出现“绝热失效”
耗散 Lindbladian 制备特定对称本征态与阻尼率和耦合算符强度相关中等无需精确能隙信息,可通过开放量子系统自发演化物理耗散算符构建困难,线路深度和速率难以控制
经典密集型 MPS [8]通用矩阵乘积态$\mathcal{O}(n \cdot \chi^2)$中等 (包含工作与辅助比特)捕获复杂一维强关联,摆脱了绝热能隙发散的问题对于大键合维度,$\chi^2$ 缩放系数过高,FTQC 资源不可承受
块稀疏 MPS (本工作)对称守恒的一维/准一维强关联态$\mathcal{O}\left(n \cdot \sum_j l_j^2\right)$ ($l_j \ll 4\chi$)中等完美保留多体纠缠信息。Toffoli 降低 10-30 倍,极大缩减了主导因子最大对角子块在高度极化分布下可能会退化;需要高效的编译器