切比雪夫多项式方法的双重应用：高效寻找多自旋系统中数千个中心本征值#

关浩宇 与 张文宪

武汉大学物理科学与技术学院，人工微纳结构教育部重点实验室，湖北武汉 430072

摘要#

计算中间谱段的大量内部本征值是量子多体系统中的一个重要问题，因为能级统计量能提供量子混沌的特征标志。我们提出了一种精确数值方法——切比雪夫多项式双重应用（DACP），可同时寻找数千个中心本征值，其中能级间距随系统尺寸呈指数级减小。为解决近简并问题，我们两次运用切比雪夫多项式：首先构建指数型半圆滤波器作为预处理步骤，然后在目标子空间中生成大量合适的基态。对伊辛自旋链和自旋玻璃碎片的数值计算证明了DACP的正确性和高效性。数值结果表明，对于伊辛自旋链，DACP比当前最先进的移位求逆方法快30倍；对于自旋玻璃碎片则快8倍。与移位求逆方法相比，DACP的计算时间受所需本征值数量的影响很小，这使其成为大规模本征值计算的强大工具。此外，对于包含多达20个自旋的自旋1/2系统，其内存消耗也保持较小的常数值（5.6 GB），因此非常适合并行计算。

第1章 引言#

能级统计是量子混沌的基本刻画手段之一[1,2]。可积系统通常表现出能级聚集和能级间距的泊松分布[2]，而混沌系统则表现出能级排斥和维格纳-戴森分布[3]。其他有用的统计工具如$\delta_{3}$统计量[4]和$\delta_{n}$统计量的功率谱也依赖于能级分布[5]。由于一般情况下无法得到解析结果，数值本征值对于精确表征这些能级统计量至关重要。此外，谱中部的单个本征态对应"无穷温度"极限，在研究多体局域化 (MBL) [6–8]中具有重要意义。数值模拟是定量理解MBL问题诸多方面的重要工具。由于关注的多体问题通常涉及巨大的希尔伯特空间，其维度随系统规模呈指数增长，完全对角化哈密顿量或求解含时薛定谔方程相当具有挑战性。因此，寻求在谱的小部分区域求解本征值问题是一种不可避免且理想的替代方案。

为了获得内部本征态，通常采用矩阵光谱学的混合策略[9–11]。该策略旨在计算谱中选定区域的本征态，并将基态求解器与谱滤波器相结合，其中滤波器被设计为将选定的内部区域变换到谱的边缘。在这些滤波器中，格林函数$( \mathscr{H} - \lambda I )^{- 1}$是一种出色的选择。应用格林函数滤波器后，能量$\lambda$附近的本征值簇被映射到非常大的正值和负值。$\lambda$附近的能级间距被放大，从而改善了收敛性。Lanczos方法[12]与此类滤波器结合使用[9,10]，并实现了格林函数的切比雪夫多项式展开[13]。特别地，移位求逆方法[14]本质上利用了这种谱变换，并广泛应用于量子多自旋系统[8,15–18]。此外，它被认为是处理MBL问题最有效的方法[8]。然而，对于大规模系统，由于需要对$( \mathscr{H} - \lambda I )^{- 1}$进行分解，该方法在计算时间和内存消耗上会迅速增长。还提出了其他谱滤波器，包括狄拉克Delta函数滤波器$\delta ( \mathcal{H} - \lambda I )$[19–23]和矩形窗函数滤波器[24]。

其他更高效的方法被提出来计算位于内部谱的大量本征对。特别地，本征值切片库 (EVSL) 在多项式滤波过程中采用狄拉克Delta函数的切比雪夫展开，并与重启的Lanczos算法相结合[21]。类似地，FILTLAN使用 (非重启) Lanczos 和多项式滤波与Krylov投影方法的组合来计算内部和极端本征值[20]。FEAST算法利用著名的柯西积分公式来表达本征投影算子，从而得到一个有理滤波器，然后对其应用子空间迭代[25]。尽管这些方法能够使用分治策略找到数千或数万个本征对，但它们在一个单一能量区间内只能计算几百或一千个本征对。注意，所有上述方法都是迭代的。计算时间大致与所需的本征值数量成正比，即更多的本征值需要更多的滤波和重新正交化[8,20,22]。

在本文中，我们提出了一种精确的数值方法DACP，用于计算能带中部的数千个本征值，这足以揭示能级统计[26,27]。对于自旋系统，这样的中部区域通常对应态密度的峰值，在此处能级近乎简并。在不分解$\mathcal{H}^{- 1}$的情况下区分这些近乎简并的本征值极具挑战性。在DACP方法中，我们首次利用切比雪夫多项式构造了一个半圆指数 (exp-semicircle) 滤波器，以快速抑制谱中不需要的部分。切比雪夫多项式的第二次应用是快速搜索一组张成特定子空间的状态，该子空间包含所有所需的本征态。结合这两个步骤，DACP方法本质上将原始的高维本征值问题转化为低维本征值问题。与多次迭代滤波不同，DACP方法仅通过单次滤波就能直接产生结果。

对于多自旋系统中的实际问题，DACP方法非常高效，因为它充分利用了切比雪夫多项式的几个优良特性，而其他方法仅利用了其中一部分。对于一大类多自旋系统，与移位求逆方法相比，DACP方法的计算速度显著提高，可达30倍。内存节省更为显著，可达100倍。此外，DACP方法区别于那些迭代滤波方法，因为当所需本征值数量在大范围内变化时，其收敛时间变化不大。

本文的组织结构如下。第2节详细描述了DACP方法的完整形式，包括指数半圆滤波、切比雪夫演化和子空间对角化。每个过程都依赖于切比雪夫多项式的一个特定性质，这使得DACP有别于其他滤波方法。在第3节中，我们使用DACP计算了两个多自旋系统——伊辛模型和自旋玻璃碎片的内部本征值，并给出了数值结果。我们在第4节中将DACP与其他方法进行了比较，并在第5节进行了总结。

2 切比雪夫多项式的双重应用#

为了获取大型自旋系统的中心本征值，我们自然受到无矩阵模式的限制，即哈密顿量的矩阵不能被显式地表示或存储。相反，我们将哈密顿量 $\mathcal{H}$ 视为一个算子，其输入和输出是态/向量。因此，我们仅对一组量子态进行操作：

其中 $k$ 是一个正整数，原始希尔伯特空间的维度最高可达 $5 \times 10^{5}$。需要注意的是，每个“矩阵-向量乘积” $\mathcal{H} \left| \psi \right.$ 作为一个基本操作，消耗大量时间但仅需少量内存。而且，如第2.3节所述，我们甚至不一定需要存储方程(1)中的态。因此，DACP对内存的要求非常小。通常，我们希望同时找到大量中心本征值。在本手稿中，我们的目标是精确找出多自旋系统能谱中部的5,000个本征值，以此作为示例。

图1：$a = 0.1$（红线）、0.05（蓝线）和0.025（黑线）时的指数半圆滤波器，其中 $\begin{array} {l c l} {k a} & {=} & {24} \end{array}$。我们设置 $E_{\mathrm{max}} ~ = ~ 1$ 和 $E_{\mathrm{min}} = - 1$。实线对应切比雪夫多项式 $T_{k} ( F ( x ) )$，其中 $\begin{array} {l} {\displaystyle{F \overline{{( x )}} = \left\lceil 2 x^{2} - ( 1 + a^{2} ) \right\rceil / ( 1 - a^{2} )}} \end{array}$，而虚线对应近似 $y = \exp ( 2 k \sqrt{a^{2} - x^{2}} ) \simeq T_{k} ( F ( x ) )$。水平的粉色点线表示滤波器对数半极大值的位置。注意，当 $x \notin [ - a , a ]$ 时，$| T_{k} ( F ( x ) ) | \leq 1$，因此滤波区域限制在 $[ - a , a ]$ 内。

DACP的思想非常直接。我们首先将一个随机初始化的态变换成一个波包，该波包基本上由至少5,000个中心本征态张成的子空间构成。这种变换是通过一个由切比雪夫多项式实现的指数半圆滤波器完成的。这里，切比雪夫多项式被用来实现指数衰减。有了这个特定的态之后，我们便生成大量线性无关的态（尽可能多），以近似张成由所需本征态构成的子空间。切比雪夫多项式再次被用来振荡（复指数）这些态，从而产生近似线性无关的态。一旦获得了所需本征态的生成集，我们就可以显式地计算出 $\mathcal{H}$ 在该子空间中的约化表示 $H$。剩余的操作限制在维度约为 $10^{4}$ 的子空间内。最后，对大小为 $10^{4} \times 10^{4}$ 的 $H$ 进行直接对角化，即可给出所需的本征值。详细步骤和讨论见以下各小节。

2.1 指数半圆滤波器#

我们利用切比雪夫多项式在区间 $[ - 1 , 1 ]$ 之外的指数增长特性，高效地构建了一个指数半圆滤波器，如图1所示。该滤波器能显著放大任何随机初始化态中所需本征态范围的分量，从而产生一个在能谱中部高度局域化的新态。我们注意到，切比雪夫滤波器也利用了相同的性质，不同之处在于它放大了能谱的低端[28,29]。类似的思想也曾应用于寻找基态的量子算法[30]中。

对于能量范围在 $[ E_{\operatorname* {m i n}} , E_{\operatorname* {m a x}} ]$ 的哈密顿量 $H$，其中 $E_{\mathrm{min}}$ 是最小值（$E_{\mathrm{min}} < 0$），$E_{\mathrm{max}}$ 是最大值（$E_{\operatorname* {m a x}} > 0$），指数半圆滤波器旨在放大对应于区间 $[ - a , a ]$ 内本征值的本征态分量，同时抑制区间 $[ E_{\operatorname* {m i n}} , - a ]$ 和 $[ a , E_{\mathrm{max}} ]$ 内的分量，其中 $a$ 是一个正实数参数。针对自旋系统，我们假设 $E_{\mathrm{min}} < - a$ 且 $a < E_{\operatorname* {m a x}}$。滤波后，生成一个主要由本征值属于区间 $[ - a , a ]$ 的本征态构成的新态。为简单起见，我们将由包含在 $[ - a , a ]$ 内的本征态张成的子空间记为 $\mathbb{L}$。因此，位于 $\mathbb{L}$ 之外的本征态的概率几乎可以忽略，正如图1和图2所示。

