scipy.sparse.linalg.eigs 和 numpy/scipy.eig 之间的不同特征值

Question

上下文：

我的目标是创建一个 Python3 程序来对大小为 N 的向量 V 进行微分运算。我这样做了，测试了它的基本运算并且它有效（微分，梯度......）。

我试图以此为基础编写更复杂的方程（Navier-Stokes、Orr-Sommerfeld 等），并试图通过计算这些方程的特征值来验证我的工作。

由于这些特征值完全出乎意料，我简化了我的问题，目前我正在尝试仅为微分矩阵计算特征值（见下文）。但是结果好像不对...

在此先感谢您的帮助，因为我没有找到解决我的问题的方法...

DM 定义：

我用切比雪夫谱法对向量进行微分运算。我使用以下 Chebyshev 包（从 Matlab 翻译成 Python）： http://dip.sun.ac.za/%7Eweideman/research/differ.html

该包允许我创建微分矩阵 DM，通过以下方式获得：

nodes, DM = chebyshev.chebdiff(N, maximal_order)

求1st,2nd,3th...阶差分，我写例如：

dVdx1 = np.dot(DM[0,:,:], V)
d2Vdx2 = np.dot(DM[1,:,:], V)
d3Vdx3 = np.dot(DM[2,:,:], V)

我测试过它并且有效。我根据该差异化过程构建了不同的运算符。我试图通过找到它们的特征值来验证它们。进展不顺利，所以我现在只尝试使用 DM。我没能找到 DM 的正确特征值。

我尝试过不同的功能：

numpy.linalg.eigvals(DM)
scipy.linalg.eig(DM)
scipy.sparse.linalg.eigs(DM)
sympy.solve( (DM-x*np.eye).det(), x) [for snall size only]

为什么我使用 scipy.sparse.LinearOperator:

我不想直接使用矩阵DM，所以我包装成一个函数来操作微分（见下面的代码）：

dVdx1 = derivative(V)

我这样做的原因来自全球项目本身。这对于更复杂的方程很有用。

创建这样的函数使我无法直接使用矩阵 DM 来查找其特征值（因为 DM 位于函数内部）。出于这个原因，我使用 scipy.sparse.LinearOperator 来包装我的方法 derivative() 并将其用作 scipy.sparse.eig().

的输入

代码和结果：

下面是计算这些特征值的代码：

import numpy as np
import scipy
import sympy

from scipy.sparse.linalg import aslinearoperator
from scipy.sparse.linalg import eigs
from scipy.sparse.linalg import LinearOperator

import chebyshev

N = 20 # should be 4, 20, 50, 100, 300
max_order = 4

option = 1
#option 1: building the differentiation matrix DM for a given order
if option == 1:
    if 0:
        # usage of package chebyshev, but I add a file with the matrix inside
        nodes, DM = chebyshev.chebdiff(N, max_order)
        order = 1
        DM = DM[order-1,:,:]
        #outfile = TemporaryFile()
        np.save('DM20', DM)
    if 1:
        # loading the matrix from the file
        # uncomment depending on N
        #DM = np.load('DM4.npy')
        DM = np.load('DM20.npy')
        #DM = np.load('DM50.npy')
        #DM = np.load('DM100.npy')
        #DM = np.load('DM300.npy')

#option 2: building a random matrix
elif option == 2:
    j = np.complex(0,1)
    np.random.seed(0)
    Real = np.random.random((N, N)) - 0.5
    Im = np.random.random((N,N)) - 0.5

    # If I want DM symmetric:
    #Real = np.dot(Real, Real.T)
    #Im = np.dot(Im, Im.T)

    DM = Real + j*Im

    # If I want DM singular:
    #DM[0,:] = DM[1,:]

# Test DM symmetric
print('Is DM symmetric ? \n', (DM.transpose() == DM).all() )        
# Test DM Hermitian
print('Is DM hermitian ? \n', (DM.transpose().real == DM.real).all() and
                                        (DM.transpose().imag == -DM.imag).all() )  

# building a linear operator which wrap matrix DM
def derivative(v):
    return np.dot(DM, v)

linop_DM = LinearOperator( (N, N), matvec = derivative)

# building a linear operator directly from a matrix DM with asLinearOperator
aslinop_DM = aslinearoperator(DM)

# comparison of LinearOperator and direct Dot Product
V = np.random.random((N))
diff_lo = linop_DM.matvec(V)
diff_mat = np.dot(DM, V)
# diff_lo and diff_mat are equals

# FINDING EIGENVALUES

#number of eigenvalues to find
k = 1
if 1:
    # SCIPY SPARSE LINALG LINEAR OPERATOR
    vals_sparse, vecs = scipy.sparse.linalg.eigs(linop_DM, k, which='SR',
                            maxiter = 10000,
                            tol = 1E-3)
    vals_sparse = np.sort(vals_sparse)
    print('\nEigenvalues (scipy.sparse.linalg Linear Operator) : \n', vals_sparse)

if 1:
    # SCIPY SPARSE ARRAY
    vals_sparse2, vecs2 = scipy.sparse.linalg.eigs(DM, k, which='SR',
                        maxiter = 10000,
                        tol = 1E-3)
    vals_sparse2 = np.sort(vals_sparse2)
    print('\nEigenvalues (scipy.sparse.linalg with matrix DM) : \n', vals_sparse2)

