酉矩阵的数值对角化

Question

为了对 酉矩阵进行数值对角化，我使用 LAPACK 例程 zgeev。

问题是：在简并的情况下，简并子空间没有正交归一化，因为该例程适用于一般矩阵。

但是，由于在我的例子中矩阵是单一的，基总是可以正交化。有没有比之后对退化子空间应用 QR 算法更好的解决方案？

Answer 1

简答：Schur decomposition!

如果方阵A 是复数，则其舒尔因式分解为A=ZTZ*，其中Z 是酉矩阵，T 是上三角矩阵。 如果A 恰好是一元的，那么T 也必须是一元的。由于T既是单位又是三角形，所以它是对角线(proof here,.or there) 让我们考虑向量Z.e_i，其中 e_i 是规范基的向量。这些向量显然形成了一个标准正交基。此外，这些向量是矩阵A 的特征向量。 因此，酉矩阵 Z 的列是酉矩阵A 的特征向量，形成一个正交基。

因此，计算酉矩阵的 Schur 分解相当于找到其特征向量的正交基之一。

ZGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices

还可以测试生成的T 以检查A 是否是单一的。

这是一段测试它的 python 代码，尽管 scipy 的 scipy.linalg.schur 使用 Lapack 的 zgees 进行 Schur 分解。我使用 hpaulj 的代码生成随机酉矩阵如How to create random orthonormal matrix in 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)