线性代数#

Simone Montangero

本章将讨论计算物理学家最常面对的一个简单但普遍存在的问题：特征值问题的求解。在量子物理中，特征值问题通常源于求解描述单粒子或多粒子量子系统的给定哈密顿量。事实上，除了少数几类存在解析解的问题（从无限深势阱到贝特假设计算），量子系统基态和激发态的搜索必须通过数值方法进行。特征值问题的数值求解等同于在一个适当基组中对表示哈密顿量的矩阵进行对角化，无论原始系统是离散的（例如一组量子比特或原子能级）还是连续的（例如自由空间或陷阱中的原子、离子，或玻色-爱因斯坦凝聚体）。后一种情况需要通过空间离散化（见附录B）或二次量子化转化为离散系统。最后，数学、物理学和计算机科学中的大量问题都可以转化为特征值问题。因此，我们首先介绍该问题，然后按复杂度顺序，介绍全局和部分特征值求解器的前沿方法。正如我们将看到的，特征值求解策略基于线性代数运算，为完整起见，我们在此简要回顾这些运算，同时也为了设定符号：全书将使用标准量子力学符号，并假设变量为复数。本章及下一章（以及附录）中的快速概述基于许多优秀的书籍和综述，例如[120–129]。

2.1 线性方程组#

线性方程组可以用算子形式表示

其中 $\hat{O}$ 是线性算子，向量 x, y 存在于希尔伯特空间 $\mathcal{H} \equiv \mathbb{C}^{N}$ 中。给定 $\hat{O}$ 和 $| y \rangle$，可以利用线性方程组等价变换的性质来求解向量 $\vert x \rangle$，这导致了高斯消去法：(1) 将方程组的行乘以因子 $f_{i , j}$，(2) 将一行替换为该行自身与另一行逐项相加的结果，这两类操作不改变方程组的解。因此，任何形如式(2.1)的线性方程组都可以转化为一个等价的（具有相同解的）方程组，其中算子具有上三角形式

给定算子 $\hat{U}$，线性方程组的解可以通过如下回代过程轻松获得：

前面的行不仅描述了高斯消去法（Gaussian elimination），还提供了求解任意 N 个线性方程组的一般算法描述，可以直接编码实现。算法的基本属性是它们的算法复杂度，定义为解决问题所需操作数量与输入规模（此处为方程数 N）的函数关系（更多关于软件复杂度和良好软件开发导论的内容，请参见附录 B）。一个简单的估计得出操作总数随 $O ( N^{3} )$ 缩放：将式 (2.2) 中矩阵下三角部分的 $N ( N - 1 ) / 2$ 个元素逐个归零需要 (N ) 次操作，总计 $\bar{O} ( N^{3} )$。式 (2.3) 描述的回代（back substitution）运算阶数较低，需要 $O \big ( N^{2} \big )$ 次操作。

注意，求解线性方程组的算法也提供了计算矩阵逆的方法：实际上，根据定义，矩阵 A 的逆满足 $\hat{O}^{- 1} \hat{O} = \hat{O} \hat{O}^{- 1} = \mathbb{1}$，因此我们可以将它们识别为

其中 $| x^{j} \rangle$ 和 $| y^{j} \rangle$ 分别构成了两个矩阵的列向量。总之，求矩阵逆的条件就简化为一组线性方程组：

这些方程可以通过构建增广矩形矩阵（由运算符 $\hat{O}$ 及其旁边的单位阵组成）的初等方法来求解。高斯消去运算符 $\hat{O}$ 的非对角项，就得到了 $\hat{O}^{- 1}$ 的计算结果。

2.1.1 LU 分解#

上一节回顾的高斯消去法只是将矩阵化为三角矩阵的一类约化方法之一。另一种在量子系统基态性质数值计算中起关键作用的约化方法是LU 分解[122]。事实上，我们可以直接证明任何运算符都可以写成左三角和右三角运算符的乘积：

其中系数 $f_{i , j}$ 正是用于将运算符 $\hat{O}$ 高斯消去为式 (2.2) 中矩阵 $\hat{U}$ 的那些系数。

正如我们将在本章后面看到的那样，当我们需要求解许多仅在向量 y 上不同的线性方程组时，这种约化的实用性就变得非常明显。实际上，在这种情况下，原问题 ${\hat{O}} | x \rangle =$ $| y \rangle$ 可以改写为两个线性方程组：

该求解仅需 $O \big ( N^{2} \big )$ 次运算，因为只需进行回代，且两个算子均为三角算子。因此，对于具有相同 $\hat{O}$ 和不同 $| y \rangle$ 的多个线性方程组，其求解算法如下：

