[矩阵计算]Davidson对角化
更新： 8 AUG 2016
花了几个礼拜写程序终于跑过Davidson对角化！至此，Davidson对角化的思路已经完全清晰。如尚有不准确之处，请务必回复指出！
一、Davidson对角化的思路
Davidson对角化是一种快速求出大规模稀疏矩阵的方法，对于求量子体系中\(\textbf{H}|C\rangle=E|C\rangle\)问题有极大的帮助。其基本原理同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.
此处总结其算法思路。
-
在某基组下给定大规模的矩阵\(A\)，要求求其最小的本征值。任取一组标准正交基，如直接从原来的基组中截取数个即可，作为初始的基组猜测，记作\(B\)。在此子空间内表示\(A\)为\(A_0\)，进行对角化求其最小本征值及对应的本征函数，分别记作\(\lambda_0\)和\(b_0\)。这里A在B中的矩阵表示\(A_0\)规模很小，可以通过QR方法等经典的对角化方法求解。
注：初始的基组构成一个子空间（Krylov子空间），Davidson对角化即通过一定的挑选不断地向这个子空间添加基函数，使其最小本征值\(\lambda_0\)趋近不变。收敛的证明略。 -
牛顿法构造新的基函数：
然后从中扣除投影在子空间B中的部分：
剩下的部分即与子空间B正交的新基，将其归一化并加入到\(B\)中。
注：实际上就是将\(b\)加入到\(B\)中并正交归一化。\(A^0\)表示矩阵\(A\)的对角元组成的对角矩阵。
- 重新计算\(A\)在新的\(B\)中的矩阵表示\(A_0\)，并求其最低本征值及对应的本征函数。重复2,3直到子空间\(B\)不变，具体说就是最低本征值\(\lambda_0\)的变化小于阈值。
注：这个判别等价于对正交化（但未归一）的\(b\)求其模长并判别，模长小于阈值时，说明新加入的本征函数对子空间的贡献极小。
以上即主要步骤，算法的详细步骤参看上文献III节。
二、Davidson对角化在计算中的简单优化
思路1. 避免重复计算\(A\)在\(B\)中的投影，每次新加入基组时对称地将\(Ab\)加入即可。对于步骤2中
可以改写为
一般来说，\(Ab_0\)是每步循环中最耗时的步骤，计算使应当单独写成函数并尽可能利用\(A\)的稀疏性与对称性减少计算。
思路2. 简化牛顿法。子空间B的规模会在循环中不断扩大，对\(A_0\)的对角化成本也会越来越高。这一点可以通过以下方法简化：
每步中仅将\(b_0\)和新的到的\(b\)作为基组构成2维子空间，将\(A\)在此子空间中表示，是一个2x2矩阵，对角化成本极低，但是可能需要增加循环步数。
或
直接将与\(b_0\)正交化，但未归一化的\(b\)与\(b_0\)相加得到新的\(b_0\)，免去对角化操作。
思路3. 减少矩阵乘法
中，$ \prod\limits_{i=1}M(I-B_iB_iT) $可以被简化为 $ (I-\sum\limits_{i=1}MB_iB_iT) \(，这里用到了\)B\(基组的正交性。这一步可能会节省很多时间。此外这一步不需要每次循环时计算，即每次在向B中扩充新基的时候减去新的\)bb^T$即可。
三、可能遇到的问题
若初始的基组\(B\)只取一个基函数，这时本征值与对应的\(A\)矩阵对角元相等，导致
中遭遇奇异矩阵。此时可以将对角矩阵的对应元加一个正的微扰即可。此处\((A^0-\lambda_0I)\)最好保持正定，特别是在使用简化的牛顿法时。牛顿法要求Hesse矩阵保持正定以求函数的极小值点。
同时，初始的子空间的估计最好能够比较接近真实的情况，否则有可能出现收敛到别的本征值上的问题，特别是简化的牛顿法。