现在我们引入具体的实现细节。注意，我们想要放大能谱的中部，但切比雪夫多项式的指数增长仅存在于两端附近。为克服这一困难，我们首先将哈密顿量平方，得到$\mathcal{H}^{2}$，其能谱范围为$[ 0 , E_{\mathrm{max}}^{2} ]$（为简单起见，假设$E_{\operatorname* {m i n}} = - E_{\operatorname* {m a x}}$）。$H$的中部能谱$[ - a , a ]$对应$\mathcal{H}^{2}$的$[ 0 , a^{2} ]$，这恰好位于低端。下一步是通过对$\mathcal{H}^{2}$的平移和归一化，将阻尼部分$[ a^{2} , E_{\mathrm{max}}^{2} ]$映射到$[ - 1 , 1 ]$中。因此，我们定义了一个算符

其中$E_{c} = ( E_{\mathrm{max}}^{2} + a^{2} ) / 2$，$E_{0} = ( E_{\mathrm{max}}^{2} - a^{2} ) / 2$。通过将$\mathcal{H}^{2}$替换为$a^{2}$或$E_{\mathrm{max}}^{2}$，并相应地得到$F ( x ) = ( x^{2} - E_{c} ) / E_{0}$，可以很容易地验证此映射的正确性。注意$\mathcal{F}$只是$H$的一个多项式表达式，因此$T_{k} ( \mathcal{F} )$也是如此。

接下来，我们探讨使用$T_{k} ( \mathcal{F} )$进行滤波的效果。由于$\mathcal{H}^{2}$在$[ 0 , a^{2} ]$内的本征值被映射到$\mathcal{F}$的$[ - 1 - 2 a^{2} / ( E_{\mathrm{max}}^{2} - a^{2} ) , - 1 ]$区间，对于能谱的低端，我们得到$T_{k} ( \mathcal{F} ) = ( - 1 )^{k} \cosh ( k \Theta )$，其中

设$\begin{array} {r} {\left| \psi \right. = \sum_{i} c_{i} \left| \phi_{i} \right. + \sum_{j} d_{j} \left| \chi_{j} \right.} \end{array}$为一个随机初始态，其中$c_{i}$和$d_{j}$为随机系数，$| \phi_{i} \rangle$为$\mathbb{L}$内部的本征态，$\left| \chi_{j} \right.$为$\mathbb{L}$外部的本征态。由$T_{k} ( \mathcal{F} )$进行的滤波为

其中 $\theta_{i}^{\mathrm{in}} = \cosh^{- 1} ( 1 + 2 ( a^{2} - E_{i}^{2} ) / ( E_{\operatorname* {m a x}}^{2} - a^{2} ) ) , \theta_{j}^{\mathrm{out}} = \cos^{- 1} ( 2 ( E_{j}^{2} - a^{2} ) / ( E_{\operatorname* {m a x}}^{2} - a^{2} ) - 1 ) , E_{i} = 0 ,$ 且$E_{j}$分别是与$| \phi_{i} \rangle$和$| \dot{\chi}_{j} \rangle$对应的本征值。在写等式(4)时，我们忽略了$( - 1 )^{k}$，因为它不影响系数的绝对值，且在最后一行是一个全局相位。当$a$很小时，即$a^{2} \ll E_{\operatorname* {m a x}}^{2}$，通过对$\cosh^{- 1} ( 1 + \varepsilon )$进行泰勒展开（其中"是一个小的正数），可以进一步推导出$\theta_{i}^{\mathrm{in}} \simeq 2 \sqrt{a^{2} - E_{i}^{2}} / E_{\mathrm{max}}$。因此，我们得到了指数-半圆滤波器

对于满足$- a \leq E_{i} \leq a$的本征态，该滤波器在$E_{i} = 0$处随着$k$增大而尖锐地达到峰值。在图1中，公式(5)的形状以虚线表示，与精确结果吻合良好。

基于初始条件$T_{0} ( \mathcal{F} ) = 1$和$T_{1} ( \mathcal{F} ) = \mathcal{F}$，可以使用递推关系式(A.3)高效地确定$k$阶切比雪夫多项式。在本文中，我们设定截断阶数$K = 12 E_{\mathrm{max}} / a$。这样的滤波器以指数速度（所有多项式中最快的速率）放大$\mathbb{L}$内本征态的分量[31]。在归一化之后，它等效地抑制了$\mathbb{L}$外部的分量，从而生成目标态：

图2：三种滤波器在区间⁻0.25, 0.25内效果的比较。我们展示了初始随机归一化态 $\textstyle | \psi \rangle = \sum_{i} c_{i} | \phi_{i} \rangle$（黑色曲线），以及经过指数半圆滤波器（$K = 240$，红色曲线）、矩形函数的切比雪夫展开（$\Theta ( x + 0.1 ) \Theta ( - x + 0.1 )$，阶数 $K = 480$，蓝色曲线）和狄拉克δ函数 $\delta ( 0 )$（$K = 480$，粉色曲线）滤波后的最终归一化态。指数半圆滤波器的参数 $a = 0.1$。x轴对应"约化能量" $E / E_{\mathrm{max}}$，y轴对应概率幅 $c_{i}$ 的绝对值。矩形函数平移了 $10^{- 10}$，δ函数平移了 $10^{- 15}$。三种滤波器使用相同数量的矩阵-向量乘积。考虑到目标幅度的放大，指数半圆滤波器的优势显而易见。

其中 $c_{i}^{\prime} = \beta \mathrm{e}^{K \theta_{i}^{\mathrm{in}}} c_{i}$，$\beta$ 是归一化常数。显然，态 $| \psi_{E} \rangle$（在能量表示中）定域在能谱中间。我们将 $| \psi_{E} \rangle$ 作为第2.2节切比雪夫演化的初始态输入。

为进一步说明指数半圆滤波器的优势，我们将其与另外两种函数的切比雪夫展开进行比较。第一种是狄拉克δ函数 $\delta ( 0 )$，曾被用于高效寻找内部本征值[21–23, 32]。第二种是矩形函数，也被用于内部本征值计算[24]。选取一个随机初始化的态 $\textstyle | \psi \rangle = \sum_{i} c_{i} | \phi_{i} \rangle$，我们在图2中绘制了经三种滤波器滤波后生成的态。所有三种滤波器使用相同数量的矩阵-向量乘积。δ滤波器的目标区间很小。然而，就我们的目的而言，重要的是目标区间（_[⁻0.1, 0.1_]）与其他区间（目标区间之外）之间的相对幅度。因此，指数半圆滤波器的效率显然是得到确认的。

2.2 切比雪夫演化#

在获得一个限制于小子空间 $\mathbb{L}$ 的态之后，我们利用切比雪夫多项式的振荡特性来高效地生成一组尽可能多的不同态，作为张成 $\mathbb{L}$ 的一组完备基。为实现这一目标，有必要将 $E_{i}$ 限制在 $[ - 1 , 1 ]$（对应于 $T_{k} ( x )$ 中的x），在此区间内切比雪夫多项式表现为类余弦函数。这一区域与指数半圆滤波器的要求相反，因此需要对H进行不同的变换。下面描述切比雪夫多项式第二次应用的具体细节。

原始哈密顿量H需要平移 $E_{c}^{\prime}$ 并按 $E_{0}^{\prime}$ 重新缩放：其中 $\begin{array} {r} {E_{c}^{\prime} = \frac{1} {2} ( E_{\mathrm{min}} + E_{\mathrm{max}} )} \end{array}$ 且 $\begin{array} {r} {E_{0}^{\prime} = \frac{1} {2} ( - E_{\mathrm{min}} + E_{\mathrm{max}} )} \end{array}$ 。与式(2)类似，我们定义一个算子 $\mathcal{G} = ( \mathcal{H} - E_{c}^{\prime} ) / E_{0}^{\prime}$ ，它显然以⁻1和1为界。再次假设 $E_{\operatorname* {m i n}} = - E_{\operatorname* {m a x}}$，得到

这样，参数a被重新缩放为 $a_{r} = a / E_{\operatorname* {m a x}} .$

让我们探讨切比雪夫演化，该演化由算子 $T_{k} ( \mathcal G )$ 控制，其中k扮演时间的角色。我们将滤波生成的态式(6)作为切比雪夫演化的初始态输入。由于 $\| \mathcal{G} \| \leq 1$，从式(A.2)可得 $T_{k} ( \mathcal{G} ) = \cos ( k \Omega )$，其中 $\Omega = \operatorname{arccos} ( \mathcal{G} )$ 。在此意义上，切比雪夫演化变为

其中 $\omega_{j} = \operatorname{arccos} ( E_{j} / E_{\operatorname* {m a x}} )$ 。注意 $\omega_{j}$ 和k都是无量纲的。当然，这不是一个物理演化，$T_{k} ( \mathcal G )$ 甚至不是一个酉算子。实际上，这一演化本质上代表了正向和反向传播的叠加。每当地多项式阶数k增加1，演化"时间"也增加1。

借助切比雪夫演化，我们能够构造张成子空间 $\mathbb{L}$ 的一组完备基。具体地，我们收集如下一组态

其中 ${\hat{X}} = \pi{\mathcal{G}} / a_{r}$，n为由关系式 $2 n + 1 \geq d$ 确定的整数，d为$\mathbb{L}$的维数。实际应用中，可能需要关系式 $2 n + 1 = 1.5 d$ 来确保式(9)的完备性。此处 $k \simeq m \pi / a_{r}$ 作为时间，其中 $m = 1 , \cdots , n$ 。类似于第2.1节，计算k阶切比雪夫多项式 $T_{k} ( \mathcal G )$。截断阶数（演化时间）为 $K^{\prime} = \lfloor n \pi / a_{r} \rfloor$ 。更多细节可见附录B。

切比雪夫多项式兼具近似三角函数和多项式性质的双重特征，在DACP方法中起着至关重要的作用。为最大限度地区分那些聚集的特征值，我们需要一组斜率尽可能陡峭的基函数，这等价于剧烈振荡。从某种意义上说，切比雪夫演化是对量子振荡的高效模拟（可能是所有多项式中最高效的）[33]。将算符 $\mathcal{G}$ 替换为实变量x，式(9)在 $E = 0$ 附近的最大斜率为 $\frac{n \pi} {a_{r}}$ ，与（最常见的）算符集合

