Как получить реальные собственные значения и собственные векторы для трехдиагональной матрицы Тёплица?

Question

Я строю матрицу 100*100 k и хочу использовать numpy.linalg.eig диагонализировать это.

k=np.zeros((100,100))
np.fill_diagonal(k,-2)
np.fill_diagonal(k[1:,:-1],1.5)
np.fill_diagonal(k[:-1,1:],0.5)

когда я пытаюсь уменьшить матрицу, такую ​​как

w,v=np.linalg.eig(k[:10,:10])

собственное значение w и собственный вектор v реальны. Но когда я пробую большую матрицу или всю полную матрицу

w,v=np.linalg.eig(k)

w а также v оказываются комплексными числами, и мнимая часть не пренебрежимо мала.

Я тоже стараюсь scipy.linalg.eigУ него есть похожая проблема.

Я хочу взять натуральный логарифм собственных значений и собственных векторов. В моей модели нет физического значения комплексных чисел.

Как я могу иметь только собственные действительные числа и собственные векторы? Если нет, то как изменить сложные собственные значения и собственные векторы на реальные с помощью python?

Answer 1

@Daniel F и @FTP были быстрее меня, смотрите их комментарии, но так как у меня здесь есть код, я могу поделиться им:

import numpy as np
from scipy import sparse

def tri_toep_eig(a, b, c, n):
    evals = a + 2*np.sqrt(b*c) * np.cos(np.pi * np.arange(1, n+1) / (n+1))
    evecs = np.sin(np.outer(np.arange(1, n+1) * np.pi / (n+1),
                            np.arange(1, n+1))) \
        * np.sqrt(b/c)**np.arange(n)[:, None]
    return evals, evecs

def tri_toep(a, b, c, n):
    return sparse.dia_matrix((np.outer((b, a, c), np.ones((n,))),
                              (-1, 0, 1)), (n, n))
def check(a, b, c, n):
    evals, evecs = tri_toep_eig(a, b, c, n)
    tt = tri_toep(a, b, c, n)
    for eva, eve in zip(evals, evecs.T):
        assert np.allclose(tt @ eve, eva * eve)

check(-2, 0.5, 1.5, 100)

Answer 2

Транспонировать исправляет собственные значения

По-видимому, LAPACK ненавидит несимметричные трехдиагональные матрицы, где большие недиагональные элементы находятся ниже диагонали. Использование транспонированной матрицы, в которой более крупные элементы находятся выше диагонали, приводит к реальным собственным значениям. (Теоретически, матрица и ее транспонирование имеют одинаковые собственные значения.) Помимо того, что они действительны, они согласуются с теоретическими значениями вплоть до сортировки и допустимых ошибок.

a, b, c, n = -2, 0.5, 1.5, 100
k = np.zeros((n, n))
np.fill_diagonal(k, a)
np.fill_diagonal(k[:-1, 1:], b)
np.fill_diagonal(k[1:, :-1], c)

theory = a - 2*np.sqrt(b*c) * np.cos(np.pi * np.arange(1, n+1) / (n+1))
computed = np.sort(np.linalg.eig(k.T)[0])
print(np.max(np.abs(theory - computed)))

Это печатает 6.183001044490766e-08наибольшее расхождение между вычисленным и теоретическим собственным значением. Без транспонирования T эта ошибка увеличивается до 0,26.

Собственные

Вы также хотели собственные векторы. Собственные векторы возвращены np.eig для транспонированной матрицы - левые собственные векторы исходной матрицы, то есть они удовлетворяют vec.dot(k) = lam*vec вместо k.dot(vec) = lam*vec, Если вы хотите получить правильные собственные векторы для исходной матрицы, используйте SciPy's eig:

from scipy import linalg as la
evals, right_evects, left_evects = (np.real(_) for _ in la.eig(k.T, left=True, right=True))

Eigensolver SciPy отличается от NumPy тем, что он возвращает собственные векторы и собственные значения с +0j прилагается; он считает их сложными, но правильно оценивает мнимую часть до 0. Я урезал эту мнимую часть выше. Обратите внимание, что порядок возврата SciPy "evals, left, right", но так как k был транспонирован, я переключился влево и вправо.

Давайте проверим эти собственные векторы:

np.max([np.linalg.norm(k.dot(vec) - ev*vec) for ev, vec in zip(evals, right_evects.T)])

возвращается 1.845213984555825e-14, неплохо. Транспонирование на массив с собственными векторами, потому что zip выбирает строки из матрицы, и нам нужны столбцы.

Бонусный контент

Итак... проблема решена? Ну, я не говорил, что могу диагонализировать вашу матрицу. Попытка инвертировать матрицу, образованную левым или правым собственным вектором, выглядит проигрышным предложением; обратное ужасно.

Кроме того, мы не должны слишком доверять вышеуказанному тесту для собственных векторов. Это дало небольшую ошибку для правильных собственных векторов... давайте попробуем это на неправильных собственных векторах, которые имеют нетривиальные мнимые части.

wrong_evals, wrong_evects = np.linalg.eig(k)
np.max([np.linalg.norm(k.dot(vec) - ev*vec) for ev, vec in zip(wrong_evals, wrong_evects.T)])

Это возвращает 1.7136373586499598e-14, Неправильные собственные векторы даже лучше, чем настоящие!