Различные собственные значения между scipy.sparse.linalg.eigs и numpy/scipy.eig
Контекст:
Моя цель - создать программу на Python3 для работы с дифференциальными операциями над вектором V размера N. Я сделал это, протестировал его на базовую операцию, и она работает (дифференцирование, градиент...).
Я попытался написать на этой основе более сложные уравнения (Навье-Стокса, Орра-Зоммерфельда,...), и я попытался подтвердить свою работу, рассчитав собственные значения этих уравнений.
Поскольку эти собственные значения были совершенно неожиданными, я упрощаю свою задачу и в настоящее время пытаюсь рассчитать собственные значения только для матрицы дифференцирования (см. Ниже). Но результаты кажутся неправильными...
Заранее спасибо за вашу помощь, потому что я не нахожу решение моей проблемы...
Определение ДМ:
Я использую чебышевский спектральный метод для управления дифференцированием векторов. Я использую следующий чебышевский пакет (перевод с Matlab на Python): http://dip.sun.ac.za/~weideman/research/differ.html
Этот пакет позволяет мне создать матрицу дифференцирования DM, полученную с помощью:
nodes, DM = chebyshev.chebdiff(N, maximal_order)
Чтобы получить 1-й, 2-й, 3-й... порядок дифференциации, я пишу, например:
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, что не так. (то же самое для симпати)
Я пробовал использовать разные случайные матрицы вместо DM (см. Вариант 2) (симметричные, несимметричные, действительные, мнимые и т. Д.), Которые имели небольшой размер N (4,5,6...) и также больше единицы (100,...). Это сработало
Изменяя такие параметры, как 'which' (LM, SM, LR...), 'tol' (10E-3, 10E-6..), 'maxiter', 'sigma' (0) в scipy.sparse... scipy.sparse.linalg.eigs всегда работал для случайных матриц, но никогда для моей матрицы DM. В лучших случаях найденные собственные значения близки к найденным Сципи, но никогда не совпадают.
Я действительно не знаю, что такого особенного в моей матрице. Я также не знаю, почему использование scipy.sparse.linagl.eig с матрицей, LinearOperator или AsLinearOperator дает разные результаты.
Я НЕ ЗНАЮ, КАК МОЖНО ВКЛЮЧИТЬ МОИ ФАЙЛЫ, СОДЕРЖАЩИЕ МАТРИЦ ДМ...
Для 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 / eigenvalues / scipy.eig /scipy.sparse.lingalg.linearOperator
Жоффруа.
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. Я подтвердлю это позже.
РЕДАКТИРОВАТЬ:
Я написал решение вопроса об обмене стека, который я публикую выше, в котором объясняется, как я решаю свою проблему. Мне кажется, что scipy.sparse.linalg.eigs действительно содержит ошибку, если матрица B не является положительно определенной и будет возвращать неправильные собственные значения.