Может ли numpy диагонализировать кососимметричную матрицу с помощью реальной арифметики?

Question

Может ли numpy диагонализировать кососимметричную матрицу с помощью реальной арифметики?

Любая кососимметричная матрица (А^Т = - А) можно превратить в Эрмитову матрицу (ИА) и диагонализирован комплексными числами. Но также возможнопривести его в блок-диагональную форму с помощью специального ортогонального преобразованияи найти его собственные значения, используя только реальную арифметику. Это реализовано где-нибудь в numpy?

2 3

python numpy matrix linear-algebra lapack

2 ответа:

Да, вы можете сделать это, вставив унитарное преобразование в середину произведения, следовательно, мы получаем

A = V * U * V^-1 = V * T' * T * U * T' * T * V^{-1}.

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

Если матрица четного размера, то все блоки являются комплексными сопряжениями. В противном случае мы получим ноль в качестве собственного значения. Собственные значения гарантированно есть нулевые вещественные части, поэтому первое, что нужно сделать, это убрать шум, а затем упорядочить так, чтобы нули находились в верхнем левом углу (произвольный выбор).
n = 5
a = np.random.rand(n,n)
a=a-np.transpose(a)
[u,v] = np.linalg.eig(a)

perm = np.argsort(np.abs(np.imag(u)))
unew = 1j*np.imag(u[perm])
Очевидно, нам также нужно переупорядочить матрицу собственных векторов,чтобы сохранить эквивалентность.
vnew = v[:,perm]
До сих пор мы не делали ничего, кроме переупорядочивания средней матрицы собственных значений в разложении собственных значений. Теперь мы переходим от сложной формы к реальной блочной диагональной форме.
Сначала мы должны знать, сколько нулей собственные значения существуют
numblocks = np.flatnonzero(unew).size // 2
num_zeros = n - (2 * numblocks)
Затем мы в основном формируем еще одно унитарное преобразование (на этот раз сложное) и придерживаемся его таким же образом
T = sp.linalg.block_diag(*[1.]*num_zeros,np.kron(1/np.sqrt(2)*np.eye(numblocks),np.array([[1.,1j],[1,-1j]])))

Eigs = np.real(T.conj().T.dot(np.diag(unew).dot(T)))
Evecs = np.real(vnew.dot(T))
Это дает вам новую действительную ценностную декомпозицию. Итак, код все в одном месте
n = 5
a = np.random.rand(n,n)
a=a-np.transpose(a)
[u,v] = np.linalg.eig(a)

perm = np.argsort(np.abs(np.imag(u)))
unew = 1j*np.imag(u[perm])
vnew = v[perm,:]
numblocks = np.flatnonzero(unew).size // 2
num_zeros = n - (2 * numblocks)
T = sp.linalg.block_diag(*[1.]*num_zeros,np.kron(1/np.sqrt(2)*np.eye(numblocks),np.array([[1.,1j],[1,-1j]])))
Eigs = np.real(T.conj().T.dot(np.diag(unew).dot(T)))
Evecs = np.real(vnew.dot(T))
print(np.allclose(Evecs.dot(Eigs.dot(np.linalg.inv(Evecs))) - a,np.zeros((n,n))))
Дает True. Заметим, что это наивный способ получения реального спектрального разложения. Есть много мест, где нужно следить за накоплением числовых ошибок.

Пример вывода
Eigs
Out[379]: 
array([[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        , -0.61882847,  0.        ,  0.        ],
       [ 0.        ,  0.61882847,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        , -1.05097581],
       [ 0.        ,  0.        ,  0.        ,  1.05097581,  0.        ]])

Evecs
Out[380]: 
array([[-0.15419078, -0.27710323, -0.39594838,  0.05427001, -0.51566173],
       [-0.22985364,  0.0834649 ,  0.23147553, -0.085043  , -0.74279915],
       [ 0.63465436,  0.49265672,  0.        ,  0.20226271, -0.38686576],
       [-0.02610706,  0.60684296, -0.17832525,  0.23822511,  0.18076858],
       [-0.14115513, -0.23511356,  0.08856671,  0.94454277,  0.        ]])

1

francis · Accepted Answer · 2016-08-20 18:59:15

Давайте взглянем на dgeev() функция библиотеки ЛАПАКА. Эта процедура вычисляет собственные значения любой реальной квадратной матрицы двойной точности. Более того, эта процедура находится прямо за функцией python numpy.linalg.eigvals() из библиотеки numpy.

Метод, используемый dgeev(), описан в документацииLAPACK . Это требует приведения матрицы A к ее вещественной форме Шура S.

Любая реальная квадратная матрица A может быть выражается следующим образом:

A=QSQ^t

Где:

Q - вещественная ортогональная матрица: QQ^t=I

S является реальной блочной верхней треугольной матрицей. Блоки на диагонали S имеют размер 1×1 или 2×2.

Действительно, если A кососимметрично, то это разложение кажется действительно близким к диагональной форме блока , полученной специальным ортогональным преобразованием из A. Более того, действительно видеть, что форма Шура S косого симметричная матрица A есть ... косо-симметрично !
Действительно, давайте вычислим транспозицию S:
S^t=(Q^tAQ)^t
S^t=Q^t(Q^tA)^t
S^t=Q^tA^tQ
S^t=Q^t(-A)Q
S^t=-Q^tAQ
S^t=-S
Следовательно, если Q является специальным ортогональным (det(Q)=1), S является блочной диагональной формой, полученной специальным ортогональным преобразованием. Кроме того, специальная ортогональная матрица P может быть вычислена путем перестановки первых двух столбцов Q, а другая форма Шура Sd матрицы A получается путем изменения знака S_{12} и S_{21}. Действительно, A=PSdP^t. Тогда Sd является блочная диагональная форма A, полученная специальным ортогональным преобразованием.

В конце концов, даже если numpy.linalg.eigvals() применительно к реальной матрице, возвращающей комплексные числа ,в этом процессе мало сложных вычислений!

Если вы просто хотите вычислить реальную форму Шура, используйте функцию scipy.linalg.schur() с аргументом output='real'.

Просто кусок кода, чтобы проверить это:
import numpy as np
import scipy.linalg as la

a=np.random.rand(4,4)
a=a-np.transpose(a)

print "a= "
print a

#eigenvalue
w, v =np.linalg.eig(a)

print "eigenvalue "
print w
print "eigenvector "
print v

# Schur decomposition
#import scipy
#print scipy.version.version

t,z=la.schur(a, output='real', lwork=None, overwrite_a=True, sort=None, check_finite=True)

print "schur form "
print t
print "orthogonal matrix "
print z