酉矩阵的数值对角化
numerical diagonalization of a unitary matrix
为了在数值上对角化 酉矩阵,我使用了 LAPACK 程序 zgeev。
问题是:在退化的情况下,退化子空间没有正交化,因为例程是针对一般矩阵的。
但是,由于在我的例子中矩阵是酉矩阵,基总是可以正交化。有没有比之后将 QR 算法应用于退化子空间更好的解决方案?
如果方阵 A
是复数,则其 Schur 分解为 A=ZTZ*
,其中 Z
是酉矩阵,T
是上三角。
如果 A
恰好是幺正的,那么 T
也必须是幺正的。由于 T
既是酉矩阵又是三角形,所以它是对角线 (proof here,.or there)
让我们考虑向量 Z.e_i
,其中 e_i 是规范基的向量。这些向量显然形成了正交基。此外,这些向量是矩阵 A
的特征向量。
因此,酉矩阵Z的列是酉矩阵A
的特征向量,构成正交基
因此,计算酉矩阵的 Schur 分解等同于找到其特征向量的正交基之一。
也可以测试结果 T
以检查 A
是否单一。
这是一段 python 测试它的代码,尽管 scipy 的 scipy.linalg.schur
使用 Lapack 的 zgees 进行 Schur 分解。我使用 hpaulj 的代码生成随机酉矩阵,如
所示
import numpy as np
import scipy.linalg
#from hpaulj,
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)
为了在数值上对角化 酉矩阵,我使用了 LAPACK 程序 zgeev。
问题是:在退化的情况下,退化子空间没有正交化,因为例程是针对一般矩阵的。
但是,由于在我的例子中矩阵是酉矩阵,基总是可以正交化。有没有比之后将 QR 算法应用于退化子空间更好的解决方案?
如果方阵 A
是复数,则其 Schur 分解为 A=ZTZ*
,其中 Z
是酉矩阵,T
是上三角。
如果 A
恰好是幺正的,那么 T
也必须是幺正的。由于 T
既是酉矩阵又是三角形,所以它是对角线 (proof here,.or there)
让我们考虑向量 Z.e_i
,其中 e_i 是规范基的向量。这些向量显然形成了正交基。此外,这些向量是矩阵 A
的特征向量。
因此,酉矩阵Z的列是酉矩阵A
的特征向量,构成正交基
因此,计算酉矩阵的 Schur 分解等同于找到其特征向量的正交基之一。
也可以测试结果 T
以检查 A
是否单一。
这是一段 python 测试它的代码,尽管 scipy 的 scipy.linalg.schur
使用 Lapack 的 zgees 进行 Schur 分解。我使用 hpaulj 的代码生成随机酉矩阵,如
import numpy as np
import scipy.linalg
#from hpaulj,
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)