给出的 $\frac{1} {a_{r}}$ 形成鲜明对比。我们发现，陡峭的斜率有助于区分近简并的特征值。此外，在附录E所述仿真条件下，使用式(10)而非式(9)的基进行切比雪夫演化时，如图8所示收敛性相当差。这些事实清晰地说明了切比雪夫多项式第二类应用的优越性。

第2.3节 子空间对角化#

结合指数半圆滤波器和切比雪夫演化计算基 $\{| \Psi_{i} \rangle : i = 1 , \cdots , 2 n + 1 \} \ ( {\mathrm{E q . ~}} ( 9 ) )$，既是DACP方法中最具挑战性的方面，也是最耗时的部分。一旦构建了合适的基，剩余任务就很简单了，即计算投影哈密顿量H的本征对。这等价于求解广义特征值问题

这里，H和S分别表示$\mathbb{L}$中的投影哈密顿量和重叠矩阵，

Λ是一个对角矩阵，其特征值位于 $[ - a , a ]$ 范围内，矩阵B将式(9)中构建的基变换为H的本征态 $\left| \phi_{j} \right.$，

所有这些矩阵的大小均为 $( 2 n + 1 ) \times ( 2 n + 1 )$。由于其规模较小，可通过LAPACK库[34]轻松对H进行对角化。

重要的是，由于切比雪夫多项式的特殊性，矩阵H和S的计算甚至可以在不显式计算和存储状态 $\left| {\Psi_{i}} \right.$ 的情况下完成。这一特性使得DACP方法在计算时间和内存方面均可得到进一步改进。此外，对于过完备基式(9)，重叠矩阵S通常奇异。因此我们采用奇异值分解。奇异值分解后，若S的特征值的绝对值低于截断条件 $\varepsilon = 10^{- 12}$ ，则将其丢弃。剩余特征值的数目有效地计入了式(9)中线性独立状态的个数。我们在附录C中给出了矩阵H和S的显式表达式（仅用 $T_{k} ( \mathcal G )$ 和 $| \psi_{E} \rangle$ 表示），以及广义特征值问题的解法。

通过子空间对角化获得的特征值可能不具有与系统真实特征值相同的精度，因此有必要进行独立检查以估计误差界。为此，如果本征态 $\left| \phi_{j} \right.$ 已知，残差范数（能量的方差）$| | r_{j} | | = \sqrt{\langle \mathcal{H}^{2} \rangle - \langle \mathcal{H} \rangle^{2}}$ ，其中 $\left. \mathcal{H}^{2} \right. = \left. \phi_{j} \right| \mathcal{H}^{2} \left| \phi_{j} \right.$ 和 $\langle \mathcal{H} \rangle = \left. \phi_{j} \middle | \mathcal{H} \middle | \phi_{j} \right.$ ，被广泛用作衡量结果精度的参数。已有研究表明，$| | \boldsymbol{r}_{j} | |$ 给出了计算特征值真实误差（绝对误差）的上界[35]。

第3章 数值结果#

我们将DACP方法应用于具有两体相互作用的量子自旋-1/2系统。这类系统是研究量子计算、固体理论和量子统计中一大类重要问题的良好模型[36–39]。大量精确的特征值有助于我们获得统计性质，区分量子混沌与可积性，并作为评估其他近似方法的基准。

一般而言，DACP方法可以处理由N个自旋中任意一对自旋之间耦合组成的自旋系统。泡利矩阵 $\sigma^{\alpha}$ 或两个耦合泡利矩阵 $\sigma^{\alpha} \otimes \sigma^{\beta}$ ，其中 $\alpha , \beta = x , y , z_{;}$ ，均由特定函数适当表示。

我们指定两个物理系统的自旋模型。一个是无序一维横向场伊辛模型[40]，其哈密顿量为

其中 $\sigma_{i}$ 是作用于自旋 i 的泡利矩阵。该系统可通过 Jordan-Wigner 变换[38]精确求解，因此是 DACP 方法的理想正确性检验器。最近邻交换相互作用常数 $J_{i , i + 1}$ 是均匀分布在 $[ - J / \sqrt{N} , J / \sqrt{N} ]$ 内的随机数，其中 $J = 10$。局部随机磁场由 $\Gamma_{i}^{z}$ 表示，它们是均匀分布在区间 $[ 0 , \Gamma ]$ 内的随机数，其中 $\Gamma = 1$。

另一个系统是自旋玻璃碎片[27]，它代表了一类全局范围相互作用的系统，需要 DMRG 方法[41]采用相对较大的键维度来处理。描述该系统的哈密顿量为

所有符号和参数与上述伊辛模型相同，只是第一个求和项遍历所有可能的自旋对，且 $J = 0.866 \Gamma$。这个系统很有趣，因为它呈现出从可积分性到量子混沌再回到可积分性的两次交叉。在极限 $J / \Gamma \to 0$ 下，基态是顺磁性的，所有自旋都沿局部磁场方向，系统是可积的[27]。在相反极限 $J / \Gamma \to \infty$ 下，基态是自旋玻璃态，并且由于存在 N 个与哈密顿量对易的算符 $(\sigma_{i}^{x})$，系统也是可积的[27]。在这两个极限之间存在一个量子混沌区域。当 $N = 20$ 时，$J = 10 \Gamma$ 大致是量子混沌到可积（自旋玻璃侧）的边界，而 $J = 0.866 \Gamma$ 表明系统处于纯量子混沌相[27]。

通过使用上界估计量（它需要很少的额外计算就能确定最大绝对值特征值 $E_{0}$ 的上界），可以估计 $E_{\mathrm{max}} = E_{0}$ 和 $E_{\mathrm{min}} = - E_{0} \ [ 42 ]$。对于此设置，我们利用了多自旋系统中态密度关于零中心呈钟形分布的对称性。由于我们需要 5,000 个中心特征值，可以设置 $n = 4 , 000$，对应维度为 8,001，足以跨越整个子空间 $\mathbb{L}$。重叠矩阵 S 通常是奇异的。DOS $\rho ( E )$ 的近似分布可以通过时间演化波函数的傅里叶变换或参考文献[43]中给出的更好估计方法高效计算。参数 a 被适当选择，以确保包含在 $[ - a , a ]$ 内的本征态数量略少于 8,000（如图 3 所示，某些收敛特征值的精度可能低于要求）。

在实际操作中，有时可能存在高度近简并的特征值，其能级间距小至 $10^{- 7} \Gamma$，而平均间距为 $10^{- 5} \Gamma$。对于切比雪夫演化来说，区分这样接近的特征值对仍然很困难（需要长两个数量级的时间）。为了应对这一挑战，我们在 DACP 方法的编程中采用了阻塞滤波技术 29，这意味着一个状态块被“同时”滤波或演化。其思想是，简并子空间中的两个或多个随机状态通常是线性无关的。对于每次数值测试，随机生成一个包含 5 个初始试验状态的状态块，并将参数 n 相应调整为 $n = 800$。

通过这些设置，我们对上述两个系统进行了数值测试，以展示 DACP 方法的精确性和效率。在本工作中，我们只考虑特征值的计算。本报告中所有的时间信息均来自 Intel(R) Xeon(R) CPU E7-8880 v4 上的顺序模式计算。

在图 3 中，我们展示了伊辛模型（$N = 19$）和自旋玻璃碎片（$N = 17$）的相对误差 $\eta$ 随系统能量 $E_{\mathrm{exact}}$ 的对数刻度图。我们将计算出的中心特征值 E 的相对误差 $\eta$ 定义为

两个系统的精确特征值已通过其它可靠方法获得。对于伊辛模型，我们利用著名的 Jordan-Wigner 变换，将原始的 $2^{N} {\times} 2^{N}$ 矩阵简化为 $2 N \times 2 N$ 矩阵，并还原出原始哈密顿量的完整谱 [38]。对于自旋玻璃碎片，我们直接使用 MATLAB 的 eigs 函数，找出最接近 $E = 0$ 的 5,000 个特征值。在我们的数值测试中，图3的参数 a 在 (a) 和 (b) 中分别为 $0.036\Gamma$ 和 $0.16\Gamma$。计算伊辛模型时，按上述设置，满足 $\eta < 10^{-6}$ 的特征值数量不足（少于 5,000 个）。我们将块大小扩展至 10，参数 n 增至 500，从而收集到足够的特征值。收敛的特征值（即满足 $\eta < 10^{-6}$ 的计算特征值）数量，对于 (a) 为 5,385 个，对于 (b) 为 5,000 个（即我们拥有的所有精确特征值），而计算的特征值总数分别为 6,232 个和 5,910 个。

图3：计算特征值的相对误差 η 的对数刻度图，(a) 为 $N = 19$ 的伊辛模型，(b) 为 $N = 17$ 的自旋玻璃碎片模型。横轴为系统能量。满足 $\eta < 10^{-6}$（黑色虚线）的特征值数量在 (a) 中为 5,385 个，在 (b) 中为 5,000 个。这种分布源于图1中滤波器的形状。

两图中 $E = 0$ 附近的尖峰是由于分母 $E_{\mathrm{exact}}$ 极小所致。最小绝对特征值在 (a) 中约为 $4.4 \times 10^{-6}\Gamma$，在 (b) 中约为 $2.9 \times 10^{-5}\Gamma$。此外，图中间部分存在一个平坦的平台，表明对于 $E = 0$ 附近的那些特征值，我们遇到了数值误差，即绝对误差达到了双精度表示的限制。在图3(a)中可以观察到其它尖峰，而在(b)中并未出现。每个这样的尖峰对应于一簇极为接近的特征值。伊辛模型的可积性意味着能级聚集和能级间距的泊松分布。我们推测切比雪夫演化对于近简并特征值效果不佳，因为它们对最终状态的相位贡献非常相似。要显著放大相位差 $\exp(i\Delta E t)$，需要 $t \gtrsim 1 / \Delta E$，其中 $\Delta E$ 为能量差。另一种解决方案是增加块大小以匹配近简并子空间的维度。忽略中央平台和尖峰，η 的分布呈现出 exp-半圆滤波器的倒置形状。我们在附录E中讨论它们的关系。

