Givens rotation QR decomposition

139 visualizaciones (últimos 30 días)
Duc Anh Le
Duc Anh Le el 11 de Feb. de 2020
Comentada: MmO el 30 de Sept. de 2021
Screen Shot 2020-02-11 at 2.08.39 PM.png
I'm trying to create a function that computes the Givens Rotation QR decomposition, following this pseudo-code.
function [Q,R] = givens(A)
[m,n] = size(A);
indexI = zeros(m,n);
indexJ = zeros(m,n);
C = zeros(m,n);
S = zeros(m,n);
for i = 1:n
for j = i+1:m
c = A(i,i)/((A(i,i))^2 + (A(j,i)^2))^0.5;
s = A(j,i)/((A(i,i))^2 + (A(j,i)^2))^0.5;
A(i,:) = c*A(i,:) + s*A(j,:);
A(j,:) = -s*A(i,:) + c*A(j,:);
indexI(j,i) = i;
indexJ(j,i) = j;
C(j,i) = c;
S(j,i) = s;
end
end
R = A;
Q = eye(m);
for i = 1:n
for j= j+1:m
Q(:,i) = c*Q(:,i) + s*Q(:,j);
Q(:,j) = -s*Q(:,i) + c*Q(:,j);
end
end
However, the R matrix, that I get, is not upper triangular. I can't seem to find the mistake here. Any help would be highly appreciated. Thanks in advance.
  2 comentarios
Benjamin Ellis
Benjamin Ellis el 10 de Mzo. de 2020
Hi! I'm in this class too. I added the lines c = C(j,i) and s = S(j,i) within the second for loop. Also the second for loop should iterate j = (i+1):m.
With these changes I got Q and R to agree with qr(A) up to a sign.
MmO
MmO el 30 de Sept. de 2021
Hello, Where you able to find the mistake? Where did you take this algorithm from? Best regards

Iniciar sesión para comentar.

Respuestas (1)

Jon
Jon el 11 de Feb. de 2020
It looks to me like your code reproduces what is in the pseudocode. Are you sure that the psuedocode that you based your code on is correct?
A few things look perhaps questionable about the psuedocode.
Why do we compute C, S, idxI idxJ and never use them?
In the second double loop you use the values c and s which are just the last values from the double loop in the first part of the algorithm. So they never change value as i and j change. Maybe that is ok but it looks strange.

Categorías

Más información sobre Creating and Concatenating Matrices en Help Center y File Exchange.

Etiquetas

Productos

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by