1. 对算子 $\hat{O}$ 进行LU分解（规模为 $O ( N^{3} )$）；
2. 执行两次回代（规模为 $O \big ( N^{2} \big )$）。

该方法不存在缺陷，因为第一步的规模与完整的高斯消去法相同。因此，它不会引入任何计算时间或内存开销，因为两个三角矩阵可以用与原始矩阵相同的存储量存储。此外，如果待求解的方程组数量 M（具有相同的 $\hat{O}$ 算子和不同的 $| y \rangle$）满足 $M \gg N$ ，则本方法更为有利，因为 $M \mathcal{O} \big ( N^{3} \big ) \gg 1 \cdot \mathcal{O} \big ( N^{3} \big ) + M \cdot \mathcal{O} \big ( N^{2} \big )$

2.2 本征值问题#

本节介绍求解作用于希尔伯特空间 $\mathcal{H} \equiv \mathbb{C}^{N}$ 上的算子 $\hat{O}$ 本征值的最常用方法。下文将专门讨论厄米算子（Hermitian operator）的情况，因为它们在量子力学问题中具有显著的重要地位。然而，以下方法可以轻松推广到更一般的情况。例如，在研究开放量子系统时，自然会遇到非厄米矩阵的处理。

我们旨在求解如下定义的本征值问题：

其中，下文中我们假定本征态是正交归一的 $\langle x_{j} | x_{k} \rangle = \delta_{j , k}$，且本征值（实数）按序排列 $| \lambda_{N} | > | \lambda_{N - 1} | > . . . | \lambda_{1} |$

2.2.1 幂法#

本征态集合 $| x_{j} \rangle$ 构成了希尔伯特空间 $\mathcal{H}$ 的一组完备基，因此任意（随机）态 $| y \rangle \in \mathcal{H}$ 可写为 $\begin{array} {r} {| y \rangle = \sum y_{j} | y_{j} \rangle} \end{array}$ ，其中 $y_{j} = \langle y_{j} | y \rangle$ 。因此，将算子 $\hat{O}$ 作用于初始随机态 $| y \rangle$ 的结果为：

也就是说，在对算子 $\hat{O}$ 进行 k 次作用后，得到的矢量态为：

将算子作用于初始态的最终结果是，根据本征值之间的比值改变相对于本征矢的投影，即将该态旋转到具有最大本征值的本征态方向：当 $k \to \infty$ 时，态的序列趋于：

相应地，比值 $| | y^{k} \rangle | / | | y^{k - 1} \rangle | \lambda_{N}^{k} / \lambda_{N}^{k - 1} = \lambda_{N}$ 。总之，可以通过多次将算子作用于初始随机矢量，并计算两个连续矢量范数之间的比值，来找到最显著的本征态和本征矢。我们稍后将展示如何推广此过程，以找到其余的本征谱和剩余的本征值。

该过程运行良好，易于实现，并且在大多数场景下都有效。然而，当系统具有近乎简并的本征值时（这构成了一类基本的物理过程，例如，正如我们将在第10章中看到的，在热力学极限下表现出临界行为的系统），其效率会显著降低。事实上，方程(2.11)中序列的收敛速度取决于比值 $\lambda_{N} / \lambda_{N - 1}$ ，如果该比值很小，实际上会阻止算法收敛。为了解决这个问题并加速算法收敛，引入了反幂法（inverse power method）。其主要思想源于方程(2.8)等价于：

因此，应用$\hat{O}^{- 1}$可以找到$\hat{O}$的最小本征值和本征矢。如前所述，收敛到最终结果的速度由比值$\lambda_{2} / \lambda_{1}$决定。显然，这两种方法除了针对谱的两端之外，没有太大区别。然而，重要的区别在于，后一种方法中，极限本征值位于收敛速度比的分母中：这在许多场景中允许我们大幅提高收敛速度。实际上，在许多情况下，可以通过某种微扰论证来估计算符的最低本征矢。另一个典型场景是，当寻找哈密顿量基态能量作为参数$g$的函数时，可以利用先前为$g - \Delta g$计算的结果。总之，可以重新定义本征值问题的偏移量（例如，为能量设置不同的零值）为

其中$q$是对最小本征值$\lambda_{1}$的估计。序列

的收敛速度现在由$( \lambda_{2} - q ) / ( \lambda_{1} - q )$给出，如果$\lambda_{1} \sim q$，则可以大幅提高。参数$q$通常被称为对角化子程序中的猜测值（guess value），提供好的猜测是区分标准结果和专业结果的关键。

