计算化学公社

标题: Davidson 算法的preconditioner是如何推导出来的？ [打印本页]

作者
Author: JohnCase 时间: 2021-8-9 13:32
标题: Davidson 算法的preconditioner是如何推导出来的？
本帖最后由 JohnCase 于 2021-9-1 04:09 编辑

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

帖子快沉了。。。

作者
Author: wzkchem5 时间: 2021-8-9 15:04
你指的是为什么preconditioner是(D-\lambda)^{-1}这个形式，还是为什么D要取矩阵对角元？

作者
Author: wangxubo 时间: 2021-8-10 00:25
因为preconditioner用什么也不算是推导出来的，而是选一个。选矩阵对角元的主要argument就是，首先这个preconditioner是近似矩阵的逆，然后ci矩阵是diagonal dominant的，所以用对角元近似矩阵，然后1/Hii自然就近似矩阵的逆了。

作者
Author: wangxubo 时间: 2021-8-10 00:26
比较详细的讨论可以看Golub的matrix computation的第十章。

作者
Author: JohnCase 时间: 2021-8-10 00:55

wzkchem5 发表于 2021-8-9 15:04
你指的是为什么preconditioner是(D-\lambda)^{-1}这个形式，还是为什么D要取矩阵对角元？

抱歉表意不清，我是想问为什么preconditioner是(D-\lambda)^{-1}这个形式。

作者
Author: JohnCase 时间: 2021-8-10 01:16

wangxubo 发表于 2021-8-10 00:25
因为preconditioner用什么也不算是推导出来的，而是选一个。选矩阵对角元的主要argument就是，首先这个prec ...

对，我就是这里不明白，这个preconditioner是为啥是近似矩阵的逆呢？如果这个近似矩阵无限近似于原矩阵的话，那么precondition这一步相当于啥也没干，只能得到当前这步的猜想解。
我正在看你推荐的书，多谢！

作者
Author: wzkchem5 时间: 2021-8-10 03:16

JohnCase 发表于 2021-8-9 18:16
对，我就是这里不明白，这个preconditioner是为啥是近似矩阵的逆呢？如果这个近似矩阵无限近似于原矩阵的 ...

没错，这个逻辑是这样的：
由微扰理论可以推出Jacobi-Davidson方法（又称Olsen方法），Jacobi-Davidson方法里面，当preconditioner取成(矩阵-近似本征值)的逆的时候，收敛速度是cubic的，否则一般是linear的。Davidson方法是Jacobi-Davidson方法的近似，忽略了一个比较小的项，所以Davidson方法大体继承了preconditioner接近(矩阵-近似本征值)的逆比较容易收敛这个特性，但是当特别接近的时候就不成立了。
参见Olsen, J. et al. Chem. Phys. Lett., 169, 463-472 (1990)

作者
Author: JohnCase 时间: 2021-8-10 09:47
本帖最后由 JohnCase 于 2021-8-10 09:54 编辑

wzkchem5 发表于 2021-8-10 03:16
没错，这个逻辑是这样的：
由微扰理论可以推出Jacobi-Davidson方法（又称Olsen方法），Jacobi-Davidson ...

多谢! Jacobi-Davidson方法里的较小的第二项直接扔掉后就是(D-\lambda)^-1*residual的形式 ✓但是Davidson本人推导preconditioner的时候用的瑞利商泰勒级数展啥的，一顿操作根本看不懂（也是微扰理论？）…… 我非常好奇他原本的的详细思路是啥？毕竟Davidson在1975年就推导的，比Olsen1990要早

作者
Author: SHENLIN 时间: 2025-7-7 09:36
占据轨道5，虚轨道44，nstate=3，为啥我写得程序一直收敛不了？
Number of occupied orbitals :    5
Number of virtual orbitals : 44
r-norm
   0.0904215215    0.0954082545    0.0955254642
Elap time for1 iteration: 0.05 min
r-norm
   0.1020421580    0.1005147596    0.3195436706
Elap time for2 iteration: 0.07 min
r-norm
   0.1147001275    0.1182726230    0.4759256149
Elap time for3 iteration: 0.07 min
r-norm
   0.0912756888    0.1447432953    1.6792727450
Elap time for4 iteration: 0.08 min
r-norm
   0.0847314365    0.1083200385    1.3451540912
Elap time for5 iteration: 0.10 min
r-norm
   0.1017993397    0.1161335163    0.9611721320
Elap time for6 iteration: 0.08 min
r-norm
   0.0977623559    0.1457940707    1.1992824782
Elap time for7 iteration: 0.10 min
r-norm
   0.0845577392    0.1163889397    0.7997297814
