Numpy的eigh和eig产生了不一致的特征值

Question

目前我正试图在NumPy中解决两个对称矩阵的广义特征值问题，我遇到了很大的麻烦，因为我期望所有的特征值都是正的，但是eigh返回了几个非常大的数字，而eig返回正确的期望值（但是，当然，非常，非常慢）。在

在本例中，请注意K与构造中预期的对称（下面是有问题的代码）：

# Calculate K matrix (<i|pHp|j> in the LGL-nodes basis)
for i in range(Ne):
    idx_s, idx_e = i*(Np-1), i*(Np-1)+Np
    K[idx_s:idx_e, idx_s:idx_e] += dmat.T.dot(diag(w*peq[idx_s:idx_e])).dot(dmat)

# Re-make matrix for efficient vector products
K = sparse.csr_matrix(K)

# Make matrix for <i|p|j> in the LGL basis as efficient diagonal sparse matrix
S = sparse.diags(peq*w_d, 0)

# Solve the generalized eigenvalue problem: Kc = lSc for hermitian matrices K and S
lQ, Q = linalg.eigh(K.todense(), S.todense())
_lQ, _Q = linalg.eig(K.todense(), S.todense())

lQ.sort()
_lQ.sort()

if not allclose(lQ, _lQ):
    print('Literally why')
    print(lQ)
    print(_lQ)
return

对于测试，dmat被定义为

^{pr2}$

所有的w[i]，w_d[i]，peq[i]本质上都是任意的正值数组。w_d和w顺序相同（~1e-1）和{}范围，顺序为（~1e-10至1e1）

我得到的一些结果是

Literally why
[ -6.25540943e+07  -4.82660391e+07  -2.62629052e+07 ...,   1.07960873e+10
   1.07967334e+10   4.26007915e+10]
[ -5.25462340e-12+0.j   4.62614812e-01+0.j   1.23357898e+00+0.j ...,
   2.17613917e+06+0.j   1.07967334e+10+0.j   4.26007915e+10+0.j]

编辑：

下面是一个自包含的代码版本，便于调试

import numpy as np
from math import *
from scipy import sparse, linalg

# Variable declarations and such (pre-computed)
Ne, Np = 256, 8
N = Ne*Np - Ne + 1
domain_size = 4/Ne


x = np.array([-0.015625  , -0.01362094, -0.00924532, -0.0032703 ,  0.0032703 ,
    0.00924532,  0.01362094,  0.015625  ])
w = np.array([ 0.00055804,  0.00329225,  0.00533004,  0.00644467,  0.00644467,
    0.00533004,  0.00329225,  0.00055804])
dmat = np.array([[ -896.        ,  1212.00631086,  -484.43454844,   275.06612251,
     -179.85209531,   124.26620323,   -83.05199285,    32.        ],
   [ -205.43460499,     0.        ,   290.78944413,  -135.17191772,
       82.83085126,   -55.64467829,    36.70818656,   -14.07728095],
   [   50.7185076 ,  -179.61445086,     0.        ,   184.03311398,
      -87.85829324,    54.08144362,   -34.37053351,    13.01021241],
   [  -23.81762789,    69.05246008,  -152.20398294,     0.        ,
      152.89115899,   -72.66291308,    42.31407046,   -15.57316561],
   [   15.57316561,   -42.31407046,    72.66291308,  -152.89115899,
        0.        ,   152.20398294,   -69.05246008,    23.81762789],
   [  -13.01021241,    34.37053351,   -54.08144362,    87.85829324,
     -184.03311398,     0.        ,   179.61445086,   -50.7185076 ],
   [   14.07728095,   -36.70818656,    55.64467829,   -82.83085126,
      135.17191772,  -290.78944413,     0.        ,   205.43460499],
   [  -32.        ,    83.05199285,  -124.26620323,   179.85209531,
     -275.06612251,   484.43454844, -1212.00631086,   896.        ]])

# More declarations
x_d = np.zeros(N)
w_d = np.zeros(N)
dmat_d = np.zeros((N, N))

for i in range(Ne):
    x_d[i*(Np-1):i*(Np-1)+Np] = x+i*domain_size
    w_d[i*(Np-1):i*(Np-1)+Np] += w
    dmat_d[i*(Np-1):i*(Np-1)+Np, i*(Np-1):i*(Np-1)+Np] += dmat

peq = (np.cos((x_d-2)*pi/4))**2

# Normalization
peq = peq/np.sum(w_d*peq)
p0 = np.maximum(peq, 1e-10)
p0 /= np.sum(p0*w_d)

# Make efficient matrix that can be built
K = sparse.lil_matrix((N, N))

# Calculate K matrix (<i|pHp|j> in the LGL-nodes basis)
for i in range(Ne):
    idx_s, idx_e = i*(Np-1), i*(Np-1)+Np
    K[idx_s:idx_e, idx_s:idx_e] += dmat.T.dot(np.diag(w*p0[idx_s:idx_e])).dot(dmat)

# Re-make matrix for efficient vector products
K = sparse.csr_matrix(K)

# Make matrix for <i|p|j> in the LGL basis as efficient diagonal sparse matrix
S = sparse.diags(p0*w_d, 0)

# Solve the generalized eigenvalue problem: Kc = lSc for hermitian matrices K and S
lQ, Q = linalg.eigh(K.todense(), S.todense())
_lQ, _Q = linalg.eig(K.todense(), S.todense())

lQ.sort()
_lQ.sort()

if not np.allclose(lQ, _lQ):
    print('Literally why')
    print(lQ)
    print(_lQ)

这真的很奇怪。在我的机器上运行所有NumPy/SciPy测试时，我没有收到任何错误。但即使运行简单的测试（使用足够大的矩阵）也可以

import numpy as np
from spicy import linalg

M = np.random.random((1000,1000))
M += M.T

np.allclose(sorted(linalg.eigh(M)[0]), sorted(linalg.eig(M)[0]))

在我的机器上失败了。尽管使用50x50矩阵运行相同的测试仍然有效——即使在重新构建SciPy/NumPy堆栈并通过所有单元测试之后也是如此。

实际上，在集群计算机上进行测试后，这似乎处处失败。我不知道为什么。

由于+=和.T作为视图而不是操作的就地行为，上述操作失败。在