量子化学中的算法(一)—大型稀疏矩阵对角化

(准备开个新坑，简单讲讲量子化学中常用的算法，主要是计算数学经典算法在量子化学方法中的应用，话不多说直接进入正题吧。)

求解矩阵的本征值及其对应的本征向量是数值代数的核心问题之一。高效、准确地求解矩阵的本征值和本征向量的算法对于量子化学程序的发展至关重要，也关乎到相应量子化学方法可以应用的范围。同时为了求解量子化学中常用的大型稀疏矩阵的本征值及其本征态，计算数学和量子化学工作者也发展了各种各样的矩阵对角化算法，本文将主要介绍量子化学中常用的 Davidson 对角化方法。

组态相互作用 (Configuration Interaction) 方法是一种常用的量子化学方法，其在 HF 方法的基础上用于修正电子相关。该方法将波函数写成如下形式：

选取合适的组态空间，求解薛定谔方程约化为如下求解矩阵本征值和本征向量问题:

$\mathbf{H}\mathbf{c}=\lambda\mathbf{c}$

其中 $H_{IJ}=\langle\Phi_I|H|\Phi_J\rangle$

一、矩阵对角化方法

(1) QR 方法

该方法是一种重要的求解矩阵全部特征值和特征向量的方法。然而该方法需要储存全部矩阵元、时间复杂度为 $O(n^3)$ 、破坏矩阵稀疏性，且量子化学问题一般只需要求解基态、低激发态或特定能量范围内的本征值，因而量子化学中一般不利用该方法去对角化矩阵。

(2) 大型稀疏矩阵对角化

哈密顿量矩阵的稀疏性来源于如下几个方面:

(1) 自旋对称性和空间对称性;

(2) 哈密顿算符是二体算符;

(3) 分子积分

其中第 (1) 条可通过自旋匹配、对称性匹配显式考虑，实际计算中的稀疏性一般来源于 (2),(3)

由于矩阵是稀疏的，从而可以快速进行矩阵向量乘这一基本操作(时间复杂度 $O(n^2)$ )。大型稀疏矩阵对角化基本思路便是从一个试探向量 $\bm{c}_0$ 出发，通过矩阵向量乘操作，同时保持矩阵的稀疏性，使得试探向量逐渐收敛到目标特征向量(一般是基态对应的特征向量)。

下简要推导该类方法基本方法–Lanczos Process 的基本思路。

求解基态本征值可对应于如下优化问题:

$\lambda_{\text{min}}=\min_{\bm{c}}\rho(\bm{c})=\min_{\bm{c}}\frac{\bm{c}^T\mathbf{H}\bm{c}}{\bm{c}^T\bm{c}}$

如果利用最速下降法求解如上优化问题，需要计算该函数的梯度。

$\left.\nabla\rho(\bm{c})\right|_{\bm{c}=\bm{c}_0}=\frac{2}{\bm{c}_0^T\bm{c}_0}(\mathbf{H}\bm{c}_0-\rho(\bm{c}_0)\bm{c}_0)$

而实际计算中并不许需要求出并做精确线搜索，因为解一定在如下子空间中 (Krylov Subspace):

$\text{span}(\bm{c}_0,\mathbf{H}\bm{c}_0)$

只需计算这两个向量之间的哈密顿矩阵元，并对角化所得到的小矩阵便相当于做了一步最速下降法。而经过 k 步迭代之后所得到的解在如下子空间中:

$\text{span}(\bm{c}_0,\mathbf{H}\bm{c}_0,\cdots,\mathbf{H}^k\bm{c}_0)$

于是大型矩阵对角化的问题，转化为子空间内矩阵对角化的问题，选取合适的初猜，经过若干次迭代之后，子空间内最小特征值可能于真实的最小特征值非常接近。(其实子空间内若干个最小特征值都可能于相应的真值非常接近)。而所求得第 j 个特征向量写成:

$\bm{c}^{(k)}_j=\sum_{i=0}^kc^{(k)}_{i,j}\mathbf{H}^i\bm{c}_0$

该方法收敛性由如下定理给出: (见《Matrix Computation》,Golub)

其中 $T_k$ 是第 k 次迭代得到的上三角矩阵， $\theta_i,\lambda_i$ 分别是子空间内第 i 个特征值和第 i 个特征值的真值。 $z_i$ 是原矩阵第 i 个特征向量。而 $cos(\phi_i)$ 描述的是初猜的质量。 $c_k()$ 是第 k 次切比雪夫多项式。

从这个定理可以得到该类方法的基本性质:

(a) 收敛速度与初始猜测的质量关系很大;

(b) 简并态收敛速度可能较慢;

然而该方法本质上是基于一阶导数，因而其收敛性质并不理想，可考虑更好的迭代修正方法。

(3) 其他思路

其实原优化问题在只变动一个分量 $c_I$ 的前提下可以精确求解:

$\left.\frac{\partial\rho}{\partial c_I}\right|_{c_I+\delta_I} = 0$

解为 $\delta_I=(\rho(\bm{c}_0)-H_{II})^{-1}q_I,\quad\bm{q}=(\mathbf{H}-\rho\bm{I})\bm{c}_0$ 。(其实 q 向量是残差向量，可用于判断收敛。)

而如果精确求解更一般的问题

$\left.\frac{\partial\rho}{\partial \bm{c}}\right|_{\bm{c}+\bm{\delta}} = \bm{0}$

将该方程展开到二阶不难发现解可以近似为:

$(\rho-H_{II})\delta_I\approx q_I+\sum_{J\neq I}H_{IJ}\delta_J+(\rho-\lambda)c_I$

而在实际计算中如果选取 $\delta_I=(\rho(\bm{c}_0)-H_{II})^{-1}q_I$ ，并不是上述方程解的好的近似，但好处是计算比较简单。而如果采取 Gauss-Seidal 迭代方法求解上述方程，计算量过大。

为了改进上述方法的不足之处，最好的方法便是结合二者的优点。也就是利用 $\delta_I=(\rho(\bm{c}_0)-H_{II})^{-1}q_I$ 改进初始猜测，而不是 $\mathbf{H}\bm{c}_0$ ,这便是量子化学中最常用的 Davidson 方法(Davidson 1975)。

二、Davidson 方法

下仅以求解基态能量为例，简述Davidson方法.

选取合适数目正交归一的向量 $\bm{c}_1,\bm{c}_2,\cdots,\bm{c}_l$ 作为初猜子空间的基组, 计算并储存如下向量和矩阵元 $\mathbf{H}\bm{c}_i, \tilde{H}_{ij} = (\bm{c}_i,\mathbf{H}c_j)$ . 对角化该矩阵，得到本征值 $\lambda^{(l)}$ 和本征向量 $\bm{\alpha}^{(l)}$ 。
构造 $\bm{q}_M=(\mathbf{H}-\lambda^{(M)}\mathbf{I})\bm{\alpha}^{(M)}$ ,其中 $\bm{\alpha}^{(M)}=\sum_{i=1}^M\alpha_{i}^{(M)}\bm{c}_i$
根据 $||\bm{q}_M||$ 判断收敛，如果收敛结束迭代；
构造 $\delta_{I,M+1}=(\lambda^{(M)}-H_{II})^{-1}q_{I,M}$ ，并与之前的基组做正交化，归一化后得到 $\bm{c}_{M+1}$
计算 $\tilde{H}_{i,M+1},i=1,2,\cdots,M+1$ 。对角化该矩阵得到新的本征值和本征向量，返回步 2 继续迭代。