在图4中，我们进一步展示了两个系统能级间距 $s$（相邻特征值之差，以平均间距为单位）的统计分布 $P(s)$。图4(a)中的统计分布由 5,000 个中心特征值（最接近 $E = 0$）构建，而图4(b)中的分布则由 2,500 个相同对称类（式(15)中的相互作用不混合自旋向上数量为偶数和奇数的状态）的中心特征值获得。横向场伊辛模型是可积的，表现出能级间距的泊松分布 $P(s) = e^{-s}$。另一方面，对于自旋玻璃碎片模型，我们在式(15)中设 $J = 0.866\Gamma$，对应于量子混沌相 [27]。事实上，图(b)中的分布遵循 Wigner-Dyson 高斯正交系综 (GOE) 统计。

图4：能级间距 $s$ 的统计分布，(a) 为 $N = 19$ 的伊辛模型，(b) 为 $N = 18$ 的自旋玻璃碎片模型。概率分布 $P(s)$ 来自 DACP 的数值结果（红色圆圈）、精确特征值（蓝色方块）、泊松统计 $P(s) = \exp(-s)$（(a)中黑色曲线）和 Wigner-Dyson GOE 统计（(b)中黑色曲线）。横轴 s 以平均能级间距为单位。

在图5中，我们展示了对于规模 $N = 15$ 的两个系统，计算时间与收敛特征值数量的关系。横轴表示每次测试中满足条件 $\eta(E) < 10^{-6}$ 的计算特征值数量。回顾一下，DACP 方法分为三个部分：exp-半圆滤波、切比雪夫演化和子空间对角化。我们将它们的时间消耗分别记为 $T_F$、$T_E$ 和 $T_D$，而求解时间为总计算时间 $T = T_F + T_E + T_D$。可以立即注意到，对于两个系统，计算时间 T 在左侧下降，在中间区域达到平台，而在右侧迅速增加。但虚线表示的 $T_E$ 在整个区域呈现出非常小的增长率，因为横轴采用对数刻度。在区域的中间部分，即 DACP 方法表现最佳的地方，求解时间 $T$ 与演化时间 $T_E$ 大致重合。

如果考虑以下四个事实，这些特征就很容易理解。(i) 注意到自旋系统的态密度 $\rho ( E )$ 通常是峰值在零处的钟形曲线，而 a 是一个典型的小参数（对于 $N = 19$，$a / E_{\mathrm{max}} \sim 10^{- 2}$），因此可以安全地将 $[ - a , a ]$ 区间内的 $\rho ( E )$ 视为常数 $\bar{\rho}$。因此，所需特征值的数量为 $R \simeq 2 \bar{\rho} a$。(ii) 将哈密顿量作用于态 $\mathcal{H} \left| \psi \right.$ 作为基本步骤，并将该基本步骤的计算时间记为 τ，则过滤时间可计算为 τ 乘以截止阶数 $K = 12 E_{\mathrm{max}} / a$ 的两倍，即 $T_{F} \simeq 24 E_{\mathrm{max}} \tau / a \propto R^{- 1}$。(iii) 类似地，切比雪夫演化时间为 $T_{E} \simeq \lfloor R \pi E_{\mathrm{max}} / a \rfloor \tau \propto R^{0}$，因为 $R \propto a$。(iv) 众所周知，完全对角化时间与矩阵大小的立方成正比，即 $T_{D} \propto R^{3}$。这些计算时间的组合清晰地解释了图5所示的行为。对于左侧区域，指数-半圆过滤占主导地位，因为 $T_{F}$ 与参数 R 成反比。然而，当进入特征值数量为几百到几千的中间区域时，演化时间 $T_{E}$ 占据了计算时间的大部分，而其他两项可以忽略不计。如图5所示，DACP 方法在计算100到3000个特征值时性能大致相同。随着 R 的持续增加，子空间对角化时间 $T_{D}$ 最终消耗了大部分计算时间。

我们指出，图5中中间区域的平台将 DACP 方法与那些迭代过滤方法区分开来[10,20,23,24,28,29]，尽管它们都使用了切比雪夫多项式。迭代方法通常需要更多的过滤和重正交化次数，因此在寻找更多特征值时收敛时间更长[8,20,22]。事实上，如文献[8]所示，可以发现对于移位-反演方法，计算时间大致与所需特征值的数量成正比。因此，DACP 方法对于大规模特征值计算非常理想。

图5：使用 DACP 方法，计算时间（CPU 秒）与满足条件 $\eta < 10^{- 6}$ 的特征值数量的关系，针对伊辛模型（带方块的红色线，左轴）和自旋玻璃碎片（带圆圈的蓝色线，右轴）。两个系统均为 $N = 15$。实线显示求解时间 $T = T_{F} + T_{E} + T_{D}$，包括指数-半圆过滤的计算时间 $T_{F}$、切比雪夫演化时间 $T_{E}$ 和子空间对角化时间 $T_{D}$。虚线显示 $T_{E}$

最后，我们注意到 DACP 方法消耗的内存相当小，因为它以矩阵-向量乘积模式工作，避免了哈密顿量的显式矩阵表示。此外，如附录 C 所示，方程(9)中的整个状态集也不保存在内存中。主要的内存消耗是存储矩阵元 $H_{i j}$ 和 $S_{i j}$，因此 DACP 方法的内存需求仅与子空间 $\mathbb{L}$ 的维度 d 有关。根据本研究的设置，在计算 $N \leq 20$ 系统中约5000个特征值时，它大约占用5.6 GB 的内存。

第4节 与其他方法的比较#

我们将 DACP 方法与其他已建立的方法在计算第3节中介绍的两种模型时的性能进行了比较。具体来说，我们选择了三种不同类型的方法，列举并讨论如下。

移位-反演方法广泛用于计算量子自旋系统能谱中部的本征对，例如文献[8,15–17]中所述。在我们的测试中，它通过 Matlab R2019b 的“eigs”函数实现，该函数使用隐式重启动 Arnoldi 方法 (ARPACK) [44]。哈密顿量 H 的矩阵完全存储在内存中，因此可以预期其计算速度快，但代价是内存使用率高。

特征值切片库 (EVSL) 提供了用于计算实对称特征值问题给定区间内特征值的例程[21]。我们选择了“ChebLanTr”例程，它基于多项式过滤，并结合了 Krylov 子空间方法和子空间迭代算法[32]。由于它允许无矩阵格式，我们实现了与 DACP 方法完全相同的矩阵-向量乘积函数。然而，在 EVSL 中，矩阵-向量乘积模式 (SetAMatvec) 仅限于应用于实向量，而量子态通常是复向量。因此，图中显示的 EVSL 的所有计算时间都是实际记录时间的两倍。

FEAST 算法是一种通用的特征值求解器，其灵感来自量子力学中的密度矩阵表示和围线积分技术[45, 46]。在我们的测试中，我们安装了 FEAST 包 v4.0 并调用了求解器函数“zfeast_hcsrev”。虽然可以使用无矩阵模式，但它要求用户提供与矩阵-向量例程相结合的直接/迭代线性系统求解器。为简单起见，我们将哈密顿量矩阵存储在稀疏 CSR 格式中。在测试过程中，我们发现当需要计算的特征值数量为几千时，它可能不收敛，因此我们将原始搜索区间划分为几个较小的区间。

关于DACP方法，为了更好地识别其扩展行为，我们去掉了子空间对角化的时间消耗，这是一个仅依赖于子空间大小的常数。需要注意的是，简单地加减一个常数不会影响扩展行为。具体而言，子空间对角化在每次测试中大约消耗 700 CPU 秒，对于 $N \ge 16$ 的系统，这样的常数可以忽略不计。尽管如此，我们在虚线中绘制了DACP方法的完整计算时间。

在比较扩展行为之前，我们发现EVSL和FEAST方法都无法在单次运行（一个搜索区间）中找到5,000个特征值。因此，必须划分搜索区间以获得更好的性能。我们假设每个较小的区间包含相似数量的特征值，其数量由它们的计算性能决定。

表1：通过对图6中曲线进行线性拟合得到的扩展常数α

图6：以计算时间（CPU秒）与系统大小 N 的关系衡量的扩展行为，针对 (a) Ising模型和 (b) 自旋玻璃碎片模型。在两个图中，我们比较了FEAST算法（粉色菱形）、shift-invert方法（蓝色圆圈）和EVSL（黑色三角形）与DACP方法（红色方块）。每个数值测试（以不同标记绘制）展示了寻找5,000个中心特征值时的计算时间。实线是通过线性拟合得到的。对于DACP方法，红色虚线表示完整求解时间 T，而红色实心方块表示演化时间 $T_{E}$

在图6中，我们比较了DACP方法与三种不同类型方法（即shift-invert、EVSL和FEAST方法）以计算时间 T（CPU时间，秒）衡量的扩展行为。图6中的每个数值测试（标记）本质上代表了寻找5,000个特征值的计算时间。在测试中，FEAST的参数设置为 $M_{0} = 1.2 R$。考虑到这两种方法的性能以及为了方便起见，我们将EVSL每个搜索区间内的特征值数量设为500，FEAST设为1,000。DACP方法和shift-invert方法都成功在单次运行中找到了5,000个特征值。EVSL方法也通过“分而治之”技术找到了5,000个特征值。然而，由于其较长的计算时间，我们将FEAST设置为寻找1,000个特征值（单个区间），总计算时间是其5倍。

如图6所示，对于 $13 \leq N \leq 17$ 的自旋系统，DACP方法在这两种模型中比其它方法快约20倍，除了与MATLAB在自旋玻璃碎片模型上的比较。由于shift-invert方法本质上是在寻找 $\mathcal{H}^{- 1}$ 的基态能量，而 $\mathcal{H}^{- 1}$ 通常不是一个稀疏矩阵，其计算时间 T 不受哈密顿量稀疏性的影响，图6(a)和(b)中的两条蓝线几乎相同。相反，由于DACP方法采用了作用在状态上的 H 的多项式组合，其计算时间 T 严重受特定哈密顿量的泡利算符数量影响。例如，当 $N = 18$ 时，自旋玻璃碎片的计算时间 T 是Ising模型的5.5倍，而泡利算符数量是6.2倍。考虑到这种影响，在交换相互作用涉及所有可能自旋对及所有三个方向的最坏情况下，DACP方法仍然优于shift-invert方法。

