[解決済み] Givens回転を用いたQR分解アルゴリズム

質問

MATLABでQR分解アルゴリズムをコーディングしているのですが、仕組みが正しいかどうか念のため確認します。 以下はメイン関数のコードです。

    function [Q,R] = QRgivens(A)
        n = length(A(:,1));
        Q = eye(n);
        R = A;

        for j = 1:(n-1)
            for i = n:(-1):(j+1)
                G = eye(n);
                [c,s] = GivensRotation( A(i-1,j),A(i,j) );
                G(i-1,(i-1):i) = [c s];
                G(i,(i-1):i)   = [-s c];
                Q = Q*G';
                R = G*R;
            end
        end
    end

サブ関数 GivensRotation を以下に示す。

    function [c,s] = GivensRotation(a,b)
        if b == 0
            c = 1;
            s = 0;
        else
            if abs(b) > abs(a)
                r = -a / b;
                s = 1 / sqrt(1 + r^2);
                c = s*r;
            else
                r = -b / a;
                c = 1 / sqrt(1 + r^2);
                s = c*r;
            end
        end
    end

調べたところ、特にMATLABでこの分解を実装する最も簡単な方法の1つであることは間違いないようです。 しかし、行列Aでテストしてみると、生成されるRは本来あるべき右三角形ではありません。 Qは直交しており、Q*R = Aなので、このアルゴリズムはいくつかのことを正しく行っていますが、正確に正しい分解を生成しているわけではありません。 おそらく、私はこの問題を長く見つめすぎているのだと思いますが、私が見落としていることについて何か洞察があれば、感謝します。

どのように解決するのですか？

は、バグが多いようです。 を見ることができます。

1. むしろ [c,s] = GivensRotation( R(i-1,j),R(i,j) ); を使うべきでしょう（注意：AではなくRです）。
2. 元の乗算 Q = Q*G'; R = G*R は正しいです。
3. r = a/b と r = b/a (マイナスなし、重要かどうかは不明)
4. G([i-1, i],[i-1, i]) = [c -s; s c];

以下はコードの例で、うまくいっているようです。

最初のファイルです。

% qrgivens.m
function [Q,R] = qrgivens(A)
  [m,n] = size(A);
  Q = eye(m);
  R = A;

  for j = 1:n
    for i = m:-1:(j+1)
      G = eye(m);
      [c,s] = givensrotation( R(i-1,j),R(i,j) );
      G([i-1, i],[i-1, i]) = [c -s; s c];
      R = G'*R;
      Q = Q*G;
    end
  end

end

と、2番目の

% givensrotation.m
function [c,s] = givensrotation(a,b)
  if b == 0
    c = 1;
    s = 0;
  else
    if abs(b) > abs(a)
      r = a / b;
      s = 1 / sqrt(1 + r^2);
      c = s*r;
    else
      r = b / a;
      c = 1 / sqrt(1 + r^2);
      s = c*r;
    end
  end

end