Geoffroy Geoffroy - 1 month ago 15
Python Question

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

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

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: http://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.