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

最佳答案

这似乎有更多的错误,
我所看到的:

  • 你应该使用 [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; 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
    

    关于matlab - 使用 Givens 旋转的 QR 分解算法,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/13438073/

    10-12 23:06
    查看更多