首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >酉矩阵的数值对角化

酉矩阵的数值对角化
EN

Stack Overflow用户
提问于 2019-01-29 21:59:50
回答 1查看 1.1K关注 0票数 3

要对正矩阵进行数值对角化,我使用LAPACK例程兹盖耶夫

问题是:在退化的情况下,退化子空间不是正规化的,因为例程是针对一般矩阵的。

但是,由于在我的例子中,矩阵是酉的,所以基总是可以正规化的。对于退化子空间,有比QR算法更好的解决方案吗?

EN

回答 1

Stack Overflow用户

回答已采纳

发布于 2019-01-30 21:38:05

简短回答:Schur分解!

如果方阵A是复的,则它的Schur分解是A,其中Z是酉的,T是上三角的。如果A碰巧是酉的,那么T也必须是酉的。由于T是幺正的和三角形的,所以它是对角的(这里有证据,或在那里),让我们考虑向量Z.e_i,其中e_i是正则基的向量。这些向量显然构成了一个正交基。此外,这些向量是矩阵A的特征向量。因此,酉矩阵Z的列是酉矩阵A 的特征向量,构成了正交基。

因此,计算酉矩阵的Schur分解等价于寻找其特征向量的正交基之一。

计算GE矩阵的特征值、Schur形式和可选的Schur向量矩阵。

生成的T也可以进行测试,以检查A是否为整体式。

下面是一段python代码来测试它,尽管use的scipy.linalg.schur使用了Lapack的zgees来进行Schur分解。我使用hpaulj代码生成随机酉矩阵,如python中随机正交矩阵的建立所示

代码语言:javascript
复制
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)
票数 3
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/54430257

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档