Algoritmo di decomposizione QR Utilizzo delle rotazioni di Givens

Question

Sto codificando un algoritmo di decomposizione QR in MATLAB, giusto per assicurarsi di avere i meccanismi corretti. Ecco il codice per la funzione principale:

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 

La funzione GivensRotation sub è il seguente:

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 

ho fatto ricerche e sono abbastanza sicuro che questo è uno dei modi più semplici per implementare questa scomposizione, specialmente in MATLAB. Ma quando lo provo su una matrice A, la R prodotta non è triangolare come dovrebbe essere. La Q è ortogonale e Q * R = A, quindi l'algoritmo sta facendo alcune cose bene, ma non produce esattamente la fattorizzazione corretta. Forse ho appena osservato il problema troppo a lungo, ma qualsiasi intuizione su ciò che ho trascurato sarebbe stata apprezzata.

Answer 1

questo sembra avere più bug, quello che vedo:

1. Si dovrebbe piuttosto usare [c, s] = GivensRotation (R (i-1, j), R (i, j)) ; (nota, R anziché A)
2. Le moltiplicazioni originali Q = Q * G '; R = G * R sono corretti.
3. r = a/b er = b/a (senza meno, non sono sicuro se è importante)
4. G ([i-1, i], [i-1, i]) = [c -s ; s c];

Ecco un esempio di codice, sembra funzionare.

primo 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 

e la seconda

% 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