0

I have following problem.

My task is to fit a polynomial to the data. I want to implenet QR algorithm using Gram-Schimdt orthogonalization process. It is built in this function:

function [ Q,R ] = QRDec( A )
n = length(A(1,:));
for i=1:n
    Q(:,i) = A(:,i);
    for j=1:i-1
        R(j,i) = (Q(:,j)')*Q(:,i);
        Q(:,i) = Q(:,i)-R(j,i)*Q(:,j);
    end
    R(i,i) = norm(Q(:,i),2);
    if R(i,i) == 0
        break;
    end
    Q(:,i)=Q(:,i)/R(i,i);
end
end

Matrices Q,R are almost the same as these Q,R which are obtained from implemented in MatLab function. The only difference is in signs. If I solve my system of equations R*x=Q*y with MatLab functions, I get exact solution. But if I use my own matrices Q and R, then I get wrong result. Can anybody tell me where is the problem in my method? I also enclose code of my script.

% clear variables
clear; clc;
N = 100;
p = ones(1,15);
d = 14;
x = linspace(0,1,N)';
y = polyval(p,x);
A = zeros(N,d+1);
for i = 1 : d+1
    A(:,i) = x.^(i-1);
end
[Qm,Rm] = QRDec(A);
[Q,R] = qr(A,0);
a_qrm = Rm\(Qm'*y);
a_qr = R\(Q'*y);
end

Do you think that such a big mistake can be caused by computation errors? I am really desperate because it seems that I have two same linear systems of equations and the solutions are different.

chip
  • 1
  • 1
  • 1
  • Could it be that your decomposition struggles with a matrix that is not square? If I run you code on [the example from the Wikipedia entry](http://en.wikipedia.org/wiki/QR_decomposition) I get the correct result. That is a `3x3` matrix. However, `size(A)` in your case gives `100x15`. – Schorsch Oct 14 '13 at 18:08
  • It is more "difficult" for MatLab to compute such system of linear system of equations. But still I don't have clue why MatLabs Q,R matrices gives the right solution and my matrices Q,R (although they are more or less the same as the previous ones) gives spoiled coefficients. It just isn't clear to me. – chip Oct 14 '13 at 19:14

1 Answers1

1

The Gram-Schmidt process in the form you implemented it is numerically unstable. In fact, your Q and Qm computed by Matlab are not the same.Furthermore your matrix is ill-conditioned, its condition number is >10^10.

This causes small errors to be magnified and could explain the effects you see.

Johanna
  • 195
  • 1
  • 10