细心的读者肯定已经注意到，在上述介绍中，有一个关键点没有被仔细处理：实际上，为了计算式(2.14)，需要算符$( \hat{O} - q \mathbb{1} )^{- 1}$。正如我们在上一节所见，这一步可以直接完成，然而，同样存在一些场景需要避免这样做。特别是，应尽可能避免写出算符$\hat{O}$的矩阵表示，因为这样做效率极低。尤其是，当算符具有稀疏矩阵的结构或表示时，实现一个计算算符作用于状态的代码比计算矩阵表示然后进行矩阵-向量乘法更高效。在标准子程序（如LAPACK）中，这些子程序通常被称为amul，并且通常是唯一需要编写的代码部分[130]。实际上，存储一个矩阵所需的内存是存储一个向量的平方倍，这使得内存限制更快地达到：这对于多体量子系统尤其如此，因为其哈密顿量通常仅包含少体相互作用，并且是具有明确张量结构的稀疏矩阵。因此，如果我们想避免显式计算算符的逆，我们可以将式(2.14)重写为

并将其识别为线性方程组：给定第k步迭代的状态，需要找到下一步的状态。要计算整个序列，我们需要针对不同的状态向量$\left| y_{k} \right.$求解相同的线性方程组，这正是上一节介绍LU分解时所研究的场景。因此，首先应计算LU分解，然后对式(2.15)中序列的每次迭代执行高斯消去。注意，尽管LU分解和求逆的复杂度均为$O ( N^{3} )$，但前者的预因子是后者的三分之一，因此仍更为方便。此外，LU分解保留了矩阵的稀疏性，而求逆则不允许利用问题的结构来提高效率。

最后，一旦计算出了第一个本征矢及其对应的本征值，可以重复该过程来计算下一个本征矢，前提是定期消除沿$\left| x_{0} \right.$方向的分量：也就是说，在初始随机状态创建后，应将其重新定义为

使得 $| y_{0} \rangle^{\prime} \bot | x_{0} \rangle$。式 (2.15) 中的序列将收敛到展开式$\begin{array} {r} {| y \rangle = \sum y_{j} | y_{j} \rangle} \end{array}$中包含的次小特征值。注意，在序列过程中应重复执行此正交化，因为数值误差可能重新引入沿$\left| x_{0} \right.$的分量，这些分量会增长并阻止算法正常工作。

2.3 三对角矩阵（Tridiagonal Matrices）#

三对角矩阵定义为仅在主对角线、上对角线和下对角线上有非零元素的矩阵，它们在物理学中扮演着重要角色，因为它们自然地出现在广泛的问题类别中，例如单体量子问题（见第 3 章）和紧束缚模型（tight binding models）。此外，正如我们将在第 2.4 节中看到的那样，大多数矩阵可以通过基变换写成这种形式，并且这种变换可以高效计算。最后，正如我们下面简要回忆的那样，求解三对角矩阵的特征值问题是非常高效的。实际上，给定一个如下形式的三对角矩阵

其特征多项式满足递推方程

因此，我们可以定义一维函数 $f_{N} ( \lambda ) = ( d_{N} - \lambda ) f_{N - 1} ( \lambda ) -$ $d_{N - 1}^{2} f_{N - 2} ( \lambda )$，其中 $f_{0} ( \lambda ) = 1$，并使用非常高效的标准方法来求其根，即原始矩阵在式 (2.2) 中的特征值 [122]。

2.4 Lanczos 方法（Lanczos Methods）#

Lanczos 对角化算法（Lanczos diagonalization algorithm）旨在找到从一般矩阵到三对角形式的变换，从而能够高效求解特征值问题 [131]。此外，它允许在不显式写出整个矩阵的情况下进行操作，从而极大地扩展了其应用领域。与幂方法（power methods）类似，它从相关希尔伯特空间（Hilbert space）中的一个随机向量 $\left| y_{1} \right.$ 开始，并将待对角化算符作用到它上得到 $| y_{2} \rangle = \hat{O} | y_{1} \rangle$。然后，它将新状态 $| y_{2} \rangle$ 分解为相对于原始状态 $\left| y_{1} \right.$ 的正交和并行分量，即：

其中，根据构造，$d_{1} = \Big \langle y_{1} \Big | \hat{O} y_{1} \Big \rangle , b_{1} | y_{2} \rangle \equiv \hat{O} | y_{1} \rangle - d_{1} | y_{1} \rangle$，并且可以通过施加 $\mid \langle b_{1} y_{2} | b_{1} y_{2} \rangle \mid^{2} = 1$ 来求出系数 $b_{1}$。得到前两个向量后，再迭代一次，但需要额外小心地将新状态分解为相对于前两个状态的正交和并行分量。然后得到：

