2015-02-12 90 views
3

爲什麼如下面所使用的scipy.sparse.linalg中的eigheigsh在求解廣義特徵值問題時給出了不正確的結果A * x = lambda * M * x如果M是非對角線?不正確的特徵值SciPy稀疏linalg.eigs,非對角線M矩陣的eigsh

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矩陣約爲30,000 x 30,000。 A總是正方形,實數和對角線,M總是正方形,實數和對稱。當M是對角線時,我得到正確的結果。但是,當解決非對角M矩陣的廣義特徵值問題時,eigseigsh都給出不正確的結果。

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 -1 * A * X =拉姆達* 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]。

回答

5

兩個eigseigsh要求Mpositive definite(見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]) 
+0

什麼Python包可以解決廣義特徵值問題與稀疏不定對稱M矩陣?看來SciPy和PySparse不能做到這一點。 – bjvanbruchem 2015-02-13 18:54:14

+2

@bjvanbruchem您可以在scipy的'eigs [h]'函數中使用'sigma'參數來使用shift-invert方法,或者在對角線上添加一個常數以使特徵值爲正。 – rubenvb 2015-10-08 09:07:32

+0

@rubenvb這確實是解決這個問題的正確方法。我已經使用scipy eigsh解決了這個問題,其值爲sigma和mode ='buckling'。 – bjvanbruchem 2015-10-26 11:43:49