4.2 修正 Gram-Schmidt の直交化法

実は前小節にあげた Gram-Schmidt の直交化法は、 浮動小数点演算システムで計算する場合、 丸め誤差の観点からは良くない (と言われている²)。 代りに次にあげる修正 Gram-Schmidt の直交化法が利用される。

修正 Gram-Schmidt の直交化法による QR 分解

   1 % ModifiedGramSchmidt.m --- ModifiedGramSchmidt()
   2 
   3 % 使用例
   4 % n=3; a=rand(n,n);
   5 % [q r]=ModifiedGramSchmidt(a)
   6 % a-q*r
   7 
   8 % A=(a1 a2 ...,an) の列ベクトルから正規直交基底
   9 % q1, q2,..., qn を並べた Q を求める。
  10 function [q,r] = ModifiedGramSchmidt(A)
  11   [n n] = size(A);
  12   %
  13   q=zeros(n,n); r=zeros(n,n);
  14   for j=1:n
  15     r(j,j)=norm(A(:,j)); % 第k列 a_k
  16     q(:,j)=A(:,j)/r(j,j);
  17     for k=j+1:n
  18       r(j,k)=q(:,j)'*A(:,k); % r_{jk} = a_k と q_j の内積
  19       A(:,k)=A(:,k)-r(j,k)*q(:,j);
  20     end
  21   end

ModifiedGramSchmidt.m の実行結果

   1 >> n=4;a=rand(n,n)
   2 a =
   3     0.6273    0.6552    0.5947    0.7764
   4     0.6991    0.8376    0.5657    0.4893
   5     0.3972    0.3716    0.7165    0.1859
   6     0.4136    0.4253    0.5113    0.7006
   7 >> [q r]=ModifiedGramSchmidt(a)
   8 q =
   9     0.5700   -0.2729   -0.7192   -0.2886
  10     0.6352    0.7226    0.2609   -0.0788
  11     0.3609   -0.5861    0.6431   -0.3356
  12     0.3759   -0.2447    0.0323    0.8932
  13 r =
  14     1.1005    1.1995    1.1492    1.0839
  15          0    0.1046   -0.2985   -0.1386
  16          0         0    0.1972   -0.2886
  17          0         0         0    0.3008
  18 >> a-q*r
  19 ans =
  20    1.0e-15 *
  21          0         0   -0.1110         0
  22          0         0         0         0
  23          0         0         0   -0.0555
  24          0         0         0         0
  25 >> 

(「良くない」のはなぜか？ 本当にそうなのか、どこかでチェックしたい。)