非对角M矩阵的不正确特征值SciPy稀疏线性

Question

如果M是非对角的，为什么下面使用的scipy.sparse.linarg中的eigh和eigsh在求解广义特征值问题A*x=lambda*M*x时给出不正确的结果？

import mkl
import numpy as np
from scipy import linalg as LA
from scipy.sparse import linalg as LAsp
from scipy.sparse import csr_matrix
A = np.diag(np.arange(1.0,7.0))
M = np.array([[ 25.1,   0. ,   0. ,  17.3,   0. ,   0. ],
       [  0. ,  33.6,  16.8,   8.4,   4.2,   2.1],
       [  0. ,  16.8,   3.6,   0. ,  11. ,   0. ],
       [ 17.3,   8.4,   0. ,   4.2,   0. ,   9.5],
       [  0. ,   4.2,  11. ,   0. ,   2.7,   8.3],
       [  0. ,   2.1,   0. ,   9.5,   8.3,   4.4]])
Asp = csr_matrix(np.matrix(A,dtype=float))
Msp = csr_matrix(np.matrix(M,dtype=float))
D, V = LA.eig(A, b=M)
eigno  = 4
Dsp0, Vsp0 = LAsp.eigs(csr_matrix(np.matrix(np.dot(np.linalg.inv(M),A))),
                         k=eigno,which='LM',return_eigenvectors=True)
Dsp1, Vsp1 = LAsp.eigs(Asp,k=eigno,M=Msp,which='LM',return_eigenvectors=True)
Dsp2, Vsp2 = LAsp.eigsh(Asp,k=eigno,M=Msp,which='LA',return_eigenvectors=True,
                          maxiter=1000)

从LA.eig和MatLab检查，这个具有测试矩阵A和M的小广义特征值问题的特征值应该是：

D = [ 0.7208+0.j,  0.3979+0.j, -0.3011+0.j, -0.3251+0.j,  0.0357+0.j,  0.0502+0.j]

我想使用稀疏矩阵，因为实际涉及的A和M矩阵大约是30000 x 30000。A总是正方形的，实的和对角线的，M总是正方形的、实的和对称的。当M是对角线时，我得到了正确的结果。然而，在求解非对角M矩阵的广义特征值问题时，eigs和eigsh都给出了不正确的结果。

Dsp1 = [-1.6526+2.3357j, -1.6526-2.3357j, -0.6243+2.7334j, -0.6243-2.7334j]
Dsp2 = [ 2.01019097,  3.09248265,  4.06799498,  7.01216316]

当我把问题转化为标准特征值形式M^-1*A*x=lambda*x时，eigs给出了正确的结果（Dsp0）。对于大矩阵，这不是一个选项，因为计算M.的逆需要太长时间

我注意到，使用或不使用mkl也会产生不同的Dsp1和Dsp2本征值。这个特征值问题可能是由我的Python安装问题引起的吗？我在Mac OS 10.10.2上运行Python 2.7.8 anaconda和SciPy 0.15.1-np19py27_p0[mkl]。

Answer 1

eigs和eigsh都要求M是正定的（有关更多详细信息，请参阅文档字符串中对M的描述）。

你的矩阵M不是正定的。注意负特征值：

In [212]: M
Out[212]: 
array([[ 25.1,   0. ,   0. ,  17.3,   0. ,   0. ],
       [  0. ,  33.6,  16.8,   8.4,   4.2,   2.1],
       [  0. ,  16.8,   3.6,   0. ,  11. ,   0. ],
       [ 17.3,   8.4,   0. ,   4.2,   0. ,   9.5],
       [  0. ,   4.2,  11. ,   0. ,   2.7,   8.3],
       [  0. ,   2.1,   0. ,   9.5,   8.3,   4.4]])
In [213]: np.linalg.eigvals(M)
Out[213]: 
array([ 45.92443169,  33.92113421, -13.12639751, -10.6991868 ,
         5.34183619,  12.23818222])