Elap time for8 iteration: 0.12 min
r-norm
   0.0913815993    0.1284229977    0.7724938236
Elap time for9 iteration: 0.12 min
r-norm
   0.1473733525    0.1621629793    0.9306342023
Elap time for10 iteration: 0.10 min
r-norm
   0.1022261208    0.1520397934    0.8568445039
Elap time for11 iteration: 0.13 min
r-norm
   0.0747816170    0.1215488567    0.5824591344
Elap time for12 iteration: 0.12 min
r-norm
   0.1038876657    0.1438944762    0.6743080739
Elap time for13 iteration: 0.13 min
r-norm
   0.1119399114    0.1492259623    0.7073477389
Elap time for14 iteration: 0.12 min
r-norm
   0.0716162756    0.0678459788    0.4724069512
Elap time for15 iteration: 0.13 min
r-norm
   0.0685578300    0.0721142019    0.4848637315
Elap time for16 iteration: 0.15 min
r-norm
   0.0506047135    0.0770437986    0.3580719103
Elap time for17 iteration: 0.15 min
r-norm
   0.0332001995    0.0589313271    0.2426525056
Elap time for18 iteration: 0.15 min
r-norm
   0.0356365333    0.0585605558    0.2491546503
Elap time for19 iteration: 0.15 min
r-norm
   0.0235190389    0.0403184234    0.1662902469
Elap time for20 iteration: 0.15 min
r-norm
   0.0168157222    0.0429556733    0.1319907180
Elap time for21 iteration: 0.17 min
r-norm
   0.0153743840    0.0338096331    0.1234901048
Elap time for22 iteration: 0.15 min
r-norm
   0.0123653356    0.0204147965    0.1167310930
Elap time for23 iteration: 0.18 min
r-norm
   0.0085262483    0.0162310911    0.2834981384
Elap time for24 iteration: 0.17 min
r-norm
   0.0063609159    0.0122873015    0.2248542586
Elap time for25 iteration: 0.18 min
r-norm
   0.0034141409    0.0105203632    0.1803029904
Elap time for26 iteration: 0.18 min
r-norm
   0.0022679864    0.0097654726    0.1318347989
Elap time for27 iteration: 0.20 min
r-norm
   0.0018805906    0.0092926372    0.1114406811
Elap time for28 iteration: 0.20 min
r-norm
   0.0016178018    0.0096152892    0.1128985775
Elap time for29 iteration: 0.20 min
r-norm
   0.0015925451    0.0102707210    0.1566133927
Elap time for30 iteration: 0.18 min
r-norm
   0.0012307279    0.0092470712    0.1726686048
Elap time for31 iteration: 0.20 min
r-norm
   0.0011781288    0.0089239939    0.1487790273
Elap time for32 iteration: 0.22 min
r-norm
   0.0010114222    0.0072834724    0.1089459243
Elap time for33 iteration: 0.22 min
r-norm
   0.0005997315    0.0058661332    0.1400052555
Elap time for34 iteration: 0.20 min
r-norm
   0.0003976856    0.0042183577    0.1415724920
Elap time for35 iteration: 0.23 min
r-norm
   0.0002629169    0.0027596067    0.1010022602
Elap time for36 iteration: 0.23 min
r-norm
   0.0001700379    0.0021657765    0.1153183994
Elap time for37 iteration: 0.23 min
r-norm
   0.0001253134    0.0013222387    0.1614236243
Elap time for38 iteration: 0.23 min
r-norm
   0.0000928468    0.0008789330    0.1389585514
Elap time for39 iteration: 0.25 min
r-norm
   0.0000759871    0.0005949575    0.0970548157
Elap time for40 iteration: 0.25 min
r-norm
   0.0000435233    0.0003580981    0.1001310626
Elap time for41 iteration: 0.25 min
r-norm
   0.0000275080    0.0002518417    0.0915286161
Elap time for42 iteration: 0.25 min
r-norm
   0.0000122202    0.0001709833    0.0891311069
Elap time for43 iteration: 0.27 min
r-norm
   0.0000071119    0.0001243516    0.0514699726
Elap time for44 iteration: 0.27 min
r-norm
   0.0000055632    0.0000882621    0.0405689223
Elap time for45 iteration: 0.27 min
r-norm
   0.0000042761    0.0000802745    0.0215980070
Elap time for46 iteration: 0.27 min
r-norm
   0.0000029298    0.0000581982    0.0124480320
Elap time for47 iteration: 0.28 min
r-norm
   0.0000024042    0.0000446184    0.0090657784
