使用 Givens 旋转的 QR 分解算法答案

【问题标题】：QR Decomposition Algorithm Using Givens Rotations使用 Givens 旋转的 QR 分解算法
【发布时间】：2014-06-26 21:56:09
【问题描述】：

我正在 MATLAB 中编写 QR 分解算法，只是为了确保我的机制正确。下面是main函数的代码：

    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 中。但是当我在矩阵 A 上测试它时，产生的 R 并不是应有的直角三角形。 Q 是正交的，并且 Q*R = A，所以该算法在做一些正确的事情，但它并没有产生完全正确的分解。也许我只是盯着这个问题太久了，但是任何关于我忽略了什么的见解都将不胜感激。

【问题讨论】：

相信我找到了错误。而不是Q = Q*G'; R = G*R;，我应该写Q = Q*G; R = G'*R 在反转矩阵上的转置时，我实际上是在错误的方向上进行了旋转，从而产生了不同的分解。

标签： matlab qr-decomposition

【解决方案1】：

这似乎有更多的错误，我所看到的：

您应该使用 [c,s] = GivensRotation( R(i-1,j),R(i,j) ); （注意，R 而不是 A）
原始乘法 Q = Q*G'; R = G*R 是正确的。
r = a/b 和 r = b/a（没有减号，不确定是否重要）
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

第二个

% 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

【讨论】：