python - 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)

为了获得一阶、二阶、三阶...阶微分，我例如写道:

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 来包装我的方法导数()并将其用作 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 找到的值中，但事实并非如此。 (同上，sympy)

我尝试过使用不同的随机矩阵而不是 DM(参见选项 2)(对称、非对称、实数、虚数等...)，其大小 N (4,5,6..)还有更大的(100，...)。这有效

通过更改 scipy.sparse 中的“which”(LM、SM、LR...)、“tol”(10E-3、10E-6..)、“maxiter”、“sigma”(0)等参数... 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，它返回相同的值。

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

λ B x= DM x

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

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

编辑:

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