if 1:
    # SCIPY SPARSE AS LINEAR OPERATOR
    vals_sparse3, vecs3 = scipy.sparse.linalg.eigs(aslinop_DM, k, which='SR',
                        maxiter = 10000,
                        tol = 1E-3)
    vals_sparse3 = np.sort(vals_sparse3)
    print('\nEigenvalues (scipy.sparse.linalg AS linear Operator) : \n', vals_sparse3)

if 0:
    # NUMPY LINALG / SAME RESULT AS SCIPY LINALG
    vals_np = np.linalg.eigvals(DM)
    vals_np = np.sort(vals_np)
    print('\nEigenvalues (numpy.linalg) : \n', vals_np)

if 1:
    # SCIPY LINALG
    vals_sp = scipy.linalg.eig(DM)
    vals_sp = np.sort(vals_sp[0])
    print('\nEigenvalues (scipy.linalg.eig) : \n', vals_sp)

if 0:
    x = sympy.Symbol('x')
    D = sympy.Matrix(DM)
    print('\ndet D (sympy):', D.det() )
    E = D - x*np.eye(DM.shape[0])
    eig_sympy = sympy.solve(E.det(), x)
    print('\nEigenvalues (sympy) : \n', eig_sympy)

这是我的结果（N=20）：

Is DM symmetric ? 
 False
Is DM hermitian ? 
 False

Eigenvalues (scipy.sparse.linalg Linear Operator) : 
 [-2.5838015+0.j]

Eigenvalues (scipy.sparse.linalg with matrix DM) : 
 [-2.58059801+0.j]

Eigenvalues (scipy.sparse.linalg AS linear Operator) : 
 [-2.36137671+0.j]

Eigenvalues (scipy.linalg.eig) : 
 [-2.92933791+0.j         -2.72062839-1.01741142j -2.72062839+1.01741142j
 -2.15314244-1.84770128j -2.15314244+1.84770128j -1.36473659-2.38021351j
 -1.36473659+2.38021351j -0.49536645-2.59716913j -0.49536645+2.59716913j
  0.38136094-2.53335888j  0.38136094+2.53335888j  0.55256471-1.68108134j
  0.55256471+1.68108134j  1.26425751-2.25101241j  1.26425751+2.25101241j
  2.03390489-1.74122287j  2.03390489+1.74122287j  2.57770573-0.95982011j
  2.57770573+0.95982011j  2.77749810+0.j        ]

scipy.sparse 返回的值应该包含在scipy/numpy 找到的值中，事实并非如此。（同情同上）

我尝试使用不同的随机矩阵代替 DM（见选项 2）（对称、非对称、实数、虚数等...），它的大小 N (4,5,6. .) 还有更大的 (100,...)。有效

通过更改参数，例如 'which'（LM、SM、LR...）、'tol'（10E-3、10E-6..）、'maxiter'、'sigma' (0) in scipy.sparse... scipy.sparse.linalg.eigs 始终适用于随机矩阵，但不适用于我的矩阵 DM。在最好的情况下，找到的特征值接近 scipy 找到的特征值，但永远不会匹配。

我真的不知道我的矩阵有什么特别之处。我也不知道为什么将 scipy.sparse.linagl.eig 与矩阵、LinearOperator 或 AsLinearOperator 一起使用会给出不同的结果。

我不知道如何包含我的包含矩阵 DM 的文件...

对于 N = 4 :

[[ 3.16666667 -4.          1.33333333 -0.5       ]
 [ 1.         -0.33333333 -1.          0.33333333]
 [-0.33333333  1.          0.33333333 -1.        ]
 [ 0.5        -1.33333333  4.         -3.16666667]]

欢迎任何想法。

版主可以用以下标记我的问题： scipy.sparse.linalg.eigs/weideman/特征值/scipy.eig/scipy.sparse.lingalg.linearOperator

若弗鲁瓦

Answer 1

我和几个同事谈过，部分解决了我的问题。我的结论是我的矩阵条件很差...

在我的项目中，我可以通过如下施加边界条件来简化我的矩阵：

DM[0,:] = 0
DM[:,0] = 0
DM[N-1,:] = 0
DM[:,N-1] = 0

生成类似于 N=4 的矩阵：

[[ 0     0               0               0]
 [ 0     -0.33333333     -1.             0]
 [ 0      1.             0.33333333      0]
 [ 0      0              0               0]]

通过使用这样的条件，我获得了 scipy.sparse.linalg.eigs 的特征值，它等于 scipy.linalg.eig 中的特征值。我也尝试使用 Matlab，它 return 具有相同的值。

为了继续我的工作，我实际上需要使用标准形式的广义特征值问题

λ B x= DM x

由于我的矩阵 B（表示拉普拉斯算子矩阵），它似乎不适用于我的情况。如果您有类似的问题，我建议您访问该问题： https://scicomp.stackexchange.com/questions/10940/solving-a-generalised-eigenvalue-problem

（我认为）矩阵 B 需要是正定的才能使用 scipy.sparse。一种解决方案是更改 B，使用 scipy.linalg.eig 或使用 Matlab。我稍后会确认。

编辑：

我在上面 post 写了一个堆栈交换问题的解决方案，它解释了我是如何解决我的问题的。如果矩阵 B 不是正定的，我发现 scipy.sparse.linalg.eigs 确实有一个错误，并且 return 的特征值会很糟糕。

scipy.sparse.linalg.eigs 和 numpy/scipy.eig 之间的不同特征值

Different eigenvalues between scipy.sparse.linalg.eigs and numpy/scipy.eig

python

scipy

eigenvalue

sparse-matrix

上下文：

DM 定义：

为什么我使用 scipy.sparse.LinearOperator:

代码和结果：

编辑：