QR Factorization Using Householder Transformations

function [Q,R]=QRfactor(A)
[m,n]=size(A);
R=A; %Start with R=A
Q=eye(m); %Set Q as the identity matrix
for k=1:m-1
x=zeros(m,1);
x(k:m,1)=R(k:m,k);
g=norm(x);
v=x; v(k)=x(k)+g;
%Orthogonal transformation matrix that eliminates one element
%below the diagonal of the matrix it is post-multiplying:
s=norm(v);
if s~=0, w=v/s; u=2*R'*w;
R=R-w*u'; %Product HR
Q=Q-2*Q*w*w'; %Product QR
end
end

for A=[-2 2 3; 1 3 5; -3 -1 2]

I got the answers Q and R different from when I use [Q,R]=qr(A). Where am I wrong with code.

[John D'Errico]

Please learn to use the code button when you post code. It takes only one click of the mouse to do so.

Answer 1

[Titus Edelhofer]

Hi Hüseyin,

I don't think something is wrong. Q*R gives A (at least for your matrix A). Having different Q and R from MATLAB's implementation does not necessarily mean something is wrong (as long as Q*R=A and Q is orthogonal, i.e. Q'*Q = identity).

Titus

[Titus Edelhofer]

BTW, your function gives the same result up to sign: if Q,R is a QR-Decomposition, then -Q, -R is one as well ;-).

[Hüseyin]

I know Q*R=A but some elements of Q and R is different (negative/positive) in my written code.. I'm trying to find how matlab computes [Q,R]=qr(A)

[John D'Errico]

You are not listening. Q and R are not unique. Your code is fine. That it produces elements with sign differences in some cases merely means that an arbitrary choice was made about sign in the MATLAB code that differs from your choice. And since the MATLAB code for QR is proprietary, you can NEVER know exactly what they did.

Answer 2

[Francesco Onorati]

The reason why there are differencies in the sign is that for numerical stability the 2-norm of each vector is taken with the opposite sign of the pivotal element of the vector itself. As you take always the norm as positive, sometimes it is in agreement with MATLAB code, sometimes it is not (here I'm supposing MATLAB uses Housolder transformation to do QR decomposition).

Answer 3

[Davide Poggiali]

You just have to change two lines

    g=-sign(x(k))*norm(x);
    v=x; v(k)=x(k)-g;

to get what you're looking for. source: wiki