三种实现矩阵QR分解的算法与程序
2017-07-08 18:46阅读：
To learn, to share, to debate, then comes
progress.
--------------------------------------------------------------------------------------------------
QR分解
若非奇异矩阵A能够分解为一个正交矩阵Q与非奇异上三角矩阵R的乘积，即：
A=QR。则称其为A的QR分解。
实现QR分解的方法有很多种，包括Givens变换法，Householder变换法，Schemit正交化法。具体原理可以参考《矩阵论》（程云鹏，西工大出版）教材，这里仅给出三种实现QR分解的程序。
1.基于Givens变换的QR分解
%****************************************
*******
%矩阵分析与应用-利用Givens变换实现矩阵的QR分解
%输入：A-NxN的矩阵
输出：Q-正交阵
R-上三角阵
17.4 by Howie
function [Q,R]=md_qrg(A)
% if size(A,1)~=size(A,2)
% disp('please input a matrix
by N x N');
% return
% end
[N,M]=size(A);
R=zeros(N);
T=eye(N);
B=A;
for j=1:N-1
Tj=eye(N+1-j);T_t=eye(N);
B=B(min(j,2):end,min(j,2):end);
b=B(:,1);
for i=2:N+1-j
T_temp=eye(N+1-j);
cs=sqrt(b(1)^2+b(i)^2);
c=b(1)/cs;
s=b(i)/cs;
T_temp(1,1)=c;T_temp(i,i)=c;
T_temp(i,1)=-s;T_temp(1,i)=s;
b=T_temp*b;
Tj=T_temp*Tj;
end
B=Tj*B;
R(j,j:end)=B(1,:);
T_t(j:end,j:end)=Tj;
%
T(j:end,j:end)=Tj*T(j:end,j:end);
T=T_t*T;
end
R(j+1,j+1:end)=B(end,end);
Q=T';
2.基于Householder变换的QR分解
%***********************************************
%矩阵分析与应用-利用Householder变换实现矩阵的QR分解
%输入：A-NxN的矩阵
输出：Q-正交阵
R-上三角阵
17.4 by Howie
function [Q,R]=md_qrh(A)
[N,M]=size(A);
B=A;
R=zeros(M);
H=eye(N);
%apply Householder transform on every column of the matrix,
the length
�creases by 1 each step
for j=1:M-1
H_t=eye(N);
B=B(min(2,j):end,min(2,j):end);
Hj=hst(B(:,1));
B=Hj*B;
R(j,j:end)=B(1,:);
H_t(j:end,j:end)=Hj;
H=H_t*H;
end
R(j+1:end,j+1:end)=B(end,end);
Q=H';
%Householder transform matrix
function [H]=hst(x)
xmod=sqrt(sum(x.*x));
alpha=-sign(x(1))*xmod;
x(1)=x(1)+alpha;
u=x/sqrt(sum(x.*x));
H=eye(length(x))-2*(u*u');
3.基于Schmidt正交化法的QR分解
%***********************************************
%矩阵分析与应用-利用Schmidt正交化实现矩阵的QR分解
%输入：A-NxN的矩阵
输出：Q-正交阵
R-上三角阵
17.4 by Howie
function [Q,R]=md_qrs(A)
[N,M]=size(A);
B=zeros(size(A));
Q=B;
C=eye(M);
%初始化
B(:,1)=A(:,1);
%第一组列向量正交化
for j=2:M
b=A(:,j);
for i=1:j-1
%正交化
kij=fkij(A(:,j),B(:,i));
%内积计算系数
b=b-kij*B(:,i);
C(i,j)=kij;
end
B(:,j)=b;
%第j列正交化完成
end
bmod=zeros(1,M);
for j=1:M
%列向量归一化
b=B(:,j);
bmod(j)=sqrt(sum(b.*b));
Q(:,j)=b/bmod(j);
end
R=diag(bmod)*C;
%以模为元素生成对角阵，产生R矩阵
function [k]=fkij(a,b)
k=sum(a.*b)/sum(b.*b);
