[矩阵计算]Davidson对角化

更新： 8 AUG 2016

花了几个礼拜写程序终于跑过Davidson对角化！至此，Davidson对角化的思路已经完全清晰。如尚有不准确之处，请务必回复指出！

一、Davidson对角化的思路

Davidson对角化是一种快速求出大规模稀疏矩阵的方法，对于求量子体系中( extbf{H}|C angle=E|C angle)问题有极大的帮助。其基本原理同Lanczos对角化很相似，此外还可以从牛顿法推出（见Molecular Electronic-Structure Theory）。其数学原理和证明不在此处列出，有深入需求者建议参看Davidson的原始论文：

E. R. Davidson.
The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real symmetric matrices.
J. Comput. Phys., 17:87-94, 1975.

此处总结其算法思路。

1. 在某基组下给定大规模的矩阵(A)，要求求其最小的本征值。任取一组标准正交基，如直接从原来的基组中截取数个即可，作为初始的基组猜测，记作(B)。在此子空间内表示(A)为(A_0)，进行对角化求其最小本征值及对应的本征函数，分别记作(lambda_0)和(b_0)。这里A在B中的矩阵表示(A_0)规模很小，可以通过QR方法等经典的对角化方法求解。
   注：初始的基组构成一个子空间（Krylov子空间），Davidson对角化即通过一定的挑选不断地向这个子空间添加基函数，使其最小本征值(lambda_0)趋近不变。收敛的证明略。

2. 牛顿法构造新的基函数：

[b=-(A^0-lambda_0I)^{-1}(A-lambda_0I)b_0 ]

然后从中扣除投影在子空间B中的部分：

[b:=prodlimits_{i=1}^M(I-B_iB_i^T) b ]

剩下的部分即与子空间B正交的新基，将其归一化并加入到(B)中。
注：实际上就是将(b)加入到(B)中并正交归一化。(A^0)表示矩阵(A)的对角元组成的对角矩阵。

3. 重新计算(A)在新的(B)中的矩阵表示(A_0)，并求其最低本征值及对应的本征函数。重复2,3直到子空间(B)不变，具体说就是最低本征值(lambda_0)的变化小于阈值。
   注：这个判别等价于对正交化（但未归一）的(b)求其模长并判别，模长小于阈值时，说明新加入的本征函数对子空间的贡献极小。

以上即主要步骤，算法的详细步骤参看上文献III节。

二、Davidson对角化在计算中的简单优化

思路1. 避免重复计算(A)在(B)中的投影，每次新加入基组时对称地将(Ab)加入即可。对于步骤2中

[b=-(A^0-lambda_0I)^{-1}(A-lambda_0I)b_0 ]

可以改写为

[b=-(A^0-lambda_0I)^{-1}(Ab_0-lambda_0b_0) ]

一般来说，(Ab_0)是每步循环中最耗时的步骤，计算使应当单独写成函数并尽可能利用(A)的稀疏性与对称性减少计算。

思路2. 简化牛顿法。子空间B的规模会在循环中不断扩大，对(A_0)的对角化成本也会越来越高。这一点可以通过以下方法简化：
每步中仅将(b_0)和新的到的(b)作为基组构成2维子空间，将(A)在此子空间中表示，是一个2x2矩阵，对角化成本极低，但是可能需要增加循环步数。
或
直接将与(b_0)正交化，但未归一化的(b)与(b_0)相加得到新的(b_0)，免去对角化操作。

思路3. 减少矩阵乘法

[b:=prodlimits_{i=1}^M(I-B_iB_i^T) b ]

中，$ prodlimits_{i=1}^M(I-B_iB_iT) $可以被简化为 $ (I-sumlimits_{i=1}^MB_iB_iT) (，这里用到了)B(基组的正交性。这一步可能会节省很多时间。此外这一步不需要每次循环时计算，即每次在向B中扩充新基的时候减去新的)bb^T$即可。

三、可能遇到的问题

若初始的基组(B)只取一个基函数，这时本征值与对应的(A)矩阵对角元相等，导致

[b=-(A^0-lambda_0I)^{-1}(A-lambda_0I)b_0 ]

中遭遇奇异矩阵。此时可以将对角矩阵的对应元加一个正的微扰即可。此处((A^0-lambda_0I))最好保持正定，特别是在使用简化的牛顿法时。牛顿法要求Hesse矩阵保持正定以求函数的极小值点。
同时，初始的子空间的估计最好能够比较接近真实的情况，否则有可能出现收敛到别的本征值上的问题，特别是简化的牛顿法。