矩陣的擬上三角化（Hessenberg矩陣）
|
%實矩陣的擬上三角分解（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矩陣將會產生不同的擬三對角矩陣