%实矩阵的拟上三角分解(Hessenberg分解) %用法:[Q,B]=hess2(A) %Q返回一个正交矩阵,B为相似于A的拟上三角矩阵 %注意:MATLAB自带了Hessenberg分解的函数hess(A) %By Castor 2010-12-25 function [Q,B]=hess2(A) tic; n=max(size(A)); s=zeros(n,1); u=zeros(n,1); c=0; I=eye(n); H=zeros(n,n); Q=I; B=I; for r=1:n-2 s=A(:,r); e=zeros(n,1); for t=1:r s(t,1)=0; end if s'*s==0 %所有分量全为零,取H=I,B保持不变 H=I; else if sign(A(r+1,r))==0 c=sqrt(s'*s); else c=-sign(A(r+1,r))*sqrt(s'*s); end e(r+1,1)=1; u=s-c*e; H=I-2*u*u'/(u'*u); end B=H*A*H; A=B; Q=Q*H; end toc; end |
测试矩阵的拟三角化:
>> A=[1 3 4;3 2 1;4 1 3]
A =
1 3 4
3 2 1
4 1 3
>> [q,b]=hess2(A)
Elapsed time is 0.354213 seconds.
q =
1.0000 0 0
0 -0.6000 -0.8000
0 -0.8000 0.6000
b =
1.0000 -5.0000 -0.0000
-5.0000 3.6000 -0.2000
-0.0000 -0.2000 1.4000
>> I=q*q'
I =
1.0000 0 0
0 1.0000 0.0000
0 0.0000 1.0000
>> q*A*q
ans =
1.0000 -5.0000 -0.0000
-5.0000 3.6000 -0.2000
-0.0000 -0.2000 1.4000
如果使用MATLAB自带的函数hess:
>> [p,t]=hess(A)
p =
0.2425 -0.9701 0
-0.9701 -0.2425 0
0 0 1.0000
t =
0.5294 2.8824 0
2.8824 2.4706 -4.1231
0 -4.1231 3.0000
>> I=p*p'
I =
1.0000 0.0000 0
0.0000 1.0000 0
0 0 1.0000
>> p*A*p
ans =
0.5294 2.8824 -0.0000
2.8824 2.4706 -4.1231
-0.0000 -4.1231 3.0000
可见Hessenberg的转换不唯一,不同的HouseHolder矩阵将会产生不同的拟三对角矩阵