QR Factorization Using Householder Transformations

123 次查看(过去 30 天)
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.
  1 个评论
John D'Errico
John D'Errico 2015-1-12
Please learn to use the code button when you post code. It takes only one click of the mouse to do so.

请先登录,再进行评论。

采纳的回答

Titus Edelhofer
Titus Edelhofer 2015-1-12
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
  3 个评论
Hüseyin
Hüseyin 2015-1-12
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
John D'Errico 2015-1-12
编辑:John D'Errico 2015-1-12
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.

请先登录,再进行评论。

更多回答(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).

Davide Poggiali
Davide Poggiali 2020-4-20
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

类别

Help CenterFile Exchange 中查找有关 Matrix Indexing 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by