我正在尝试编写一个程序,该程序获取mxn矩阵并对其进行QR分解
它还没有完成,但我遇到了一个问题。我试着用维基百科上的例子运行我的程序http://en.wikipedia.org/wiki/QR_decomposition
A=[12,-51,4;6,167,-68;-4,24,-41]
他们称之为Q1,Q2…我称之为Qtemp。每次我计算Qtemp时,我都会把它打印出来,以便看到我得到和维基百科相同的结果。第一季度我会,但第二季度我不会。
他们的Q2和我的Q2价值相同,但符号不同。他们有a+的地方我有a-,他们有a-的地方,我有a+。
这是我的代码:
Q=eye(m);
R=A;
for i=1:min(m-1,n)
ei=zeros(n,1);
ei(i,1)=1;
x=A(:,i);
for j=1:i-1
x(j,1)=0;
end
u=x-norm(x)*ei;
v=u/norm(u);
Qtemp=eye(m)-2*(v*v');
A=Qtemp*A;
disp(Qtemp);
end
我只是复制了他们的算法并将其翻译成代码,但第二个Qtemp的输出仍然很糟糕。
看起来并没有在每次迭代中减少块的大小。一切似乎都是相同的m
和n
的函数(您没有在代码中定义它们)。请参阅维基百科页面上他们定义A′并使用它来构建Q2(仅较低的三分之二)。以下是我的一些代码,适用于执行3乘3矩阵的QR分解,这可能会有所帮助。特别注意,第二块仅在A(:,2)
和q(2:3,:)
:上工作
function [q,r]=qr3(A)
u = A(:,1);
u(1) = u(1)-(1-2*(u(1)<0))*norm(u); % Flip < to > to match sign convention of qr
u = u/norm(u);
u(~isfinite(u)) = sqrt(3)/3;
q = -2*(u*u');
q([1 5 9]) = q([1 5 9])+1;
u = q(2:3,:)*A(:,2);
u(1) = u(1)-(1-2*(u(1)<0))*norm(u); % Flip < to > to match sign convention of qr
u = u/norm(u);
u(~isfinite(u)) = sqrt(2)/2;
q(:,2:3) = q(:,2:3)*[1-2*u(1)^2 -2*u(1)*u(2);
-2*u(1)*u(2) 1-2*u(2)^2];
r = triu(q'*A);
上面的代码和维基百科上详细介绍的方法使用了与Matlab的qr
函数不同的符号约定。请参阅我在代码中的注释,了解如何翻转符号。