Elap time for48 iteration: 0.28 min
r-norm
   0.0000020247    0.0000367220    0.0086109354
Elap time for49 iteration: 0.30 min
r-norm
   0.0000018733    0.0000243191    0.0080043443
Elap time for50 iteration: 0.30 min
Warning: Algorithm did not converge!!

作者
Author: SHENLIN 时间: 2025-7-7 09:39
和上面一样的,nstate=5需要迭代43圈才收敛，nstate=10需要22圈才收敛。到底是啥问题？？
Number of occupied orbitals :    5
Number of virtual orbitals : 44
r-norm
   0.0859670869    0.0961734206    0.0961779972    0.0823855538    0.0823923430
Elap time for1 iteration: 0.07 min
r-norm
   0.0939422432    0.1136267943    0.1133913454    0.0932864648    0.0919924146
Elap time for2 iteration: 0.10 min
r-norm
   0.1104706139    0.1402461731    0.1362228944    0.1047233933    0.0992314979
Elap time for3 iteration: 0.10 min
r-norm
   0.0839918056    0.1276662383    0.1345923188    0.1114826929    0.1167705100
Elap time for4 iteration: 0.10 min
r-norm
   0.0879364097    0.1367220350    0.1345198104    0.0817144452    0.1023588033
Elap time for5 iteration: 0.12 min
r-norm
   0.1333545255    0.1256084830    0.1509244491    0.0733161709    0.0816386380
Elap time for6 iteration: 0.12 min
r-norm
   0.0918918996    0.1194373967    0.1129820572    0.0778894700    0.0745357280
Elap time for7 iteration: 0.13 min
r-norm
   0.0945563448    0.1496321715    0.1360978203    0.0954277892    0.1042596870
Elap time for8 iteration: 0.13 min
r-norm
   0.1382295514    0.1113698975    0.1474383490    0.0743684072    0.0854812935
Elap time for9 iteration: 0.13 min
r-norm
   0.1025797830    0.1000518348    0.0937715458    0.0662613920    0.0638135975
Elap time for10 iteration: 0.15 min
r-norm
   0.0931504528    0.0922992344    0.1080575882    0.0620406054    0.0748651619
Elap time for11 iteration: 0.17 min
r-norm
   0.0926317499    0.0772896051    0.0809489506    0.0613332612    0.0641595639
Elap time for12 iteration: 0.15 min
r-norm
   0.0665347374    0.0629086255    0.0770907103    0.0425396491    0.0594336014
Elap time for13 iteration: 0.18 min
r-norm
   0.0523427534    0.0622265817    0.0625065519    0.0356882962    0.0386663151
Elap time for14 iteration: 0.17 min
r-norm
   0.0499933112    0.0456126482    0.0626748593    0.0299150655    0.0414691036
Elap time for15 iteration: 0.18 min
r-norm
   0.0278982335    0.0292917812    0.0391304847    0.0252510147    0.0324943526
Elap time for16 iteration: 0.20 min
r-norm
   0.0221077037    0.0167310702    0.0230958364    0.0192126420    0.0215620310
Elap time for17 iteration: 0.22 min
r-norm
   0.0137119275    0.0121285101    0.0149942035    0.0181570962    0.0159096701
Elap time for18 iteration: 0.20 min
r-norm
   0.0101229064    0.0065642979    0.0114052417    0.0246803063    0.0145937636
Elap time for19 iteration: 0.22 min
r-norm
   0.0053748381    0.0045956244    0.0060060404    0.0340680875    0.0173046677
Elap time for20 iteration: 0.23 min
r-norm
   0.0032055613    0.0030581182    0.0037756580    0.1187479188    0.0225333528
Elap time for21 iteration: 0.23 min
r-norm
   0.0019937965    0.0019969316    0.0021297509    0.1485934236    0.0355366962
Elap time for22 iteration: 0.23 min
r-norm
   0.0015196083    0.0011215612    0.0015955654    0.0864943371    0.0502523481
Elap time for23 iteration: 0.25 min
r-norm
   0.0011999264    0.0007936677    0.0009122443    0.0704919186    0.0892320727
Elap time for24 iteration: 0.25 min
r-norm
   0.0008931291    0.0005064775    0.0006463799    0.0716211999    0.1176367125
Elap time for25 iteration: 0.27 min
r-norm
   0.0006871657    0.0002780361    0.0004348996    0.0506645645    0.1007364757
Elap time for26 iteration: 0.27 min
r-norm
   0.0005234747    0.0001451260    0.0003245010    0.0397565107    0.0870168418
Elap time for27 iteration: 0.27 min
r-norm
   0.0003975630    0.0000903843    0.0002054682    0.0370669757    0.0694277818
