Algorytm rozkładu QR przy użyciu rotacji Givens

Question

Kodowanie algorytmu dekompozycji QR w MATLAB, aby upewnić się, że mam poprawną mechanikę. Oto kod dla funkcji 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 

Funkcja sub GivensRotation podano poniżej:

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 

Zrobiłem badania i jestem całkiem pewny, że jest to jeden z najprostszych sposobów na do wdrożenia tego rozkładu, szczególnie w MATLAB. Ale kiedy testuję to na macierzy A, wytworzony R nie jest trójkątny, jak powinien. Q jest ortogonalne, a Q * R = A, więc algorytm robi pewne rzeczy dobrze, ale nie produkuje dokładnie poprawnej faktoryzacji. Być może właśnie patrzyłem na problem zbyt długo, ale doceniłbym każdy wgląd w to, co przeoczyłem.

Answer 1

ten wydaje się mieć więcej błędów, co widzę:

1. Powinieneś raczej użyć [c, s] = GivensRotation (R (i-1, j), R (i, j)) ; (uwaga, R zamiast A)
2. Oryginalne multiplikacje Q = Q * G '; R = G * R są prawidłowe.
3. r = a/b i r = b/a (bez minus, nie wiem czy to ważne)
4. G ([i-1, i], [i-1, i]) = [c -s ; s c];

Oto przykład kodu, wydaje się działać.

pierwszy plik,

% 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 

a drugi

% 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