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);