5
$\begingroup$

For the qr factorization using classic Gram-Schmidt algorithm, I found the 2 different implementations below. The first one uses the for loop to compute the upper triangular matrix R, but the second one uses the matrix-vector multiplication. Since they are mathematically equivalent, but they produce different results. I cannot find why it is so, can anyone give me some indications? Furthermore, I often find that whenever I understand an algorithm theoretically, if I implemented it, I will have some trouble because of different numerical results. This is why I always ask questions about how, actually, to implement the algorithm with specific code. Because even if I know the algorithm, I still cannot generate the correct results. Thanks.

%   test the classical Gram-Schmidt algorithm
clc;clear;format compact;
A = hilb(7);
%%  method 1
R(1,1) = norm(A(:,1));
Q(:,1) = A(:,1)/R(1,1);
m = size(A,2);
for j=2:m
    for i=1:j-1
        R(i,j) = Q(:,i)'*A(:,j);%   the difference between 2 methods from below
    end
    q_hat = A(:,j)-Q(:,1:j-1)*R(1:j-1,j);
    R(j,j) = norm(q_hat);
    Q(:,j) = q_hat/R(j,j);
end
%%  method 2
R1(1,1) = norm(A(:,1));
Q1(:,1)=A(:,1)/R1(1,1);
for j=2:m
    R1(1:j-1,j) = Q1(:,1:j-1)'*A(:,j);  %   the difference between 2 methods from above
    q_hat = A(:,j)-Q1(:,1:j-1)*R1(1:j-1,j);
    R1(j,j) = norm(q_hat);
    Q1(:,j) = q_hat/R(j,j);
end
%%  compare
norm(A-Q*R)
norm(Q'*Q-eye(m))
norm(A-Q1*R1)
norm(Q1'*Q1-eye(m))

My results are as follows:

ans =
     0
ans =
    0.3369
ans =
   8.9238e-10
ans =
    0.1890
$\endgroup$
5
  • $\begingroup$ I haven't really looked at your code yet, but running it I get results different from yours (ans = 3.9252e-17; ans = 0.0867; ans = 3.9252e-17; ans = 0.0867), and they seem to match almost exactly. Maybe you forgot to clear your workspace or did something else before running it that may have affected the results? $\endgroup$ – cliesens Jan 17 '20 at 0:14
  • $\begingroup$ ...though after running it with larger dimensions for A, the results seem to get more and more different the bigger A is. Unless I'm missing something, both methods do exactly the same thing, but the second one vectorizes the inner for loop in the first example. $\endgroup$ – cliesens Jan 17 '20 at 0:22
  • $\begingroup$ @cliesens, Thanks for your reply, but I have added the command clear to clear all the variables in the workspace, my matlab is 2018b, so I am afraid that your implementation may be something wrong. could you have a check? $\endgroup$ – sunshine Jan 17 '20 at 2:51
  • 2
    $\begingroup$ i'm really not sure about this, but i think it might have something to do with the instability of the classic Gram-Schmidt algorithm. You can see that both Q and Q1 are pretty far from orthogonal. That instability might lead the very tiny differences that come from the different order of operations to add up to large errors. $\endgroup$ – Thijs Steel Jan 17 '20 at 16:31
  • $\begingroup$ @cliesens, Thanks very much, I get it. I have used the modified Gram-Schmidt algorithm and it works well. So, it is indeed because of the instability of the classic Gram-Schmidt algorithm. Thanks very much. hope to communicate with you often. $\endgroup$ – sunshine Jan 18 '20 at 3:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.