5.2 MATLAB での実験

まず補題 5.2 (1) にあるように、 与えられた ($n$ 次) 実対称行列 $A$ と $p$ , $q\in \{1,2,\dots,n\}$ , $p\ne q$ に対して、 $B:=U^T A U$ , $U:=G(p,q,\theta)$ を計算する関数は次のようになる ($\theta$ を与える代りに $\cos\theta$ , $\sin\theta$ を与えるようにしてある)。

jacobi0.m

   1 % jacobi0.m
   2 % Jacobi 法の基礎ステップ --- x_p,x_q 平面の角θの回転 (Givens 変換)
   3 %  c=cosθ, s=sinθ
   4 function b=jacobi0(a,p,q,c,s)
   5   b=a;
   6   b(p,:)=a(p,:)*c-a(q,:)*s;
   7   b(q,:)=a(p,:)*s+a(q,:)*c;
   8   b(:,p)=b(p,:)';
   9   b(:,q)=b(q,:)';
  10   b(p,p)=a(p,p)*c*c-2*a(p,q)*s*c+a(q,q)*s*s;
  11   b(q,q)=a(p,p)+a(q,q)-b(p,p);
  12   b(p,q)=(a(p,p)-a(q,q))*c*s+a(p,q)*(c*c-s*s);
  13   b(q,p)=b(p,q);

$b_{pq}=b_{qp}=0$ とするためには、 補題 5.2 (2) にあるように $\theta$ を 選ばねばならないが、 これについては Rutishauser の計算式を一部利用すると、 次のような関数を得る。

jacobi1.m

   1 % jacobi1.m
   2 % Jacobi 法の1ステップ --- (p,q)成分を叩く (0 になるよう Given 変換をする)
   3 function b=jacobi1(a,p,q)
   4   % cosθ,sinθを得る --- Rutishauser の計算式
   5   z=(a(q,q)-a(p,p))/(2*a(p,q));   % cot 2θ
   6   t=sign(z)/(abs(z)+sqrt(1+z^2)); % tanθ
   7   c=1/sqrt(1+t^2);                % cosθ
   8   s=c*t;                          % sinθ
   9   % u=s/(1+c); % tan(θ/2)
  10   %
  11   b=jacobi0(a,p,q,c,s);
  12   % 丸め誤差により 0 になっていないのを強制的にクリア
  13   b(p,q)=0;
  14   b(q,p)=0;

特別巡回 Jacobi 法は次のようになる。

jacobi.m

   1 % jacobi.m
   2 % 巡回 Jacobi 法
   3 function lambda=jacobi(a,iter)
   4   [n,n]=size(a);
   5   for k=1:iter
   6     for p=1:n-1
   7       for q=p+1:n
   8         a=jacobi1(a,p,q);
   9       end
  10     end
  11     %a
  12   end
  13   lambda=sort(diag(a));

次の実験プログラムは、 与えられた自然数 $n$ に対して、 $n$ 次実対称行列 $A$ を作り、 繰り返しの数を変えながら巡回 Jacobi 法による近似固有値を求めて、 それを信頼できる値 (MATLAB の eig の結果) と比較している。

jacobiexp.m

   1 % jacobiexp.m
   2 % 
   3 function jacobiexp(n)
   4   % n次実対称行列を生成
   5   a=rand(n,n);
   6   a=(a+a')/2;
   7   % チェック用にMATLAB の関数で固有値を求めておく
   8   eigen=eig(a);
   9   % タイマーをリセット
  10   tic
  11   for i=1:100
  12     % jacobi() で i 回反復したときの近似固有値の精度を求める
  13     error=norm(jacobi(a,i)-eigen);
  14     % 結果を表示
  15     fprintf('%d %e\n', i, error);
  16     % 精度が十分ならば繰り返しをやめる
  17     if (error < 1e-13)
  18       break;
  19     end
  20   end
  21   % 経過時間を表示
  22   toc

実行結果

>> jacobiexp(5)  
1 1.741441e-01
2 7.220988e-02
3 3.632924e-04
4 2.167315e-13
5 2.914633e-15
elapsed_time =
    0.1253
>> jacobiexp(10)
1 2.753837e-01
2 5.652731e-02
3 8.028047e-05
4 3.017277e-09
5 8.183436e-15
elapsed_time =
    0.5129
>> jacobiexp(20) 
1 7.215286e-01
2 1.455863e-01
3 1.153771e-02
4 2.517012e-05
5 7.633347e-11
6 3.977614e-14
elapsed_time =
    6.6447
>>