📚 课程目录第4章精确对角化与稀疏矩阵方法

4. 精确对角化与稀疏矩阵方法

精确对角化（Exact Diagonalization, ED）是量子多体动力学中最直接、最可控的数值方法之一。它的核心思想并不复杂：把哈密顿量在某个有限维希尔伯特空间中写成矩阵，然后求本征值、本征矢，或者直接计算时间演化。真正决定ED能否用于实际问题的，不是“对角化”这一步，而是基矢如何构造、矩阵如何存储、矩阵-向量乘法如何实现，以及时间演化如何避免显式矩阵指数。本节围绕这四个关键环节展开。

一、基矢构造：从占据数到整数编码

以自旋-1/2链或费米子格点模型为例，设格点数为 L，每个格点有两个局域态。若固定粒子数或总磁化，则只需保留满足守恒量的基矢。最常用的编码方式是用一个整数表示一个占据组态：第 i 位为1表示该格点被占据，为0表示空。例如 L=4、两个粒子时，组态 |1010⟩ 对应整数 0b1010 = 10。

构造基矢的典型流程是：先枚举所有 2^L 个整数，筛选出满足粒子数守恒的组态，再按升序或降序排列，并建立“整数→索引”的哈希映射。这个映射是后续所有矩阵元计算的基础。对于 L=20、半填充的情形，基矢维度约为 C(20,10)=184756，完全在单机内存可承受范围内；而 L=30 时维度约1.55亿，必须依赖稀疏存储和并行计算。

基矢排序方式会影响稀疏矩阵的非零元分布，但对Lanczos收敛速度影响很小；真正影响效率的是矩阵-向量乘法的实现是否利用了位运算和局域性。

二、稀疏矩阵存储与矩阵-向量乘法

多体哈密顿量通常是局域的：每一项只涉及少数几个格点。因此，哈密顿量矩阵是高度稀疏的。以最近邻海森堡模型为例，每个基矢只与 O(L) 个其他基矢相连，非零元总数约为 O(L·D)，其中 D 是基矢维度。存储时常用CSR（Compressed Sparse Row）格式：三个数组分别保存非零元值、列索引和行起始位置。

更高效的做法是不显式存储矩阵，而是实现一个函数 apply_H(v)，输入向量 v，输出 H·v。对每个基矢 |s⟩，遍历哈密顿量中的所有局域项，计算它作用后的新组态 |s'⟩，通过哈希映射找到 s' 的索引，累加贡献。这样内存只需存储基矢表和哈密顿量参数，而非 O(L·D) 个非零元。对于 D=10^6、L=20 的典型问题，显式CSR约需数百MB，而on-the-fly乘法只需几十MB。

// 伪代码：on-the-fly 矩阵-向量乘法
for i in 0..D-1:
    s = basis[i]
    for each bond (j,k) in lattice:
        // 对角项
        out[i] += Jz * sz(s,j) * sz(s,k) * v[i]
        // 交换项
        if s[j] != s[k]:
            s2 = flip(s, j, k)
            idx = hash[s2]
            out[idx] += 0.5 * Jxy * v[i]

三、Lanczos算法：求基态与低激发态

Lanczos算法是求解大型稀疏厄米矩阵极端本征值的最有效方法之一。它通过构造Krylov子空间 K_m = span{v, Hv, H²v, ..., H^{m-1}v}，将原矩阵投影为一个 m×m 的三对角矩阵 T_m，然后对 T_m 做稠密对角化。当 m 远小于 D 时，计算量大幅降低。

算法迭代过程如下：给定归一化初始向量 v_0，令 β_0=0，v_{-1}=0。对 n=0,1,2,...：

w = H·v_n - β_n·v_{n-1}
α_n = ⟨v_n, w⟩
w = w - α_n·v_n
β_{n+1} = ||w||
v_{n+1} = w / β_{n+1}

每一步只需一次矩阵-向量乘法、两次向量内积和若干向量更新。三对角矩阵 T_m 的对角元为 α_n，次对角元为 β_n。对 T_m 做本征分解，最小本征值即为基态能量的近似，对应的Ritz向量可通过 V_m y 还原。

实际使用中需要注意：Lanczos在有限精度下会失去正交性，出现“幽灵本征值”。解决办法是定期重正交化，或采用完全重正交化的Lanczos。对于基态计算，通常 m=100~300 即可收敛到机器精度；对于低激发态，可用块Lanczos或位移-反迭代。

四、时间演化算符的Krylov子空间实现

计算 |ψ(t)⟩ = e^{-iHt}|ψ(0)⟩ 时，显式计算矩阵指数 e^{-iHt} 对大型稀疏矩阵不可行。Krylov子空间方法提供了高效近似：在每一步时间演化中，以当前态 |ψ(t)⟩ 为初始向量，构造小规模Krylov子空间 K_m，在该子空间内计算 e^{-iHt} 的作用。

具体地，用Arnoldi过程构造正交基 V_m = [v_1, v_2, ..., v_m] 和上Hessenberg矩阵 H_m，满足 H V_m ≈ V_m H_m。则

|ψ(t+Δt)⟩ ≈ V_m · exp(-i H_m Δt) · e_1 · ||ψ(t)||

其中 e_1 是第一个单位向量。由于 m 通常取20~50，exp(-i H_m Δt) 可直接用稠密矩阵指数计算。误差随 m 指数收敛，且与 Δt 和哈密顿量谱宽有关。实际实现中，Δt 不宜过大，否则需要更大的 m；通常取 Δt·||H|| ~ 0.1~1 量级。

Krylov时间演化与Lanczos求基态共享同一套矩阵-向量乘法内核，区别在于：Lanczos利用三对角结构求极端本征值，而时间演化需要完整的Arnoldi正交化来保证 exp(-iHΔt) 的精度。

五、方法对比与选择依据

方法	适用目标	内存需求	典型维度	关键限制
稠密对角化	全谱、小系统	O(D²)	D ≤ 10⁴	内存爆炸
Lanczos	基态、低激发态	O(D·m)	D ≤ 10⁸	重正交化开销
Krylov时间演化	实时动力学	O(D·m)	D ≤ 10⁷	步长与m的平衡

选择方法时，先明确目标：若只需基态能量和关联函数，Lanczos是首选；若要模拟量子淬火后的实时演化，Krylov时间演化更合适；若系统维度很小且需要全谱信息，稠密对角化仍然简单可靠。无论哪种方法，基矢构造和稀疏矩阵-向量乘法都是共同的底层基础，值得优先优化。

最后给出一个可操作的参数选择参考：对于 L=16 的海森堡链，基矢维度 D=12870，Lanczos取 m=100 即可收敛到 10^{-12} 精度；Krylov时间演化取 m=30、Δt=0.05/J，单步误差约 10^{-8}。若系统增大到 L=24，D≈2.7×10^6，需使用并行矩阵-向量乘法，并将Lanczos的 m 提高到200以上以应对更密集的能谱。这些数值属于典型经验范围，实际使用时建议先用小系统验证收敛性，再逐步放大。