Thuật toán phân tách QR bằng cách sử dụng Givens Rotations

Question

Tôi đang mã hóa thuật toán phân tách QR trong MATLAB, chỉ để đảm bảo rằng tôi có cơ chế chính xác. Dưới đây là mã cho các chức năng chính:

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 

Chức năng GivensRotation phụ được đưa ra dưới đây:

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 

tôi đã thực hiện nghiên cứu và tôi khá chắc chắn đây là một trong những cách đơn giản nhất để thực hiện phân tách này, đặc biệt là trong MATLAB. Nhưng khi tôi thử nghiệm nó trên một ma trận A, R được tạo ra không phải là hình tam giác bên phải. Q là trực giao, và Q * R = A, vì vậy thuật toán đang làm một số điều đúng, nhưng nó không tạo ra chính xác yếu tố chính xác. Có lẽ tôi đã chỉ nhìn chằm chằm vào vấn đề quá lâu, nhưng bất kỳ cái nhìn sâu sắc như những gì tôi đã bỏ qua sẽ được đánh giá cao.

Answer 1

này dường như có nhiều lỗi, những gì tôi thấy:

1. Bạn khá nên sử dụng [c, s] = GivensRotation (R (i-1, j), R (i, j)) ; (lưu ý, R thay vì A)
2. Các phép nhân gốc Q = Q * G '; R = G * R là chính xác.
3. r = a/b và r = b/a (không có dấu trừ, không chắc chắn nếu nó quan trọng)
4. G ([i-1, i], [i-1, i]) = [c -s ; s c];

Dưới đây là mã ví dụ, có vẻ như hoạt động.

tập đầu tiên,

% 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 

và thứ hai

% 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