QR Decomposition Algorithm Using Givens Rotations

Question

I am coding a QR decomposition algorithm in MATLAB, just to make sure I have the mechanics correct. Here is the code for the main function:

    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

The sub function GivensRotation is given below:

    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

I've done research and I'm pretty sure this is one of the most straightforward ways to implement this decomposition, in MATLAB especially. But when I test it on a matrix A, the R produced is not right triangular as it should be. The Q is orthogonal, and Q*R = A, so the algorithm is doing some things right, but it is not producing exactly the correct factorization. Perhaps I've just been staring at the problem too long, but any insight as to what I've overlooked would be appreciated.

Answer 1

this seems to have more bugs, what I see:

1. You should rather use [c,s] = GivensRotation( R(i-1,j),R(i,j) ); (note, R instead of A)
2. The original multiplications Q = Q*G'; R = G*R are correct.
3. r = a/b and r = b/a (without minus, not sure if it's important)
4. G([i-1, i],[i-1, i]) = [c -s; s c];

Here is an example code, seems to work.

First file,

% 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

and the second

% 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