B. MATLAB で QR変換して遊んだときのガラクタ箱

(今度講義する羽目になったらちゃんと説明を書こう。)

% 『三重対角行列の QR 変換は三重対角ではないが、Hessenberg 形である。』
% つまり『Hessenberg の QR 変換は Hessenberg』ということかな？
% 三重対角行列のサンプルを作るには diag() が利用できる。
n=4; d=rand(n,1); u=rand(n-1,1); l=rand(n-1,1);
a=diag(d,0)+diag(u,1)+diag(l,-1)
for i=1:100
  [q r]=qr(a);
  a=r*q
end

% ちなみに関数にしておくと、
% function a = randtrid(n)
%   a=diag(rand(n,1),0)+diag(rand(n-1,1),1)+diag(rand(n-1,1),-1);

% 『実対称行列の QR 変換は実対称行列である。』
% 実対称行列のサンプルを作るには、三角部分の切り出し tril(), triu() の
% 利用も考えられるが、良く知られた命題
% 『任意の行列 A について、行列の「対称部分」 (A+A^T)/2 は対称行列になる。』
% を使うのが簡単だろう。
n=4;a=rand(n,n);a=(a+a')/2
for i=1:100
  [q r]=qr(a);
  a=r*q
end

% 『実対称三重対角行列の QR 変換は三重対角行列である。』
% (∵対称なHessenberg行列は三重対角であるから)
n=4; d=rand(n,1); u=rand(n-1,1); a=diag(d,0)+diag(u,1)+diag(u,-1)
for i=1:100
  [q r]=qr(a);
  a=r*q
end

% 最終的な結果と、eig() を使って得た近似固有値との比較
% a を変換してしまう前に近似固有値を計算して保存しておく必要がある。
n=4; d=rand(n,1); u=rand(n-1,1); a=diag(d,0)+diag(u,1)+diag(u,-1)
e=eig(a);
for i=1:100
  [q r]=qr(a);
  a=r*q
end
a
e

% 毎回、対角成分を eig() で計算した近似固有値と並べて見る
% 遅くなるので、あまりループは使いたくないが
%   d=zeros(n,1);
%   for j=1:n
%     d(j)=a(j,j)
%   end
% で対角成分を収めたベクトル d が得られる。
% e と d を並べて表示するには、
%   [e d]
% とするのが簡単。比較するためには、両方とも大きさの順に並べておくのがよい。
% それには sort() を使えばよい。
n=4;a=rand(n,n);a=(a+a')/2
e=sort(eig(a));
for i=1:100
  [q r]=qr(a);
  a=r*q
  d=zeros(n,1);
  for j=1:n
    d(j)=a(j,j);
  end
  d=sort(d);
  format long
  [e d]
  norm(d-e)
  format short
end

% 収束の様子を対角線の下側の成分が小さくなることで確認してみよう。
% 上三角部分 (対角線込み) を取り出す triu() を用いて
n=4;a=rand(n,n);a=(a+a')/2
e=sort(eig(a));
for i=1:100
  [q r]=qr(a);
  a=r*q
  norm(a-triu(a))
end
% とするのが考えられる。
% 元の行列が対称の場合は、非対称部分の大きさを見るのも良いだろう。
% 対角部分の切り出しは、a .* eye(n,n) で出来るので、d-a のノルムを見てみる。
n=4;a=rand(n,n);a=(a+a')/2
e=sort(eig(a));
for i=1:100
  [q r]=qr(a);
  a=r*q
  d=a .* eye(n,n);
  norm(d-a)
end

% 結果が大きいとき、それを読むために、less のようなぺージャーが呼び出
% されて、読み終わった後、画面から消えてしまうが、それをファイルに残し
% たい場合は、(END) が出ているときに、s と打ってからLOGファイルの名前
% を入力する。

%
function ret = rand_hessenberg(n)
  ret = triu(rand(n,n))w + diag(rand(n-1,1),-1);
%
function QR = qr_iteration(a,n)
  QR = a;
  for i=1:n
    [q r]=qr(QR);
    QR=r*q;
  end
%
a=rand_hessenberg(10)
a=qr_iteration(a,1)
% 後は繰り返し、というのもいいかも
-----------------------------------------------------------------
function ret = rand_h(n)
  ret = triu(rand(n,n)) + diag(rand(n-1,1),-1);

function ret = rand_t(n)
  ret=diag(rand(n,1),0)+diag(rand(n-1,1),1)+diag(rand(n-1,1),-1);

function ret = rand_st(n)
  u=rand(n-1,1);
  ret=diag(rand(n,1),0)+diag(u,1)+diag(u,-1);

function ret = rand_s(n)
  a=rand(n,n);
  ret = (a+a')/2;

function QR = qr_iteration(a,n)
  QR = a;
  for i=1:n
    [q r]=qr(QR);
    QR=r*q;
  end

function n = norm_l(a)
  n = norm(a-triu(a));