※ 引述《snow3804 (snow3804)》之銘言：
: matlab有hess(A)可以將矩陣A轉換成Hessenberg矩陣
: 但我想知道這中間是怎麼計算的
: 搜尋電腦裡有hess.m檔只是打開卻只看到註解而已
: http://tinyurl.com/m6v7p96
: 有沒有辦法可以知道是怎麼實作的
推 iHakka: matlab是用lapack去實作的，可以追回去看看 08/13 19:27
推 infernodimon: lapack裡有 印象中是用householder 08/14 19:50
我講一下原由,我在這裡查到Hessenberg的matlab程式碼
http://www.auburn.edu/~tamtiny/lecture%2010.pdf
function [H,Q]=houshess(A)
n=max(size(A));
Q=eye(n);
H=A;
for k=1:(n-2)
[v,beta]=vhouse(H(k+1:n,k));
I=eye(k);
N=zeros(k,n-k);
m=length(v);
R=eye(m)-beta*v*v';
HH=H(k+1:n,k:n);
H(k+1:n,k:n)=R*H(k+1:n,k:n);
HH=H(1:n,k+1:n);
H(1:n,k+1:n)=H(1:n,k+1:n)*R;
P=[I, N; N', R];
Q=Q*P;
end
return
function [v,beta]=vhouse(x)
n=length(x);
x=x/norm(x);
s=x(2:n)'*x(2:n);
v=[1; x(2:n)];
if (s==0), beta=0;
else
mu=sqrt(x(1)^2+s);
if (x(1) <= 0)
v(1)=x(1)-mu;
else
v(1)=-s/(x(1)+mu);
end
beta=2*v(1)^2/(s+v(1)^2);
v=v/v(1);
end
return
但執行的結果和matla內建的hess()結果不一樣
http://www.mathworks.com/help/matlab/ref/hess.html
A=[-149 -50 -154;537 180 546;-27 -9 -25];
H=houshess(A)
H =
-149.0000 -42.2037 156.3165
537.6783 152.5511 -554.9272
0 0.0728 2.4489
H=hess(A)
H =
-149.0000 42.2037 -156.3165
-537.6783 152.5511 -554.9272
0 0.0728 2.4489
看得出來其實只是有些地方差個負號而已
當然這兩個都是Hessenberg矩陣
但我還是想知道程式碼是差在哪裡造成正負號不同
lapack我努力去找到程式碼,只是lapack是用fortran寫的,但我看不懂程式碼
http://www.netlib.org/lapack/explore-3.1.1-html/sgehrd.f.html