此外，这四种方法之间的扩展性是可比的。为了定量比较，我们通过拟合数值结果提取扩展常数α，其中 α 定义为 $T = T_{0} \exp ( \alpha N )$ 。α的值如表1所示。确实，对于自旋玻璃碎片模型，扩展常数非常接近，而EVSL和FEAST在Ising模型上具有较小的α。假设 α 对于更大的 N 保持不变，那么DACP方法在时间成本上的优势可能会在相当大规模的系统中持续存在。具体而言，对于Ising模型，FEAST和DACP之间的扩展线在 $N \simeq 28$ 处交叉，而对于玻璃模型，EVSL和DACP之间的扩展线在 $N \simeq 80$ 处交叉。此外，正如文献8中shift-invert方法的数值测试所示，寻找 $\mathcal{H}^{- 1}$ 的因式分解时间主导了其他计算步骤。仅考虑因式分解部分，文献[8]中的执行时间显示出扩展常数 $\alpha \simeq 1.66$，这表明对于更大自旋的系统，扩展行为更差。因此，DACP方法的时间效率得到了确认。

最后，在表2中，我们比较了这些方法在 $N \geq 16$ 的Ising系统中的内存消耗。由于子空间矩阵元素构成了DACP方法的主要存储，其内存需求在测试期间保持恒定。更精确地说，$N = 18$ 系统的量子态占用约4 MB内存，且内存中仅保留约10个量子态。主要项是两个矩阵S和H的完整存储，这随着N的增加保持恒定。相反，对于shift-invert方法，$\mathcal{H}^{- 1}$ 的矩阵大小随着系统尺寸N的增加而快速增加（正比于 $4^{N} )$ ，需要大量内存。例如，在我们的测试中，$N = 17$ 的系统消耗了573 GB内存，这已经是DACP方法的100倍。由于EVSL方法也以矩阵自由模式运行，并且一次性仅计算500个特征值（子空间维度 d $'$ 1, 000），它所需的内存最少。使用稀疏矩阵格式，FEAST方法节省了大量内存，但仍显著多于矩阵自由方法。

表2：四种方法在Ising模型上的内存消耗。

5 结论#

我们提出了DACP方法，通过使用两次Chebyshev多项式，高效计算大规模（至少5,000个）中心特征值，并具有高精度。它探索了Chebyshev多项式的几个卓越性质，以有效滤除不需要的幅度，并为谱中间部分构建合适的基态。特别地，所提出的方法专门用于求解谱的中间部分（$E = 0$ 附近的特征值）。这一限制使得可以将Chebyshev多项式的不同性质结合在一起。最终，两个关键特征使DACP区别于多项式滤波方法：其计算时间对所需特征值数量的依赖性较弱，而其内存开销与系统尺寸无关。此外，如附录D所示，DACP方法比应用了exp-semicircle滤波器的Lanczos方法更稳定、更高效。

针对Ising模型和自旋玻璃碎片的数值测试证实了DACP方法的准确性和效率。与广泛使用的shift-invert方法相比，DACP方法显著提高了计算速度，对于Ising模型提高了高达30倍，而对于自旋玻璃碎片，速度提升较小但仍然显著（8倍）。此外，内存需求大幅降低，对于 $N = 17$ 的自旋系统降低了高达100倍，对于更大的系统则降低更多。作为一种计算中心特征值的强大工具，DACP方法可能在许多物理问题中找到潜在应用，例如凝聚态物理中的多体局域化[7,15,17,47]，以及量子混沌中的能级统计[1,16,27,48]。

致谢#

我们感谢X.-H. Deng的讨论。本文中的数值计算是在武汉大学超级计算中心的超级计算系统上完成的。

资助信息 本工作得到了国家自然科学基金（NSFC）资助号91836101、U1930201和11574239的支持。

A Chebyshev多项式#

在本文中，我们限制自己仅允许对H进行多项式组合的操作。我们利用Chebyshev多项式来完成上述任务。它是将 $\mathcal{H}^{k}$ 的组合与 ${\mathcal{H}}$ 的近似指数函数（无论是 exp(kH) 还是 exp(ikH)）联系起来的桥梁。由于其几个显著的性质，Chebyshev多项式可能会充分发挥这类方法的潜力。

图7：$k = 3$ (红线)、$k = 4$ (蓝线)、$k = 5$ (黑线) 的 Chebyshev多项式 $T_{k} ( x )$。第一行((a-c)实线)说明了 $T_{k} ( x )$ 在区间 $[-1, 1]$ 内的振荡，而第二行((d-f)虚线)显示了Chebyshev多项式在 $[-1, 1]$ 外的快速增长。

第一类k阶Chebyshev多项式定义为

初始条件为 $T_{0} ( x ) = 1$ 和 $T_{1} ( x ) = x ~ [ 49 ]$。它是一个包含两种不同表达式的分段函数，同时作为多项式函数，它既是连续的也是光滑的。为简单起见，当 $x \in [ - 1 , 1 ]$ 时，我们设 $\theta = \cos^{- 1} ( x )$ (cos $\theta = x )$；当 $x \in [ 1 , \infty )$ 时，设 $\theta = \cosh^{- 1} ( x )$，相应的 $\theta$ 范围分别为 $\theta \in [ 0 , \pi ]$ 和 $\theta \in [ 0 , \infty )$。用变量 $\theta$ 表示，方程(A.1)变为

我们很容易观察到，在区间 $[-1, 1]$ 内，$T_{k} ( x )$ 是一个类似正弦或余弦的振荡函数，值域介于 $-1$ 和 1 之间，如图 7(a-c) 所示，而在区间 $[-1, 1]$ 之外，它增长得极快，如图 7(d-f) 所示。

注意 $\cosh ( k \theta ) = ( \mathrm{e}^{k \theta} + \mathrm{e}^{- k \theta} ) / 2$。因此可以预见切比雪夫多项式在区间 $[ - 1 , 1 ]$ 之外呈指数增长。实际上，已知在满足可比条件的情况下，所有次数 $\leq k$ 的多项式中，切比雪夫多项式 $T_{k} ( x )$ 在区间 $[-1, 1]$ 之外增长最快 [31]。

与这些性质相关的一个实用性质是：$T_{k + 1} ( x )$ 可以通过三项递推高效地确定

切比雪夫多项式的所有这些性质使其成为一个强大的工具，并对 DACP 方法非常有用。

B 方程(9)的详细推导#

在此我们推导出通过切比雪夫演化构造集合方程(9)的显式表达式。我们重点关注 $a \ll E_{\operatorname* {m a x}}$ 的情况，这对于大型（$N \ge 15$）自旋系统是相当合理的。

利用当 $x$ 很小时，$\arccos ( x ) = \pi / 2 - x + o ( x )$ 的事实，我们在方程(B.1)的第三行得到 $\omega_{j} = \operatorname{arccos} ( E_{j} / E_{\mathrm{max}} ) \simeq \pi / 2 - E_{j} / E_{\mathrm{max}}$。因此，我们得到表达式

然后，我们进行一个截断阶数为 $K^{\prime} = \lfloor n \pi / a_{r} \rfloor$ 的切比雪夫演化，当 $k = \lfloor m \pi / a_{r} \rfloor$ 且 $m = 1 , \cdots , n$ 时，记录 $T_{k - 1} ( \mathcal{G} ) | \psi_{E} \rangle$ 和 $T_{k} ( \mathcal{G} ) \left| \psi_{E} \right.$。演化结束后，状态集合方程(9)便自动生成。

C 矩阵元的评估与广义特征值问题的求解#

如附录 B 所示，我们使用切比雪夫多项式重写基组方程(9)，即 $\left\{\left| \Psi_{i} \right. : i = 1 , \cdots , 2 n + 1 \right\}$：

其中 $k_{m} = \lfloor m \pi / a_{r} \rfloor , m = 1 , \cdots , n , \lvert \Psi_{1} \rangle = \lvert \psi_{E} \rangle , \lvert \Psi_{2} \rangle = T_{k_{1} - 1} ( \mathcal{G} ) \lvert \psi_{E} \rangle , \lvert \Psi_{3} \rangle = T_{k_{1}} ( \mathcal{G} ) \lvert \psi_{E} \rangle$ 等等。为简单起见，可以进一步假设 $k_{m}$ 是偶数。

矩阵元 $S_{i j} = \left. \Psi_{i} | \Psi_{j} \right. = \left. \psi_{E} | T_{x} ( \mathcal{G} ) T_{y} ( \mathcal{G} ) | \psi_{E} \right.$，其中 $x$ 和 $y$ 分别由 $i$ 和 $j$ 直接确定：

而 $y$ 与 $j$ 的关系相同。通过利用关系式 [49]

甚至可以在不记录切比雪夫演化过程中任何状态的情况下找到矩阵元。相反，我们只需在适当的时候，即当 $k = k_{m} - 2 , k = k_{m} - 1 , k = k_{m}$ 和 $k = k_{m} + 1$ 时，记录两个数 $\left. \psi_{E} | \tau_{k} ( \mathcal{G} ) | \psi_{E} \right.$ 和 $\langle \psi_{E} | \mathcal{H} T_{k} ( \mathcal{G} ) | \psi_{E} \rangle$。

最后，我们得到矩阵元的显式表达式

还需注意，由于需要 $T_{x + y}$，且 $x$ 和 $y$ 都可能达到最大值 $k_n$，在此模式下截断阶数 $K$ 加倍为 $K = 2 k_n$。当采用块技术时，情况变得稍复杂一些。此时需要记录交叉项（如 $2 \times 2$ 矩阵）而非单个数值。在此情形下，$i, j$ 与 $x, y$ 的关系保持不变。假设我们有两条不同的演化态 $| \psi_1 \rangle$ 和 $| \psi_2 \rangle$，则需要同时（相同 $k$ 下）记录以下三个数值：

最终版本的矩阵 $S$ 由四个矩阵组成：

其中 $S_{i j}^{a b} = \left. \psi_{a} \right| T_{x} ( \mathcal{G} ) T_{y} ( \mathcal{G} ) \left| \psi_{b} \right.$

对于广义本征值问题式(11)，首先将厄米矩阵 $S$ 对角化为：