Elap time for28 iteration: 0.32 min
r-norm
   0.0003009416    0.0000650028    0.0001465787    0.0258479447    0.0710442528
Elap time for29 iteration: 0.32 min
r-norm
   0.0001913904    0.0000389681    0.0001114010    0.0177343223    0.0602410773
Elap time for30 iteration: 0.35 min
r-norm
   0.0001331266    0.0000345184    0.0000624007    0.0101300747    0.0347834598
Elap time for31 iteration: 0.30 min
r-norm
   0.0000798987    0.0000438127    0.0000215644    0.0083319747    0.0251518196
Elap time for32 iteration: 0.32 min
r-norm
   0.0000452280    0.0000222379    0.0000115394    0.0077254885    0.0145465833
Elap time for33 iteration: 0.32 min
r-norm
   0.0000248855    0.0000115116    0.0000058604    0.0059863711    0.0076289659
Elap time for34 iteration: 0.33 min
r-norm
   0.0000072587    0.0000049618    0.0000030861    0.0035746366    0.0040836744
Elap time for35 iteration: 0.33 min
r-norm
   0.0000037596    0.0000030474    0.0000014237    0.0014577607    0.0024615593
Elap time for36 iteration: 0.35 min
r-norm
   0.0000020340    0.0000008507    0.0000005975    0.0006395652    0.0010532152
Elap time for37 iteration: 0.47 min
r-norm
   0.0000009556    0.0000002606    0.0000003213    0.0004556286    0.0002249363
Elap time for38 iteration: 0.47 min
r-norm
   0.0000002619    0.0000001297    0.0000000842    0.0001885796    0.0001356042
Elap time for39 iteration: 0.37 min
r-norm
   0.0000000769    0.0000000345    0.0000000754    0.0000684111    0.0000568973
Elap time for40 iteration: 0.37 min
r-norm
   0.0000000686    0.0000000138    0.0000000510    0.0000336428    0.0000179266
Elap time for41 iteration: 0.37 min
r-norm
   0.0000000663    0.0000000107    0.0000000477    0.0000185220    0.0000102552
Elap time for42 iteration: 0.50 min
r-norm
   0.0000000432    0.0000000099    0.0000000336    0.0000037480    0.0000024045
Converged in 43 iterations.

作者
Author: wzkchem5 时间: 2025-7-7 10:38

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

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

作者
Author: SHENLIN 时间: 2025-7-31 09:21
本帖最后由 SHENLIN 于 2025-7-31 10:44 编辑

谢谢老师指导，确实是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

作者
Author: hebrewsnabla 时间: 2025-7-31 12:21
本帖最后由 hebrewsnabla 于 2025-7-31 12:23 编辑

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

你的max space是多少？

作者
Author: SHENLIN 时间: 2025-7-31 13:48
   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这样写就是错误的，一直收敛不了。这两者不是等价的吗？好奇怪

作者
Author: SHENLIN 时间: 2025-7-31 16:11
本帖最后由 SHENLIN 于 2025-7-31 16:16 编辑

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

你的max space是多少 ...

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

作者
Author: hebrewsnabla 时间: 2025-7-31 16:25

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

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

作者
Author: JohnCase 时间: 2025-8-7 18:19

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

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

作者
Author: JohnCase 时间: 2025-8-7 18:29

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

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

作者
Author: SHENLIN 时间: 2025-9-2 10:36
本帖最后由 SHENLIN 于 2025-9-2 10:56 编辑

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

是这个维度，就是对角线元素。我还是没搞明白之前preconditiner是这样写得 (, 下载次数 Times of downloads: 2) (, 下载次数 Times of downloads: 2)
我看MRCC程序，把分母改成下面的这样
(, 下载次数 Times of downloads: 2)
就没有问题了。我觉得是等价的啊

作者
Author: SHENLIN 时间: 2025-9-2 10:50
本帖最后由 SHENLIN 于 2025-9-2 10:57 编辑

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

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

作者
Author: wzkchem5 时间: 2025-9-2 10:58

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

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

作者
Author: JohnCase 时间: 2025-9-3 12:37
本帖最后由 JohnCase 于 2025-9-3 14:23 编辑

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

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

lowest等于size(residues, 2)吗

作者
Author: SHENLIN 时间: 2025-9-3 16:54
本帖最后由 SHENLIN 于 2025-9-3 17:03 编辑

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

lowest等于size(residues, 2)吗

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

作者
Author: JohnCase 时间: 2025-9-3 18:18

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

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

欢迎光临计算化学公社 (http://bbs.keinsci.com/)