A. MATLAB で Householder 法プログラムを書いて遊ぶ

ガラクタ箱だけど、後で使おう。

naive_householder.m

function H=naive_householder(A)
[N N]=size(A);
for k=1:N-2
 b=A(k+1:N,k);
 B=A(k+1:N,k+1:N);
 s=-sign(b(1))*norm(b);
 v=b;
 v(1)=v(1)-s;
 v=v/norm(v);
 Q=eye(N-k,N-k)-2*v*v';
 B=Q*B*Q;
 b=zeros(N-k,1);
 b(1)=s;
 A(k+1:N,k)=b;
 A(k,k+1:N)=b';
 A(k+1:N,k+1:N)=B;
end
H=A;

householder.m

function H=householder(A)
[N N]=size(A);
for k=1:N-2
 b=A(k+1:N,k);
 b1=b(1);
 B=A(k+1:N,k+1:N);
 s=-sign(b1)*norm(b);
 w=b;
 w(1)=w(1)-s;
 norm_w=2*(s*s-b1*s);
 p=B*w/(norm_w/2);
 alpha=w'*p/norm_w;
 q=p-alpha*w;
 B=B-w*q'-q*w';
 b=zeros(N-k,1);
 b(1)=s;
 A(k+1:N,k)=b;
 A(k,k+1:N)=b';
 A(k+1:N,k+1:N)=B;
end
H=A;

householder_exp2.m

%format long
A=[1 1 1 1;
   1 2 2 2;
   1 2 3 3;
   1 2 3 4]
householder(A)

householder_exp2.m の実行結果

>> householder_exp2

A =

     1     1     1     1
     1     2     2     2
     1     2     3     3
     1     2     3     4


ans =

    1.0000   -1.7321         0         0
   -1.7321    7.6667    1.2472         0
         0    1.2472    0.9762   -0.1237
         0         0   -0.1237    0.3571

>> 

少し大きめの行列を作って、 それを householder() で三重対角化してみた。 ついでに固有値を計算して元の行列と変化がないことを確かめた (まあ相似変換が正しいことの状況証拠)。 そのうち三重対角行列の固有値の計算をするプログラムを書こう。

householder_exp3.m

format short
N=10
A=zeros(N,N);
for i=1:N
   t=i*ones(N-i+1,1);
   A(i:N,i)=t;
   A(i,i:N)=t';
end
A
H=householder(A)
e1=sort(eig(A))
e2=sort(eig(H))
norm(e1-e2)

householder_exp3.m の実行結果