其中 $V$ 是 $S$ 的本征向量矩阵，$V V^{\dagger} = I$，$\Lambda_{s}$ 是相应的本征值矩阵。由于 $S$ 通常是奇异的，我们通过剔除本征值的绝对值低于截断值 $\varepsilon = 10^{-12}$ 的列来压缩 $(2n+1) \times (2n+1)$ 的矩阵 $V$。记保留的本征值个数为 $m$，则压缩后的本征向量矩阵 $\widetilde{V}$ 的阶数为 $(2n+1) \times m$，且

下一步是形成压缩后的哈密顿矩阵 $\widetilde{H}$。由于

记变换矩阵 $U = \widetilde{V} \widetilde{\Lambda}_{s}^{- \frac 12}$，则压缩后的 $m \times m$ 哈密顿矩阵为

然后直接对角化阶数 $m \times m$ 的厄米矩阵 $\widetilde{H}$（$m$ 的取值范围为 $10^3$ 到 $10^4$）：

其中 $\widetilde{Y}$ 是 $\widetilde{H}$ 的本征向量矩阵，$\widetilde{\Lambda}$ 是由原始哈密顿量 $H$ 中落在区间 $[-a, a]$ 内的那些目标本征值组成的对角矩阵。投影哈密顿量 $H$ 的本征态可通过基本矩阵代数获得：

记式(C.3)为一个 $2^N \times (2n+1)$ 的矩阵 $A$，其中 $| \chi_i \rangle$ 是第 $i$ 列，则原始哈密顿量 $H$ 中落在 $[-a, a]$ 内的本征态为

我们最终得到特征值 $\widetilde{\Lambda}$ 和本征态 $\Phi$。

D 与混合方法的比较#

由于Lanczos方法广泛用于计算最低（或最高）特征值，而指数半圆滤波器提供了一种将中心特征值转换为最高特征值的方法，因此很自然会想到将两者结合起来高效计算中心特征值，这比DACP方法简单得多。但有一些理由倾向于采用DACP方法中使用的复杂技术，而不是标准的Krylov空间构造，特别是在需要大规模计算时。

首先，已知Lanczos方法在实际条件下数值稳定性不佳。具体来说，尽管Lanczos算法在理论上表现出完美特性，但在实践中它会失去许多设计特征，例如全局正交性以及Lanczos递归状态间的线性独立性[35,50,51]。这些缺陷有效地限制了可计算特征值的数量。此外，由于误差累积，广义特征值问题的出现（用于处理非正交基态）似乎是不可避免的[52,53]。相比之下，切比雪夫递归在理想和实际计算中都拥有许多有趣的共同性质[50]。特别地，它在 $x \in [ -1, 1 ]$ 上准确且稳定，允许在切比雪夫阶域中进行数万步的传播而不会积累显著误差[33]。

其次，考虑到重正交化步骤，Lanczos算法的效率会下降。重正交化步骤在Lanczos方法和Arnoldi方法中都是必要部分。当需要大量（例如数千个）特征值时，总成本实际上由重正交化主导[28,29,54]。假设所需特征值总数为R，希尔伯特空间维度为D，则Lanczos类方法的规模为 $D R^{2}$，因为每个生成的Ritz向量都需要与已有向量进行重正交化（另见文献[29,44]的讨论）。相比之下，我们在第3节中已表明 $T_{F} \propto D R^{-1}$，$T_{E} \propto D R^{0}$，$T_{D} \propto D^{0} R^{3}$，其中每一项都优于 $D R^{2} \ ( R \ll D )$。尽管已经提出了部分重正交化方案，但它们会导致计算成本和内存流量增加[20,55]，并且在精度要求严格时不能保证成功[55]。此外，DACP方法和Lanczos算法在空间复杂度上不同，前者为 $\operatorname* {m a x} ( D R^{0} , R^{2} )$，而后者一旦需要重正交化就要求 $D R$[44]。因此，DACP方法在时间和空间复杂度上都优于Lanczos算法。

最后，忽略重正交化的需求，我们发现两种方法在时间和空间复杂度上具有可比性。回忆一下，切比雪夫演化时间作为DACP方法中的主导项，变为 $T_{E} \simeq \pi R E_{\operatorname* {m a x}} \tau / a$，其中 τ 是矩阵-向量乘积的时间。然后我们讨论混合版本（Lanczos方法与指数半圆滤波器结合）的时间复杂度。如文献22的补充材料图1以及文献28的讨论所示，Lanczos步数 m 至少应为所需特征值数量R的两倍，即 $m \geq 2R$。由于需要2R次迭代，且每次迭代需要应用一次滤波器（$T_{F} \simeq 24 E_{\operatorname* {m a x}} \tau / a$），计算时间则为 $2R \cdot 24 E_{\mathrm{max}} \tau / a = 48 R E_{\mathrm{max}} \tau / a$。因此，DACP方法与Lanczos结合指数半圆滤波器具有相同的时间复杂度。

图8：DACP方法第二步中不同基的相对误差分布 $\eta ( E )$ 比较，包括 $\{T_{k} ( \mathcal{H} ) \}$（红色方块）、$\{\mathcal{H}^{k} \}$（蓝色方块）以及经过Schmidt正交化的 $\{\mathcal{H}^{k} \}$（蓝色叉号）。

此外，我们还使用DACP模拟（由于已知精确特征值，计算DACP三步产生的误差）测试了式(10)中的基组，其他条件保持不变。结果如图8所示。这样的基组在条件 $\eta ( E ) \leq 10^{-6}$ 下不会收敛到任何特征值。与DACP方法相比，基组 $\{\mathcal{H}^{k} \}$ 收敛最好的特征值位于区间的两端。总之，DACP方法比Lanczos方法结合指数半圆滤波器更稳定、更高效。

E 误差分析#

在本节中，我们推导相对误差分布的形状。我们使用两种不同方法得到误差结果。第一种方法大致遵循DACP方法的流程，只是用直接计算结果替代了那些计算困难的步骤。第二种方法是在若干近似下的分析推导。

我们首先阐述两种方法共有的假设。

1. 我们从初始态 $\begin{array} {r} {\left| \psi \right. = \beta \sum_{i} {c_{i} \left| E_{i} \right.}} \end{array}$ 出发，其中 $\beta$ 是归一化常数，$c_{i}$ 是实值概率幅，$\left| E_{i} \right.$ 是能量为 $E_{i}$ 的本征态，并且 $- a \leq E_{i} \leq a$。即，当 $| E_{i} | > a$ 时 $c_{i} = 0$。

2. 精确能级已知且等间距，即 $E_{i + 1} - E_{i} = E_{i} - E_{i - 1}$。

3. 最大抽象能量 $E_{\mathrm{max}} = 1$。

4. 我们按照切比雪夫演化和DACP方法的子空间对角化计算特征值。

在这些设定下，我们采用DACP方法的模拟程序（表示部分解析结果）来高效生成接近DACP方法的误差。具体来说，我们不执行DACP方法中的K阶切比雪夫演化，而是显式写出式(9)中态的表达，并构建式(12)中相应的子空间矩阵。结果与DACP几乎完全一致，除了式(B.2)两边存在微小差异。然而，这种模拟捕捉了切比雪夫演化和子空间对角化过程中产生的主要误差，而无需执行困难的计算任务（假设已知精确特征值）。此外，它也很方便测试由式(10)构建的基态的有效性。

接下来，我们介绍分析推导方法。模拟方法和分析方法之间的区别本质上在于DACP的第三步。前者执行子空间矩阵的完全对角化。为简化而抓住本质，这里我们只考虑两个不同的特征值 $\{E_{1} = 0 , E_{2} \}$ ^}，结合两个构造态 $\{| \psi_{1} \rangle , | \psi_{2} \rangle = \sin ( \pi \mathcal{H} ) | \psi_{1} \rangle \}$，其中

$\varepsilon = c_{2} / c_{1} \ll 1_{:}$，且 $\beta_{1} \simeq 1$。在式$(5)$的分布情况下，$c_{1} \simeq \exp ( 24 )$ $\varepsilon \ll 1$ 意味着 $E_{2}$ 远离 $E_{1} = 0$。显式表达 $| \psi_{2} \rangle$，我们得到

数值舍入误差（记为δ）在归一化过程中引入，即像 $\varepsilon^{2}$ 这样的小数被修改为 $\varepsilon^{2} + \delta$。具体来说，重叠矩阵S的表达式为

假设 $10^{- 12} \le \varepsilon^{2} \ll 1$，对式(E.1)进行对角化得到

相应地，式(C.11)中矩阵H的特征值为 $\left\{E_{1} , E_{2} \varepsilon^{2} \sin^{2} ( \pi E_{2} ) \right\}$，然后乘以 $\Lambda_{s}^{- \frac{1} {2}}$ 两次（回忆 $\widetilde{H} = \widetilde{\Lambda}_{s}^{- \frac 12} \widetilde{V}^{\dagger} H \widetilde{V} \widetilde{\Lambda}_{s}^{- \frac 12} )$ ）。受舍入误差影响，最小数 $E_{2} \varepsilon^{2} \sin^{2} ( \pi E_{2} )$ 变为 $E_{2} \varepsilon^{2} \sin^{2} ( \pi E_{2} ) + \delta $，产生相对误差

推广式(E.2)并使用 $c_{1}^{2} \simeq \exp ( 48 )$，最终得到

在推导式(E.3)的右端时，我们使用了双精度数有16位有效数字以及我们图9(a)模拟中的 $\beta^{2} \simeq 4 \times 10^{- 23}$，其中 $\beta$ 是满足 $\begin{array} {r} {\beta^{2} \sum_{i} c_{i}^{2} = 1} \end{array}$ 的归一化常数。当概率幅分布变化时，$\beta$ 相应改变，给出 $\eta ( E )$ 的不同表达式。

另一种方式，我们也可以考虑一个更复杂的情况，即构造的态集为 $\{| \psi_{1} \rangle , | \psi_{2} \rangle = \cos ( \pi \mathcal{H} ) | \psi_{1} \rangle \}$ ^}。于是我们有

矩阵

给出对角矩阵

类似的论证同样导出式(E.3)。

在图9中，我们比较了上述两种方法的相对误差分布。具体地，还展示了DACP方法的结果（同图3(b)）。前两种方法的初态直接设为 $\begin{array} {r} {| \psi \rangle = \beta \sum_{i} c_{i} | \phi_{i} \rangle} \end{array}$，其中

