Davidson 算法的preconditioner是如何推导出来的？

JohnCase

求解大型稀疏对角矩阵特征值的Davidson算法，使用矩阵对角元素作为preconditioner, 以此计算出下一循环中的subspace basis，即 basis_new = (D - \lambda)^-1 *residual
Davidson在原文献里给出了这一步的证明（见附件），但是没有给出公式（3）具体是如何推导出来的。
我查阅了关于Davidson算法的很多资料，但是唯独没有人复述Davidosn本人的推导过程，只是在陈述如何code。。。当然，由Jacobi-Davidson近似一下就可以得出Davidson preconditioner。
总觉得这样用着不太安心, 这用了快五十年的算法究竟是怎么推理出来的？
求助老师们帮忙解答疑惑 :-)


帖子快沉了。。。

JohnCase

SHENLIN 发表于 2025-9-3 16:54
是等于这个，就是要算的态的个数。为啥最好取这个，我看高斯前两圈取的都是nstate*3

感觉这是两段代码最不一样的地方

SHENLIN

JohnCase 发表于 2025-9-3 12:37
初猜数量是多少？算10个态的话，最好是设置前18个态作为初猜。

lowest等于size(residues, 2)吗

是等于这个，就是要算的态的个数。为啥最好取这个，我看高斯前两圈取的都是nstate*3

JohnCase

SHENLIN 发表于 2025-9-2 10:50
再次感谢老师，不然这个bug不知道要找到什么时候。请问还有什么加速收敛的策略吗？我看每一圈计算norm的 ...

初猜数量是多少？算10个态的话，最好是设置前18个态作为初猜。

lowest等于size(residues, 2)吗

wzkchem5

SHENLIN 发表于 2025-9-2 10:50
再次感谢老师，不然这个bug不知道要找到什么时候？请问还有什么加速收敛的策略吗？我看每一圈计算norm的 ...

大的特征值收敛慢是正常的，比如你算10个特征值，第10个特征值和第11个特征值的距离一般很接近，导致需要迭代很多轮才能把第10个特征值和第11个特征值之间的矩阵非对角元消干净。
迭代次数少一般是因为第一步求解比较多的特征值，收敛后再把高能量的解扔掉。比如第一步求解20个根，但从第二步开始只计算前10个根的residual vector，后10个根仍然放在Krylov子空间里，但不产生新的vector，也就是Krylov子空间按照20->30->40...的规律增加

SHENLIN

wzkchem5 发表于 2025-7-7 10:38
感觉大体还是有收敛趋势的，可能preconditioner实现得有问题？可以对角化一个对角矩阵，看是不是一步收敛 ...

再次感谢老师，不然这个bug不知道要找到什么时候。请问还有什么加速收敛的策略吗？我看每一圈计算norm的时候，前面小的特征值收敛比较快，后面大的收敛圈数要多得多，请问这里有什么改进的方法吗？另外我用高斯迭代的圈数好像要少很多，它是用了什么策略吗？

SHENLIN

JohnCase 发表于 2025-8-7 18:29
diag(i)是什么，i是nocc*vir这个维度吗

是这个维度，就是对角线元素。我还是没搞明白  之前preconditiner是这样写得
我看MRCC程序，把分母改成下面的这样

就没有问题了。我觉得是等价的啊

JohnCase

SHENLIN 发表于 2025-7-31 13:48
do j = 1, lowest
            do i=1,nocc
                do a=1,nvir

diag(i)是什么，i是nocc*vir这个维度吗

JohnCase

SHENLIN 发表于 2025-7-31 16:11
好的，谢谢。一定要设定max space吗？直接用最大迭代圈数控制可以吗？

讲道理不需要设定max space。或者说设定max space不能大于A矩阵大小

hebrewsnabla

SHENLIN 发表于 2025-7-31 16:11
好的，谢谢。一定要设定max space吗？直接用最大迭代圈数控制可以吗？

我不清楚你的程序有没有trial vector的space，也就是个数的概念，如果没有的话，大概是求n个根就只有n个trial vector？那效率应该是比较低的。
如果有space的概念，每一步会扩大space，但是没有max space（超过了就重来），也可能会有收敛问题，但也不一定。

SHENLIN

hebrewsnabla 发表于 2025-7-31 12:21
可以参考一下这个实现 https://github.com/pyscf/pyscf/b ... nalg_helper.py#L291

你的max space是多少 ...

好的，谢谢。一定要设定max space吗？直接用最大迭代圈数控制可以吗？

SHENLIN

       do j = 1, lowest
            do i=1,nocc
                do a=1,nvir
                ia=(i-1)*nvir+a
                denom = -eorb(a+nocc)+eorb(i)+eigenvalues(j)           
                ! 处理分母过小的情况
                if (abs(denom) < t) then
                    denom = sign(t, denom)
                end if           
                ! 计算修正项
                correction(ia, j) = residues(ia, j) / denom
                end do
            end do
        end do   我这样写preconditioner就是正确的，能正常收敛
     do j = 1, size(residues, 2)
        do i = 1, size(correction, 1)
            D = eigenvalues(j) - diag(i)
            
            ! 处理分母过小的情况
            if (abs(D) < t) then
                D = sign(t, D)
            end if
            
            ! 计算修正项
            correction(i, j) = residues(i, j) / D
        end do
    end do这样写就是错误的，一直收敛不了。这两者不是等价的吗？好奇怪

hebrewsnabla

可以参考一下这个实现 https://github.com/pyscf/pyscf/b ... nalg_helper.py#L291

你的max space是多少？

SHENLIN

谢谢老师指导，确实是diag的原因，因为我如果把矩阵元都求出来，然后取对角线元素，上面的例子16步左右能收敛。但是我在写CIS的过程中，为了节约存储，不是不求出H(ia,jb),而是直接求sigma向量HV吗？所以我的diag取的都是HF轨道能级差项，忽略了后面电子积分项（我看很多程序好像都是这么做的）。初猜和preconditioner用的都是这个，需要怎么调整呢？
do i=1,nocc
     do a=1,nvir
           ia=(i-1)*nvir+a
            diag(ia)=eorb(a+nocc)-eorb(i)
       enddo
enddo

wzkchem5

SHENLIN 发表于 2025-7-7 09:36
占据轨道5，虚轨道44，nstate=3，为啥我写得程序一直收敛不了？
Number of occupied orbitals   :     5
...

感觉大体还是有收敛趋势的，可能preconditioner实现得有问题？可以对角化一个对角矩阵，看是不是一步收敛，这样验证preconditioner有没有实现对