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

Question

Context:

My goal is to create a Python3 program to operate differential operations on a vector V of size N. I did so, test it for basic operation and it works (differentiation, gradient...).

I tried to write with that basis more complex equations (Navier-Stokes, Orr-Sommerfeld,...) and I tried to validate my work by calculating the eigenvalues of these equations.

As these eigenvalues were completely unexpected, I simplify my problem and I am currently trying to calculate the eigenvalues only for the differentiation matrix (see below). But the results seem wrong...

Thanks in advance for your help, because I do not find any solution to my problem...

Definition of DM:

I use Chebyshev spectral method to operate the differentiation of vectors. I use the following Chebyshev package (translated from Matlab to Python): http://dip.sun.ac.za/%7Eweideman/research/differ.html

That package allow me to create a differentiation matrix DM, obtained with:

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

To obtain the 1st, 2nd, 3th... order differentiation, I write for example:

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

I tested that and it works. I've build different operators based on that differentiation process. I've tried to validate them by finding their eigenvalues. It didn't go well so I am just trying right now with DM only. I do not manage to find the right eigenvalues of DM.

I've tried with different functions:

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]

Why I use scipy.sparse.LinearOperator:

I do not want to directly use the matrix DM, so I wrapped into a function which operates the differentiation (see code below) like that:

dVdx1 = derivative(V)

The reason why I do that comes from the global project itself. This is useful for more complicated equations.

Creating such a function prevents me from using directly the matrix DM to find its eigenvalues (because DM stay inside the function). For that reason, I use a scipy.sparse.LinearOperator to wrap my method derivative() and use it as an input of scipy.sparse.eig().

Code and results:

Here is the code to compute these eigenvalues:

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)

Here are my results (for 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        ]

The values returned by scipy.sparse should be included in the ones found by scipy/numpy, which is not the case. (idem for sympy)

I've tried with different random matrices instead of DM (see option 2) (symmetric, non-symmetric, real, imaginary, etc...), which had small size N (4,5,6..) and also bigger ones (100,...). That worked

By changing parameters like 'which' (LM, SM, LR...), 'tol' (10E-3, 10E-6..), 'maxiter', 'sigma' (0) in scipy.sparse... scipy.sparse.linalg.eigs always worked for random matrices but never for my matrix DM. In best cases, found eigenvalues are close to the ones found by scipy, but never match.

I really do not know what is so particular in my matrix. I also dont know why using scipy.sparse.linagl.eig with a matrix, a LinearOperator or a AsLinearOperator gives different results.

I DO NOT KNOW HOW I COULD INCLUDE MY FILES CONTAINING MATRICES DM...

For 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]]

Every idea is welcome.

May a moderator could tag my question with : scipy.sparse.linalg.eigs / weideman / eigenvalues / scipy.eig /scipy.sparse.lingalg.linearOperator

Geoffroy.

Answer 1

I spoke with a few colleague and solve partly my problem. My conclusion is that my matrix is simply very ill conditioned...

In my project, I can simplify my matrix by imposing boundary condition as follow:

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

which produces a matrix similar to that for N=4:

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

By using such condition, I obtain eigenvalues for scipy.sparse.linalg.eigs which are equal to the one in scipy.linalg.eig. I also tried using Matlab, and it return the same values.

To continue my work, I actually needed to use the generalized eigenvalue problem in the standard form

λ B x= DM x

It seems that it does not work in my case because of my matrix B (which represents a Laplacian operator matrix). If you have a similar problem, I advise you to visit that question: https://scicomp.stackexchange.com/questions/10940/solving-a-generalised-eigenvalue-problem

(I think that) the matrix B needs to be positive definite to use scipy.sparse. A solution would be to change B, to use scipy.linalg.eig or to use Matlab. I will confirm that later.

EDIT:

I wrote a solution to the stack exchange question I post above which explains how I solve my problem. I appears that scipy.sparse.linalg.eigs has indeed a bug if matrix B is not positive definite, and will return bad eigenvalues.