在DACP方法中，经过指数半圆过滤后也制备了类似的态。显然，$c_{i}$的分布直接与相对误差 $\eta ( E )$ 相关，特别是对于远离 $E = 0$ 的特征值。因此，为减少误差并发现更多特征值，切比雪夫演化的最佳初态 $| \psi_{E} \rangle$ 应为窗函数

其中R是 $[ - a , a ]$ 内特征空间的维度。

参考文献#

[1] F. Haake, Quantum signatures of chaos, In Quantum coherence in mesoscopic systems, Springer, Boston, MA (1991), doi:10.1007 978-1-4899-3698-1_38.

[2] T. A. Brody, J. Flores, J. B. French, P. A. Mello, A. Pandey and S. S. M. Wong, Randommatrix physics: spectrum and strength fluctuations, Rev. Mod. Phys. 53, 385 (1981), doi:10.1103_/RevModPhys.53.385.

[3] O. Bohigas, M. J. Giannoni and C. Schmit, Characterization of chaotic quantum spectra and universality of level fluctuation laws, Phys. Rev. Lett. 52, 1 (1984), doi:10.1103_/PhysRevLett.52.1.

[4] M. L. Mehta, Random matrices, Academic Press, Cambridge, Massachusetts, US (1991).

[5] A. Relaño, J. M. G. Gómez, R. A. Molina, J. Retamosa and E. Faleiro, Quantum chaos and 1/ f noise, Phys. Rev. Lett. 89, 244102 (2002), doi:10.1103_/PhysRevLett.89.244102.

[6] D. M. Basko, I. L. Aleiner and B. L. Altshuler, Metal-insulator transition in a weakly interacting many-electron system with localized single-particle states, Ann. Phys. 321, 1126 (2006), doi:10.1016_/j.aop.2005.11.014.

7 R. Nandkishore and D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Annu. Rev. Condens. Matter Phys. 6, 15 (2015), doi:10.1146 annurev-conmatphys-031214-014726.

[8] F. Pietracaprina, N. Macé, D. J. Luitz and F. Alet, Shift-invert diagonalization of large many-body localizing spin chains, SciPost Phys. 5, 045 (2018), doi:10.21468_/SciPostPhys.5.5.045.

[9] T. Ericsson and A. Ruhe, The spectral transformation Lánczos method for the numerical solution of large sparse generalized symmetric eigenvalue problems, Math. Comput. 35, 1251 (1980), doi:10.1090_/S0025-5718-1980-0583502-2.

[10] R. E. Wyatt, Matrix spectroscopy: Computation of interior eigenstates of large matrices using layered iteration, Phys. Rev. E 51, 3643 (1995), doi:10.1103_/PhysRevE.51.3643.

[11] T. J. Minehardt, J. D. Adcock and R. E. Wyatt, Enhanced matrix spectroscopy: The preconditioned Green-function block Lánczos algorithm, Phys. Rev. E 56, 4837 (1997), doi:10.1103 PhysRevE.56.4837.

[12] C. Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, J. Res. Natl. Bur. Stan. 45, 255 (1950), doi:10.6028_/jres.045.026.

[13] V. A. Mandelshtam and H. S. Taylor, A simple recursion polynomial expansion of the Green’s function with absorbing boundary conditions. Application to the reactive scattering, J. Chem. Phys. 103, 2903 (1995), doi:10.1063_/1.470477.

14 Z. Bai, J. Demmel, J. Dongarra, A. Ruhe and H. van der Vorst, Templates for the solution of algebraic eigenvalue problems: A practical guide, SIAM, Philadelphia (2000), doi:10.1137 1.9780898719581.

[15] D. J. Luitz, N. Laflorencie and F. Alet, Many-body localization edge in the random-field Heisenberg chain, Phys. Rev. B 91, 081103 (2015), doi:10.1103 PhysRevB.91.081103.

16 P. Sierant and J. Zakrzewski, Model of level statistics for disordered interacting quantum many-body systems, Phys. Rev. B 101, 104201 (2020), doi:10.1103_/PhysRevB.101.104201.

[17] M. Hopjan and F. Heidrich-Meisner, Many-body localization from a one-particle perspective in the disordered one-dimensional Bose-Hubbard model, Phys. Rev. A 101, 063617 (2020), doi:10.1103_/PhysRevA.101.063617.

[18] J. Šuntajs, J. Bonˇca, T. Prosen and L. Vidmar, Ergodicity breaking transition in finite disordered spin chains, Phys. Rev. B 102, 064207 (2020), doi:10.1103_/PhysRevB.102.064207.

[19] R. Li, Y. Xi, E. Vecharynski, C. Yang and Y. Saad, A thick-restart Lanczos algorithm with polynomial filtering for hermitian eigenvalue problems, SIAM J. Sci. Comput. 38, A2512 (2016), doi:10.1137_/15M1054493.

[20] H.-r. Fang and Y. Saad, A filtered Lanczos procedure for extreme and interior eigenvalue problems, SIAM J. Sci. Comput. 34, A2220 (2012), doi:10.1137 110836535.

[21] R. Li, Y. Xi, L. Erlandson and Y. Saad, The eigenvalues slicing library (EVSL): Algorithms, implementation, and software, SIAM J. Sci. Comput. 41, C393 (2019), doi:10.1137_/18M1170935.

[22] P. Sierant, M. Lewenstein and J. Zakrzewski, Polynomially filtered exact diagonalization approach to many-body localization, Phys. Rev. Lett. 125, 156601 (2020), doi:10.1103_/PhysRevLett.125.156601.

[23] H. Guan and W. Zhang, Delta-Davidson method for interior eigenproblem in many-spin systems, Chin. Phys. B 30, 030205 (2021), doi:10.1088_/1674-1056_/abd74a.

[24] A. Pieper, M. Kreutzer, A. Alvermann, M. Galgon, H. Fehske, G. Hager, B. Lang and G. Wellein, High-performance implementation of Chebyshev filter diagonalization for interior eigenvalue computations, J. Comput. Phys. 325, 226 (2016), doi:10.1016_/j.jcp.2016.08.027.

[25] P. T. P. Tang and E. Polizzi, FEAST as a subspace iteration eigensolver accelerated by approximate spectral projection, SIAM J. Matrix Anal. Appl. 35, 354 (2014), doi:10.1137_/13090866X.

26 B. I. Shklovskii, B. Shapiro, B. R. Sears, P. Lambrianides and H. B. Shore, Statistics of spectra of disordered systems near the metal-insulator transition, Phys. Rev. B 47, 11487 (1993), doi:10.1103 PhysRevB.47.11487.

[27] B. Georgeot and D. L. Shepelyansky, Integrability and quantum chaos in spin glass shards, Phys. Rev. Lett. 81, 5129 (1998), doi:10.1103_/PhysRevLett.81.5129.

[28] Y. Zhou and Y. Saad, A Chebyshev-Davidson algorithm for large symmetric eigenproblems, SIAM J. Matrix Anal. Appl. 29, 954 (2007), doi:10.1137_/050630404.

29 Y. Zhou, A block Chebyshev-Davidson method with inner-outer restart for large eigenvalue problems, J. Comput. Phys. 229, 9188 (2010), doi:10.1016_/j.jcp.2010.08.032.

[30] D. Poulin and P. Wocjan, Preparing ground states of quantum many-body systems on a quantum computer, Phys. Rev. Lett. 102, 130503 (2009), doi:10.1103 PhysRevLett.102.130503.

[31] T. J. Rivlin, An introduction to the approximation of functions, Blaisdell Publishing Company, Waltham, Massachusetts (1969), doi:https:_//doi.org_/10.1137_/1012069.

32 R. Li, Y. Xi, E. Vecharynski, C. Yang and Y. Saad, A thick-restart Lanczos algorithm with polynomial filtering for Hermitian eigenvalue problems, SIAM J. Sci. Comput. 38, A2512 (2016), doi:10.1137 15M1054493.

[33] R. Chen and H. Guo, The Chebyshev propagator for quantum systems, Comput. Phys. Commun. 119, 19 (1999), doi:10.1016 S0010-4655(98)00179-9.

[34] E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. D. Croz, A. Greenbaum, S. Hammarling, A. McKenney, S. Ostrouchov and D. Sorensen, LAPLACK user’s guide—release 2.0, SIAM, Philadelphia (1994), doi:10.1137_/1.9780898719604.

[35] Y. Saad, Numerical methods for large eigenvalue problems: Revised edition, SIAM, Philadelphia, (2011), doi:10.1137_/1.9781611970739.

[36] M. A. Nielsen and I. L. Chuang, Quantum computation and quantum information, Cambridge University Press, Cambridge, UK, ISBN 9781107002173 (2011), doi:10.1017 CBO9780511976667.

[37] N. N. Bogolubov and N. N. Bogolubov Jr., Introduction to quantum statistical mechanics, World Scientific, Singapore, ISBN 9789814295192 (2009), doi:10.1142_/7623.

[38] S. Sachdev, Quantum phase transitions, Cambridge University Press, Cambridge, UK, ISBN 9780511973765 (2011), doi:10.1017_/CBO9780511973765.

[39] A. Dutta, G. Aeppli, B. K. Chakrabarti, U. Divakaran, T. F. Rosenbaum and D. Sen, Quantum phase transitions in transverse field spin models: From statistical physics to quantum information, Cambridge University Press, Cambridge, UK, ISBN 9781107706057 (2015), doi:10.1017_/CBO9781107706057.

[40] D. S. Fisher, Critical behavior of random transverse-field Ising spin chains, Phys. Rev. B 51, 6411 (1995), doi:10.1103 PhysRevB.51.6411.

[41] U. Schollwöck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005), doi:10.1103_/RevModPhys.77.259.

[42] Y. Zhou and R.-C. Li, Bounding the spectrum of large Hermitian matrices, Linear Algebra Appl. 435, 480 (2011), doi:10.1016_/j.laa.2010.06.034.

[43] A. Hams and H. De Raedt, Fast algorithm for finding the eigenvalue distribution of very large matrices, Phys. Rev. E 62, 4365 (2000), doi:10.1103_/PhysRevE.62.4365.

[44] R. B. Lehoucq, D. C. Sorensen and C. Yang, ARPACK user’s guide: Solution of large scale eigenvalue problems with implicitly restarted Arnoldi methods, SIAM, Philadelphia (1998), doi:10.1137 1.9780898719628.

