Valores propios incorrectos SciPy linalg.eigs disperso, eigsh para matriz M no diagonal

Por quéeigh yeigsh de scipy.sparse.linalg como se usa a continuación da resultados incorrectos al resolver el problema de valor propio generalizado A * x = lambda * M * x, si M no es diagonal?

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)

Desde LA.eig y comprobando con MatLab, los valores propios de este pequeño problema de valor propio generalizado con las matrices de prueba A y M deberían ser:

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]

Quiero usar matrices dispersas porque las matrices A y M reales involucradas son alrededor de 30,000 x 30,000. A es siempre cuadrado, real y diagonal, M siempre es cuadrado, real y simétrico. Cuando M es diagonal obtengo los resultados correctos. Sin embargo, amboseigs yeigsh dar resultados incorrectos al resolver el problema del valor propio generalizado para una matriz M no diagonal.

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]

Cuando convierto el problema a la forma de valor propio estándar M ^ -1 * A * x = lambda * x,eigs da el resultado correcto (Dsp0). Para matrices grandes, esta no es una opción porque lleva demasiado tiempo calcular el inverso de M.

Me di cuenta de que usandomkl o no produce diferentes valores propios de Dsp1 y Dsp2 también. ¿Podría este problema de valor propio ser causado por un problema con mi instalación de Python? Estoy ejecutando Python 2.7.8 anaconda con SciPy 0.15.1 - np19py27_p0 [mkl] en Mac OS 10.10.2.

Respuestas a la pregunta(1)

Su respuesta a la pregunta