численная диагонализация унитарной матрицы

Question

Для числовой диагонали унитарной матрицы я использую процедуру LAPACK zgeev .

Проблема в том, что в случае вырождений вырожденное подпространство не ортонормируется, поскольку процедура предназначена для общих матриц.

Однако, поскольку в моем случае матрицы унитарные, базис всегда можно ортонормировать .Есть ли лучшее решение, чем применение QR-алгоритма к вырожденному подпространству?

Answer 1

Краткий ответ: Schur decomposition!

Если квадратная матрица A является сложной, то ее факторизация Шура равна A=ZTZ*, где Z является унитарным, а T - верхней треугольной.Если A окажется унитарным, T также должен быть унитарным.Поскольку T является унитарным и треугольным, он имеет диагональ (доказательство здесь, . или там) Давайте рассмотрим векторы Z.e_i, где e_i - векторы канонического базиса,Эти векторы, очевидно, образуют ортонормированный базис.Более того, эти векторы являются собственными векторами матрицы A. Следовательно, столбцы унитарной матрицы Z являются собственными векторами унитарной матрицы A и образуют ортонормированный базис.

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

ZGEESX вычисляет собственные значения, форму Шура и, необязательно, матрицу векторов Шура для матриц GE

Полученный T также может быть проверен, чтобы проверить, что A является унитарным.

Вот фрагмент кода Python, проверяющего его, хотя scipy.linalg.schur scipy использует zgees Лапака для разложения Шура.Я использовал код hpaulj для генерации случайной унитарной матрицы, как показано в Как создать случайную ортонормированную матрицу в python numpy

import numpy as np
import scipy.linalg

#from hpaulj, https://stackoverflow.com/questions/38426349/how-to-create-random-orthonormal-matrix-in-python-numpy
def rvs(dim=3):
     random_state = np.random
     H = np.eye(dim)
     D = np.ones((dim,))
     for n in range(1, dim):
         x = random_state.normal(size=(dim-n+1,))
         D[n-1] = np.sign(x[0])
         x[0] -= D[n-1]*np.sqrt((x*x).sum())
         # Householder transformation
         Hx = (np.eye(dim-n+1) - 2.*np.outer(x, x)/(x*x).sum())
         mat = np.eye(dim)
         mat[n-1:, n-1:] = Hx
         H = np.dot(H, mat)
         # Fix the last sign such that the determinant is 1
     D[-1] = (-1)**(1-(dim % 2))*D.prod()
     # Equivalent to np.dot(np.diag(D), H) but faster, apparently
     H = (D*H.T).T
     return H

n=42
A= rvs(n)
A = A.astype(complex)
T,Z=scipy.linalg.schur(A,output='complex',lwork=None,overwrite_a=False,sort=None,check_finite=True)

#print T
normT=np.linalg.norm(T,ord=None) #2-norm
eigenvalues=[]
for i in range(n):
    eigenvalues.append(T[i,i])
    T[i,i]=0.
normTu=np.linalg.norm(T,ord=None)
print 'must be very low if A is unitary: ',normTu/normT

#print Z
for i in range(n):
    v=Z[:,i]
    w=A.dot(v)-eigenvalues[i]*v
    print i,'must be very low if column i of Z is eigenvector of A: ',np.linalg.norm(w,ord=None)/np.linalg.norm(v,ord=None)