45 E. Polizzi, FEAST eigenvalue solver v4.0 user guide, arXiv:2002.04807.

[46] E. Polizzi, Density-matrix-based algorithm for solving eigenvalue problems, Phys. Rev. B 79, 115112 (2009), doi:10.1103 PhysRevB.79.115112.

[47] D. A. Abanin, E. Altman, I. Bloch and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019), doi:10.1103_/RevModPhys.91.021001.

[48] D. Mondal, S. Sinha and S. Sinha, Chaos and quantum scars in a coupled top model, Phys. Rev. E 102, 020101 (2020), doi:10.1103 PhysRevE.102.020101.

[49] J. C. Mason and D. C. Handscomb, Chebyshev polynomials, Chapman and Hall Publisher, Boca Raton, ISBN 9780849303555 (2003).

[50] R. Chen and H. Guo, On the convergence scaling laws of Lanczos and Chebyshev recursion methods, J. Chem. Phys. 119, 5762 (2003), doi:10.1063_/1.1598958.

[51] J. K. Cullum and R. A. Willoughby, Lanczos algorithms for large symmetric eigenvalue computations: Vol. I: Theory, SIAM, Philadelphia (1985), doi:10.1137_/1.9780898719192.

[52] B. N. Parlett, The symmetric eigenvalue problem, SIAM, Philadelphia (1987), doi:10.1137 1.9781611971163.

[53] M. R. Wall and D. Neuhauser, Extraction, through filter-diagonalization, of general quantum eigenvalues or classical normal mode frequencies from a small number of residues or a short-time segment of a signal. I. Theory and application to a quantum-dynamics model, J. Chem. Phys. 102, 8011 (1995), doi:10.1063_/1.468999.

[54] Y. Saad, Y. Zhou, C. Bekas, M. L. Tiago and J. Chelikowsky, Diagonalization methods in PARSEC, Phys. Status Solidi B 243, 2188 (2006), doi:10.1002_/pssb.200666816.

[55] C. Bekas, E. Kokiopoulou and Y. Saad, Computation of large invariant subspaces using polynomial filtered Lanczos iterations with applications in density functional theory, SIAM J. Matrix Anal. Appl. 30, 397 (2008), doi:10.1137 060675435.

阅读笔记#

一句话概括#

本文提出DACP（切比雪夫多项式双重应用）方法，通过两次利用切比雪夫多项式——首先构造指数半圆滤波器将谱中部本征值放大到谱端，再通过切比雪夫演化生成线性无关基态——将高维本征值问题转化为低维子空间对角化问题，在伊辛自旋链上相比移位求逆方法提速达30倍，内存降低100倍。

核心论证链#

1. 指数半圆滤波器（切比雪夫第一次应用）：构建算符 $\mathcal{F} = (\mathcal{H}^2 - E_c)/E_0$ 将中心能谱 $[-a,a]$ 映射到 $\mathcal{F}$ 的 $[-1,1]$ 之外，利用切比雪夫多项式在 $x\notin[-1,1]$ 区间指数增长的特性（式(4)），以 $K = 12E_{\mathrm{max}}/a$ 阶截断将谱中部本征态分量指数放大、外部分量抑制，生成局域在子空间 $\mathbb{L}$ 的目标态 $|\psi_E\rangle$。

2. 切比雪夫演化（切比雪夫第二次应用）：对目标态 $|\psi_E\rangle$ 施加算符 $\mathcal{G} = \mathcal{H}/E_{\mathrm{max}}$（保证谱在 $[-1,1]$ 内），利用 $T_k(\mathcal{G}) = \cos(k\arccos(\mathcal{G}))$ 的振荡特性，按 $k_m = \lfloor m\pi/a_r \rfloor$ 步长记录 $T_{k}(\mathcal{G})|\psi_E\rangle$，生成一组 $\sin(\cdot)/\cos(\cdot)$ 形式的线性无关基态（式(9)），其斜率 $\sim n\pi/a_r$ 远优于幂级数基的 $1/a_r$。

3. 子空间对角化：利用切比雪夫乘积恒等式（式(C.4)）直接计算投影哈密顿量 $H_{ij} = \langle\Psi_i|\mathcal{H}|\Psi_j\rangle$ 和重叠矩阵 $S_{ij}$，无需显式存储基态向量；再通过奇异值分解处理 $S$ 的奇异性（截断条件 $\varepsilon=10^{-12}$），最后对压缩后 $m\times m$ 的 $\widetilde{H}$ 直接对角化得到目标本征值。

4. 计算复杂度分析：滤波时间 $T_F \propto (E_{\mathrm{max}}/a)\tau \propto R^{-1}$，演化时间 $T_E \propto R^0$（常数），对角化时间 $T_D \propto R^3$；当所需本征值数 $R$ 在 $10^2$ 到 $10^3$ 量级时，$T_E$ 主导，总时间与 $R$ 几乎无关——这与其他迭代方法（计算时间 $\propto R$）形成本质区别。

5. 数值验证与对比：对伊辛模型（可积，泊松统计）和自旋玻璃碎片（混沌，Wigner-Dyson GOE统计）验证了DACP的精度（$\eta<10^{-6}$）；通过与移位求逆、EVSL、FEAST三种方法的扩展行为比较，确认DACP在 $N=13$–$17$ 自旋系统中快20倍的性能优势，且内存恒为5.6 GB（移位求逆在 $N=17$ 需573 GB）。

实验参数详解#

批判性思考#

1. 原子性假设的有效性：指数半圆滤波器在 $a$ 很小时才退化到精确的半圆形式 式(5)的近似 $T_k(\mathcal{F}) \simeq \exp(2k\sqrt{a^2-\mathcal{H}^2}/E_{\mathrm{max}})$ 依赖于泰勒展开 $\cosh^{-1}(1+\varepsilon) \simeq 2\sqrt{a^2-E_i^2}/E_{\mathrm{max}}$，这要求 $a^2 \ll E_{\mathrm{max}}^2$。当系统能谱不对称或 $a$ 增大时，滤波器形状偏离半圆形，可能导致目标区间外分量被放大的风险增加，论文未量化这一偏差对更大 $a$ 值的实际影响。

2. 近简并本征值问题的本质局限未被充分解决 论文承认能级间距小至 $10^{-7}\Gamma$ 时切比雪夫演化难以区分（需要长两个数量级的演化时间），采取的阻塞滤波技术（block size=5）只是暴力扩展子空间维度而非从根本上解决区分问题。对于高度简并系统（如具有对称性的自旋系统），简并度可能超过块大小，这时误差尖峰（如图3(a)中观测到的）会显著增多。

3. 与 Lanczos + 滤波器混合方法的比较可能存在时序不对等 附录D中声称混合方法的计算时间为 $2R \cdot 24E_{\mathrm{max}}\tau/a$，但这一估计假设每次Lanczos迭代都需要重新应用完整的指数半圆滤波器（$K=12E_{\mathrm{max}}/a$ 次矩阵-向量乘积）。实际上，Lanczos迭代过程中预计算的滤波器系数可以复用，且实际应用中的部分重正交化策略可能使复杂度 $DR^2$ 中的系数远小于论文的估计。

4. 自旋玻璃碎片模型中DACP的时间优势显著减弱 伊辛模型DACP比移位求逆快30倍，但自旋玻璃碎片仅快8倍。论文归因于泡利算符数量（6.2倍差异），但这也暗示DACP对哈密顿量的稀疏性/算子结构高度敏感。对于全连通模型（所有自旋对耦合），每个矩阵-向量乘积的成本从伊辛模型的 $O(N)$ 增加到 $O(N^2)$，这会直接放大 $T_F$ 和 $T_E$。

局限性#

• 适用性局限于谱中点附近：DACP方法专门针对 $E=0$ 附近的谱中部设计，无法直接用于谱边缘或任意能量区间的本征值计算。
• 对称性假设：方法隐含假设 $E_{\mathrm{min}} = -E_{\mathrm{max}}$（能谱对称），对于能谱不对称的系统需要额外平移处理，论文未验证此情形下的性能。
• $a$ 参数依赖先验的态密度信息：参数 $a$ 需根据目标本征值数目和态密度 $\rho(E)$ 手动选取，若 $\rho(E)$ 估计不准，可能导致滤波后的子空间 $\mathbb{L}$ 维度不足或过度。
• 切比雪夫递推的数值稳定性在超大阶数下未验证：论文中滤波器阶数 $K = 12E_{\mathrm{max}}/a$，当 $a$ 很小时 $K$ 可达数千甚至上万，虽然附录A声称切比雪夫递推稳定，但实际舍入误差在极端阶数下的累积效应未量化。
• 子空间对角化时间 $T_D \propto R^3$ 在高 $R$ 时成为瓶颈：当所需本征值数超过约5,000时，$T_D$ 将主导总计算时间，此时DACP不再具有"计算时间与 $R$ 无关"的优势。
• 误差分析限于简化模型：附录E的误差分析假设等间距能级和双态近似，无法完全刻画实际多体系统复杂谱分布下的误差行为。

关键公式速查#

• $$\mathcal{F} = (\mathcal{H}^2 - E_c)/E_0$$— 平移缩放算符，将谱中部映射到切比雪夫指数增长区，式(2)
• $$T_k(\mathcal{F}) \simeq \exp\!\left(\frac{2k}{E_{\mathrm{max}}}\sqrt{a^2 - \mathcal{H}^2}\right)$$— 指数半圆滤波器近似，式(5)
• $$\mathcal{G} = \mathcal{H}/E_{\mathrm{max}}$$— 归一化算符，将谱映射到 $[-1,1]$ 用于切比雪夫振荡，式(7)
• $$T_k(\mathcal{G}) = \cos(k\arccos(\mathcal{G}))$$— 切比雪夫多项式的三角函数形式，用于振荡演化，式(A.2)
• $$T_{k+1}(x) = 2x\,T_k(x) - T_{k-1}(x)$$— 三项递推关系，高效计算切比雪夫多项式，式(A.3)
• $$\widetilde{H} = U^{\dagger} H U,\quad U = \widetilde{V}\widetilde{\Lambda}_s^{-1/2}$$— 广义特征值问题的重正化变换，式(C.10)–(C.11)

术语对照#