4.3 Givens 変換による QR 分解

(説明を書く暇がないのでとりあえず算法だけでも)

Givens 変換による QR 分解

   1 % QR_Givens.m --- QR_Givens()
   2 % 注意: このプログラムは原理の説明用であり、計算の効率は非常に低い
   3 
   4 % 使用例
   5 % n=3; a=rand(n,n); [q r]=QR_Givens(a)
   6 % norm(a-q*r)
   7 % norm(q*q')
   8 
   9 % A=(a1 a2 ...,an) の列ベクトルから正規直交基底
  10 % q1, q2,..., qn を並べた Q を求める。
  11 function [Q,R] = QR_Givens(A)
  12   [n n] = size(A);
  13   R=A; Q=eye(n,n);
  14   %
  15   for q=1:n-1
  16     for p=q+1:n
  17       d=sqrt(R(p,q)*R(p,q)+R(q,q)*R(q,q));
  18       c=R(q,q)/d;
  19       s=R(p,q)/d;
  20       G=Givens(n,p,q,c,s); % xp, xq 平面でθの回転 (c=cosθ, s=sinθ)
  21       R=G'*R;
  22       Q=Q*G;
  23     end
  24   end

QR_Givens.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]=QR_Givens(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.0000    0.1046   -0.2985   -0.1386
  16     0.0000    0.0000    0.1972   -0.2886
  17     0.0000   -0.0000    0.0000    0.3008
  18 >> a-q*r
  19 ans =
  20    1.0e-15 *
  21          0    0.1110    0.1110    0.1110
  22          0    0.1110    0.1110    0.0555
  23    -0.0555         0         0   -0.0278
  24    -0.0555    0.0555    0.1110         0
  25 >>