其中，与之前一样，$b_{1} = \Big \langle y_{1} \Big | \hat{O} y_{2} \Big \rangle , d_{2} = \Big \langle y_{2} \Big | \hat{O} y_{2} \Big \rangle$，并且可以通过施加 $b_{2} | y_{3} \rangle$ 的归一化条件来计算系数 $b_{2}$。

在算法的前两次迭代之后，我们得到了系数 $d_{1} , d_{2} , b_{1}$ 和 $b_{2}$ 的值，以及状态 $| y_{1} \rangle , | y_{2} \rangle$ 和 $\left| y_{3} \right.$ 。该算法的主要亮点在第三次迭代中展现出来：再次将算符 $\hat{O}$ 应用于状态 $\left| y_{3} \right.$ ，并将其分解为相对于所有先前计算的状态的并行和正交分量，我们得到：

第一个系数如前定义 $r = \left. y_{1} \middle | \hat{O} y_{3} \right.$ ，由于 $\hat{O}$ 是厄米算子（Hermitian operator），且构造时保证了正交关系 $\vert y_{3} \rangle \bot \vert y_{2} \rangle , \vert y_{1} \rangle$，因此该系数为零：$y_{1} \mid \hat{O} y_{3} = \hat{O} y_{1} \mid y_{3} = ( d_{1} y_{1} \mid + b_{1}^{*} y_{2} \mid ) \mid y_{3} = 0$。一般情况下可以证明，对于任意连续迭代，有

综上，我们将前 $N$ 个向量 $\vert y_{k} \rangle$ 构成的矩阵定义为

方程 (2.22) 的关系意味着 $T_{N} = Y^{\dagger} {\hat{O}} Y$，其中 $T_{N}$ 是如方程 (2.17) 定义的三对角矩阵（tridiagonal matrix）。

注意，要计算三对角矩阵 $T_{N}$ 的系数，无需显式写出矩阵 $Y$，也无需写出原始算子 $\hat{O}$ 的完整矩阵表示——实际上，矩阵通常是稀疏的（sparse）或具有某种可利用的结构。事实上，$T_{N}$ 中的系数可以通过方程 (2.22) 获得：每一步只需存储向量 $\vert y_{k - 1} \rangle , \vert y_{k} \rangle , \vert y_{k + 1} \rangle$，以及两个系数向量 $a_{k}$ 和 $b_{k}$。一旦这些系数计算完毕，就可以应用 2.3 节回顾的方法，找出算子 $\hat{O}$ 在其三对角形式 $T_{N}$ 下的特征值和特征向量。因此，通过 Lanczos 方法（Lanczos method），即使是无法存入内存的巨大矩阵也可以对角化。我们请感兴趣的读者查阅关于 Lanczos 方法及其所基于的 Krylov 子空间（Krylov subspace）的广泛文献，例如 [132]。

2.5 习题#

1. 考虑一个大小为 N 的随机厄米矩阵 A。 (a) 编写一个程序，初始化 A 并执行 LU 分解。其复杂度随 N 如何缩放？ (b) 将 A 对角化，并按升序存储 N 个特征值 $\lambda_{i}$。 (c) 计算特征值之间的归一化间距 $s_{i} = \Delta \lambda_{i} / \bar{\Delta \lambda}$，其中

   $$ \Delta \lambda_{i} = \lambda_{i + 1} - \lambda_{i} , $$

   且 $\bar{\Delta \lambda}$ 是 $\Delta \lambda_{i}$ 的平均值。 (d) 计算在 $\lambda_{i}$ 附近不同层级数（即 $N / 100 , N / 50 , N / 10 \ldots N$）上的平均间距 $\bar{\Delta \lambda}$，并比较不同选择下下一个练习的结果。

2. 研究 $P ( s )$，即上一习题中定义的 $s_{i}$ 的分布，从不同大小的随机矩阵中积累 $s_{i}$ 的值。 (a) 计算随机厄米矩阵的 $P ( s )$。 (b) 计算对角线上为随机实数的对角矩阵的 $P ( s )$。 (c) 用以下函数拟合相应的分布：

   $$ P ( s ) = a s^{\alpha} \exp \bigl ( - b s^{\beta} \bigr ) $$

   并将其与随机矩阵理论（Random Matrix Theory, RMT）[133] 的预测进行比较。

3. 考虑方程 (4.4) 定义的伊辛哈密顿量（Ising Hamiltonian）。 (a) 编写一个程序，对不同的 N 计算 $2^{N} \times 2^{N}$ 矩阵 H。 (b) 对不同的 $N = 1 , \ldots , N$ 和 $\lambda \in [ 0 : 3 ]$ 将 H 对角化。 (c) 针对不同的 N，绘制前 k 个能级随 $\lambda$ 变化的图。讨论谱的结构。