Givens rotation QR decomposition

138 views (last 30 days)
Duc Anh Le
Duc Anh Le on 11 Feb 2020
Commented: MmO on 30 Sep 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 Comments
Benjamin Ellis
Benjamin Ellis on 10 Mar 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 on 30 Sep 2021
Hello, Where you able to find the mistake? Where did you take this algorithm from? Best regards

Sign in to comment.

Answers (1)

Jon
Jon on 11 Feb 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.

Categories

Find more on Creating and Concatenating Matrices in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!