计算化学公社

标题: 请教Python求本征值的最佳方法 [打印本页]

作者
Author: granvia 时间: 2020-4-6 10:54
标题: 请教Python求本征值的最佳方法
用Numpy编写Lowdin正交化的脚本时发现，使用eig函数得到的本征值精度不是特别高，而采用svd函数获取的本征值精度要高好几个数量级，并且计算效率也快些。此外，我对比了一下Numpy与Matlab的eig结果，发现Numpy的精度和可靠性也远不如Matlab。因此，有两个问题想请教一下：

1. 有没有利用svd求一般实数本征值问题的简单代码实现？因为SVD分解得到的对角阵都是非负的，因此无法直接给出非半正定矩阵的本征值。当然可以通过U和V矩阵的相应矢量的正负关系来判断，但总觉得有些笨拙。所以想问下有没有更简单的实现方法？

2. 在Python中，有没有其它高效可靠地求本征值的方法（针对于一般矩阵，不一定是稀疏矩阵）？

作者
Author: snljty 时间: 2020-4-6 11:17
我记得numpy.linalg就是基于BLAS和LAPACK库啊，之前测试过的numpy.linalg.eigh的结果和LAPACK库的dsyev函数结果完全一样

作者
Author: liyuanhe211 时间: 2020-4-6 17:51
能不能给出一个例子？

作者
Author: granvia 时间: 2020-4-8 15:58

liyuanhe211 发表于 2020-4-6 17:51
能不能给出一个例子？

比如说，由于数值误差，某实对称阵A不完全对称，eig函数给出的本征值和本征矢量就会是复数。如果采用如下技巧，事先人为将A对称化：A = ( A + A.T ) / 2，则Numpy的eig()有时仍会给出复数解（尽管虚部很小）。但在Matlab却没上述问题。

作者
Author: granvia 时间: 2020-4-8 16:00

snljty 发表于 2020-4-6 11:17
我记得numpy.linalg就是基于BLAS和LAPACK库啊，之前测试过的numpy.linalg.eigh的结果和LAPACK库的dsyev函数 ...

我在网上找了下，好像是这样的。也有人说SciPy解本征值比NumPy更靠谱些，不知道大家有没有经验？

作者
Author: snljty 时间: 2020-4-8 17:02

granvia 发表于 2020-4-8 16:00
我在网上找了下，好像是这样的。也有人说SciPy解本征值比NumPy更靠谱些，不知道大家有没有经验？

如果已知是实对称矩阵，用eigh比eig更好。
我感觉至少scipy的linalg更方便(是一个子包，要import scipy.linalg as spla这样)，还可以直接求解广义特征值问题等。不像numpy还得手动做个正交化。

作者
Author: granvia 时间: 2020-4-8 17:15

snljty 发表于 2020-4-8 17:02
如果已知是实对称矩阵，用eigh比eig更好。
我感觉至少scipy的linalg更方便(是一个子包，要import scipy. ...

哦哦，多谢。我试试

作者
Author: granvia 时间: 2020-4-9 10:37

snljty 发表于 2020-4-8 17:02
如果已知是实对称矩阵，用eigh比eig更好。
我感觉至少scipy的linalg更方便(是一个子包，要import scipy. ...

确实，使用eigh()得到的Lowdin正交化结果与svd的一致，而与eig()的不同。

作者
Author: snljty 时间: 2020-4-9 11:32

granvia 发表于 2020-4-9 10:37
确实，使用eigh()得到的Lowdin正交化结果与svd的一致，而与eig()的不同。

eig好像更接近LAPACK的dgeev函数甚至zgeev之类的函数，我也记不清了，得翻翻Numpy，Scipy和LAPACK的手册。不过针对厄米矩阵，专门开发的函数不管从效率还是结果上应该都是更好的。

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